还剩7页未读,继续阅读
文本内容:
实验四 求微分方程的解
一、问题背景与实验目的实际应用问题通过数学建模所归纳而得到的方程,绝大多数都是微分方程,真正能得到代数方程的机会很少.另一方面,能够求解的微分方程也是十分有限的,特别是高阶方程和偏微分方程(组).这就要求我们必须研究微分方程(组)的解法,既要研究微分方程(组)的解析解法(精确解),更要研究微分方程(组)的数值解法(近似解).对微分方程(组)的解析解法精确解,__tlab有专门的函数可以用,本实验将作一定的介绍.本实验将主要研究微分方程组的数值解法(近似解),重点介绍Euler折线法.
二、相关函数(命令)及简介1.dsolveequ1equ2…__tlab求微分方程的解析解.equ
1、equ
2、…为方程(或条件).写方程(或条件)时用Dy表示y关于自变量的一阶导数,用用D2y表示y关于自变量的二阶导数,依此类推.2.simplifys对表达式s使用__ple的化简规则进行化简.例如symsxsimplifysinx^2+cosx^2ans=13.[rhow]=______s由于__tlab提供了多种化简规则,______命令就是对表达式s用各种规则进行化简,然后用r返回最简形式,how返回形成这种形式所用的规则.例如symsx[rhow]=______cosx^2-sinx^2r=cos2*xhow=combine4.[TY]=solverodefuntspany0求微分方程的数值解.说明1其中的solver为命令ode
45、ode
23、ode
113、ode15s、ode23s、ode23t、ode23tb之一.2odefun是显式常微分方程3在积分区间tspan=上,从到,用初始条件求解.4要获得问题在其他指定时间点上的解,则令tspan=(要求是单调的).5因为没有一种算法可以有效地解决所有的ODE问题,为此,__tlab提供了多种求解器Solver,对于不同的ODE问题,采用不同的Solver.求解器SolverODE类型特点说明ode45非刚性单步算法;
4、5阶Runge-Kutta方程;累计截断误差达大部分场合的首选算法ode23非刚性单步算法;
2、3阶Runge-Kutta方程;累计截断误差达使用于精度较低的情形ode113非刚性多步法;Adams算法;高低精度均可到计算时间比ode45短ode23t适度刚性采用梯形算法适度刚性情形ode15s刚性多步法;Gears反向数值微分;精度中等若ode45失效时,可尝试使用ode23s刚性单步法;2阶Rosebrock算法;低精度当精度较低时,计算时间比ode15s短ode23tb刚性梯形算法;低精度当精度较低时,计算时间比ode15s短6要特别的是ode
23、ode45是极其常用的用来求解非刚性的标准形式的一阶常微分方程组的初值问题的解的__tlab的常用程序,其中ode23采用龙格-库塔2阶算法,用3阶公式作误差估计来调节步长,具有低等的精度.ode45则采用龙格-库塔4阶算法,用5阶公式作误差估计来调节步长,具有中等的精度.5.ezplotxy[tmint__x]符号函数的作图命令.xy为关于参数t的符号函数,[tmint__x]为t的取值范围.6.inline建立一个内联函数.格式inlineexprvar1var2…,注意括号里的表达式要加引号.例Q=dblquadinliney*sinxpi2*pi0pi
三、实验内容
1.几个可以直接用__tlab求微分方程精确解的例子例1求解微分方程,并加以验证.求解本问题的__tlab程序为symsxy%line1y=dsolveDy+2*x*y=x*exp-x^2x%line2diffyx+2*x*y-x*exp-x^2%line3simplifydiffyx+2*x*y-x*exp-x^2%line4说明1行line1是用命令定义xy为符号变量.这里可以不写,但为确保正确性,建议写上;2行line2是用命令求出的微分方程的解1/2*exp-x^2*x^2+exp-x^2*C13行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*C14行line4用simplify函数对上式进行化简,结果为0,表明的确是微分方程的解.例2求微分方程在初始条件下的特解,并画出解函数的图形.求解本问题的__tlab程序为symsxyy=dsolvex*Dy+y-expx=0y1=2*exp1xezploty微分方程的特解为y=1/x*expx+1/x*exp1__tlab格式,即,解函数的图形如图1图1例3求微分方程组在初始条件下的特解,并画出解函数的图形.求解本问题的__tlab程序为symsxyt[xy]=dsolveDx+5*x+y=exptDy-x-3*y=0x0=1y0=0t______x;______y;ezplotxy[
01.3];axisauto微分方程的特解式子特别长以及解函数的图形均略.
2.用ode
23、ode45等求解非刚性的标准形式的一阶常微分方程组的初值问题的数值解近似解.例4求解微分方程初值问题的数值解,求解范围为区间[
00.5].fun=inline-2*y+2*x^2+2*xxy;[xy]=ode23fun[
00.5]1;x;y;plotxyo-xans=
0.
00000.
04000.
09000.
14000.
19000.
24000.
29000.
34000.
39000.
44000.
49000.5000yans=
1.
00000.
92470.
84340.
77540.
71990.
67640.
64400.
62220.
61050.
60840.
61540.6179图形结果为图2.图2例5求解描述振荡器的经典的VerderPol微分方程分析令则先编写函数文件verderpol.m functionxprime=verderpoltxglobalmu;xprime=[x2;mu*1-x1^2*x2-x1];再编写命令文件vdp
1.m globalmu;mu=7;y0=[1;0][tx]=ode45verderpol
[040]y0;x1=x:1;x2=x:2;plottx1图形结果为图3.图
33.用Euler折线法求解前面讲到过,能够求解的微分方程也是十分有限的.下面介绍用Euler折线法求微分方程的数值解近似解的方法.Euler折线法求解的基本思想是将微分方程初值问题化成一个代数方程,即差分方程,主要步骤是用差商替代微商于是记,从而,则有例6用Euler折线法求解微分方程初值问题的数值解步长h取
0.4,求解范围为区间
[02].解本问题的差分方程为相应的__tlab程序见附录1.数据结果为
01.
00000.
40001.
40000.
80002.
12331.
20003.___
51.
60004.
45932.
00006.3074图形结果见图4图4特别说明本问题可进一步利用四阶Runge-Kutta法求解,读者可将两个结果在一个图中显示,并和精确值比较,看看哪个更“精确”?(相应的__tlab程序参见附录2).
四、自己动手
1.求微分方程的通解.
2.求微分方程的通解.
3.求微分方程组在初始条件下的特解,并画出解函数的图形.
4.分别用ode
23、ode45求上述第3题中的微分方程初值问题的数值解近似解,求解区间为.利用画图来比较两种求解器之间的差异.
5.用Euler折线法求解微分方程初值问题的数值解步长h取
0.1,求解范围为区间
[02].
6.用四阶Runge-Kutta法求解微分方程初值问题的数值解步长h取
0.1,求解范围为区间
[03].四阶Runge-Kutta法的迭代公式为Euler折线法实为一阶Runge-Kutta法相应的__tlab程序参见附录2.试用该方法求解第5题中的初值问题.
7.用ode45方法求上述第6题的常微分方程初值问题的数值解近似解,从而利用画图来比较两者间的差异.
五、附录附录1fulu
1.mclearf=symy+2*x/y^2;a=0;b=2;h=
0.4;n=b-a/h+1;x=0;y=1;szj=[xy];fori=1:n-1y=y+h*subsf{xy}{xy};x=x+h;szj=[szj;xy];endszjplotszj:1szj:2附录2fulu
2.mclearf=symy-expx*cosx;a=0;b=3;h=
0.1;n=b-a/h+1;x=0;y=1;szj=[xy];fori=1:n-1l1=subsf{xy}{xy};l2=subsf{xy}{x+h/2y+l1*h/2};l3=subsf{xy}{x+h/2y+l2*h/2};l4=subsf{xy}{x+hy+l3*h};y=y+h*l1+2*l2+2*l3+l4/6;x=x+h;szj=[szj;xy];endszjplotszj:1szj:2。