用MATLAB解常微分方程.doc

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

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

3、r = cos(2*x)how = combine4T,Y = solver(odefun,tspan,y0) 求微分方程的数值解说明：(1) 其中的 solver 为命令 ode45、ode23、ode113、ode15s、ode23s、ode23t、ode23tb 之一(2) odefun 是显式常微分方程： 0)(,ytfd(3) 在积分区间 tspan= 上，从 到 ，用初始条件 求解,0ftf 0y(4) 要获得问题在其他指定时间点 上的解，则令 tspan= ,210t（要求是单调的） ,210ftt(5) 因为没有一种算法可以有效地解决所有的 ODE 问题，为此，Matlab 提

4、供了多种求解器 Solver，对于不同的 ODE 问题，采用不同的 Solver求解器SolverODE 类型 特点 说明ode45 非刚性 单步算法；4、5 阶 Runge-Kutta 方程；累计截断误差达3)(x大部分场合的首选算法ode23 非刚性 单步算法；2、3 阶 Runge-Kutta 方程；累计截断误差达)(x使用于精度较低的情形ode113 非刚性 多步法；Adams 算法；高低精度均可到 6310计算时间比 ode45 短ode23t 适度刚性 采用梯形算法 适度刚性情形ode15s 刚性 多步法；Gears 反向数值微分；精度中等若 ode45 失效时，可尝试使用ode2

5、3s 刚性 单步法；2 阶 Rosebrock 算法；低精度当精度较低时，计算时间比 ode15s 短ode23tb 刚性 梯形算法；低精度 当精度较低时，计算时间比 ode15s 短(6) 要特别的是：ode23、ode45 是极其常用的用来求解非刚性的标准形式的一阶常微分方程(组)的初值问题的解的 Matlab 的常用程序，其中：ode23 采用龙格-库塔 2 阶算法，用 3 阶公式作误差估计来调节步长，具有低等的精度ode45 则采用龙格-库塔 4 阶算法，用 5 阶公式作误差估计来调节步长，具有中等的精度5ezplot(x,y,tmin,tmax) ：符号函数的作图命令 x,y 为关于

6、参数 t 的符号函数，tmin,tmax 为 t 的取值范围6inline()：建立一个内联函数格式：inline(expr, var1, var2,) ，注意括号里的表达式要加引号例：Q = dblquad(inline(y*sin(x), pi, 2*pi, 0, pi)三、实验内容1. 几个可以直接用 Matlab 求微分方程精确解的例子：例 1：求解微分方程 ，并加以验证2xeydx求解本问题的 Matlab 程序为：syms x y %line1y=dsolve(Dy+2*x*y=x*exp(-x2),x) %line2diff(y,x)+2*x*y-x*exp(-x2) %line

7、3simplify(diff(y,x)+2*x*y-x*exp(-x2) %line4说明：(1) 行 line1 是用命令定义 x,y 为符号变量这里可以不写，但为确保正确性，建议写上；(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， 表明的确

8、是微分方程的解)(xy例 2：求微分方程 在初始条件 下的特解，并画出解函0xeyey2)1(数的图形求解本问题的 Matlab 程序为：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 格式)，即 ，xey解函数的图形如图 1：-6 -4 -2 0 2 4 6-30-20-1001020304050x1/x exp(x)+1/x exp(1)图 1例 3：求微分方程组 在初始条件 下的特解，035yxdtet 0|,1|0ttyx并画出解函数

9、的图形求解本问题的 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 等求解非刚性的标准形式的一阶常微分方程(组)的初值问题的数值解( 近似解)例 4：求解微分方程初值问题 的数值解，求解范围为1)0(2yxdx区间0, 0.5fun=inline(-2*y+2*x2+2*x,x,y);x,y=ode23(fun,0,0.

10、5,1);x;y;plot(x,y,o-) xans =0.0000 0.0400 0.0900 0.1400 0.1900 0.24000.2900 0.3400 0.3900 0.4400 0.4900 0.5000 yans =1.0000 0.9247 0.8434 0.7754 0.7199 0.67640.6440 0.6222 0.6105 0.6084 0.6154 0.6179图形结果为图 20 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.50.60.650.70.750.80.850.90.951图 2例 5：求解描述振荡器的经典的

11、 Ver der Pol 微分方程 .7,0)(,1)0(,)1(22 yydtdty分析：令 则,121txx.)(, 12121 xtxt 先编写函数文件 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=x(:,1);x2=x(:,2);plot(t,x1)图形结果为图 30 5 10 15 20 25 30 35

12、40-2.5-2-1.5-1-0.500.511.522.5图 33. 用 Euler 折线法求解前面讲到过，能够求解的微分方程也是十分有限的下面介绍用 Euler 折线法求微分方程的数值解(近似解)的方法Euler 折线法求解的基本思想是将微分方程初值问题 0)(,yxfd化成一个代数方程，即差分方程，主要步骤是用差商 替代微商hxy)(,于是：dxy)(),(,0xyxyfhkkk记 ，从而 ，则有,1kkkyhx)(1hkk,2).,()10 nfykk 例 6：用 Euler 折线法求解微分方程初值问题1)0(,2yxd的数值解(步长 h 取 0.4)，求解范围为区间0,2 解：本问题

13、的差分方程为1,20).),( ),(,4.0,1, 210 nkyxxfyxhfyxykk （ 其 中 ：相应的 Matlab 程序见附录 1数据结果为：0 1.00000.4000 1.40000.8000 2.12331.2000 3.11451.6000 4.45932.0000 6.3074图形结果见图 4：0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 21234567图 4特别说明：本问题可进一步利用四阶 Runge-Kutta 法求解，读者可将两个结果在一个图中显示，并和精确值比较，看看哪个更“ 精确”？（相应的 Matlab 程序参见附录 2） 四、自

14、己动手1. 求微分方程 的通解0sin2)1(2xyx2. 求微分方程 的通解ey53. 求微分方程组 0yxdt在初始条件 下的特解，并画出解函数 的图形0|,1|0ttyx ()yfx4. 分别用 ode23、ode45 求上述第 3 题中的微分方程初值问题的数值解(近似解 )，求解区间为 利用画图来比较两种求解器之间的差异,2t5. 用 Euler 折线法求解微分方程初值问题 1)0(,23yx的数值解(步长 h 取 0.1)，求解范围为区间0,2 6. 用四阶 Runge-Kutta 法求解微分方程初值问题 1)0(,cosyxe的数值解(步长 h 取 0.1)，求解范围为区间0,3

15、四阶 Runge-Kutta 法的迭代公式为(Euler 折线法实为一阶 Runge-Kutta 法)：1,20),(2),()2(6,)(343121 43110 nkhLyxfLhfyxLLhxykkkkk 相应的 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;for i=1:n-1y=y+h*sub

16、s(f,x,y,x,y);x=x+h;szj=szj;x,y;endszjplot(szj(:,1),szj(:,2)附录 2：(fulu2.m)clearf=sym(y-exp(x)*cos(x);a=0;b=3;h=0.1;n=(b-a)/h+1;x=0;y=1;szj=x,y;for i=1:n-1l1=subs(f,x,y,x,y);l2=subs(f,x,y,x+h/2,y+l1*h/2);l3=subs(f,x,y,x+h/2,y+l2*h/2);l4=subs(f,x,y,x+h,y+l3*h);y=y+h*(l1+2*l2+2*l3+l4)/6;x=x+h;szj=szj;x,y;endszjplot(szj(:,1),szj(:,2)