用MATLAB解常微分方程.pdf
实验四实验四 求微分方程的解求微分方程的解一、问题背景与实验目的一、问题背景与实验目的实际应用问题通过数学建模所归纳而得到的方程，绝大多数都是微分方程，真正能得到代数方程的时机很少另一方面，能够求解的微分方程也是十分有限的，特别是高阶方程和偏微分方程组 这就要求我们必须研究微分方程组的解法， 既要研究微分方程 组 的解析解法 精确解 ， 更要研究微分方程 组的数值解法近似解 对微分方程组的解析解法(精确解)，Matlab 有专门的函数可以用，本实验将作一定的介绍本实验将主要研究微分方程(组)的数值解法 近似解 ， 重点介绍 Euler 折线法二、相关函数命令及简介二、相关函数命令及简介1dsolvedsolve(equ1,equ2,)：Matlab 求微分方程的解析解equ1、equ2、为方程或条件 写方程或条件时用 DyDy 表示表示 y y 关于自变量的一阶导数关于自变量的一阶导数，用用用用 D2yD2y 表示表示 y y 关于自变量的二阶导数，依此类推关于自变量的二阶导数，依此类推2simplify(ssimplify(s) )：对表达式 s 使用 maple 的化简规则进行化简例如：syms xsimplify(sin(x)2 + cos(x)2)ans=13r,how=simple(s)r,how=simple(s)：由于 Matlab 提供了多种化简规则，simple 命令就是对表达式 s 用各种规则进行化简，然后用 r 返回最简形式，how 返回形成这种形式所用的规则例如：syms xr,how=simple(cos(x)2-sin(x)2)r = cos(2*x)how = combine4T,Y = solver(odefun,tspan,y0) 求微分方程的数值解说明：(1) 其中的 solversolver 为命令为命令 ode45ode45、ode23ode23、ode113ode113、ode15sode15s、ode23sode23s、ode23tode23t、ode23tbode23tb 之一之一dy f (t, y)(2) odefun 是显式常微分方程：dty(t0) y0(3) 在积分区间 tspan=t0,tf上，从t0到tf，用初始条件y0求解(4) 要 获 得 问 题 在 其 他 指 定 时 间 点t0,t1,t2,上 的 解 ， 则 令tspan=t0,t1,t2,tf要求是单调的(5) 因为没有一种算法可以有效地解决所有的 ODE 问题，为此，Matlab 提供了多种求解器 Solver，对于不同的 ODE 问题，采用不同的 Solver求解器求解器ODEODE 类型类型特点特点说明说明SolverSolverode45非刚性单步算法；4、5 阶 Runge-Kutta大部分场合的首选算法方程；累计截断误差达(x)3单步算法；2、3 阶 Runge-Kutta使用于精度较低的情形方程；累计截断误差达(x)3多步法；Adams 算法；高低精计算时间比 ode45 短度均可到103106ode23tode15sode23s适度刚性刚性刚性采用梯形算法多步法；Gears 反向数值微分；精度中等单步法；2 阶 Rosebrock 算法；低精度梯形算法；低精度适度刚性情形假设 ode45 失效时，可尝试使用当精度较低时， 计算时间比 ode15s 短当精度较低时， 计算时间比 ode15s 短ode23非刚性ode113非刚性ode23tb刚性(6) 要特别的是：ode23、ode45 是极其常用的用来求解非刚性的标准形式的一阶常微分方程(组)的初值问题的解的 Matlab 的常用程序，其中：ode23 采用龙格-库塔 2 阶算法，用 3 阶公式作误差估计来调节步长，具有低等的精度ode45 则采用龙格-库塔 4 阶算法，用 5 阶公式作误差估计来调节步长，具有中等的精度5ezplotezplot(x,y,tmin,tmax)：符号函数的作图命令x,y 为关于参数 t 的符号函数，tmin,tmax 为 t 的取值范围6inline()inline()：建立一个内联函数：建立一个内联函数格式：inline(expr, var1, var2,) ，注意括号里的表达式要加引号例：Q = dblquad(inline(y*sin(x), pi, 2*pi, 0, pi)三、实验内容三、实验内容1.几个可以直接用 Matlab 求微分方程精确解的例子：2dy例 1：求解微分方程 2xy xex，并加以验证dx求解本问题的 Matlab 程序为：syms x y%line1y=dsolve(Dy+2*x*y=x*exp(-x2),x)%line2diff(y,x)+2*x*y-x*exp(-x2)%line3simplify(diff(y,x)+2*x*y-x*exp(-x2)%line4说明：(1) 行 line1 是用命令定义定义 x,yx,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，说明y y(x)确实是微分方程的解例 2：求微分方程xyy ex 0在初始条件y(1) 2e下的特解，并画出解函数的图形求解本问题的 Matlab 程序为：syms x yy=dsolve(x*Dy+y-exp(x)=0,y(1)=2*exp(1),x)ezplot(y)ezplot(y)e ex微分方程的特解为：y=1/x*exp(x)+1/x* exp (1) (Matlab 格式)，即y ，x解函数的图形如图 1：1/x exp(x)+1/x exp(1)50403020100-10-20-30-6-4-20 x246图 1dxt5x y edt例 3：求微分方程组在初始条件x |t01, y |t0 0下的特解，dy x 3y 0dt并画出解函数的图形求解本问题的 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 等求解非刚性的标准形式的一阶常微分方程(组)的初值问题的数值解(近似解)dy2 2y 2x 2x例 4：求解微分方程初值问题dx的数值解，求解范围为区y(0) 1间0, 0.5fun=inline(-2*y+2*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图形结果为图 210.950.90.850.80.750.70.650.600.050.10.150.20.250.30.350.40.450.5图 2例 5：求解描述振荡器的经典的 Ver der Pol 微分方程d2y2dy(1 y ) y 0, y(0) 1, y(0) 0, 7.2dtdt分析：令x1 y,x2dx1dxdx2,则1 x2,2(1 x1)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=x(:,1);x2=x(:,2);plot(t,x1)图形结果为图 32.521.510.50-0.5-1-1.5-2-2.50510152025303540图 33. 用 Euler 折线法求解前面讲到过，能够求解的微分方程也是十分有限的下面介绍用 Euler 折线法求微分方程的数值解(近似解)的方法Euler 折线法求解的基本思想是将微分方程初值问题dy f (x, y),dxy(x0) y0化成一个代数方程， 即差分方程， 主要步骤是用差商于是：dyy(x h) y(x)替代微商,dxh记xk1y(xk h) y(xk) f (xk, y(xk),hy0 y(x0) xk h, yk y(xk)，从而yk1 y(xk h)，则有y0 y(x0),k 0,1,2,n 1xk1 xk h,y y hf (x , y ).kkkk1例 6：用 Euler 折线法求解微分方程初值问题2xdy y ,2ydxy(0) 1的数值解(步长 h 取 0.4)，求解范围为区间0,2解：本问题的差分方程为x0 0, y01,h 0.4,xk1 xk h,k 0,1,2,n 12xyk1 yk hf (xk, yk)（其中：f (x, y) y 2).y相应的 Matlab 程序见附录 1数据结果为：01.00000.40001.40000.80002.12331.20003.11451.60004.45932.00006.3074图形结果见图 4：765432100.20.40.60.811.21.41.61.82图 4特别说明：本问题可进一步利用四阶 Runge-Kutta 法求解，读者可将两个结果在一个图中显示， 并和精确值比较， 看看哪个更 “精确” ？ 相应的 Matlab 程序参见附录 2 四、自己动手四、自己动手1. 求微分方程(x21)y2xy sin x 0的通解2. 求微分方程y2y5y exsin x的通解3. 求微分方程组dx x y 0dtdy x y 0dt在初始条件x |t01, y |t0 0下的特解，并画出解函数y f (x)的图形4. 分别用 ode23、ode45 求上述第 3 题中的微分方程初值问题的数值解(近似解)，求解区间为t0,2利用画图来比较两种求解器之间的差异5. 用 Euler 折线法求解微分方程初值问题12x2y y 3,yy(0) 1的数值解(步长 h 取 0.1)，求解范围为区间0,26. 用四阶 Runge-Kutta 法求解微分方程初值问题y y excosx,y(0) 1的数值解(步长 h 取 0.1)，求解范围为区间0,3四阶Runge-Kutta 法的迭代公式为 (Euler 折线法实为一阶Runge-Kutta法)：y0 y(x0),xk1 xk h,hyk1 yk(L1 2L2 2L3 L4)6k 0,1,2,n 1L1 f (xk, yk)hhL2 f (xk, ykL1)22hhL3 f (xk, ykL2)22L4 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;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)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)
