精华内容
下载资源
问答
  • matlab求解微分方程的解析
    2022-03-21 10:08:38
    简 介:本文将研究微分方程的解析解算法,介绍在MATLAB 环境中如何用微分方程求解函数直接得出线性微分方程组的解析解,并对一阶简单的非线性微分方程的解析解求解进行探讨,从而得出结论,一般非线性微分方程是没有解析解的。

    关键词 微分方程解析解MATLAB

    §01


    假设已知常系数线性微分方程的一般描述方法为
    d n y ( t ) d t n + a 1 d n − 1 y ( t ) d t n − 1 + a 2 d n − 2 y ( t ) d t n − 2 + ⋯ + a n − 1 d y ( t ) d t + a n y ( t )   = b 1 d m u ( t ) d t m + b 2 d m − 1 u ( t ) d t m − 1 + ⋯ + b m d u ( t ) d t + b m + 1 u ( t ) \frac{\mathrm{d}^{n} y(t)}{\mathrm{d} t^{n}}+a_{1} \frac{\mathrm{d}^{n-1} y(t)}{\mathrm{d} t^{n-1}}+a_{2} \frac{\mathrm{d}^{n-2} y(t)}{\mathrm{d} t^{n-2}}+\cdots+a_{n-1} \frac{\mathrm{d} y(t)}{\mathrm{d} t}+a_{n} y(t) \\ \ \\ =b_{1} \frac{\mathrm{d}^{m} u(t)}{\mathrm{d} t^{m}}+b_{2} \frac{\mathrm{d}^{m-1} u(t)}{\mathrm{d} t^{m-1}}+\cdots+b_{m} \frac{\mathrm{d} u(t)}{\mathrm{d} t}+b_{m+1} u(t) dtndny(t)+a1dtn1dn1y(t)+a2dtn2dn2y(t)++an1dtdy(t)+any(t) =b1dtmdmu(t)+b2dtm1dm1u(t)++bmdtdu(t)+bm+1u(t)

    其中, a i a_i ai, b i b_i bi均为常数,对零初值问题有 L [ d m y ( t ) / d t m ] = s m L [ y ( t ) ] \mathscr{L}\left[\mathrm{d}^{m} y(t) / \mathrm{d} t^{m}\right]=s^{m} \mathscr{L}[y(t)] L[dmy(t)/dtm]=smL[y(t)],可以对应得出下面的多项式代数方程

    s n + a 1 s n − 1 + a 2 s n − 2 + ⋯ + a n − 1 s + a n = 0 s^{n}+a_{1} s^{n-1}+a_{2} s^{n-2}+\cdots+a_{n-1} s+a_{n}=0 sn+a1sn1+a2sn2++an1s+an=0

    假设代数方程的特征根 s i s_i si均可以求出,且假设它们均相异,则可以得出原微分方程的解析解一般形式为

    y ( t ) = C 1 e r 1 t + C 2 e r 2 t + ⋯ + C n e r n t + γ ( t ) y(t)=C_{1} \mathrm{e}^{r_{1} t}+C_{2} \mathrm{e}^{r_{2} t}+\cdots+C_{n} \mathrm{e}^{r_{n} t}+\gamma(t) y(t)=C1er1t+C2er2t++Cnernt+γ(t)

    其中, C i C_i Ci为待定系数,而 γ ( t ) \gamma(t) γ(t)是满足 u ( t ) u(t) u(t)输入的一个特解。 s i s_i si有重根的情况也有相应的解析解形式。

    从得出的代数方程看,由著名的Abel-Ruffini定理可知,4次及以下的多项式代数方程是能求出根的解析解的,故可以得出结论,低阶线性微分方程有一般意义下的解析解,结合多项式方程的数值解法可以得出一般高次多项式代数方程的数值解法,即高阶线性微分方程的准解析解方法。本文将介绍用 MATLAB 语言及其符号运算工具箱求解线性常系数微分方程解析解的方法。

    §02 数调用格式


    y = dsolve(fun1, fun2, ... , funm)  % 默认自变量t
    y = dsolve(fun1, fun2, ... , funm, 'x')   % 指明自变量x
    
    • 可以同时求解多个方程、已知条件
    • 注意自变量设置,否则可能得出无用的结果

    §03 用举例


    例1:微分方程解析解

    题目: 假设输入信号为 u ( t ) = e − 5 t cos ⁡ ( 2 t + 1 ) + 5 u(t)=\mathrm{e}^{-5 t} \cos (2 t+1)+5 u(t)=e5tcos(2t+1)+5,试求下面微分方程的通解:

    y ( 4 ) ( t ) + 10 y ( 3 ) ( t ) + 35 y ¨ ( t ) + 50 y ˙ ( t ) + 24 y ( t ) = 5 u ¨ ( t ) + 4 u ˙ ( t ) + 2 u ( t ) y^{(4)}(t)+10 y^{(3)}(t)+35 \ddot{y}(t)+50 \dot{y}(t)+24 y(t)=5 \ddot{u}(t)+4 \dot{u}(t)+2 u(t) y(4)(t)+10y(3)(t)+35y¨(t)+50y˙(t)+24y(t)=5u¨(t)+4u˙(t)+2u(t)

    求解:

    syms y(t);
    u = exp(-5*t)*cos(2*t+1)+5;
    uu = 5*diff(u,t,2) + 4*diff(u,t) + 2*u;
    f(t) = diff(y,t,4) + 10*diff(y,t,3) + 35*diff(y,t,2) + 50*diff(y,t) + 24*y;
    ySol = dsolve(f(t)==uu);
    ySol = simplify(ySol)
    

    运行程序可得:

    sy =
     
    C1*exp(-4*t) - (547*exp(-5*t)*sin(2*t + 1))/520 - (343*exp(-5*t)*cos(2*t + 1))/520 + C2*exp(-3*t) + C3*exp(-2*t) + C4*exp(-t) + 5/12
    

    化简可得 解的数学形式:

    y ( t ) = 5 12 − 343 520 e − 5 t cos ⁡ ( 2 t + 1 ) − 547 520 e − 5 t sin ⁡ ( 2 t + 1 ) + C 1 e − 4 t + C 2 e − 3 t + C 3 e − 2 t + C 4 e − t y(t)=\frac{5}{12}-\frac{343}{520} \mathrm{e}^{-5 t} \cos (2 t+1)-\frac{547}{520} \mathrm{e}^{-5 t} \sin (2 t+1)+C_{1} \mathrm{e}^{-4 t}+C_{2} \mathrm{e}^{-3 t}+C_{3} \mathrm{e}^{-2 t}+C_{4} \mathrm{e}^{-t} y(t)=125520343e5tcos(2t+1)520547e5tsin(2t+1)+C1e4t+C2e3t+C3e2t+C4et

    解析解的检验:

    diff(ySol,4) + 10*diff(ySol,3) + 35*diff(ySol,2) + 50*diff(ySol) + 24*ySol - uu;
    simplify(ans)
    

    运行程序可得:

    ans =
     
    0
    

    已知初值条件:
    y ( 0 ) = 3 ,    y ˙ ( 0 ) = 2 ,    y ¨ ( 0 ) = y ( 3 ) ( 0 ) = 0 y(0)=3, ~~\dot{y}(0)=2, ~~\ddot{y}(0)=y^{(3)}(0)=0 y(0)=3,  y˙(0)=2,  y¨(0)=y(3)(0)=0

    syms y(t);
    u = exp(-5*t)*cos(2*t+1)+5;
    uu = 5*diff(u,t,2) + 4*diff(u,t) + 2*u;
    f(t) = diff(y,t,4) + 10*diff(y,t,3) + 35*diff(y,t,2) + 50*diff(y,t) + 24*y;
    Dy = diff(y,t); D2y = diff(y,t,2); D3y = diff(y,t,3);
    cond = [y(0)==3, Dy(0)==2, D2y(0)==0, D3y(0)==0];
    ySol= dsolve(f(t)==uu, cond);
    ySol = simplify(ySol) 
    ezplot(ySol,[0,5])
    
    ▲ 图1 已知初值条件微分方程解的图像

    复杂边界条件:

    y ( 0 ) = 1 / 2 ,    y ˙ ( π ) = 1 ,    y ¨ ( 2 π ) = 0 ,    y ˙ ( 2 π ) = 1 / 5 y(0)=1 / 2, ~~\dot{y}(\pi)=1, ~~\ddot{y}(2 \pi)=0, ~~\dot{y}(2 \pi)=1 / 5 y(0)=1/2,  y˙(π)=1,  y¨(2π)=0,  y˙(2π)=1/5

    syms y(t);
    u = exp(-5*t)*cos(2*t+1)+5;
    uu = 5*diff(u,t,2) + 4*diff(u,t) + 2*u;
    f(t) = diff(y,t,4) + 10*diff(y,t,3) + 35*diff(y,t,2) + 50*diff(y,t) + 24*y;
    Dy = diff(y,t); D2y = diff(y,t,2);
    cond = [y(0)==1/2, Dy(pi)==1, D2y(2*pi)==0, Dy(2*pi)==1/5];
    ySol= dsolve(f(t)==uu, cond);
    ySol = vpa(ySol);
    ySol = simplify(ySol) 
    ezplot(ySol,[0,5])
    
    ▲ 图2 复杂边界条件微分方程解的图像

    例2:求解微分方程组

    题目: 试求解线性微分方程组的解析解

    { x ′ ′ ( t ) + 2 x ′ ( t ) = x ( t ) + 2 y ( t ) − e − t y ′ ( t ) = 4 x ( t ) + 3 y ( t ) + 4 e − t \left\{\begin{array}{l}x^{\prime \prime}(t)+2 x^{\prime}(t)=x(t)+2 y(t)-\mathrm{e}^{-t}\\ \\ y^{\prime}(t)=4 x(t)+3 y(t)+4 \mathrm{e}^{-t}\end{array}\right. x(t)+2x(t)=x(t)+2y(t)ety(t)=4x(t)+3y(t)+4et

    syms x(t) y(t)
    eqn1 = diff(x,t,2) + 2*diff(x,t) == x + 2*y(t) - exp(-t);
    eqn2 = diff(y,t) == 4*x + 3*y(t) + 4*exp(-t);
    [xSol, ySol] = dsolve([eqn1, eqn2])
    

    例3:变系数微分方程

    题目: 变系数微分方程

    ( 2 x + 3 ) 3 y ′ ′ ′ + 3 ( 2 x + 3 ) y ′ − 6 y = 0 (2 x+3)^{3} y^{\prime \prime \prime}+3(2 x+3) y^{\prime}-6 y=0 (2x+3)3y+3(2x+3)y6y=0

    syms y(x)
    eqn = (2*x+3)^3*diff(y,x,3) + 3*(2*x+3)*diff(y,x) - 6*y == 0;
    ySol = dsolve(eqn)
    

    例4:多维微分方程组

    题目:

    { x ′ ′ − x + y + z = 0 x + y ′ ′ − y + z = 0 x + y + z ′ ′ − z = 0 \left\{\begin{array}{l}x^{\prime \prime}-x+y+z=0 \\ \\x+y^{\prime \prime}-y+z=0 \\ \\ x+y+z^{\prime \prime}-z=0\end{array}\right. xx+y+z=0x+yy+z=0x+y+zz=0

    初值: x ( 0 ) = 1 ,    y ( 0 ) = z ( 0 ) = x ′ ( 0 ) = y ′ ( 0 ) = z ′ ( 0 ) = 0 x(0)=1,~~ y(0)=z(0)=x^{\prime}(0)=y^{\prime}(0)=z^{\prime}(0)=0 x(0)=1,  y(0)=z(0)=x(0)=y(0)=z(0)=0

    syms x(t) y(t) z(t)
    eqn1 = diff(x,t,2) - x + y + z == 0;
    eqn2 = x + diff(y,t,2) - y + z == 0;
    eqn3 = x + y + diff(z,t,2) - z == 0;
    Dx = diff(x, t); Dy = diff(y, t); Dz = diff(z, t); 
    cond = [x(0)==1, y(0)==0, z(0)==0, Dx(0)==0, Dy(0)==0, Dz(0)==0];
    [xSol, ySol, zSol] = dsolve([eqn1, eqn2, eqn3], cond)
    
    更多相关内容
  • MATLAB求解微分方程微分方程组方法介绍和例子。Matlab
  • MATLAB算法-求解微分方程数值和解析.ppt
  • matlab求解微分方程组代码 计算方法(calculation methodb) 该项目是《计算方法》一书中提到的经典方法和算法的matlab程序实现,包含代码详解和运行过程。 :grinning_face_with_big_eyes: 1.简介 2.线性方程组的...
  • TMU_BME_2013 Topic: 如何使用 MATLAB微分方程组 a.What ? 微分方程 指描述未知函数的导数与自变 量之间的关系的方程未知函数是一元函 数的微分方程称作 常微分方程 未知函数 是多元函数的微分方程称作 偏...
  • Matlab 求解微分方程(ODE)
  • 简 介:前面介绍了微分方程的解析方法,同时也指出很多非线性微分方程是不存在解析解法的,需要使用数值解法对之进行研究。本文着重讨论基于 MATLAB/Simulink语言的各类微分方程的数值方法。关键词:微分方程、...
    简 介:前面介绍了微分方程的解析解方法,同时也指出很多非线性微分方程是不存在解析解法的,需要使用数值解法对之进行研究。本文着重讨论基于 MATLAB/Simulink语言的各类微分方程的数值解方法。

    关键词 微分方程、数值解、MATLAB

    §01


    一般微分方程的数值解法很大一类是关于微分方程初值问题的数值解法,这类问题需要用一阶显式的微分方程组来描述为

    x ˙ ( t ) = f ( t , x ( t ) ) \dot{\boldsymbol{x}}(t)=\boldsymbol{f}(t, \boldsymbol{x}(t)) x˙(t)=f(t,x(t))

    其中, x T ( t ) = [ x 1 ( t ) , x 2 ( t ) , ⋯   , x n ( t ) ] \boldsymbol{x}^{\mathrm{T}}(t)=\left[x_{1}(t), x_{2}(t), \cdots, x_{n}(t)\right] xT(t)=[x1(t),x2(t),,xn(t)]称为状态向量, f T ( ⋅ ) = [ f 1 ( ⋅ ) , f 2 ( ⋅ ) , ⋯   , f n ( ⋅ ) ] \boldsymbol{f}^{\mathrm{T}}(\cdot)=\left[f_{1}(\cdot), f_{2}(\cdot), \cdots, f_{n}(\cdot)\right] fT()=[f1(),f2(),,fn()]可以是任意非线性函数。所谓初值问题是指,若已知初始状态 x 0 = [ x 1 ( 0 ) , ⋯   , x n ( 0 ) ] T \boldsymbol{x}_{0}=\left[x_{1}(0), \cdots, x_{n}(0)\right]^{\mathrm{T}} x0=[x1(0),,xn(0)]T,用数值求解方法求出在某个时间区间 t ∈ [ 0 , t f ] t \in\left[0, t_{\mathrm{f}}\right] t[0,tf]内各个时刻状态变量 x ( t ) \boldsymbol{x}(t) x(t) 的数值,这里 t f t_{\mathrm{f}} tf 又称为终止时间。

    其他类型微分方程求解必须先转换:

    1. 单个高阶常微分方程处理方法

    一个高阶常微分方程的一般形式为:

    y ( n ) = f ( t , y , y ′ , ⋯   , y ( n − 1 ) ) y^{(n)}=f\left(t, y, y^{\prime}, \cdots, y^{(n-1)}\right) y(n)=f(t,y,y,,y(n1))

    输出变量的各阶导数初始值为:

    y ( 0 ) ,    y ′ ( 0 ) ,    ⋯   ,    y ( n − 1 ) ( 0 ) y(0), ~~y^{\prime}(0), ~~\cdots,~~ y^{(n-1)}(0) y(0),  y(0),  ,  y(n1)(0)

    选择一组状态变量:

    x 1 = y ,    x 2 = y ′ ,    ⋯   ,    x n = y ( n − 1 ) x_{1}=y, ~~x_{2}=y^{\prime}, ~~\cdots,~~ x_{n}=y^{(n-1)} x1=y,  x2=y,  ,  xn=y(n1)

    原高阶常微分方程模型变换为一阶微分方程组:

    { x 1 ′ = x 2 x 2 ′ = x 3 ⋮ x n ′ = f ( t , x 1 , x 2 , ⋯   , x n ) \left\{\begin{aligned} x_{1}^{\prime} &=x_{2} \\ x_{2}^{\prime} &=x_{3} \\ & \vdots \\ x_{n}^{\prime} &=f\left(t, x_{1}, x_{2}, \cdots, x_{n}\right) \end{aligned}\right. x1x2xn=x2=x3=f(t,x1,x2,,xn)

    初值为:

    x 1 ( 0 ) = y ( 0 ) ,    x 2 ( 0 ) = y ′ ( 0 ) ,    ⋯   ,    x n ( 0 ) = y ( n − 1 ) ( 0 ) x_{1}(0)=y(0), ~~x_{2}(0)=y^{\prime}(0), ~~\cdots, ~~x_{n}(0)=y^{(n-1)}(0) x1(0)=y(0),  x2(0)=y(0),  ,  xn(0)=y(n1)(0)

    2. 高阶常微分方程组的变换方法

    多元高阶常微分方程组的处理

    { x ( m ) = f ( t , x , x ′ , ⋯   , x ( m − 1 ) , y , ⋯   , y ( n − 1 ) ) y ( n ) = g ( t , x , x ′ , ⋯   , x ( m − 1 ) , y , ⋯   , y ( n − 1 ) ) \left\{\begin{array}{l}x^{(m)}=f\left(t, x, x^{\prime}, \cdots, x^{(m-1)}, y, \cdots, y^{(n-1)}\right) \\ \\ y^{(n)}=g\left(t, x, x^{\prime}, \cdots, x^{(m-1)}, y, \cdots, y^{(n-1)}\right)\end{array}\right. x(m)=f(t,x,x,,x(m1),y,,y(n1))y(n)=g(t,x,x,,x(m1),y,,y(n1))

    状态变量的选择不唯一,建议:选择如下状态变量

    x 1 = x ,   x 2 = x ′ ,   ⋯   ,   x m = x ( m − 1 ) ,   x m + 1 = y ,   x m + 2 = y ′ ,   ⋯   ,   x m + n = y ( n − 1 ) x_{1}=x, ~x_{2}=x^{\prime}, ~\cdots, ~x_{m}=x^{(m-1)}, \\ \ \\ x_{m+1}=y, ~x_{m+2}=y^{\prime},~ \cdots, ~x_{m+n}=y^{(n-1)} x1=x, x2=x, , xm=x(m1), xm+1=y, xm+2=y, , xm+n=y(n1)

    新的状态方程

    { x 1 ′ = x 2 ⋮ x m ′ = f ( t , x 1 , x 2 , ⋯   , x m + n ) x m + 1 ′ = x m + 2 ⋮ x m + n ′ = g ( t , x 1 , x 2 , ⋯   , x m + n ) \left\{\begin{aligned} x_{1}^{\prime}=& x_{2} \\ & \vdots \\ x_{m}^{\prime}=& f\left(t, x_{1}, x_{2}, \cdots, x_{m+n}\right) \\ x_{m+1}^{\prime} &=x_{m+2} \\ & \vdots \\ x_{m+n}^{\prime} &=g\left(t, x_{1}, x_{2}, \cdots, x_{m+n}\right) \end{aligned}\right. x1=xm=xm+1xm+nx2f(t,x1,x2,,xm+n)=xm+2=g(t,x1,x2,,xm+n)

    可以描述该方程,然后用 ode45 等求解。

    §02 分方程求解的误差与步长问题


    1. 不能无限制地减小步长 h h h的值的两条原因

    • 减慢计算速度
    • 增加累积误差

    2. 在对微分方程求解过程中应采取的三个措施

    • 选择适当的步长
    • 改进近似算法精度
    • 采用变步长方法

    §03 数调用格式(ode45)


    [t, x] = ode45(fun,[t0, tf], x0)  % 直接求解
    [t, x] = ode45(fun,[t0, tf], x0, options)  % 带有控制选项
    [t, x] = ode45(fun,[t0, tf], x0, options, p1, p2, ...)  % 带有附加参数
    

    当然,还存在其他求解函数,如ode15sode23等,不同的算法(函数)适合于不同类型的微分方程。

    §04 分方程的MATLAB描述


    描述需要求解的微分方程组

       f u n c t i o n     x d = f u n ( t , x ) \rm{function}~~~ \boldsymbol{x}_{d}= fun (t, \boldsymbol{x}) function   xd=fun(t,x)

       f u n c t i o n     x d = f u n ( t , x , p 1 , p 2 , ⋯   ) \rm{function} ~~~ \boldsymbol{x}_{d}= fun \left(t, \mathrm{x}, p_{1}, p_{2}, \cdots\right) function   xd=fun(t,x,p1,p2,)

    修改控制变量

    options = odeset('RelTol', 1e-7) ;
    options = odeset;  options.RelTol = 1e-7;
    

    §05 用举例:Lorenz方程


    例1:无参数

    题目: 求解下列Lorenz模型

    { x ˙ 1 ( t ) = − β x 1 ( t ) + x 2 ( t ) x 3 ( t ) x ˙ 2 ( t ) = − ρ x 2 ( t ) + ρ x 3 ( t ) x ˙ 3 ( t ) = − x 1 ( t ) x 2 ( t ) + σ x 2 ( t ) − x 3 ( t ) \left\{\begin{array}{l}\dot{x}_{1}(t)=-\beta x_{1}(t)+x_{2}(t) x_{3}(t) \\ \\ \dot{x}_{2}(t)=-\rho x_{2}(t)+\rho x_{3}(t) \\ \\ \dot{x}_{3}(t)=-x_{1}(t) x_{2}(t)+\sigma x_{2}(t)-x_{3}(t)\end{array}\right. x˙1(t)=βx1(t)+x2(t)x3(t)x˙2(t)=ρx2(t)+ρx3(t)x˙3(t)=x1(t)x2(t)+σx2(t)x3(t)

    式中参数为 β = 8 / 3 ,   ρ = 10 ,   σ = 28 \beta=8 / 3, ~\rho=10, ~\sigma=28 β=8/3, ρ=10, σ=28

    初始条件为 x 1 ( 0 ) = x 2 ( 0 ) = 0 ,    x 3 ( 0 ) = ϵ ,    i . e . ,   ϵ = 1 0 − 10 x_{1}(0)=x_{2}(0)=0, ~~x_{3}(0)=\epsilon,~~ \rm{i.e.}, ~\epsilon=10^{-10} x1(0)=x2(0)=0,  x3(0)=ϵ,  i.e., ϵ=1010

    求解:

    1. 写出标准型

    x ′ ( t ) = [ − β x 1 ( t ) + x 2 ( t ) x 3 ( t ) − ρ x 2 ( t ) + ρ x 3 ( t ) − x 1 ( t ) x 2 ( t ) + σ x 2 ( t ) − x 3 ( t ) ] , x ( 0 ) = [ 0 0 ϵ ] \boldsymbol{x}^{\prime}(t)=\left[\begin{array}{c}-\beta x_{1}(t)+x_{2}(t) x_{3}(t) \\ \\ -\rho x_{2}(t)+\rho x_{3}(t) \\ \\ -x_{1}(t) x_{2}(t)+\sigma x_{2}(t)-x_{3}(t)\end{array}\right], \quad \boldsymbol{x}(0)=\left[\begin{array}{l}0 \\ \\ 0 \\ \\ \epsilon\end{array}\right] x(t)=βx1(t)+x2(t)x3(t)ρx2(t)+ρx3(t)x1(t)x2(t)+σx2(t)x3(t),x(0)=00ϵ

    2. 微分方程的MATLAB描述

    • M-文件描述
    function y = lorenzeq(t, x)
    	y = [-8/3*x(1)+x(2)*x(3); -10*x(2)+10*x(3); -x(1)*x(2)+28*x(2)-x(3)];
    
    • 匿名函数
    f = @(t,x)[-8/3*x(1)+x(2)*x(3); -10*x(2)+10*x(3); -x(1)*x(2)+28*x(2)-x(3)];
    

    3. 微分方程求解

    • 匿名函数
    t_final = 100; 
    x0 = [0; 0; 1e-10]; 
    f = @(t,x)[-8/3*x(1)+x(2)*x(3); -10*x(2)+10*x(3); -x(1)*x(2)+28*x(2)-x(3)];
    [t,x] = ode45(f, [0,t_final], x0); 
    plot(t,x), figure; plot3(x(:,1),x(:,2),x(:,3));
    
    • M-文件描述
    t_final = 100; 
    x0 = [0; 0; 1e-10]; 
    f = @(t,x)[-8/3*x(1)+x(2)*x(3); -10*x(2)+10*x(3); -x(1)*x(2)+28*x(2)-x(3)];
    [t,x] = ode45(@lorenzeq, [0,t_final], x0); 
    plot(t,x), figure; plot3(x(:,1),x(:,2),x(:,3));
    
    function y = lorenzeq(t, x)
    	y = [-8/3*x(1)+x(2)*x(3); -10*x(2)+10*x(3); -x(1)*x(2)+28*x(2)-x(3)];
    end
    
    ▲ 图 1

    ▲ 图 2

    例2:带有附加参数

    引入附加参数的目的: 如果参数 β \beta β, ρ \rho ρ, σ \sigma σ 改变,不需要修改原函数。

    重新求解Lorenz方程

    方式一

    f = @(t,x,beta,rho,sigma)[-beta*x(1) + x(2)*x(3); -rho*x(2) + rho*x(3); -x(1)*x(2) + sigma*x(2) - x(3)];
    t_final=100; 
    x0 = [0; 0; 1e-10]; 
    b1 = 8/3; r1 = 10; s1 = 28; 
    [t,x]=ode45(f, [0,t_final], x0, [], b1, r1, s1);
    

    方式二

    beta = 2; rho = 5; sigma = 20; 
    f = @(t,x)[-beta*x(1) + x(2)*x(3); -rho*x(2) + rho*x(3); -x(1)*x(2) + sigma*x(2) - x(3)];
    t_final = 100;
    x0 = [0; 0; 1e-10]; 
    [t2,x2] = ode45(f, [0,t_final], x0);
    
    展开全文
  • MATLAB求解微分方程

    千次阅读 2021-05-23 20:48:31
    然而能够求得解析微分方程十分有限,绝大多数微分方程需要利用数值方法来近似求解。下面介绍如何用 Matlab 来计算微分方程(组)的数值。 Euler折线法 微分方程的基本数值解法——Euler折线法。 步骤: 分割...

    求微分方程的解

    自牛顿发明微积分以来,微分方程在描述事物运动规律上已发挥了重要的作用。实际应用问题通过数学建模所得到的方程,绝大多数是微分方程。由于实际应用的需要,人们必须求解微分方程。然而能够求得解析解的微分方程十分有限,绝大多数微分方程需要利用数值方法来近似求解。下面介绍如何用 Matlab 来计算微分方程(组)的数值解。

    Euler折线法

    微分方程的基本数值解法——Euler折线法。

    步骤:
    分割求解区间,差商代替微商,解代数方程

    例子:
    在这里插入图片描述解:

    等距剖分: a = x 0 < x 1 < x 2 < ⋯ < x n − 1 < x n = b a=x_0<x_1<x_2< \cdots <x_{n-1}<x_n=b a=x0<x1<x2<<xn1<xn=b

    步长: h = x k + 1 − x k = ( b − a ) / n , k = 0 , 1 , 2 , . . . , n − 1 h=x_{k+1}-x_k=(b-a)/n,k=0,1,2,...,n-1 h=xk+1xk=(ba)/n,k=0,1,2,...,n1

    y ( x k + 1 ) ≈ y ( x k ) + h y ′ ( x k ) y(x_{k+1})\approx y(x_k)+hy\prime(x_k) y(xk+1)y(xk)+hy(xk)

    取步长 h = (2 - 0)/n = 2/n,得差分方程
    在这里插入图片描述
    当 h=0.4,即 n=5 时,代码如下:

    clear
    a=0;  b=2;
    h=0.4;
    n=(b-a)/h+1; 
    x(1)=0;  y(1)=1;
    for i=1:n-1        
        y(i+1)=y(i)+h*(y(i)+2*x(i)/y(i)^2);
        x(i+1)=x(i)+h;
    end
    y1=((5/3)*exp(3*x)-2*x-2/3).^(1/3);
    plot(x,y1,x,y);
    h = legend('解析解','Euler解');
    

    在这里插入图片描述

    解得:
    y = ( 5 3 e 3 x − 2 x − 2 3 ) 1 / 3 y=(\frac{5}{3}e^{3x}-2x-\frac{2}{3})^{1/3} y=(35e3x2x32)1/3

    Runge-Kutta方法

    龙格-库塔法(Runge-Kutta methods)是用于非线性常微分方程的解的重要的一类隐式或显式迭代法。

    同样是上面的例子:
    在这里插入图片描述
    在这里插入图片描述

    clear;
    a=0; b=2; h=0.4;
    n=(b-a)/h+1; 
    x(1)=0; y(1)=1;  
    for i=1:n-1
        L1=y(i)+2*x(i)/y(i)^2;
        t=x(i)+h/2;t1=y(i)+h*L1/2;
        L2=t1+2*t/t1^2;t1=y(i)+h*L2/2;
        L3=t1+2*t/t1^2;t1=y(i)+h*L3/2;
        L4=t1+2*(x(i)+h)/t1^2;
        y(i+1)=y(i)+h*(L1+2*L2+2*L3+L4)/6; 
        x(i+1)=x(i)+h;
    end
    plot(x,y)
    

    在这里插入图片描述

    MATLAB自带的ODE求解器

    ODE (常微分方程)

    没有一种算法可以有效地解决所有的 ODE 问题,因此MATLAB 提供了多种ODE求解器,对于不同的ODE,可以调用不同的求解器。
    在这里插入图片描述

    基本语法
    [T,Y] = solver(odefun,tspan,y0)

    其中 y0 为初值条件,tspan为求解区间;Matlab在数值求解时自动对求解区间进行分割,T (向量) 中返回的是分割点的值(自变量),Y (向量) 中返回的是解函数在这些分割点上的函数值。
    solver 为Matlab的ODE求解器(可以是 ode45、ode23、ode113、ode15s、ode23s、ode23t、ode23tb)

    展开全文
  • 一个微分方程可以通过多种方法求解。 我已经介绍了 simulink 方法来解决微分方程。 que 显示在屏幕截图中。 非常欢迎查询和评论。 :)
  • 本文是自己写的关于怎样利用MATLAB求解微分方程数值的,文中从Euler法讲起,最后总结了常用的odeXX的用法及其原理,其中包含各个函数怎样使用的MATLAB代码
  • matlab求解微分方程

    2022-05-09 16:54:12
    Matlab求常微分方程组的解析_y=520(2sinM-sin2M)的博客-CSDN博客_matlab求解微分方程的解析 Matlab求常微分方程组的解析_y=520(2sinM-sin2M)的博客-CSDN博客_matlab求解微分方程的解析 Matlab...
    展开全文
  • matlab解微分方程

    2018-09-17 21:51:49
    代码为2018年全国数学建模竞赛a题第一问matlab源程序,可动态生成三层隔热服距离与温度的关系图,以及三层隔热服的温度分布图。主要内容建立在一维非稳态热传导以及偏微分方程的解法程序
  • 微分方程在实际应用中十分广泛,涉及领域众多,但对于微分方程的数值的计算仍然有很大挑战。本文着重对微分方程数值解求解的常用的一类基础方法——欧拉法进行在MAILAB的应用下的一个简单介绍。
  • 结合MATLAB微分方程数值工具箱介绍偏微分方程求解,分GUI和MATLAB函数两种实现方式进行介绍。
  • MATLAB四阶龙格库塔法 求解微分方程数值 部分源码 clear;clc;close all h=0.2; t=0:h:3; x(1)=1; %使用Runge-Kutta方法,计算微分方程的数值
  • Matlab微分方程的数值解法常用程序-偏微分方程的数值解法_程序.rar 包括解决一些微分方程的常用程序,希望对大家有用,欢迎下载!!:)
  • 龙格库塔法是用于非线性常微分方程的重要的一类隐式或显式迭代法。
  • //使用命名空间XSLSF//数组xArray存放x的值;ti为当前有效值的个数;tmax为ti对应的时间;tmin为起始时间。xt(t:k:xArray,ti,tmax,tmin)={k=(t-tmin)/(tmax-tmin)*ti-1,if[k<0, k=0], if[k>=ti, k=ti-1],...
  • 1_基于MATLAB求解微分方程数值和解析.ppt数学建模
  • Matlab求解微分方程(组)及偏微分方程(组)
  • 利用MATLAB求解一阶微分方程,输入信号为单位阶跃脉冲信号,代码绘制了系统脉冲响应的计算结果。
  • 微分方程MATLAB求解

    2018-09-07 14:24:12
    微分方程
  • matlab解微分方程

    2021-05-24 20:08:37
    研究了使用matlab如何解微分方程
  • 如何用matlab求解非线性微分方程组(基于龙格库塔的数值微分算法)?.docx
  • matlab求解微分方程

    2022-04-07 22:30:25
    实例 微分方程: 初始条件: %构建 syms y(t) s= dsolve('-Dy=y +1','y(0) =1'); %作图 t = 0:0.1:5 y = eval(subs(s)) plot(t,y) 在命令窗口行键入结果 's',显示解析式:
  • MATLAB使用欧拉Euler法求解微分方程组 部分源码 clear;clc c=2/3; %设置c的值 x(1)=0.1; %设置x初值为0.1 y(1)=0.3; %设置y初值为0.3 h=0.05; %设置步长为0.05
  • Matlab 解微分方程;一微分方程的解析解;输入命令: y=dsolve'D2y+4*Dy+29*y=0'y(0)=0,Dy(0)=15'x;二微分方程的数值解;二建立数值解法的一些途径;2使用数值积分;3使用泰勒公式; 1在解n个未知函数的方程组时x0和x均...
  • MATLAB求解微分方程

    千次阅读 2020-02-09 15:48:28
    微分方程的解析微分方程(组)的解析命令: dsolve('方程1','方程2',...'方程n','初始条件','自变量') 注:字母D表示求微分,D2、D3等表示求高阶微分,D后所跟的字母为因变量,自变量可以指定或由系统规则...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 11,928
精华内容 4,771
关键字:

利用matlab求解微分方程的解

matlab 订阅