用MATLAB解常微分方程

７２实验四求微分方程的解一、问题背景与实验目的实际应用问题通过数学建模所归纳而得到的方程，绝大多数都是微分方程，真正能得到代数方程的机会很少．另一方面，能够求解的微分方程也是十分有限的，特别是高阶方程和偏微分方程（组）．这就要求我们必须研究微分方程（组）的解法，既要研究微分方程（组）的解析解法（精确解），更要研究微分方程（组）的数值解法（近似解）．对微分方程（组）的解析解法(精确解)，Matlab有专门的函数可以用，本实验将作一定的介绍．本实验将主要研究微分方程(组)的数值解法（近似解），重点介绍Euler折线法．二、相关函数（命令）及简介1．dsolve('equ1','equ2',…)：Matlab求微分方程的解析解．equ1、equ2、…为方程（或条件）．写方程（或条件）时用Dy表示y关于自变量的一阶导数，用用D2y表示y关于自变量的二阶导数，依此类推．2．simplify(s)：对表达式s使用maple的化简规则进行化简．例如：symsxsimplify(sin(x)^2+cos(x)^2)ans=13．[r,how]=simple(s)：由于Matlab提供了多种化简规则，simple命令就是对表达式s用各种规则进行化简，然后用r返回最简形式，how返回形成这种形式所用的规则．例如：symsx[r,how]=simple(cos(x)^2-sin(x)^2)r=cos(2*x)how=combine4．[T,Y]=solver(odefun,tspan,y0)求微分方程的数值解．说明：(1)其中的solver为命令ode45、ode23、ode113、ode15s、ode23s、ode23t、ode23tb之一．(2)odefun是显式常微分方程：00)(),(ytyytfdtdy(3)在积分区间tspan=],[0ftt上，从0t到ft，用初始条件0y求解．７３(4)要获得问题在其他指定时间点,210,,ttt上的解，则令tspan=],,,[,210ftttt（要求是单调的）．(5)因为没有一种算法可以有效地解决所有的ODE问题，为此，Matlab提供了多种求解器Solver，对于不同的ODE问题，采用不同的Solver．求解器SolverODE类型特点说明ode45非刚性单步算法；4、5阶Runge-Kutta方程；累计截断误差达3)(x大部分场合的首选算法ode23非刚性单步算法；2、3阶Runge-Kutta方程；累计截断误差达3)(x使用于精度较低的情形ode113非刚性多步法；Adams算法；高低精度均可到6310~10计算时间比ode45短ode23t适度刚性采用梯形算法适度刚性情形ode15s刚性多步法；Gear's反向数值微分；精度中等若ode45失效时，可尝试使用ode23s刚性单步法；2阶Rosebrock算法；低精度当精度较低时，计算时间比ode15s短ode23tb刚性梯形算法；低精度当精度较低时，计算时间比ode15s短(6)要特别的是：ode23、ode45是极其常用的用来求解非刚性的标准形式的一阶常微分方程(组)的初值问题的解的Matlab的常用程序，其中：ode23采用龙格-库塔2阶算法，用3阶公式作误差估计来调节步长，具有低等的精度．ode45则采用龙格-库塔4阶算法，用5阶公式作误差估计来调节步长，具有中等的精度．5．ezplot(x,y,[tmin,tmax])：符号函数的作图命令．x,y为关于参数t的符号函数，[tmin,tmax]为t的取值范围．6．inline()：建立一个内联函数．格式：inline('expr','var1','var2',…)，注意括号里的表达式要加引号．例：Q=dblquad(inline('y*sin(x)'),pi,2*pi,0,pi)７４三、实验内容1.几个可以直接用Matlab求微分方程精确解的例子：例1：求解微分方程22xxexydxdy，并加以验证．求解本问题的Matlab程序为：symsxy%line1y=dsolve('Dy+2*x*y=x*exp(-x^2)','x')%line2diff(y,x)+2*x*y-x*exp(-x^2)%line3simplify(diff(y,x)+2*x*y-x*exp(-x^2))%line4说明：(1)行line1是用命令定义x,y为符号变量．这里可以不写，但为确保正确性，建议写上；(2)行line2是用命令求出的微分方程的解：1/2*exp(-x^2)*x^2+exp(-x^2)*C1(3)行line3使用所求得的解．这里是将解代入原微分方程，结果应该为0，但这里给出：-x^3*exp(-x^2)-2*x*exp(-x^2)*C1+2*x*(1/2*exp(-x^2)*x^2+exp(-x^2)*C1)(4)行line4用simplify()函数对上式进行化简，结果为0，表明)(xyy的确是微分方程的解．例2：求微分方程0'xeyxy在初始条件ey2)1(下的特解，并画出解函数的图形．求解本问题的Matlab程序为：symsxyy=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格式)，即xeeyx，解函数的图形如图1：７５-6-4-20246-30-20-1001020304050x1/xexp(x)+1/xexp(1)图1例3：求微分方程组035yxdtdyeyxdtdxt在初始条件0|,1|00ttyx下的特解，并画出解函数的图形．求解本问题的Matlab程序为：symsxyt[x,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]);axisauto微分方程的特解(式子特别长)以及解函数的图形均略．2.用ode23、ode45等求解非刚性的标准形式的一阶常微分方程(组)的初值问题的数值解(近似解)．例4：求解微分方程初值问题1)0(2222yxxydxdy的数值解，求解范围为区间[0,0.5]．fun=inline('-2*y+2*x^2+2*x','x','y');[x,y]=ode23(fun,[0,0.5],1);x';y';plot(x,y,'o-')x'ans=0.00000.04000.09000.14000.19000.24000.29000.34000.39000.44000.49000.5000y'７６ans=1.00000.92470.84340.77540.71990.67640.64400.62220.61050.60840.61540.6179图形结果为图2．00.050.10.150.20.250.30.350.40.450.50.60.650.70.750.80.850.90.951图2例5：求解描述振荡器的经典的VerderPol微分方程.7,0)0(',1)0(,0)1(222yyydtdyydtyd分析：令,,121dtdxxyx则.)1(,1221221xxxdtdxxdtdx先编写函数文件verderpol.m：functionxprime=verderpol(t,x)globalmu;xprime=[x(2);mu*(1-x(1)^2)*x(2)-x(1)];再编写命令文件vdp1.m：globalmu;mu=7;y0=[1;0][t,x]=ode45('verderpol',[0,40],y0);x1=x(:,1);x2=x(:,2);plot(t,x1)图形结果为图3．７７0510152025303540-2.5-2-1.5-1-0.500.511.522.5图33.用Euler折线法求解前面讲到过，能够求解的微分方程也是十分有限的．下面介绍用Euler折线法求微分方程的数值解(近似解)的方法．Euler折线法求解的基本思想是将微分方程初值问题00)(),,(yxyyxfdxdy化成一个代数方程，即差分方程，主要步骤是用差商hxyhxy)()(替代微商dxdy,于是：)()),(,()()(00xyyxyxfhxyhxykkkk记)(,1kkkkxyyhxx，从而)(1hxyykk，则有1,,2,1,0).,(,),(1100nkyxhfyyhxxxyykkkkkk例6：用Euler折线法求解微分方程初值问题1)0(,22yyxydxdy的数值解(步长h取0.4)，求解范围为区间[0,2]．解：本问题的差分方程为７８1,,2,1,0).2),(),(,,4.0,1,021100nkyxyyxfyxhfyyhxxhyxkkkkkk（其中：相应的Matlab程序见附录1．数据结果为：01.00000.40001.40000.80002.12331.20003.11451.60004.45932.00006.3074图形结果见图4：00.20.40.60.811.21.41.61.821234567图4特别说明：本问题可进一步利用四阶Runge-Kutta法求解，读者可将两个结果在一个图中显示，并和精确值比较，看看哪个更“精确”？（相应的Matlab程序参见附录2）．四、自己动手1.求微分方程0sin2')1(2xxyyx的通解．2.求微分方程xeyyyxsin5'2''的通解．3.求微分方程组00yxdtdyyxdtdx７９在初始条件0|,1|00ttyx下的特解，并画出解函数()yfx的图形．4.分别用ode23、ode45求上述第3题中的微分方程初值问题的数值解(近似解)，求解区间为[0,2]t．利用画图来比较两种求解器之间的差异．5.用Euler折线法求解微分方程初值问题1)0(,12'32yyxyy的数值解(步长h取0.1)，求解范围为区间[0,2]．6.用四阶Runge-Kutta法求解微分方程初值问题1)0(,cos'yxeyyx的数值解(步长h取0.1)，求解范围为区间[0,3]．四阶Runge-Kutta法的迭代公式为(Euler折线法实为一阶Runge-Kutta法)：1,,2,1,0),()2,2()2,2(),()22(6,),(342312143211100nkhLyhxfLLhyhxfLLhyhxfLyxfLLLLLhyyhxxxyykkkkkkkkkkkk相应的Matlab程序参见附录2．试用该方法求解第5题中的初值问题．7.用ode45方法求上述第6题的常微分方程初值问题的数值解(近似解)，从而利用画图来比较两者间的差异．五、附录附录1：(fulu1.m)clearf=sym('y+2*x/y^2');a=0;b=2;h=0.4;n=(b-a)/h+1;x=0;y=1;８０szj=[x,y];fori=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];fori=1:n-1l1=subs(f,{'x','y'},{x,y});l2=subs(f,{'x','y'},{x+h/2,y+l1*