用MATLAB解常微分方程

1、实验四 求微分方程的解一、问题背景与实验目的实际应用问题通过数学建模所归纳而得到的方程 ， 绝大多数都是微分方 程 ，真正能得到代数方程的机会很少 另一方面 ， 能够求解的微分方程也是十 分有限的 ，特别是高阶方程和偏微分方程 （组） 这就要求我们必须研究微分 方程（组）的解法 ，既要研究微分方程 （组）的解析解法 （精确解 ），更要研 究微分方程 （组）的数值解法 （近似解 ）对微分方程 （组）的解析解法 （精确解 ），Matlab 有专门的函数可以用 ，本 实验将作一定的介绍 本实验将主要研究微分方程 （组）的数值解法 （近似解 ），重点介绍 Euler 折 线法二、相关函数 （命令 ）及

2、简介1 dsolve （equ1,equ2, ） ： Matlab 求微分方程的解析解 equ1 、 equ2 、为方程 （或条件 ）写方程 （或条件 ）时用 Dy 表示 y 关于自变量的 一阶导数 ，用用 D2y 表示 y 关于自变量的二阶导数 ，依此类推 2 simplify（s ）：对表达式 s 使用 maple 的化简规则进行化简 例如：syms xsimplify（sin（x）2 + cos（x）2）ans=13 r,how=simple（s） ：由于 Matlab 提供了多种化简规则 ， simple 命令专业 word 可编辑就是对表达式 s 用各种规则进行化简 ，然后用 r 返

3、回最简形式 ，how 返回形 成这种形式所用的规则 例如：syms xr,how=simple(cos(x)2-sin(x)2)r = cos(2*x)how = combine4 T,Y = solver( odefun,tspan,y 0) 求微分方程的数值解说明：ode45 、 ode23 、 ode113 、 ode15s 、ddyt f (t,y)y(t0) y0(1) 其 中 的 solver 为 命 令ode23s 、ode23t 、 ode23tb 之一 (2) odefun 是显式常微分方程 ：(3) 在积分区间 tspan=t0,tf上，从t0到tf ，用初始条件 y0求解

4、(4) 要获得问题 在其他 指定 时间 点 t0,t1,t2, 上的解，则 令 tspan = t0,t1,t2, ,tf(要求是单调的 )(5) 因为没有一种算法可以有效地解决所有的ODE 问题 ，为此 ， Matlab提供了多种求解器 Solver，对于不同的 ODE 问题，采用不同的 Solver 求解器ODE 类型特点说明专业 word 可编辑Solverode45非刚性单 步 算 法 ； 4 、 5 阶大部分场合的首选算Runge-Kutta 方程； 累计截断法误差达 ( x)3ode23非刚性单 步 算 法 ； 2 、 3 阶使用于精度较低的情Runge-Kutta 方程； 累计截

5、断形误差达 ( x)3ode113非刚性多步法 ；Adams 算法；高低精计 算 时 间 比 ode45度均可到 10 3 10 6短ode23t适度刚性采用梯形算法适度刚性情形ode15s刚性多 步 法 ； Gears 反 向 数 值 微若 ode45 失效 时 ，分；精度中等可尝试使用ode23s刚性单步法 ； 2 阶 Rosebrock 算当精度较低时 ， 计算法；低精度时间比 ode15s 短ode23t刚性梯形算法 ；低精度当精度较低时 ， 计算b时间比 ode15s 短(6) 要特别的是 ：ode23 、ode45 是极其常用的用来求解非刚性的标准形式 的一阶常微分方程 (组)的初

6、值问题的解的 Matlab 的常用程序 ，其中 ：ode23 采用龙格 -库塔 2 阶算法 ，用 3 阶公式作误差估计来调节步长 ，具 有低等的精度 ode45 则采用龙格 -库塔 4 阶算法 ，用 5 阶公式作误差估计来调节步长 ，专业 word 可编辑具有中等的精度 5 ezplot (x,y,tmin,tmax) ：符号函数的作图命令 x,y 为关于参数 t 的 符号函数 ，tmin,tmax 为 t 的取值范围 6 inline() ：建立一个内联函数 格式 ：inline(expr, var1, var2, ) ，注 意括号里的表达式要加引号 三、实验内容1. 几个可以直接用 Mat

