计算机应用基础-5-常微分方程求解

1、第五章第五章 常微分方程的求解常微分方程的求解5.1 常微分方程基本特征常微分方程基本特征5.2 常微分方程的初值问题求解常微分方程的初值问题求解5.3 常微分方程的边值问题求解常微分方程的边值问题求解5.4 高阶常微分方程的求解高阶常微分方程的求解5 5.1 常微分方程的数值解常微分方程的数值解一般微分方程式描述系统内部变量的变化率如何受系一般微分方程式描述系统内部变量的变化率如何受系统内部变量和外部激励，如输入的影响。统内部变量和外部激励，如输入的影响。当常微分方程式能够解析求解时，可用当常微分方程式能够解析求解时，可用MATLAB的的符号工具箱中的功能找到精确解。符号工具箱中的功能找到精

2、确解。在微分方程难以获得解析解的情况下，可以方便地在在微分方程难以获得解析解的情况下，可以方便地在数值上求解。数值上求解。微分方程数值解是数值计算的基本内容，由于微分方微分方程数值解是数值计算的基本内容，由于微分方程的多样性，有不同的解法。程的多样性，有不同的解法。MATLAB给出的给出的7种解种解法分别由法分别由7个不同的函数来完成。个不同的函数来完成。 一、刚性问题一、刚性问题(Stiff Problem)Stiff Problem)1、刚性方程、刚性方程对于一个常微分方程组，如果其对于一个常微分方程组，如果其Jacobian矩阵的特征矩阵的特征值相差十分悬殊，那么这个方程组就称为刚性方程

3、组。值相差十分悬殊，那么这个方程组就称为刚性方程组。对于刚性方程组，为保持解法的稳定，步长选取很对于刚性方程组，为保持解法的稳定，步长选取很困难。有些解法不能用来解刚性方程组，而有些解困难。有些解法不能用来解刚性方程组，而有些解法对稳定性的要求不严格，可以用来解决刚性问题。法对稳定性的要求不严格，可以用来解决刚性问题。5 5.1 常微分方程的数值解常微分方程的数值解MATLAB为解常微分方程初值问题提供的命令包为解常微分方程初值问题提供的命令包括：微分方程结算命令、被解算命令、调用的常括：微分方程结算命令、被解算命令、调用的常微分方程文件格式命令、积分算法参数选项微分方程文件格式命令、积分算法

4、参数选项(Options)处理命令以及输出处理命令等。处理命令以及输出处理命令等。这些命令的含义特点和使用范围如表所示。这些命令的含义特点和使用范围如表所示。 5 5.1 常微分方程的数值解常微分方程的数值解5 5.1 常微分方程的数值解常微分方程的数值解5.1 常微分方程基本特征常微分方程基本特征常微分方程（组）包括初值问题和边值问题两种：常微分方程（组）包括初值问题和边值问题两种： （1） 初值问题：初值问题： 初始条件已知，形式如下初始条件已知，形式如下 0yayyxa,y, xfdxdy对于初值问题方程组，由多个常微分方程及其初始条对于初值问题方程组，由多个常微分方程及其初始条件构成件

5、构成（2） 边值问题：已知首末两个端点值，形式如下：边值问题：已知首末两个端点值，形式如下：; 0y, bg, bx; 0y, ag, axbxa,y, xfdxdy边值问题求解：边值问题求解： 打靶法：将边值问题转化成初值问题打靶法：将边值问题转化成初值问题 松弛法：包括有限差分法和配置法松弛法：包括有限差分法和配置法5.1 常微分方程基本特征常微分方程基本特征一一 欧拉公式欧拉公式 对于常微分方程初值问题对于常微分方程初值问题,在求解区间在求解区间a,b上作等距分割，步长上作等距分割，步长 记记xn = xn-1+h , n =1,2,m。用差商近似导数计算常微分方程。用差商近似导数计算常

