第十五章 常微分方程的解法
建立微分方程只是解决问题的第一步，通常需要求出方程的解来说明实际现象，并加以检验。如果能得到解析形式的解固然是便于分析和应用的，但是我们知道，只有线性常系数微分方程，并且自由项是某些特殊类型的函数时，才可以肯定得到这样的解，而绝大多数变系数方程、非线性方程都是所谓“解不出来”的，即使看起来非常简单的方程如
dy
=y2+x2，于是对于用微分方程解决实际问题来说，数值解法就是一个十dx
分重要的手段。
§1 常微分方程的离散化
下面主要讨论一阶常微分方程的初值问题，其一般形式是
⎧dy
⎪=f(x,y)
⎨dx
⎪⎩y(a)=y0
a≤x≤b
(1)
在下面的讨论中，我们总假定函数f(x,y)连续，且关于y满足李普希兹(Lipschitz)条
件，即存在常数L，使得
|f(x,y)−f(x,|≤L|y−|
这样，由常微分方程理论知，初值问题(1)的解必定存在唯一。
所谓数值解法，就是求问题（1）的解y(x)在若干点
a=x0
处的近似值yn(n=1,2,L,N)的方法，yn(n=1,2,L,N)称为问题（1）的数值解，
hn=xn+1−xn称为由xn到xn+1的步长。今后如无特别说明，我们总取步长为常量h。
建立数值解法，首先要将微分方程离散化，一般采用以下几种方法： （i）用差商近似导数
y(xn+1)−y(xn)
代替y'(xn)代入（1）中的微分方程，则得
h
y(xn+1)−y(xn)
≈f(xn,y(xn))(n=0,1,L)
h
若用向前差商化简得
y(xn+1)≈y(xn)+hf(xn,y(xn))
如果用y(xn)的近似值yn代入上式右端，所得结果作为y(xn+1)的近似值，记为yn+1，
则有
yn+1=yn+hf(xn,yn)
(n=0,1,L) （2）
这样，问题（1）的近似解可通过求解下述问题 ⎨
⎧yn+1=yn+hf(xn,yn)(n=0,1,L)
（3）
⎩y0=y(a)
得到，按式（3）由初值y0可逐次算出y1,y2,L。式（3）是个离散化的问题，称为差分方程初值问题。
-179-
需要说明的是，用不同的差商近似导数，将得到不同的计算公式。 （ii）用数值积分方法
将问题（1）的解表成积分形式，用数值积分方法离散化。例如，对微分方程两端积分，得
y(xn+1)−y(xn)=∫
xn+1
xn
f(x,y(x))dx(n=0,1,L) (4)
右边的积分用矩形公式或梯形公式计算。
（iii）Taylor多项式近似
将函数y(x)在xn处展开，取一次Taylor多项式近似，则得
y(xn+1)≈y(xn)+hy'(xn)=y(xn)+hf(xn,y(xn))
再将y(xn)的近似值yn代入上式右端，所得结果作为y(xn+1)的近似值yn+1，得到离
散化的计算公式
yn+1=yn+hf(xn,yn)
以上三种方法都是将微分方程离散化的常用方法，每一类方法又可导出不同形式的计算公式。其中的Taylor展开法，不仅可以得到求数值解的公式，而且容易估计截断误差。
§2 欧拉（Euler）方法
2.1 Euler方法
Euler 方法就是用差分方程初值问题（3）的解来近似微分方程初值问题（1）的解，即由公式（3）依次算出y(xn)的近似值yn(n=1,2,L)。这组公式求问题（1）的数值解称为向前Euler公式。
如果在微分方程离散化时，用向后差商代替导数，即y'(xn+1)≈则得计算公式
⎨
y(xn+1)−y(xn)
，
h
⎧yn+1=yn+hf(xn+1,yn+1)(n=0,1,L)
（5）
=yy(a)⎩0
用这组公式求问题（1）的数值解称为向后Euler公式。
向后Euler法与Euler法形式上相似，但实际计算时却复杂得多。向前Euler公式是显式的，可直接求解。向后Euler公式的右端含有yn+1，因此是隐式公式，一般要用迭代法求解，迭代公式通常为
(0)⎧⎪yn+1=yn+hf(xn,yn) ⎨(k+1)(k)⎪⎩yn+1=yn+hf(xn+1,yn+1)
(k=0,1,2,L)
（6）
2.2 Euler方法的误差估计
对于向前Euler公式（3）我们看到，当n=1,2,L时公式右端的yn都是近似的，所以用它计算的yn+1会有累积误差，分析累积误差比较复杂，这里先讨论比较简单的所谓局部截断误差。
假定用（3）式时右端的yn没有误差，即yn=y(xn)，那么由此算出
yn+1=y(xn)+hf(xn,y(xn)) （7）
-180-
局部截断误差指的是，按（7）式计算由xn到xn+1这一步的计算值yn+1与精确值y(xn+1)之差y(xn+1)−yn+1。为了估计它，由Taylor展开得到的精确值y(xn+1)是
h2
y(xn+1)=y(xn)+hy'(xn)+y''(xn)+O(h3) （8）
2
（7）、（8）两式相减（注意到y'=f(x,y)）得
y(xn+1)−yn+1
2
h2
=y''(xn)+O(h3)≈O(h2) （9） 2
p+1
即局部截断误差是h阶的，而数值算法的精度定义为：
若一种算法的局部截断误差为O(h)，则称该算法具有p阶精度。
显然p越大，方法的精度越高。式（9）说明，向前Euler方法是一阶方法，因此它的精度不高。
§3 改进的Euler方法
3.1 梯形公式
利用数值积分方法将微分方程离散化时，若用梯形公式计算式(4)中之右端积分，即
h
f(x,y(x))dx≈f(xn,y(xn))+f(xn+1,y(xn+1))] ∫xn
2
并用yn,yn+1代替y(xn),y(xn+1)，则得计算公式
h
yn+1=yn+[f(xn,yn)+f(xn+1,yn+1)]
2
xn+1
这就是求解初值问题（1）的梯形公式。
直观上容易看出，用梯形公式计算数值积分要比矩形公式好。梯形公式为二阶方法。 梯形公式也是隐式格式，一般需用迭代法求解，迭代公式为
(0)⎧yn
+1=yn+hf(xn,yn)⎪
h⎪(k+1)(k)
yy[f(xn,yn)+f(xn+1,yn=+⎨n+1n+1)] (10)
2⎪
(k=0,1,2,L)⎪⎩
由于函数f(x,y)关于y满足Lipschitz条件，容易看出
hL(k)(k+1)(k)(k−1)
|yn−y|≤|yn+1−yn+1n+1+1|
2
hL
如果实际计算时精度要求不太高，用公式（10）求解时，每步可以只迭代一次，由此导出一种新的方法—改进Euler法。
3.2 改进Euler法
按式（5）计算问题（1）的数值解时，如果每步只迭代一次，相当于将Euler公式与梯形公式结合使用：先用Euler公式求yn+1的一个初步近似值n+1，称为预测值，然后用梯形公式校正求得近似值yn+1，即
-181-
预测⎧n+1=yn+hf(xn,yn)
⎪
⎨ （11） h
yy[f(x,y)f(x,) ] 校正=++n+1nnnn+1n+1⎪2⎩
式（11）称为由Euler公式和梯形公式得到的预测—校正系统，也叫改进Euler法。
为便于编制程序上机，式（11）常改写成
⎧
⎪yp=yn+hf(xn,yn)⎪
⎨yq=yn+hf(xn+h,yp) （12）
⎪
⎪yn+1=1(yp+yq)⎩2
改进Euler法是二阶方法。
§4 龙格—库塔（Runge—Kutta）方法
回到Euler方法的基本思想—用差商代替导数—上来。实际上，按照微分中值定理应有
y(xn+1)−y(xn)
=y'(xn+θh),0
h
注意到方程y'=f(x,y)就有
y(xn+1)=y(xn)+hf(xn+θh,y(xn+θh)) （13）
不妨记=f(xn+θh,y(xn+θh))，称为区间[xn,xn+1]上的平均斜率。可见给出一种斜率，（13）式就对应地导出一种算法。
向前Euler公式简单地取f(xn,yn)为，精度自然很低。改进的Euler公式可理解为取f(xn,yn)，f(xn+1,n+1)的平均值，其中n+1=yn+hf(xn,yn)，这种处理提高了精度。
如上分析启示我们，在区间[xn,xn+1]内多取几个点，将它们的斜率加权平均作为
，就有可能构造出精度更高的计算公式。这就是龙格—库塔方法的基本思想。
4.1 首先不妨在区间[xn,xn+1]内仍取2个点，仿照（13）式用以下形式试一下
⎧yn+1=yn+h(λ1k1+λ2k2)⎪
（14） ⎨k1=f(xn,yn)
⎪k=f(x+αh,y+βhk),0
nn1⎩2
其中λ1,λ2,α,β为待定系数，看看如何确定它们使（14）式的精度尽量高。为此我们分析局部截断误差y(xn+1)−yn+1，因为yn=y(xn)，所以（14）可以化为
-182-
⎧yn+1=y(xn)+h(λ1k1+λ2k2)⎪
⎪k1=f(xn,y(xn))=y'(xn)⎪
（15） ⎨k2=f(xn+αh,y(xn)+βhk1)
⎪ =f(x,y(x))+αhf(x,y(x))
nnxnn
⎪
2
⎪ βhkf(x,y(x))O(h)++nn1y⎩
其中k2在点(xn,y(xn))作了Taylor展开。（15）式又可表为
yn+1=y(xn)+(λ1+λ2)hy'(xn)+λ2αh2(fx+
注意到
β
ffy)+O(h3) α
h2
y(xn+1)=y(xn)+hy'(xn)+y''(xn)+O(h3)
2
3
中y'=f，y''=fx+ffy，可见为使误差y(xn+1)−yn+1=O(h)，只须令
λ1+λ2=1，λ2α=
1β
，=1 （16） 2α
1
，α=β=1，即为改2
待定系数满足（16）的（15）式称为2阶龙格—库塔公式。由于（16）式有4个未知数而只有3个方程，所以解不唯一。不难发现，若令λ1=λ2=
进的Euler公式。可以证明，在[xn,xn+1]内只取2点的龙格—库塔公式精度最高为2阶。
4.2 4阶龙格—库塔公式
要进一步提高精度，必须取更多的点，如取4点构造如下形式的公式：
⎧yn+1=yn+h(λ1k1+λ2k2+λ3k3+λ4k4)⎪k=f(x,y)
nn1⎪⎪
（17） ⎨k2=f(xn+α1h,yn+β1hk1)
⎪k=f(x+αh,y+βhk+βhk)
nn22132
⎪3⎪⎩k4=f(xn+α3h,yn+β4hk1+β5hk2+β6hk3)
其中待定系数λi,αi,βi共13个，经过与推导2阶龙格—库塔公式类似、但更复杂的计
算，得到使局部误差y(xn+1)−yn+1=O(h)的11个方程。取既满足这些方程、又较简单的一组λi,αi,βi，可得
5
-183-
h⎧
=y(k1+2k2+2k3+k4)+n1⎪6⎪
⎪k1=f(xn,yn)⎪hkh⎪
⎨k2=f(xn+,yn+1 （18）
22⎪
hk2h⎪
=++kf(x,ynn⎪3
22
⎪⎪⎩k4=f(xn+h,yn+hk3)
这就是常用的4阶龙格—库塔方法（简称RK方法）。
§5 线性多步法
以上所介绍的各种数值解法都是单步法，这是因为它们在计算yn+1时，都只用到前一步的值yn，单步法的一般形式是
(n=0,1,L,N−1) （19）
其中ϕ(x,y,h)称为增量函数，例如Euler方法的增量函数为f(x,y)，改进Euler法的
增量函数为
yn+1=yn+hϕ(xn,yn,h)
ϕ(x,y,h)=[f(x,y)+f(x+h,y+hf(x,y))]
如何通过较多地利用前面的已知信息，如yn,yn−1,L,yn−r，来构造高精度的算法计算yn+1，这就是多步法的基本思想。经常使用的是线性多步法。
让我们试验一下r=1，即利用yn,yn−1计算yn+1的情况。 从用数值积分方法离散化方程的（4）式 y(xn+1)−y(xn)=
12
∫
xn+1
xn
f(x,y(x))dx
出发，记f(xn,yn)=fn，f(xn−1,yn−1)=fn−1，式中被积函数f(x,y(x))用二节点
(xn−1,fn−1)，(xn,fn)的插值公式得到（因x≥xn)，所以是外插。
x−xn−1x−xn
f(x,y(x))=fn+fn−1
xn−xn−1xn−1−xn
=
1
[(x−xn−1)fn−(x−xn)fn−1]h
h3h
fn−fn−1 22
（20）
此式在区间[xn,xn+1]上积分可得 于是得到
yn+1=yn+
∫
xn+1
xn
f(x,y(x))dx=
h
(3fn−fn−1) （21） 2
注意到插值公式（20）的误差项含因子(x−xn−1)(x−xn)，在区间[xn,xn+1]上积分后
-184-
得出h，故公式（21）的局部截断误差为O(h)，精度比向前Euler公式提高1阶。
若取r=2,3,L可以用类似的方法推导公式，如对于r=3有
33
h
(55fn−59fn−1+37fn−2−9fn−3) （22） 245
其局部截断误差为O(h)。
如果将上面代替被积函数f(x,y(x))用的插值公式由外插改为内插，可进一步减小误差。内插法用的是yn+1,yn,L,yn−r+1，取r=1时得到的是梯形公式，取r=3时
yn+1=yn+可得
h
(9fn+1+19fn−5fn−1+fn−2) （23） 24
5
与（22）式相比，虽然其局部截断误差仍为O(h)，但因它的各项系数（绝对值）大为减小，误差还是小了。当然，（23）式右端的fn+1未知，需要如同向后Euler公式一
yn+1=yn+
样，用迭代或校正的办法处理。
§6 一阶微分方程组与高阶微分方程的数值解法
6.1 一阶微分方程组的数值解法 设有一阶微分方程组的初值问题
⎧y'i=fi(x,y1,y2,L,ym)
⎨ (i=1,2,L,m) （24）
y(a)y=i0⎩i
TTT
若记y=(y1,y2,L,ym)，y0=(y10,y20,L,ym0)，f=(f1,f2,L,fm)，则初值
问题（24）可写成如下向量形式 ⎨
⎧y'=f(x,y)
（25）
⎩y(a)=y0
如果向量函数f(x,y)在区域D：
a≤x≤b,
y∈Rm
m
连续，且关于y满足Lipschitz条件，即存在L>0，使得对∀x∈[a,b]，y1,y2∈R，都有
f(x,y1)−f(x,y2)≤Ly1−y2
那么问题（25）在[a,b]上存在唯一解y=y(x)。
问题（25）与（1）形式上完全相同，故对初值问题（1）所建立的各种数值解法可全部用于求解问题（25）。
6.2 高阶微分方程的数值解法
高阶微分方程的初值问题可以通过变量代换化为一阶微分方程组初值问题。 设有m阶常微分方程初值问题
⎧y(m)=f(x,y,y',L,y(m−1))a≤x≤b
（26） ⎨(1)(m−1)(m−1)
(a)=y0⎩y(a)=y0,y'(a)=y0,L,y
(m−1)
引入新变量y1=y,y2=y',L,ym=y，问题（26）就化为一阶微分方程初值问题
-185-
⎧y'1=y2 y1(a)=y0⎪(1)
==' ()yyyay320⎪2
⎪
（27） ⎨M M
⎪y'=y y(a)=y(m−2)
0mm−1⎪m−1
(m−1)⎪⎩y'm=f(x,y1,L,ym) ym(a)=y0
然后用6.1中的数值方法求解问题（27），就可以得到问题（26）的数值解。
最后需要指出的是，在化学工程及自动控制等领域中，所涉及的常微分方程组初值问题常常是所谓的“刚性”问题。具体地说，对一阶线性微分方程组
dy
=Ay+Φ(x) （28） dx
m
其中y,Φ∈R,A为m阶方阵。若矩阵A的特征值λi(i=1,2,L,m)满足关系
(i=1,2,L,m)
max|Reλi|>>min|Reλi|
1≤i≤m
1≤i≤m
Reλi
则称方程组（28）为刚性方程组或Stiff方程组，称数
s=max|Reλi|/min|Reλi|
1≤i≤m
1≤i≤m
为刚性比。对刚性方程组，用前面所介绍的方法求解，都会遇到本质上的困难，这是由数值方法本身的稳定性限制所决定的。理论上的分析表明，求解刚性问题所选用的数值方法最好是对步长h不作任何限制。
§7 Matlab解法
7.1 Matlab数值解
7.1.1 非刚性常微分方程的解法
Matlab的工具箱提供了几个解非刚性常微分方程的功能函数，如ode45，ode23，ode113，其中ode45采用四五阶RK方法，是解非刚性常微分方程的首选方法，ode23采用二三阶RK方法，ode113采用的是多步法，效率一般比ode45高。
Matlab的工具箱中没有Euler方法的功能函数。 （I）对简单的一阶方程的初值问题
⎨改进的Euler公式为
⎧y'=f(x,y)
=y(x)y00⎩
⎧
⎪yp=yn+hf(xn,yn)⎪
⎨yq=yn+hf(xn+h,yp)
⎪
⎪yn+1=1(yp+yq)⎩2
我们自己编写改进的Euler 方法函数eulerpro.m如下：
function [x,y]=eulerpro(fun,x0,xfinal,y0,n); if nargin
-186-
h=(xfinal-x0)/n; x(1)=x0;y(1)=y0; for i=1:n
x(i+1)=x(i)+h;
y1=y(i)+h*feval(fun,x(i),y(i)); y2=y(i)+h*feval(fun,x(i+1),y1); y(i+1)=(y1+y2)/2; end
例1 用改进的Euler方法求解
y'=−2y+2x+2x，(0≤x≤0.5)，y(0)=1 解 编写函数文件doty.m如下： function f=doty(x,y); f=-2*y+2*x^2+2*x;
在Matlab命令窗口输入：
[x,y]=eulerpro('doty',0,0.5,1,10) 即可求得数值解。
（II）ode23，ode45，ode113的使用 Matlab的函数形式是
[t,y]=solver('F',tspan，y0)
这里solver为ode45，ode23，ode113，输入参数 F 是用M文件定义的微分方程y'=f(x,y)右端的函数。tspan=[t0,tfinal]是求解区间，y0是初值。
例2 用RK方法求解
2
y'=−2y+2x2+2x，(0≤x≤0.5)，y(0)=1
解 同样地编写函数文件doty.m如下： function f=doty(x,y); f=-2*y+2*x^2+2*x;
在Matlab命令窗口输入：
[x,y]=ode45('doty',0,0.5,1) 即可求得数值解。
7.1.2 刚性常微分方程的解法
Matlab的工具箱提供了几个解刚性常微分方程的功能函数，如ode15s，ode23s，ode23t，ode23tb，这些函数的使用同上述非刚性微分方程的功能函数。 7.1.3 高阶微分方程y=f(t,y,y',L,y)解法 例3 考虑初值问题
y'''−3y''−y'y=0y(0)=0y'(0)=1y''(0)=−1
解 （i）如果设y1=y,y2=y',y3=y''，那么
(n)
(n−1)
⎧y1'=y2 y1(0)=0⎪
⎨y'2=y3 y2(0)=1
⎪y'=3y+yy y(0)=−1
3213⎩3
初值问题可以写成Y'=F(t,Y),Y(0)=Y0的形式，其中
Y=[y1;y2;y3]。
（ii）把一阶方程组写成接受两个参数t和y，返回一个列向量的M文件F.m：
function dy=F(t,y);
-187-
dy=[y(2);y(3);3*y(3)+y(2)*y(1)];
注意：尽管不一定用到参数t和y，M—文件必须接受此两参数。这里向量dy必须是列向量。
（iii）用Matlab解决此问题的函数形式为
[T,Y]=solver('F',tspan,y0)
这里solver为ode45、ode23、ode113，输入参数F是用M文件定义的常微分方程组，tspan=[t0 tfinal]是求解区间，y0是初值列向量。在Matlab命令窗口输入
[T,Y]=ode45('F',[0 1],[0;1;-1])
就得到上述常微分方程的数值解。这里Y和时刻T是一一对应的，Y(:,1)是初值问题的解，Y(:,2)是解的导数，Y(:,3)是解的二阶导数。
例4 求 van der Pol 方程
y''−μ(1−y)y'+y=0 的数值解，这里μ>0是一参数。
解 （i）化成常微分方程组。设y1=y,y2=y'，则有
2
⎧y'1=y2
⎨2
=−−y'μ(1y)yy⎩2121
（ii）书写M文件（对于μ=1）vdp1.m:
function dy=vdp1(t,y);
dy=[y(2);(1-y(1)^2)*y(2)-y(1)];
（iii）调用Matlab函数。对于初值y(0)=2,y'(0)=0，解为 [T,Y]=ode45('vdp1',[0 20],[2;0]);
（iv）观察结果。利用图形输出解的结果： plot(T,Y(:,1),'-',T,Y(:,2),'--')
title('Solution of van der Pol Equation,mu=1'); xlabel('time t');
ylabel('solution y');
解 (i)书写M文件vdp1000.m:
-188-
function dy=vdp1000(t,y);
dy=[y(2);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中按如下规定重新表达：
符号D表示对变量的求导。Dy表示对变量y求一阶导数，当需要求变量的n阶导数时，用Dn表示，D4y表示对变量y求4阶导数。
由此，常微分方程y''+2y'=y在Matlab中，将写成'D2y+2*Dy=y'。
7.2.1 求解常微分方程的通解
无初边值条件的常微分方程的解就是该方程的通解。其使用格式为：
dsolve('diff_equation')
dsolve(' diff_equation'，'var')
式中diff_equation 为待解的常微分方程，第1种格式将以变量t为自变量进行求解，第2种格式则需定义自变量var。
例6 试解常微分方程
x+y+(x−2y)y'=0
解 编写程序如下：
syms x y
diff_equ='x^2+y+(x-2*y)*Dy=0';
dsolve(diff_equ,'x')
7.2.2 求解常微分方程的初边值问题
求解带有初边值条件的常微分方程的使用格式为：
dsolve('diff_equation'，'condition1，condition2，…'，'var')
其中condition1，condition2，… 即为微分方程的初边值条件。
例7 试求微分方程
y'''−y''=x，y(1)=8,y'(1)=7,y''(2)=4
的解。
解 编写程序如下：
y=dsolve('D3y-D2y=x','y(1)=8,Dy(1)=7,D2y(2)=4','x')
7.2.3 求解常微分方程组
求解常微分方程组的命令格式为：
dsolve('diff_equ1，diff_equ2，…'，'var')
dsolve('diff_equ1，diff_equ2，…'，'condition1，condition2，…'，'var')
第1种格式用于求解方程组的通解，第2种格式可以加上初边值条件，用于具体求解。
例8 试求常微分方程组： 2
⎧f''+3g=sinx ⎨g'+f'=cosx⎩
-189-
的通解和在初边值条件为f'(2)=0,f(3)=3,g(5)=1的解。
解 编写程序如下：
clc,clear
equ1='D2f+3*g=sin(x)';
equ2='Dg+Df=cos(x)';
[general_f,general_g]=dsolve(equ1,equ2,'x')
[f,g]=dsolve(equ1,equ2,'Df(2)=0,f(3)=3,g(5)=1','x')
7.2.4 求解线性常微分方程组
（i）一阶齐次线性微分方程组
⎡x1⎤⎡a11La1n⎤
X'=AX，X=⎢M⎥，A=⎢MMM⎥ ⎥⎢⎥⎢⎢⎢⎦⎦⎣xn⎥⎣an1Lann⎥
At这里的’表示对t求导数。e是它的基解矩阵。X'=AX，X(t0)=X0的解为
X(t)=eA(t−t0)X0。
例9 试解初值问题
⎡213⎤⎡1⎤⎥X，X(0)=⎢2⎥ X'=⎢02−1⎥⎢⎢⎥⎢⎢⎣1⎥⎦⎦⎣002⎥
解 编写程序如下：
syms t
a=[2,1,3;0,2,-1;0,0,2];
x0=[1;2;1];
x=expm(a*t)*x0
（ii）非齐次线性方程组
由参数变易法可求得初值问题
X'=AX+f(t)，X(t0)=X0
的解为
X(t)=eA(t−t0)X0+∫eA(t−s)f(s)ds. t0t
例10 试解初值问题
⎡100⎤⎡0⎤⎡0⎤
X'=⎢21−2⎥X+⎢0⎥，X(0)=⎢1⎥。 ⎢⎥⎢⎥⎢⎥t⎢⎢⎢⎣321⎥⎦⎣ecos2t⎥⎦⎣1⎥⎦
解 编写程序如下：
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)
-190-
习 题 十 五
1. 用欧拉方法和龙格—库塔方法求微分方程数值解，画出解的图形，对结果进行分析比较。
（i）xy''+xy'+(x−n)y=0,y⎜2222⎛π⎞⎛π⎞⎟=2,y'⎜⎟=−（Bessel 方程，令π⎝2⎠⎝2⎠
21n=)，精确解y=sinx。 x2
（ii）y''+ycosx=0,y(0)=1,y'(0)=0，幂级数解 122496558 y=1−x+x−x+x−L 2!4!6!8!
2. 一只小船渡过宽为d的河流，目标是起点A正对着的另一岸B点。已知河水流速v1与船在静水中的速度v2之比为k。
（i）建立小船航线的方程，求其解析解。
（ii）设d=100m，v1=1m/s，v2=2m/s，用数值解法求渡河所需时间、任意时刻小船的位置及航行曲线，作图，并与解析解比较。
-191-