7、lab 求微分方程精确解的例子 ： 例 1：求解微分方程 dy 2xy xe x ，并加以验证 例：Q = dblquad(inline(y*sin(x), pi, 2*pi, 0, pi)dx求解本问题的 Matlab 程序为 ：syms x y%line1y=dsolve(Dy+2*x*y=x*exp(-x2),x) %line2 diff(y,x)+2*x*y-x*exp(-x2)%line3simplify(diff(y,x)+2*x*y-x*exp(-x2)%line4说明：(1) 行 line1 是用命令定义 x,y 为符号变量这里可以不写 ，但为确保正确 性，建议写上；专业 wo

8、rd 可编辑(2) 行 line2 是用命令求出的微分方程的解 ： 1/2*exp(-x2)*x2+exp(-x2)*C1(3) 行 line3 使用所求得的解 这里是 将解代入原微分方程 ，结果应该为 0， 但这里给出 ：-x3*exp(-x2)-2*x*exp(-x2)*C1+2*x*(1/2*exp(-x2)*x2+exp(-x2 )*C1)(4) 行 line4 用 simplify() 函数对上式进行化简 ，结果为 0， 表明 y y(x) 的确是微分方程的解 例 2：求微分方程 xy y ex 0 在初始条件 y(1) 2e 下的特解 ，并画出解函 数的图形 求解本问题的 Matl

9、ab 程序为 ：syms x yy=dsolve(x*Dy+y-exp(x)=0,y(1)=2*exp(1),x)ezplot(y)微分方程的特解为 ： y=1/x*exp(x)+1/x* exp (1) (Matlab 格式)，即 e exy e e ， 解函数的图形如图 1：x专业 word 可编辑1/x exp(x)+1/x exp(1)-6 -4 -2 0 2 4 6x50403020100-10-20-30图1dx 5x y et例 3：求微分方程组 dt 在初始条件 x|t 0 1,y|t 0 0 下的特dy x 3y 0 t 0 t 0 dt解 ，并画出解函数的图形 求解本问题的

10、 Matlab 程序为 ：syms x y tx,y=dsolve(Dx+5*x+y=exp(t),Dy-x-3*y=0,x(0)=1,y(0)=0,t)simple(x);simple(y);ezplot(x,y,0,1.3);axis auto微分方程的特解 (式子特别长 )以及解函数的图形均略 2. 用 ode23 、ode45 等求解非刚性的标准形式的一阶常微分方程 (组)的初值 问题的数值解 (近似解 )dy 2y 2x2 2x例4：求解微分方程初值问题 dx 2y 2x 2x的数值解 ，求解范围为区y(0) 1专业 word 可编辑间0, 0.5fun=inline(-2*y+2*

11、x2+2*x,x,y);x,y=ode23(fun,0,0.5,1);x;y;plot(x,y,o-) xans =0.00000.04000.09000.14000.19000.24000.29000.34000.39000.44000.49000.5000 yans =1.00000.92470.84340.77540.71990.67640.64400.62220.61050.60840.61540.6179图形结果为图 20.950.90.850.80.750.70.650.60 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5图2专业 word

12、 可编辑例 5： 求解描述振荡器的经典的 Ver der Pol 微分方程d 2y(1y2) dy y 0, y(0) 1, y(0) 0, 7.dt 2dt分析：令 x1 y,x2 dx1 ,则 dx1 x2,dx2(1 x12)x2 x1.dtdtdt先编写函数文件 verderpol.m ：function xprime = verderpol(t,x)global mu;xprime = x(2);mu*(1-x(1)2)*x(2)-x(1);再编写命令文件 vdp1.m ：global mu;mu = 7;y0=1;0t,x = ode45(verderpol,0,40,y0);x1

13、=x(:,1);x2=x(:,2);plot(t,x1)图形结果为图 3专业 word 可编辑2.521.510.50-0.5-1-1.5-2-2.505 10 15 20 25 30 35图33. 用 Euler 折线法求解前面讲到过 ， 能够求解的微分方程也是十分有限的面介绍用 Euler 折线法求微分方程的数值解 （近似解 ）的方法 Euler 折线法求解的基本思想是将微分方程初值问题ddyx f (x,y),dxy(x hh) y(x) 替代微商y(x0) y0化成一个代数方程 ， 即差分方程 ， 主要步骤是用差商ddyx ,于是：y(xk h) y(xk )hf(xk , y(xk

14、),y0 y(x0 )记 xk 1 xk h,yk y(xk)，从而 yk 1 y(xk h) ，则有y0 y(x0 ),xk 1 xk h, k 0,1,2, ,n 1 yk 1 yk hf (xk ,yk).例 6：用 Euler 折线法求解微分方程初值问题专业 word 可编辑ddyx y 2yx2 ,y(0) 1的数值解 （步长 h 取 0.4），求解范围为区间 0,2 解：本问题的差分方程为x0 0,y0 1,h 0.4,xk 1 xk h,yk 1 yk hf(xk,yk)(其中： f(x,y)k 0,1,2, ,n 1 2yx2).y相应的 Matlab 程序见附录 1数据结果为

15、 ：01.00000.40001.40000.80002.12331.20003.11451.60004.45932.00006.3074图形结果见图 4：654320 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2专业 word 可编辑1图4特别说明 ：本问题可进一步利用四阶 Runge-Kutta 法求解 ， 读者可将两个 结果在一个图中显示 ，并和精确值比较 ，看看哪个更 “精确 ”？(相应的 Matlab程序参见附录 2)四、自己动手1. 求微分方程 (x2 1)y 2xy sin x 0的通解2. 求微分方程 y 2y 5y exsin x的通解3. 求微分方

16、程组dx x y 0 dtdy x y 0 dt在初始条件 x|t 0 1, y|t 0 0 下的特解 ，并画出解函数 y f (x) 的图形4. 分别用 ode23 、ode45 求上述第 3 题中的微分方程初值问题的数值解(近似解 )，求解区间为 t 0,2 利用画图来比较两种求解器之间的差异5. 用 Euler 折线法求解微分方程初值问题y y 12yx32yy(0) 1的数值解 (步长 h 取 0.1)，求解范围为区间 0,2 6. 用四阶 Runge-Kutta 法求解微分方程初值问题xy y e cosx, y(0) 1的数值解 (步长 h 取 0.1)，求解范围为区间 0,3 专

17、业 word 可编辑四阶 Runge-Kutta 法的迭代公式为 (Euler 折线法实为一阶 Runge-Kutta法)：y0 y(x0 ),xk 1 xk h,hyk 1 yk (L1 2L2 2L3 L4)6L1 f (xk ,yk )k 0,1,2, ,n 1hhL2 f (xk 2,yk 2 L1) hhL3 f (xk h2, yk h2L2 )L4 f (xk h,yk hL3 )相应的 Matlab 程序参见附录 2试用该方法求解第 5 题中的初值问题 7. 用 ode45 方法求上述第 6 题的常微分方程初值问题的数值解( 近似解 )，从而利用画图来比较两者间的差异 五、附录附录 1： (fulu1.m)clearf=sym(y+2*x/y2);a=0;b=2;h=0.4;n=(b-a)/h+1;x=0;y=1;szj=x,y;专业 word 可编辑for i=1:n-1y=y+h*subs(f,x,y,x,y);x=x+h;szj=szj;x,y;endszjplot(szj(:,1),szj(:,2)附录 2： (fulu2.m)clea