精华内容
下载资源
问答
  • 2021-07-10 00:03:20


    一、ode45?

    引用以下百度百科的内容:ode45,常微分方程的数值求解。MATLAB提供了求常微分方程数值解的函数。当难以求得微分方程的解析解时,可以求其数值解,Matlab中求微分方程数值解的函数有七个:ode45,ode23,ode113,ode15s,ode23s,ode23t,ode23tb 。
    ode是Matlab专门用于解微分方程的功能函数。该求解器有变步长(variable-step)和定步长(fixed-step)两种类型。不同类型有着不同的求解器,其中ode45求解器属于变步长的一种,采用Runge-Kutta算法;其他采用相同算法的变步长求解器还有ode23。
    ode45表示采用四阶-五阶Runge-Kutta算法,它用4阶方法提供候选解,5阶方法控制误差,是一种自适应步长(变步长)的常微分方程数值解法,其整体截断误差为(Δx)^5。解决的是Nonstiff(非刚性)常微分方程。
    ode45是解决数值解问题的首选方法,若长时间没结果,应该就是刚性的,可换用ode15s试试。
    这就是重点啦,ode45是使用了龙格库塔方法四阶五阶方法的一种求解常微分方程的方法,也是目前最常用的一种数值解计算方案。

    二、案例解析

    1.求解一阶常微分方程

    考虑常微分方程为:
    220sin( ω ( t ) \omega(t) ω(t)) = Rs * y + Ls * dy + j * ω \omega ω * y
    除y和dy之外所有变量都已知,将方程进行变形,得到dy的表达式:
    dy = ( 220sin( ω ( t ) \omega(t) ω(t)) - Rs * y - j * ω \omega ω * y ) / Ls
    然后用odefun定义一个dy的表达式,格式如下:

    odefun=@(t,y) (220*sin(wt) - Rs * y - j * w * y) / Ls; %定义dy表达式
    tspan=[1 4]; %求解区间
    y0=-2; %初值
    [t,y]=ode45(odefun,tspan,y0);
    plot(t,y) %作图
    xlabel('t')
    ylabel('dy')
    

    2.求解二阶及高阶常微分方程

    需要求解的高阶常微分方程:y’’ = -ty + e^ty’ + 3sin2t
    求解的关键是将高阶转为一阶,odefun的书写.
    F(y,y’,y’’…y(n-1),t)=0用变量替换,y1=y,y2=y’…注意odefun方程定义为行向量
    dxdy=[y(1),y(2)…]

    function Testode45
    tspan=[3.9 4.0]; %求解区间
    y0=[2 8]; %初值
    [t,x]=ode45(@odefun,tspan,y0);
    plot(t,x(:,1),'-o',t,x(:,2),'-*')
    legend('y1','y2')
    title('y'' ''=-t*y + e^t*y'' +3sin2t')
    xlabel('t')
    ylabel('y')
    function y=odefun(t,x)
    y=zeros(2,1); % 列向量
    y(1)=x(2);
    y(2)=-t*x(1)+exp(t)*x(2)+3*sin(2*t);
    end
    end
    

    更多相关内容
  • 朋友问题: 有微分方程如下:m d 2 y d t 2 + d y d t e x p ( t ) − y 2 = 5 m \frac{d^2y}{dt^2} + \frac{dy}{dt} exp(t) - y^2 = 5mdt2d2y​+dtdy​exp(t)−y2=5其中,y ( t = 0 ) = 1 y(t=0)=1y(t=0)=1,d y / ...

    朋友问题: 有微分方程如下:

    m d 2 y d t 2 + d y d t e x p ( t ) − y 2 = 5 m \frac{d^2y}{dt^2} + \frac{dy}{dt} exp(t) - y^2 = 5mdt2d2y​+dtdy​exp(t)−y2=5

    其中,y ( t = 0 ) = 1 y(t=0)=1y(t=0)=1,d y / d t ( t = 0 ) = − 10 dy/dt (t=0) = -10dy/dt(t=0)=−10。

    请在区间[ 0 , 5 ] [0, 5][0,5]内绘制两个子图,分别为d y / d t dy/dtdy/dt与y yy,每个子图涵盖m = 1 m=1m=1与m = 2 m=2m=2两种情况。

    所以本题的核心问题在于:用数值计算的方法求解该方程,得到各点,绘制点图。

    使用 matlab 自带的 ode45 ,方程组用句柄表示。

    ode45 参见教程:如何使用ODE45

    首先把题目方程转换,转换为 ode45 能理解的方式。

    先声明变量:

    y 1 = y y 2 = y ′ \begin{aligned} y_1 & = y \\ y_2 & = y' \\ \end{aligned}y1​y2​​=y=y′​

    于是整理方程:

    y 1 ′ = y 2 y 2 ′ = 1 m ( 5 − y 1 ′ e y 1 + y 1 2 ) \begin{aligned} y_1' & = y_2 \\ y_2' & = \frac{1}{m} (5 - y_1' e^{y_1} + y_1^2) \end{aligned}y1′​y2′​​=y2​=m1​(5−y1′​ey1​+y12​)​

    于是我们知道,ode45中要有2个变量,且将其右边的式子分别表示出来,即:

    dy = @(t, y)[y(2); (5 - y(2)*exp(y(1)) + y(1)^2)/m];

    其中:

    y(2)代表 y 1 ′ = y 2 y_1' = y_2y1′​=y2​;

    (5 - y(1)*exp(y(2)) + y(2)^2)/m代表 y 2 ′ = 1 m ( 5 − y 1 ′ e y 1 + y 1 2 ) y_2' = \frac{1}{m} (5 - y_1' e^{y_1} + y_1^2)y2′​=m1​(5−y1′​ey1​+y12​)。

    接着,规定初值:y ( t = 0 ) = 1 y(t=0)=1y(t=0)=1,d y / d t ( t = 0 ) = − 10 dy/dt (t=0) = -10dy/dt(t=0)=−10。

    y10 = 1;

    y20 = -10;

    规定自变量 t tt 范围:

    tspan = [0, 5];

    输入 ode45 则为:

    [t, y] = ode45(dy, tspan, [y10, y20]);

    整个题目的代码为:

    % 表示该方程组

    m = 1;

    dy = @(t, y)[y(2); (5 - y(2)*exp(y(1)) + y(1)^2)/m];

    y10 = 1;

    y20 = -10;

    tspan = [0, 5];

    % m = 1

    [t_m_1, y_m_1] = ode45(dy, tspan, [y10, y20]);

    % m = 2

    m = 2;

    dy = @(t, y)[y(2); (5 - y(2)*exp(y(1)) + y(1)^2)/m];

    [t_m_2, y_m_2] = ode45(dy, tspan, [y10, y20]);

    % plot

    subplot(1, 2, 1);

    plot(t_m_1, y_m_1(:, 2));

    hold on

    plot(t_m_2, y_m_2(:, 2));

    title('dy/dt')

    legend('m=1','m=2')

    subplot(1, 2, 2);

    plot(t_m_1, y_m_1(:, 1));

    hold on

    plot(t_m_2, y_m_2(:, 1));

    title('y')

    legend('m=1','m=2')

    ae93d99a1ef4c4496bcb76a96d246247.png

    顺便学了 ode45 ,不错。

    展开全文
  • 详见mathwork 函数句柄的创建 1、 直接@ ...[t,y] = ode45(odefun,tspan,y0) t:自变量 y(t):因变量 odefun:函数句柄(自动y’=) tspan:t的范围,如[0,10] y0:初值 1、一阶方程 a=@...

    详见mathworks

    龙格库塔方法

    在这里插入图片描述
    写成矩阵(状态方程)的形式更简洁一点(其实这两种方法结果是一样的,如果C是[1,0,0]的话,就很明显了)
    在这里插入图片描述
    例如:求系统在0-5s内的单位阶跃相应,已知传递函数:
    在这里插入图片描述
    状态方程:
    在这里插入图片描述

    clear
    %状态方程
    A=[0,1
        -100,-10];
    B=[0;
        1];
    C=[100,0];
    %初值
    x=[0;0];y=0;t0=0;tf=5;  
    h=0.05;%步长
    t=t0;%从t0开始迭代,储存时间
    N=round((tf-t0)/h);     %迭代步数
    r=1;%单位阶跃输入
    for i=1:N   
        k1=A * x+B*r;
        k2=A * (x+h*k1/2)+B*r;
        k3=A * (x+h*k2/2)+B*r;
        k4=A * (x+h*k3)+B*r;
        x=x+h*(k1+2*k2+2*k3+k4)/6;                    %采用四阶龙格库塔法
        y=[y,C*x];                                    %输出值
        t=[t,t(i)+h];
    end
    plot(t,y);
    

    结果:
    在这里插入图片描述

    • matlab中传函转化状态空间方程的函数

    构造传函:

    M=100;%分子
    N=[1,10,100];%分母
    G=tf(M,N);%传函
    

    转换:

        g=ss(G);%状态方程
        A=g.A;%各个系数矩阵
        B=g.B;
        C=g.C;
    

    ODE45

    • 函数句柄的创建

    1、 直接@

    A=@sin(x);
    A(pi)
    

    2、创建匿名函数

    %定义函数f(x)=x^2
    f = @(x) x^2;
    f(2)
    
    • ODE45使用
      [t,y] = ode45(odefun,tspan,y0)
      t:自变量
      y(t):因变量
      odefun:函数句柄(自动y’=)
      tspan:t的范围,如[0,10]
      y0:初值

    1、一阶方程
    在这里插入图片描述

    a=@(t,y) 2*t;
    [t,y] = ode45(a, [0 5], 0);
    plot(t,y)
    

    2、二阶方程
    在这里插入图片描述
    先转化为两个一阶微分:
    在这里插入图片描述
    创建m文件:

    function dydt = vdp1(t,y)
    dydt = [y(2); 
    		(1-y(1)^2)*y(2)-y(1)];%括号不可省,不是变量,y是2x1列向量
    

    调用ode45

    [t,y] = ode45(@vdp1,[0 20],[2; 0]);
    plot(t,y)
    

    输出y两列,分别为y(1)和y(2)

    展开全文
  • 对于形如以下的常微分方程: 采用四阶龙格库塔法的计算公式: 四阶 龙格库塔法精度为4,属于单步递推法。 单步递推法的基本思想是从(x(i),y(i))点出发,以某一斜率沿直线达到(x(i+1),y(i+1))点。 二、...

    一、计算公式

    对于形如以下的常微分方程:

    \frac{dy}{dx}=f(x,y)

    y(a)=y0

    采用四阶龙格库塔法的计算公式

    y_{i+1}=y_i+h/6*(K_1+2K_2+2*K_3+K_4)

    K_1=f(x_i,y_i)

    K_2=f(x_i+h/2,y_i+h/2*K_1)

    K_3=f(x_i+h/2,y_i+h/2*K_2)

    K_4=f(x_i+h,y_i+h*K_3)

    四阶 龙格库塔法精度为4,属于单步递推法

    单步递推法的基本思想是从(x(i),y(i))点出发,以某一斜率沿直线达到(x(i+1),y(i+1))点。

    二、实例计算

    从上述定义可以看出,龙格库塔实质上是求一阶微分方程,但是如果将一阶导看作变量,则二阶导也不过是这个变量的一阶导而已。

    接下来就看实例吧:

    对于下述二阶方程:

    {\ddot{q}}+2*\varepsilon *{\dot{q}}+q=Fm+Fah*cos(\Omega *t)

    1、基本思想:

    令位移为    q=u(1)

    q的一阶导,即位移的一阶导(速度)为   \dot{q}=u(2)

    q的二阶导   \dot{u(2))}=-2* \varepsilon *u(2)-u(1)+Fm+Fah*cos(\Omega *t)

    求解位移u(1)的龙格库塔计算方法如下:

    KK1=u(2);
    KK2=u(2)+h/2*KK1;
    KK3=u(2)+h/2*KK2;
    KK4=u(2)+h*KK3;
    u(1)=u(1)+h/6*(KK1+2*KK2+2*KK3+KK4);

    求解速度u(2)的龙格库塔计算方法如下:

    K1=-2*eptheton*u(2)-u(1)+Fm+Fah*wh*wh*sin(wh*tao+fav_h);
    K2=-2*eptheton*(u(2)+h/2*K1)-(u(1)+h/2)+Fm+Fah*wh*wh*sin(wh*tao+fav_h);
    K3=-2*eptheton*(u(2)+h/2*K2)-(u(1)+h/2)+Fm+Fah*wh*wh*sin(wh*tao+fav_h);
    K4=-2*eptheton*(u(2)+h*K3)-(u(1)+h)+Fm+Fah*wh*wh*sin(wh*tao+fav_h);
    u(2)=u(2)+h/6*(K1+2*K2+2*K3+K4);

    2、编程实现

    参数设置:

    h=0.001;       % 步长
    t0=0:h:400;    % 总时长
    w=5;
    ep=0.02;
    Fm=0.1;
    Fah=0.05;
    u(1)=0;u(2)=0;  % 初值

    总的程序实现

    h=0.001;
    t0=0:h:400;
    w=5;
    ep=0.02;
    Fm=0.1;
    Fah=0.05;
    u(1)=0;u(2)=0;
    for i=1:length(t0)   % 进行多次迭代
        tao=t0(i);
        u=RK(u,tao,h,ep,w,Fm,Fah);     
        Result(i,:)=u;   % 将每次迭代的位移和速度保存
    end
    figure(1)
    subplot(2,1,1)
    plot(t0,Result(:,1))      % 绘制位移图
    xlabel('Time')
    ylabel('displacement')
    subplot(2,1,2)
    plot(t0,Result(:,2))      % 绘制速度图
    xlabel('Time')
    ylabel('velocity')
    
    function u=RK(u,tao,h,ep,w,Fm,Fah)
    KK1=u(2);
    KK2=u(2)+h/2*KK1;
    KK3=u(2)+h/2*KK2;
    KK4=u(2)+h*KK3;
    u(1)=u(1)+h/6*(KK1+2*KK2+2*KK3+KK4);
    K1=-2*ep*u(2)-u(1)+Fm+Fah*cos(w*tao);
    K2=-2*ep*(u(2)+h/2*K1)-u(1)-h/2+Fm+Fah*cos(w*tao);
    K3=-2*ep*(u(2)+h/2*K2)-u(1)-h/2+Fm+Fah*cos(w*tao);
    K4=-2*ep*(u(2)+h*K3)-u(1)-h+Fm+Fah*cos(w*tao);
    u(2)=u(2)+h/6*(K1+2*K2+2*K3+K4);
    end

    结果图如下所示

    值得特别注意的是

    u=RK(u,tao,h,ep,w,Fm,Fah);

    输入的u与输出的u一定要符号一致,从而确保下一次迭代的初始值是上一次的值。确保迭代能一直进行下去。

    三、辅助验证

    接下来用MATLAB自带的ode45函数来进行验证。

    之前已经写过ode45函数的用法,在此不再进行介绍。

    源程序如下:

    t0=0:0.001:400;
    w=5;
    ep=0.02;
    Fm=0.1;
    Fah=0.05;
    y0=[0 0];
    [t,u]=ode45(@(t,u) odefun(t,u,w,ep,Fm,Fah),t0,y0);
    figure(1)
    subplot(2,1,1)
    plot(t,u(:,1))
    xlabel('Time')
    ylabel('displacement')
    subplot(2,1,2)
    plot(t,u(:,2))
    xlabel('Time')
    ylabel('velocity')
    function du=odefun(t,u,w,ep,Fm,Fah)
    du=[u(2);
        -2*ep*u(2)-u(1)+Fm+Fah*cos(w*t)]; 
    end

    运算结果图如下所示

    两中方法计算的结果是一样的。

    创作不易,如有帮助,记得点个赞呐。

    展开全文
  • matlab 多次求解偏微分方程 ode45

    千次阅读 2019-04-27 20:39:47
    师兄和我讨论了一个问题,就是在matlab中求解偏微分方程, 其中,偏微分方程中有的常数是一直变化的,要求很多次,而不是一个固定的常数求一次就行了。 其中,A1和A2是要求解的因变量,x是自变量,其他为...
  • 用MATLAB按二阶龙格库塔法求解微分方程组,急用于毕业设计!问题补充:今天才发现自己之前做的一点都不对,17号就交论文了,我傻了,急死了!求各位大侠帮帮忙。谢谢!要求解的微分方程如图所示。初始条件为t=0,x=y=...
  • 第四十九篇 二阶微分方程的打靶法 边界值问题 当我们试图用自变量的不同值所提供的信息来解二阶或更高阶的常微分方程时,我们必须使用与前面描述的不同的数值方法。 前面提到的初值问题经常涉及到“时间”作为自...
  • 实际上不管是多少个自由度的减摇,其实都需要解运动方程,当初做单自由度横摇的时候就是用四阶龙格库塔求解二阶微分方程,具体代码见: https://blog.csdn.net/Mezikov/article/details/103898387 0 写在前面 代码...
  • 记录一下MATLAB中ode45函数求解非刚性微分方程
  • 这是一个如何将二阶微分方程改写成两个一阶微分方程以便利用MATLAB的诸如ODE45等求解器求解的例子。下面的方程组包含了一个一阶与一个二阶微分方程: x'= - y*exp(-t/5) y' * exp(-t/5) 1; (1) y''= -2*...
  • Matlab学习——求解微分方程(组)

    千次阅读 2021-04-20 05:05:43
    Matlab学习——求解微分方程(组)发布时间:2018-05-29 17:18,浏览次数:738, 标签:Matlab介绍:1.在 Matlab 中,用大写字母 D 表示导数,Dy 表示 y 关于自变量的一阶导数,D2y 表示 y 关于自变量的二阶导数,...
  • ode45求解 clc clear 问题 % u'=-3u+6x+5 % u(0)=3 % 解析解:u=2e^(-3x)+2x+1 欧拉法 h=0.01;%步长 x=0:h:1; u=zeros(length(x),1); u(1)=3; fori=1:length(x)-1 du=-3*u(i)+6*x(i)+5; u(i+1)=u(i)+h*du; ...
  • nbspmatlabMatlab解微分方程(ODE+PDE).pdf33页本文档一共被下载:次,您可全文免费在线阅读后下载本文档。 下载提示1.本站不保证该用户上传的文档完整性,不预览、不比对内容而直接下载产生的反悔问题本站不予受理...
  • 在python里如何用4阶龙格库塔法解微分方程组呢?(在matlab里是写一个函数文件,然后ode45调用并赋初值) 大佬可否根据实例在python写一段程序我参考学习一下 (最后还需要绘制曲线t...
  • TECHNOLOGY INFORMATION 学 术 论 坛 传统的解微分方程组的方法有近似分析解法﹑表解法和图解法。这些方法有一定的局限性。MATLAB 是一种基于矩阵的数学软件包, 该软件包包括了一个数值程序扩展库, 并且有高级编程...
  • 第四十五篇 常微分方程编程基础 在工程分析的许多问题中,物理定律是用变量的导数表达的而不是变量本身来表示的,因此需要解决微分方程微分方程的解本质上涉及积分,有时可以用解析法得到。例如,最简单的微分...
  • Maple笔记2--常微分方程求解

    千次阅读 2021-02-07 02:24:03
    来源:网络论坛转载(VB资料库)常微分方程求解微分方程求解是数学研究与应用的一个重点和难点. Maple能够显式或隐式地解析地求解许多微分方程求解.在常微分方程求解器dsolve中使用了一些传统的技术例如laplace变换和...
  • 文章目录常微分方程求解过程中ode45与C语言自编程计算结果比较前言一、四阶Runge-Kutta算法原理二、C语言代码三、实例比较实例1实例2总结 前言 在求解常微分方程时,MATLAB中的ode求解器是必不可少的利器,几乎所有...
  • 数学建模方法:微分方程预测学习笔记
  • matlab微分方程代码概率Nordsieck方法 关于 该软件版本补充了: Michael Schober, Simo Särkkä, Philipp Hennig: "A probabilistic model for the numerical solution of initial value problems", 2017. 在Matlab...
  • 微分方程基础

    千次阅读 2021-05-25 09:31:28
    这与之前描述的常微分方程(ODE)相反,后者只涉及一个自变量。 工程和科学中的许多现象都是用偏微分方程描述的。例如,一个因变量,如压力或温度,可以作为时间(t)和空间(x, y, z)的函数变化。 求解偏微分方程的两种...
  • 欧拉法(Euler)求解常微分方程的Matlab程序及案例

    万次阅读 多人点赞 2021-02-27 23:07:56
    概念1) 常微分方程2) 一阶常微分方程3) 微分方程的解析解法和数值解法2. 算法2.1 显式欧拉法2.2 隐式欧拉法2.3 两步欧拉法2.4 改进欧拉法2.5 四种欧拉方法的对比3. 程序4. 案例5. 联系作者 1. 概念 1) 常微分方程 ...
  • 3、 数学建模实例 微分方程Matlab求解 求微分方程的数值解 (一)常微分方程数值解的定义 (二)建立数值解法的一些途径 (三)用Matlab软件求常微分方程的数值解 返 回 1、目标跟踪问题一:导弹追踪问题 2、目标跟踪问题...
  • 第四十七篇 线性常微分方程的theta-法 线性方程的θ-方法 迄今为止所描述的初值问题数值解的方法,只要它们是一阶导的,并以标准形式排列,适用于线性或非线性方程组。 如果微分方程是线性的,则可能采用的另一种...
  • 实例简介】和Matlab应用有关的,具体介绍常微分方程的使用和解法,原理性介绍,帮助理解。局部截断误差指的是,按()式计算由到这一步的计算值与精确值之差+。为了估计它,由展开得到的精确值是()、()两式相减(注意到...
  • SciML可以用AI求解偏微分方程,可以说十分黑科技了。 官网提供了notebook教程,可通过如下方式开启, >] pkg> add https://github.com/SciML/SciMLTutorials.jl # 按退格键退出pkg julia> using ...
  • 微分方程组的龙格-库塔公式求解 matlab 版 南京大学 王寻 1. 一阶常微分方程组 考虑方程组          00 00 zxz,z , y, xg' z yxy,z , y, xf' y 其经典四阶龙格-库塔格式如下: 对于 n...
  • ode45(ode45用法举例)

    千次阅读 2021-05-08 04:06:20
    ode45是用4阶方法提供候选解,5阶方法控制误差,是一种自适应步长的方法。而我们平时用的4阶和5阶龙格库塔法的公式中步长是给定的。具体算法和原理你可以看.ode45的初始条件是否必须是在x=0处没有必要 只要是选取的...
  • 用MATLAB方法进行变参量微分方程处理1 变参数微分方程数值求解例子2 求function dydt=fun(t,y,u,v)r=u+2;s=v-2;dydt=[r+y(2); s*y(1)-2*s*y(2)];u=[1;5;15;20;25];v=[6;12;18;24;30];tspan=0:1:4;y0=[0 2];yy=y0;for...
  • 参考:《高等应用数学问题的MATLAB求解(第四版)》 本文涵盖常微分方程的解析法、数值法(Euler、Runge-Kutta法等)

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 382
精华内容 152
关键字:

二阶微分方程ode45实例