6、微分方程。 做做 的在的在x = x0处的一阶向前差商得：处的一阶向前差商得： 于是得到：于是得到：mabh)(xyhxyxyxy)()()( 010)(,()( 000 xyxfxy)(,()()(0001xyxfhxyxy5.2 常微分方程的初值问题求解常微分方程的初值问题求解故故y(x1)的近似值的近似值y1可按可按 求得。类似地，由求得。类似地，由 得到计算近似值的向前欧拉公式得到计算近似值的向前欧拉公式)(,()( )()()( 1nnnnnnxyxfxyhxyxyxy以及),(),(00010001yxhfyyyxfhyy或),(1nnnnyxhfyy5.2 常微分方程的初值问题求

7、解常微分方程的初值问题求解 由由yn直接算出直接算出yn+1值的计算格式称为显式格式，值的计算格式称为显式格式，向前欧拉公式是显式格式。而欧拉方法的几何意义是向前欧拉公式是显式格式。而欧拉方法的几何意义是以以y1作为斜率，通过点作为斜率，通过点(x0, y0)做一条直线，它与直线做一条直线，它与直线x=x1的交点就是的交点就是y1。依此类推，是以。依此类推，是以yn+1作为斜率，经作为斜率，经过点的直线与直线过点的直线与直线x=xn+1的交点。故欧拉法也称为欧的交点。故欧拉法也称为欧拉折线法，如下图所示。拉折线法，如下图所示。 图4.1 欧拉折线法5.2 常微分方程的初值问题求解常微分方程的初

8、值问题求解 龙格龙格- -库塔法是求解常微分方程较常用的一种方法，它通库塔法是求解常微分方程较常用的一种方法，它通过巧妙的线性组合，在显式格式的情况下获得理想的计算精度，过巧妙的线性组合，在显式格式的情况下获得理想的计算精度，大大提高了计算速度。该方法的推导过程不是本课程要研究的大大提高了计算速度。该方法的推导过程不是本课程要研究的对象，本课程主要研究其实际应用对象，本课程主要研究其实际应用, ,，故直接给出各类龙格，故直接给出各类龙格- -库库塔公式。塔公式。 1 1、 二阶龙格二阶龙格- -库塔库塔 ),(),()(2121211hkyhxfkyxfkkkhyynnnnnn)2,2(),(

9、12121khyhxfkyxfkhkyynnnnnn二二 龙格龙格-库塔方法库塔方法 5.2 常微分方程的初值问题求解常微分方程的初值问题求解二阶龙格二阶龙格-库塔公式的局部截断误差为库塔公式的局部截断误差为O(h3), 是二阶是二阶精度的计算公式。也可建立高阶的龙格精度的计算公式。也可建立高阶的龙格-库塔公式，库塔公式，如三阶龙格如三阶龙格-库塔公式、四阶龙格库塔公式、四阶龙格-库塔公式。较常用库塔公式。较常用的是四阶龙格的是四阶龙格-库塔公式，四阶龙格库塔公式，四阶龙格-库塔公式的局部库塔公式的局部截断误差界为截断误差界为O(h5), 是四阶精度的计算公式。是四阶精度的计算公式。5.2 常

10、微分方程的初值问题求解常微分方程的初值问题求解)2,()2,2(),()4(6) 1 (2131213211hkhkyhxfkkhyhxfkyxfkkkkhyynnnnnnnn2312131132,3231,31,34 )2(hkyhxfkhkyhxfkyxfkkkhyynnnnnnnn23121321143,4321,21,4329 ) 3(hkyhxfkhkyhxfkyxfkkkkhyynnnnnnnn5.2 常微分方程的初值问题求解常微分方程的初值问题求解2 三阶龙格三阶龙格-库塔库塔3、四阶龙格、四阶龙格库塔公式库塔公式 342312143211,21,2121,21,226 ) 1

11、(hkyhxfkhkyhxfkhkyhxfkyxfkkkkkhyynnnnnnnnnn321421312143211,31,3231,31,338 )2(hkhkhkyhxfkhkhkyhxfkhkyhxfkyxfkkkkkhyynnnnnnnnnn5.2 常微分方程的初值问题求解常微分方程的初值问题求解龙格龙格-库塔方法的调用格式：库塔方法的调用格式：t, y=ode45(ODEfun, tspan, y0, options, p1)自定义的函数自定义的函数求解器求解器积分限或离积分限或离散的向量散的向量初值初值Options：积分参数值，不选为默认值：积分参数值，不选为默认值p1: 传递给

12、传递给ODEfun的参数的参数输出：输出： t: 时间变量，因变量，列向量时间变量，因变量，列向量 y: 状态变量，求解状态变量，求解5.2 常微分方程的初值问题求解常微分方程的初值问题求解例例5-1 求解初值问题：求解初值问题： 10y1x0yx2yy计算结果：计算结果： ode45 ode23 精确值精确值 x=1 1.73205081676653 1.73215487941780 1.73205085.2 常微分方程的初值问题求解常微分方程的初值问题求解求解程序求解程序”xODE.m” 龙格龙格-库塔方法库塔方法步长的选择步长的选择 怎样选取合适的步长，这在实际计算中是很重怎样选取合适的

13、步长，这在实际计算中是很重要的。由于步长愈小，每步计算的截断误差就愈小；要的。由于步长愈小，每步计算的截断误差就愈小；但在一定的求解范围内，需要完成的步数就愈多，但在一定的求解范围内，需要完成的步数就愈多，这不但引起计算量的增加，而且计算步数的增加往这不但引起计算量的增加，而且计算步数的增加往往造成舍入误差的严重积累，反而降低了计算精度。往造成舍入误差的严重积累，反而降低了计算精度。上面介绍的龙格上面介绍的龙格-库塔方法是对定步长（即步长库塔方法是对定步长（即步长h为为常数）而言的，但为了保证精确度，一种有效的措常数）而言的，但为了保证精确度，一种有效的措施是在计算过程中自动进行步长调整，即变

14、步长技施是在计算过程中自动进行步长调整，即变步长技巧。巧。 5.2 常微分方程的初值问题求解常微分方程的初值问题求解 下面以四阶龙格下面以四阶龙格-库塔方法为例，说明如何自动选择步库塔方法为例，说明如何自动选择步长，使计算结果满足给定精度的要求。长，使计算结果满足给定精度的要求。 设从节点设从节点xn出发，先以出发，先以h为步长，利用四阶龙格为步长，利用四阶龙格-库塔库塔公式方法经过一步计算得公式方法经过一步计算得y(xn+1)的近似值，记为的近似值，记为 ，由，由于公式的局部截断误差是于公式的局部截断误差是y(h5)，故有，故有 当当h不大时，不大时，c可近似地看作常数。然后将步长可近似地看

15、作常数。然后将步长h对折，对折，即取即取h/2为步长，从出发经过两步计算求为步长，从出发经过两步计算求y(xn+1)的近似值，的近似值，记为记为 ，每一步计算的局部截断误差为，每一步计算的局部截断误差为c(h/2)5，于是就，于是就有有 )(1hny5)h(1n1nchy)x(y)2/(1hny5)2/(11)2/()(hcyxyhnn5.2 常微分方程的初值问题求解常微分方程的初值问题求解 把它与上式相比，可得把它与上式相比，可得 经整理可得经整理可得 这表明以这表明以 作为作为y(xn+1) 的近似值，其误差可用先后两次计算的近似值，其误差可用先后两次计算结果之差来表示，因而，只需考察结果

16、之差来表示，因而，只需考察 是否是否成立。若成立，则可将成立。若成立，则可将 作为作为y(xn+1)的近似值；若不成立，的近似值；若不成立，则将步长再次对折进行计算，直到不等式成立为止，并取最后则将步长再次对折进行计算，直到不等式成立为止，并取最后的的 作为计算结果。以上方法就是计算过程中自动选择步长作为计算结果。以上方法就是计算过程中自动选择步长的方法，也称为变步长方法。从表面上看，为了选择适当的步的方法，也称为变步长方法。从表面上看，为了选择适当的步长，每一步的计算量增加了，但从总体考虑，工作量往往还是长，每一步的计算量增加了，但从总体考虑，工作量往往还是减减少的。 161)()()(11

17、)2/(11hnnhnnyxyyxy)()()(11151)2/(11hnnhnnyxyyxy)2/(1hny)2/(11)(hnnyxy)2/(1hny)2/(1hny5.2 常微分方程的初值问题求解常微分方程的初值问题求解龙格龙格-库塔方法的主要优点是计算精确度较高，能满足通常库塔方法的主要优点是计算精确度较高，能满足通常的计算要求；容易编制程序；每次计算的计算要求；容易编制程序；每次计算y(xn+1) ，只用到前，只用到前一步的计算结果一步的计算结果yn，因此，在已知初始值，因此，在已知初始值y0的条件下，就可的条件下，就可自动地进行计算，是单步法，而且计算过程中可随时改变自动地进行计算

18、，是单步法，而且计算过程中可随时改变步长。它的缺点是每前进一步需要计算多次步长。它的缺点是每前进一步需要计算多次f(x,y)的值，因的值，因此，计算工作量较大，且其截断误差难以估计。此，计算工作量较大，且其截断误差难以估计。 在实际应用上，一般当要求更高精确度时，采用的办在实际应用上，一般当要求更高精确度时，采用的办法是缩小步长，而不是采用更高阶的公式，因为高阶公式法是缩小步长，而不是采用更高阶的公式，因为高阶公式的计算太复杂，一般选用标准四阶龙格的计算太复杂，一般选用标准四阶龙格-库塔方法即可。库塔方法即可。 5.2 常微分方程的初值问题求解常微分方程的初值问题求解下面用一个例子说明：下面用

19、一个例子说明： 用标准四阶龙格用标准四阶龙格- -库塔方法计算得库塔方法计算得 由积分法算得由积分法算得y(2)=2.23607,y(2)=2.23607,当当h=0.5h=0.5时相差时相差0.000130.00013，而，而h=0.25h=0.25误误差小于差小于0.0000010.000001，但当，但当 h=0.125h=0.125时虽然足够精确但计算次数却比时虽然足够精确但计算次数却比h=0.25h=0.25多了一倍。因此合理选择步长既能保证精度又能减少计算量多了一倍。因此合理选择步长既能保证精度又能减少计算量。 1)0( yyxdxdy5.2 常微分方程的初值问题求解常微分方程的初

20、值问题求解 边值问题求解之一边值问题求解之一-打靶算法打靶算法5.3 常微分方程的边值问题求解常微分方程的边值问题求解二阶微分方程的边值问题的二阶微分方程的边值问题的 求解：求解：在区间在区间a;b上求解该方程的上求解该方程的 解，且在这两个边界点上满足：解，且在这两个边界点上满足：上式即为边界条件上式即为边界条件 ),( yyxFxy bbybyaayaybbaa5.3 常微分方程的边值问题求解常微分方程的边值问题求解对于线性微分方程对于线性微分方程其中其中 p(x), q(x),f(x)均为给定的函数，边界条件均为给定的函数，边界条件可简化为：可简化为： y(a)= a y(b)= b x

21、fxyxqxyxpxy bbybyaayaybbaa5.3 常微分方程的边值问题求解常微分方程的边值问题求解打靶法的步骤：打靶法的步骤：（1）求出下面方程初值问题的数值解）求出下面方程初值问题的数值解 y1(b)(2)求出下面方程初值问题的数值解求出下面方程初值问题的数值解 y1(b)(3)求出下面方程初值问题的数值解求出下面方程初值问题的数值解 yp(b) 0, 1, 01111 1ayayxyxqxyxpxy 1, 0, 02222 2ayayxyxqxyxpxy 1, 0, 0 ayayxyxqxyxpxyppppp（4） 若若y2(b) =0,则计算则计算 bybybympab215.

22、3 常微分方程的边值问题求解常微分方程的边值问题求解 (5) 计算下面初值问题的计算下面初值问题的 数值解，则数值解，则y(x)即为原边值问即为原边值问题的数值解题的数值解 mayayxfxyxqxyxpxya11 ,),(得出一阶微分方程组模型：得出一阶微分方程组模型： 取变量取变量 x1=y, x2=y得到：得到： xfxxpxxqdxdxxdxdx21221,function t, y=shooting(f1, f2, tspan, x0f, varargin)t0=tspan(1); tfinal=tspan(2);ga=x0f(1); gb=x0f(2);t, y1=ode45(f1

23、, tspan, 1; 0, varargin);t, y2=ode45(f1, tspan, 0; 1, varargin);t, yp=ode45(f2, tspan, 0; 0, varargin);m=(gb-ga*y1(end, 1)-yp(end,1)/y2(end, 1);t, y=ode45(f2, tspan, ga; m, varargin);将上述文件存为将上述文件存为“shooting.m”5.3 常微分方程的边值问题求解常微分方程的边值问题求解编写求解通用函数编写求解通用函数5.3 常微分方程的边值问题求解常微分方程的边值问题求解其中：其中： tspan为初始和终止计

24、算时间构成的向量为初始和终止计算时间构成的向量 为边界值为边界值,0bafx 12111221xfxxpxxqdxdxxdxdx 2111221xxpxxqdxdxxdxdx辅助函数辅助函数f1( )辅助函数辅助函数f2( )【例例5-25-2】function xdot=c7fun1(t, x)xdot=x(2); -2*x(1)+3*x(2);function xdot=c7fun2(t, x)xdot=x(2); t-2*x(1)+3*x(2);编写两个函数，分别存为编写两个函数，分别存为m文件文件5.3 常微分方程的边值问题求解常微分方程的边值问题求解 段方程的数值解在求1 , 021

25、, 10,23 yyxyyyC7fun1.mC7fun2.m求解：求解： t, y=shooting(c7fun1, c7fun2, 0, 1, 1, 2) ; plot (t,y) ME_5_105.3 常微分方程的边值问题求解常微分方程的边值问题求解边值问题求解之二边值问题求解之二- -有限差分算法有限差分算法5.3 常微分方程的边值问题求解常微分方程的边值问题求解 有限差分法是用差分格式近似替代有限差分法是用差分格式近似替代ODE方程中的方程中的导数项，从而将导数项，从而将ODE离散为线性或非线性方程组而求离散为线性或非线性方程组而求解，差分格式在偏微分方程中讲解。解，差分格式在偏微分方

26、程中讲解。122112211112222211nnnnnnnnnbbbbtwvtwvtwyt5.3 常微分方程的边值问题求解常微分方程的边值问题求解其中其中n为中间计算点的个数为中间计算点的个数 i, i=1,2,n-1 为微分方程的数值解为微分方程的数值解 yi+1= i, y1= a和和yn= b其它参数的计算其它参数的计算 hiaxnabhhwbbwbbnixfhbxphwxphvxqhtibnnnaiiiiiiii1,1/,1, 2 , 1,2121,211111122为计算步长其中 Function x,y=fdiff(funcs, tspan, x0f,n)t0=tspan(1);

27、tfinal=tspan(2);ga=x0f(1);gb=x0f(2);for i=1:n x(i)=t0+h*(i-1);endx0=x(1:n-1);t=-2+h2*feval(funcs,x0,2);tmp= feval(funcs,x0,1);v=1+h*tem/2;w=1-h*tem/2;b= h2*feval(funcs,x0,3);b(1)=b(1)-w(1)*ga;b(n-1)=b(n-1)-v(n-1)*gb;b=b;A=diag(t);for i=1:n-2 A(i,i+1)=v(i); A(i+1,i)=w(i+1);endy=inv(A)*b;x=x, tfinal;

28、y=ga; y; gb;5.3 常微分方程的边值问题求解常微分方程的边值问题求解边值问题求解之三边值问题求解之三- -内部函数内部函数求解函数：求解函数：bvpinit( ), bvp4c( ), 和和 deval( )sol = bvp4c(odefun, bcfun, solinit)Odefun: A function that evaluates the differential equations f(x,y). It can have the formdydx = odefun(x,y)dydx = odefun(x,y,p1,p2,.)dydx = odefun(x,y,para

29、meters)dydx = odefun(x,y,parameters,p1,p2,.)5.3 常微分方程的边值问题求解常微分方程的边值问题求解Bcfun: A function that computes the residual in the boundary conditions bc(y(a), y(b). It can have the formres = bcfun(ya,yb)res = bcfun(ya,yb,p1,p2,.)res = bcfun(ya,yb,parameters)res = bcfun(ya,yb,parameters,p1,p2,.) 10y0y00y0y

30、x2cosq2y 5.3 常微分方程的边值问题求解常微分方程的边值问题求解求解方程，其中求解方程，其中q=5ME_5_185.4 5.4 微分方程的解析解方法微分方程的解析解方法边值问题求解之四边值问题求解之四- -符号算法符号算法 y=dsolve(f1,f2,fm) y=dsolve(f1,f2,fm,x) % 指明自变量指明自变量求求 x(t)=x(t)(1-x2(t)的解析解的解析解Syms x t X=dsolve(Dx=x*(1-x2)例例5-8X = 1/(1+exp(-2*t)*C1)(1/2) -1/(1+exp(-2*t)*C1)(1/2)【例例5-55-5】输入信号输入信

31、号syms t yu=exp(-5*t)*cos(2*t+1)+5;uu=5*diff(u,t,2)+4*diff(u,t)+2*u uu =87*exp(-5*t)*cos(2*t+1)+92*exp(-5*t)*sin(2*t+1)+105.4 5.4 微分方程的解析解方法微分方程的解析解方法 512cos5tetut tutututytytytyty24 5)(24503510234 求以下微分方程的通解：求以下微分方程的通解：步骤步骤1 步骤步骤2 Syms t yy=dsolve(D4y+10*D3y+35*D2y+50*Dy+24*y=, . 87*exp(-5*t)*cos(2*

32、t+1)+92*exp(-5*t)*sin(2*t+1)+10)y = -87/52*cos(1)*cos(t)2*exp(-5*t)-87/260*cos(1)*sin(t)*cos(t)*exp(-5*t)+87/104*exp(-5*t)*cos(1)+23/26*exp(-5*t)*sin(1)-87/520*sin(1)*cos(2*t)*exp(-5*t)+23/130*cos(1)*cos(2*t)*exp(-5*t)-23/65*sin(1)*sin(t)*cos(t)*exp(-5*t)+5/12+87/104*sin(1)*sin(2*t)*exp(-5*t)-23/26*

33、cos(1)*sin(2*t)*exp(-5*t)-23/13*sin(1)*cos(t)2*exp(-5*t)+C1*exp(-4*t)+C2*exp(-3*t)+C3*exp(-2*t)+C4*exp(-t)5.4 5.4 微分方程的解析解方法微分方程的解析解方法syms t y y=dsolve(D4y+10*D3y+35*D2y+50*Dy+24*y=, . 87*exp(-5*t)*cos(2*t+1)+92*exp(-5*t)*sin(2*t+1)+10, y(0)=3, . Dy(0)=2, D2y(0)=0, D3y(0)=0)y(t)= -87/52*cos(1)*cos(t

34、)2*exp(-5*t)-87/520*sin(1)*cos(2*t)*exp(-5*t)+23/130*cos(1)*cos(2*t)*exp(-5*t)-23/65*sin(1)*sin(t)*cos(t)*exp(-5*t)-87/260*cos(1)*sin(t)*cos(t)*exp(-5*t)+87/104*sin(1)*sin(2*t)*exp(-5*t)-23/26*cos(1)*sin(2*t)*exp(-5*t)-23/13*sin(1)*cos(t)2*exp(-5*t)+5/12+87/104*exp(-5*t)*cos(1)+23/26*exp(-5*t)*sin(1

35、)+(-271/30*cos(1)+41/15*sin(1)-25/4)*exp(-4*t)+(179/8*cos(1)+5/8*sin(1)+73/3)*exp(-3*t)+(-445/26*cos(1)-51/13*sin(1)-69/2)*exp(-2*t)+(133/30*cos(1)+97/60*sin(1)+19)*exp(-t)5.4 5.4 微分方程的解析解方法微分方程的解析解方法假设假设 y(0)=3, y(0)=2, y(0)=y(0)=0 x,y=dsolve(D2x+2*Dx=x+2*y-exp(-t), Dy=4*x+3*y+4*exp(-t)x = -1/6*(-2

36、4*C1*exp(-t)+6*C1*exp(1+6(1/2)*t)-3*C1*6(1/2)*exp(1+6(1/2)*t)+6*C1*exp(-(-1+6(1/2)*t)+3*C1*6(1/2)*exp(-(-1+6(1/2)*t)+5*C2*6(1/2)*exp(-(-1+6(1/2)*t)+12*C2*exp(-(-1+6(1/2)*t)-5*C2*6(1/2)*exp(1+6(1/2)*t)+12*C2*exp(1+6(1/2)*t)-24*C2*exp(-t)-2*C3*6(1/2)*exp(-(-1+6(1/2)*t)-6*C3*exp(-(-1+6(1/2)*t)+2*C3*6(1

37、/2)*exp(1+6(1/2)*t)-6*C3*exp(1+6(1/2)*t)+12*C3*exp(-t)-150*exp(-t)+72*exp(-t)*t)/(-2+6(1/2)/(2+6(1/2)5.4 5.4 微分方程的解析解方法微分方程的解析解方法 tt e4ty3tx4tyety2txtx2xx求解求解y = 2/3*(3*C1*exp(1+6(1/2)*t)-6*C1*exp(-t)+3*C1*exp(-(-1+6(1/2)*t)+C2*6(1/2)*exp(-(-1+6(1/2)*t)+3*C2*exp(-(-1+6(1/2)*t)-C2*6(1/2)*exp(1+6(1/2)

38、*t)+3*C2*exp(1+6(1/2)*t)-6*C2*exp(-t)+3*C3*exp(-t)-C3*6(1/2)*exp(-(-1+6(1/2)*t)+C3*6(1/2)*exp(1+6(1/2)*t)+18*exp(-t)*t-36*exp(-t)/(-2+6(1/2)/(2+6(1/2)5.4 5.4 微分方程的解析解方法微分方程的解析解方法5.5 高阶常微分方程求解高阶常微分方程求解一一 单个高阶常微分方程处理方法单个高阶常微分方程处理方法一个高阶微分方程的一般形式为：一个高阶微分方程的一般形式为：输出变量输出变量y(t)的各阶导数初始值为：的各阶导数初始值为：选出一组状态变量：

39、选出一组状态变量： 1,nnyyytfy 0,0,01nyyy121,nnyxyxyx5.5高阶常微分方程求解高阶常微分方程求解原高阶常微分方程变换为一阶微分方程组原高阶常微分方程变换为一阶微分方程组初值为：初值为：n21nn1nn3 22211x,x,x, tfyxyxxyxyxxyxyx 00,00,00121nnyxyxyx【例例5-35-3】5.5 高阶常微分方程求解高阶常微分方程求解function y=vdp_eq(t, x, flag, mu)y=x(2) ; -mu*(x(1)2-1)*x(2)-x(1);建立函数建立函数“vdp_eq.m”x0=-0.2; -0.7; t_f

40、inal=20; mu=1; t1, y1=ode45(vdp_eq, 0, t_final, x0, , mu); mu=2; t2, y2=ode45(vdp_eq, 0, t_final, x0, , mu); plot(t1, y1, t2, y2, :) figure; plot(y1( :, 1), y1( :, 2), y2( :, 1), y2( :, 2), :)已知：已知：y(0)=-0.2, y(0)=-0.7, 用数值方法求用数值方法求Van der pol 方程的解：方程的解：01 2yyyyME_5_16.m5.5 高阶常微分方程求解高阶常微分方程求解5.5 高阶常

41、微分方程求解高阶常微分方程求解二二 高阶常微分方程组高阶常微分方程组5.5 高阶常微分方程求解高阶常微分方程求解 1111,nmnnmmyyyxxxtgyyyyxxxtfx121121,nymmmmmxxyxyxxxxxxxnmnmmmmmnmmxxxtgxxxxxxxxtfxxxxx,213221213221选择一组状态变量：选择一组状态变量：原方程变为：原方程变为：【例例5-35-3】5.5 高阶常微分方程求解高阶常微分方程求解Apollo卫星的运动轨迹卫星的运动轨迹(x,y) 满足一下方程：满足一下方程： 0493575. 10, 00, 2 . 1, 2 . 10,-145.82/1*2*2222221*3231 32*31 yyxxyxryxrryryyxyrxrxxyx初值，其中求上述方程的解并绘制求上述方程的解并绘制Apollo位置的位置的(x,y)