1、 -293- 第十五章 常微分方程的解法 建立微分方程只是解决问题的第一步，通常需要求出方程的解来说明实际现象，并加以检验。如果能得到解析形式的解固然是便于分析和应用的，但是我们知道，只有线性常系数微分方程，并且自由项是某些特殊类型的函数时，才可以得到这样的解，而绝大多数变系数方程、非线性方程都是所谓“解不出来”的，即使看起来非常简单的方程如22xydxdy+= 。于是对于用微分方程解决实际问题来说，数值解法就是一个十分重要的手段。 1 常微分方程的离散化 下面主要讨论一阶常微分方程的初值问题，其一般形式是 =0)(),(yaybxayxfdxdy（ 1） 在下面的讨论中，我们总假定函数 ),
2、( yxf 连续，且关于 y 满足李普希兹 (Lipschitz)条件，即存在常数 L，使得 |),(),(| yyLyxfyxf 这样，由常微分方程理论知，初值问题 (1)的解必定存在唯一。 所谓数值解法，就是求问题（ 1）的解 )(xy 在若干点 bxxxxaN=L ，使得 , bax ，mRyy 21, ，都有 2121),(),( yyLyxfyxf 那么问题（ 25）在 , ba 上存在唯一解 )(xyy = 。 问题（ 25）与（ 1）形式上完全相同，故对初值问题（ 1）所建立的各种数值解法可全部用于求解问题（ 25） 。 6.2 高阶微分方程的数值解法 高阶微分方程的初值问题可以
3、通过变量代换化为一阶微分方程组初值问题。 设有 m 阶常微分方程初值问题 =)1(0)1()1(00)1()()(,)(,)(),(mmmmyayyayyaybxayyyxfyLL（ 26） 引入新变量)1(21,=mmyyyyyy L ，问题（ 26）就化为一阶微分方程组初值问题 =)1(01)2(011)1(02320121)( ),()( )( )( mmmmmmmmyayyyxfyyayyyyayyyyayyyLMM （ 27） 然后用 6.1 中的数值方法求解问题（ 27） ，就可以得到问题（ 26）的数值解。 最后需要指出的是，在化学工程及自动控制等领域中，所涉及的常微分方程组初值
4、问题常常是所谓的“刚性”问题。具体地说，对一阶线性微分方程组 )(xAydxdy+= （ 28） 其中 ARym, 为 m 阶方阵。若矩阵 A的特征值 ),2,1( miiL= 满足关系 ),2,1(0Re miiL= 则称方程组（ 28）为刚性方程组或 Stiff 方程组，称数 |Re|min/|Re|max11imiimis = 为刚性比。对刚性方程组，用前面所介绍的方法求解，都会遇到本质上的困难，这是由数值方法本身的稳定性限制所决定的。理论上的分析表明，求解刚性问题所选用的数值方法最好是对步长 h不作任何限制。 7 初值问题的 Matlab 解法和符号解 7.1 Matlab 数值解 7
5、.1.1 非刚性常微分方程的解法 Matlab 的工具箱提供了几个解非刚性常微分方程的功能函数，如 ode45， ode23，ode113，其中 ode45 采用四五阶 RK 方法，是解非刚性常微分方程的首选方法， ode23采用二三阶 RK 方法， ode113 采用的是多步法，效率一般比 ode45 高。 Matlab 的工具箱中没有 Euler 方法的功能函数。 （ I）对简单的一阶方程的初值问题 =00)(),(yxyyxfy改进的 Euler 公式为 +=+=+=+)(21),(),(1 qpnpnnqnnnpyyyyhxhfyyyxhfyy我们自己编写改进的 Euler 方法函数
6、eulerpro.m 如下： function x,y=eulerpro(fun,x0,xfinal,y0,n); if nargin 是一参数。 解 （i）化成常微分方程组。设 ,21yyyy = ，则有 =1221221)1(yyyyyy（ ii）书写 M 文件（对于 1= ） vdp1.m: function dy=vdp1(t,y); dy=y(2);(1-y(1)2)*y(2)-y(1); （ iii）调用 Matlab函数。对于初值 0)0(,2)0( = yy ，解为 T,Y=ode45(vdp1,0 20,2;0); （ iv）观察结果。利用图形输出解的结果： plot(T,Y
7、(:,1),-,T,Y(:,2),-) title(Solution of van der Pol Equation,mu=1); xlabel(time t); ylabel(solution y); legend(y1,y2); 0 2 4 6 8 10 12 14 16 18 20-3-2-10123Solution of van der Pol Equation,mu=1time tsolutionyy1y2图 1 例 5 van der Pol 方程， 1000= （刚性） 解 (i)书写M 文件vdp1000.m: function dy=vdp1000(t,y); dy=y(2)
8、;1000*(1-y(1)2)*y(2)-y(1); （ii）观察结果 t,y=ode15s(vdp1000,0 3000,2;0); plot(t,y(:,1),o) title(Solution of van der Pol Equation,mu=1000); xlabel(time t); ylabel(solution y(:,1); 7.2 常微分方程的解析解 在 Matlab 中，符号运算工具箱提供了功能强大的求解常微分方程的符号运算命令dsolve。常微分方程在 Matlab 中按如下规定重新表达： -303- 符号 D 表示对变量的求导。 Dy 表示对变量 y 求一阶导数，当
9、需要求变量的 n 阶导数时，用 Dn 表示， D4y 表示对变量 y 求 4 阶导数。 由此，常微分方程 yyy =+ 2 在 Matlab 中，将写成 D2y+2*Dy=y。 7.2.1 求解常微分方程的通解 无初边值条件的常微分方程的解就是该方程的通解。其使用格式为： dsolve(diff_equation) dsolve( diff_equation， var) 式中 diff_equation 为待解的常微分方程， 第 1种格式将以变量 t为自变量进行求解，第 2 种格式则需定义自变量 var。 例 6 试解常微分方程 0)2(2=+ yyxyx 解 编写程序如下： syms x y
10、 diff_equ=x2+y+(x-2*y)*Dy=0; dsolve(diff_equ,x) 7.2.2 求解常微分方程的初边值问题 求解带有初边值条件的常微分方程的使用格式为： dsolve(diff_equation， condition1， condition2， ， var) 其中 condition1， condition2， 即为微分方程的初边值条件。 例 7 试求微分方程 xyy = ， 4)2(,7)1(,8)1( = yyy 的解。 解 编写程序如下： y=dsolve(D3y-D2y=x,y(1)=8,Dy(1)=7,D2y(2)=4,x) 7.2.3 求解常微分方程组
11、求解常微分方程组的命令格式为： dsolve(diff_equ1， diff_equ2， ， var) dsolve(diff_equ1， diff_equ2， ， condition1， condition2， ， var) 第 1 种格式用于求解方程组的通解，第 2 种格式可以加上初边值条件，用于具体求解。 例 8 试求常微分方程组： =+=+xfgxgfcossin3的通解和在初边值条件为 1)5(,3)3(,0)2( = gff 的解。 解 编写程序如下： clc,clear equ1=D2f+3*g=sin(x); equ2=Dg+Df=cos(x); general_f,gener
12、al_g=dsolve(equ1,equ2,x) f,g=dsolve(equ1,equ2,Df(2)=0,f(3)=3,g(5)=1,x) 7.2.4 求解线性常微分方程组 （i）一阶齐次线性微分方程组 -304- AXX = ，=nxxX M1，=nnnnaaaaALMMML1111这里的表示对 t求导数。Ate 是它的基解矩阵。 AXX = ，00)( XtX = 的解为0)(0)( XetXttA = 。 例 9 试解初值问题 XX=200120312 ，=121)0(X 解 编写程序如下： syms t a=2,1,3;0,2,-1;0,0,2; x0=1;2;1; x=expm(a
13、*t)*x0 （ii）非齐次线性方程组 由参数变易法可求得初值问题 )( tfAXX += ，00)( XtX = 的解为 +=ttstAttAdssfeXetX00)()()(0)(. 例 10 试解初值问题 +=teXXt2cos00123212001 ，=110)0(X 。 解 编写程序如下： clc,clear syms t s a=1,0,0;2,1,-2;3,2,1;ft=0;0;exp(t)*cos(2*t); x0=0;1;1; x=expm(a*t)*x0+int(expm(a*(t-s)*subs(ft,s),s,0,t); x=simple(x) 8 边值问题的 Matl
14、ab解法 Matlab中用 bvp4c命令求解常微分方程的两点边值问题，微分方程的标准形式为 ),( yxfy = ， 0)(),( =byaybc 或者是 ),( pyxfy = ， 0),(),( =pbyaybc 其中 p 是有关的参数。这里 fy, 可以为向量函数，求解的区间为 , ba ， bc为边界条件。 一般地说，边值问题在计算上比初值问题困难得多，特别地，由于边值问题的解可能是多值的， bvp4c需要提供猜测的初始值。下面我们首先给出一个简单的例子。 -305- 例 11 考察描述在水平面上一个小水滴横截面形状的标量方程 0)(1)(1()(2/3222=+ xhdxdxhxh
15、dxd， 0)1()1( = hh 这里 )(xh 表示 x处水滴的高度。设 )()(1xhxy = ，dxxdhxy)()(2= ，把上述微分方程写成两个一阶微分方程组 )()(21xyxydxd= 2/32212)(1)(1)()( xyxyxydxd+= 上述微分方程组可以由如下函数表示 function yprime=drop(x,y); yprime=y(2);(y(1)-1)*(1+y(2)2)(3/2); 边界条件通过残差函数指定，边界条件通过如下函数表示 function res=dropbc(ya,yb); res=ya(1);yb(1); 我们使用211)( xxy = 和
16、 )11.0/()(22xxxy += 作为初始猜测解，由如下函数定义 function yinit=dropinit(x); yinit=sqrt(1-x.2);-x./(0.1+sqrt(1-x.2); 利用如下的程序就可以求微分方程的边值问题并画出图 2。 solinit=bvpinit(linspace(-1,1,20),dropinit); sol=bvp4c(drop,dropbc,solinit); fill(sol.x,sol.y(1,:),0.7,0.7,0.7) axis(-1,1,0,1) xlabel(x,FontSize,12) ylabel(h,Rotation,0
17、,FontSize,12) -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 100.10.20.30.40.50.60.70.80.91xh图 2 这里调用函数 bvpinit，计算区间 1,1 上等间距的 20个点的数据，然后调用函数bvp4c，得到数值解的结构 sol，填充命令 fill填充1yx 平面上的解曲线。 一般地， bvp4c的调用格式如下： -306- sol=bvp4c(odefun,bcfun,solinit,options,p1,p2, ); 函数 odefun的格式为 yprime=odefun(x,y,p1,p2, ); 函数 bcf
18、un的格式为 res=bcfun(ya,yb, p1,p2, ); 初始猜测结构 solinit有两个域： solinit.x提供初始猜测的 x值，排列顺序从左到右排列，其中 solinit.x(1)和 solinit.x(end)分别为 a和 b 。对应地，solinit.y(:,i)给出点 solinit.x(i)处初始猜测解。 输出参数 sol是包含数值解的一个结构，其中 sol.x给出了计算数值解的 x点，sol.x(i)处的数值解由 sol.y(:,i)给出，类似地， sol.x(i)处数值解的一阶导数值由 sol.yp(:,i)给出。 我们可以把上面的所有函数都放在一个文件中，程序
19、如下： function sol=example11; solinit=bvpinit(linspace(-1,1,20),dropinit); sol=bvp4c(drop,dropbc,solinit); fill(sol.x,sol.y(1,:),0.7,0.7,0.7) axis(-1,1,0,1) xlabel(x,FontSize,12) ylabel(h,Rotation,0,FontSize,12) function yprime=drop(x,y); yprime=y(2);(y(1)-1)*(1+y(2)2)(3/2); function res=dropbc(ya,yb)
20、; res=ya(1);yb(1); function yinit=dropinit(x); yinit=sqrt(1-x.2);-x./(0.1+sqrt(1-x.2); 例 12 描述 0=x 处固定， 1=x 处有弹性支持，沿着 x轴平衡位置以均匀角速度旋转的绳的位移方程 0)()(22=+ xyxydxd 具有边界条件 0)0( =y ， 1)(0=xxydxd， 0)()(1=+=xxydxdxy 这个边值问题是一个特征问题，我们必须找到参数 的值使得方程的解存在。如果我们提供了参数 的猜测值和对应解的猜测值，我们也可以利用函数 bvp4c求解特征问题。上述微分方程可以写成下面的微分
21、方程组： )()(21xyxydxd= )()(12xyxydxd= 编写程序如下（下面的所有程序放在一个文件中）： function sol=skiprun; -307- solinit=bvpinit(linspace(0,1,10),skipinit,5); sol=bvp4c(skip,skipbc,solinit); plot(sol.x,sol.y(1,:),-,sol.x,sol.yp(1,:),-,LineWidth,2) xlabel(x,FontSize,12) legend(y_1,y_2,0) function yprime=skip(x,y,mu); yprime=y
22、(2);-mu*y(1); function res=skipbc(ya,yb,mu); res=ya(1);ya(2)-1;yb(1)+yb(2); function yinit=skipinit(x); yinit=sin(x);cos(x); 图 3给出上述边值问题解的图象。 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1-0.500.51xy1y2图 3 例 13 微分方程组为 vuwuu /)(5.0 = )(5.0 uwv = zuwwyww /)(5.0)(10009.0( = )(5.0 uwz = )(100 wyy = 边界条件为 1)0()
23、0()0( = wvu ， 10)0( =z ； )1()1( yw = 。 我们使用如下猜测解 1)( =xu 1)( =xv 191.85.4)(2+= xxxw 10)( =xz 91.095.4)(2+= xxxy 我们编写如下程序（所有程序放在一个文件中）： function sol=example13 solinit=bvpinit(linspace(0,1,5),exinit); sol=bvp4c(exode,exbc,solinit); plot(sol.x,sol.y) -308- function yprime=exode(x,y) yprime=0.5*y(1)*(y(
24、3)-y(1)/y(2) -0.5*(y(3)-y(1) (0.9-1000*(y(3)-y(5)-0.5*y(3)*(y(3)-y(1)/y(4) 0.5*(y(3)-y(1) 100*(y(3)-y(5); function res=exbc(ya,yb) res=ya(1)-1;ya(2)-1;ya(3)-1;ya(4)+10;yb(3)-yb(5); function yinit=exinit(x) yinit=1;1;-4.5*x2+8.91*x+1;-10;-4.5*x2+9*x+0.91; 习 题 十 五 1. 用欧拉方法和龙格库塔方法求微分方程数值解，画出解的图形，对结果进行分
25、析比较。 （ i） 22,22,0)(222=+ yyynxxyyx （ Bessel 方程，令)21=n ，精确解xxy2sin= 。 （ ii） 0)0(,1)0(,0cos =+ yyxyy ，幂级数解 L+=8642!855!69!42!211 xxxxy 2. 一只小船渡过宽为 d 的河流，目标是起点 A正对着的另一岸 B 点。已知河水流速1v 与船在静水中的速度2v 之比为 k 。 （ i）建立小船航线的方程，求其解析解。 （ ii）设 100=d m， 11=v m/s， 22=v m/s，用数值解法求渡河所需时间、任意时刻小船的位置及航行曲线，作图，并与解析解比较。 3. Lorenz方程是一个三阶的非线性系统，它是由描述大气动力系统的Navier-Stokes偏微分方程演化而来的。自由系统如下： +=xyzzxzyxyxyx& )(当系统参数 , 在一定范围内，系统就出现混沌，如 10= ， 28= ，3/8= 时，出现混沌现象。求在初始条件 17,13,5)0(),0(),0( =zyx 时，方程组的数值解，并画出解的图形。
