常微分方程的数值解法分析解析

第八章 常微分方程的数值解法
一．内容要点
考虑一阶常微分方程初值问题:⎪⎩⎪⎨⎧==0
0)()
,(y x y y x f dx dy
微分方程的数值解：设微分方程的解y (x )的存在区间是[a,b ]，在[a,b ]内取一系列节
点a= x 0< x 1<…< x n =b ,其中h k =x k+1-x k ；(一般采用等距节点，h=(b-a)/n 称为步长)。

在每个节点x k 求解函数y(x)的近似值:y k ≈y(x k )，这样y 0 , y 1 ,...,y n 称为微分方程的数值解。

用数值方法，求得f(x k )的近似值y k ，再用插值或拟合方法就求得y(x)的近似函数。

(一)常微分方程处置问题解得存在唯一性定理
对于常微分方程初值问题:⎪⎩⎪⎨⎧==0
0)()
,(y x y y x f dx dy
如果：
(1) 在B y y A x x 00≤-≤≤,的矩形内),(y x f 是一个二元连续函数。

(2) ),(y x f 对于y 满足利普希茨条件，即
2121y y L y x f y x f -≤-),(),(则在C x x 0≤≤上方程⎪⎩⎪⎨⎧==0
0)()
,(y x y y x f dx
dy
的解存在且唯一，这里C=min((A-x 0),x 0+B/L),L 是利普希茨常数。

定义：任何一个一步方法可以写为),,(h y x h y y k k k 1k Φ+=+，其中),,(h y x k k Φ称为算法的增量函数。

收敛性定理：若一步方法满足： (1)是p 解的.
(2) 增量函数),,(h y x k k Φ对于y 满足利普希茨条件.
(3) 初始值y 0是精确的。

则),()()(p h O x y kh y =-kh =x -x 0,也就是有
0x y y lim k x x kh 0h 0
=--=→)(
(一)、主要算法 1．局部截断误差
局部截断误差：当y(x k )是精确解时，由y(x k )按照数值方法计算出来的1~
+k y 的误差y (x k+1)- 1~
+k y 称为局部截断误差。

注意：y k+1和1~
+k y 的区别。

因而局部截断误差与误差e k +1=y (x k +1) -y k +1不同。

如果局部截断误差是O (h p+1),我们就说该数值方法具有p 阶精度。

1. 显式欧拉格式：⎩⎨⎧=+=+001)()
,(y x y y x hf y y k k k k k =1,2,⋯n-1.
显式欧拉格式的特点：
(1)、单步方法； (2)、显式格式；
(3)、局部截断误差)(2h O 因而是一阶精度 。

隐式欧拉格式⎩
⎨⎧=+=+++00111)()
,(y x y y x hf y y k k k k
2. 两步欧拉格式：⎩⎨⎧=+=-+001)()
,(y x y y x hf 2y y k k 1k k k =1,2,⋯n-1.
两步欧拉格式的局部截断误差)(3h O ，因而是二阶精度 3. 梯形方法:
=+1k y []⎪⎩
⎪⎨⎧
=+⋅+=+++00)(),(),(y x y y x f y x f 2h y y k k 1k 1k k 1
k ； 3.改进的欧拉方法：
预测值: ),(^
k k k 1k y x f h y y ⋅+=+
校正值: =+1k y ⎥⎦
⎤
⎢⎣⎡+⋅+++),(),(^k k 1k 1k k y x f y x f 2h y 。

或改写为:⎪⎪⎩
⎪
⎪⎨⎧
+=+=+=++)
(21
),(),(11c p k p k k c k k k p y y y y x f h y y y x f h y y
4、梯形方法与改进欧拉方法的截断误差是O(h 3),具有二阶精度。

5、龙格-库塔法的思想
1). 二阶龙格-库塔法计算公式： ⎪⎩⎪
⎨⎧λ+λ-+=++==+)
)((),(),(211
121K K 1h y y hK p y h p x f K y x f K k k k k k k
当:2
1
=λp 时，得一簇龙格-库塔公式，其局部截断误差均为O (h 3)都是二阶精度。

特别取1p ==
λ,21
,就是改进欧拉公式。

取1p =λ=,2
1
,得二阶龙格—库塔法为：
⎪
⎩⎪⎨⎧+=++==+2
1121)2,2(),(hK y y K h y h x f K y x f K k k k
k k k ，称为二阶中点格式。

2)、经典龙格－库塔格式(也称为标准龙格－库塔方法)：
⎪⎪⎪⎪
⎩
⎪
⎪
⎪⎪⎨⎧++++=++=++=++==+)22(6),()2,2()2,2(),(432113423121K K K K h y y hK y h x f K K h y h x f K K h y h x f K y x f K k
k k k k
k k
k k k 四阶龙格－库塔方法的截断误差为()5h O ，具有四阶精度。

一般一阶常微分方程初值问题用四阶龙格－库塔方法计算，其精度均满足实际问题精度要求。

3).变步长龙格－库塔方法： 从节点x k 出发，以步长h 据四阶龙格－库塔方法求出一
个近似值)
(1h k y +,然后以步长h/2求出一个近似值)2/(1h k y +，得误差事后估计式：
)()
(1)2/(1)
2/(1
h k h k h k 1k y y 15
1y y ++++-≈- 根据)
(1)2/(1h k h k y y ++-=∆来选取步长h 。

4). RKF 格式：变步长龙格－库塔方法，因频繁加倍或折半步长会浪费计算量。

Felhberg 改进了传统龙格－库塔方法，得到RKF 格式，较好解决了步长的确定，而且提高了精度与稳定性，为Matlab 等许多数值计算软件采用。

4/5阶RKF 格式：由4阶龙格－库塔方法与5阶龙格－库塔方法结合而成。

⎪
⎭⎫
⎝⎛+-+++=⎪
⎭
⎫
⎝⎛-+++=++654311^
5431155250956430285611282566561351651410121972565140821625K K K K K h y y K K K K h y y n n n n
⎪⎪
⎪
⎪⎪⎪⎪⎩
⎪⎪⎪⎪
⎪⎪⎪⎨
⎧
-+-+-+=-+-++=+-++=++
+=++==)401141041859256535442278,2()
410484551336808216439,()219772962197720021971932,1312()32
9
323,83()
4,4(),(443216432153214213121K K K K K y h x f K K K K K y h x f K K K K y h x f K K K y h x f K K h y h x f K y x f K n n n n n n
n n n n n n
Felhberg 得到的最佳步长hs ，其中h 为当前步长，4
/111^
⎪⎪⎪⎪⎭
⎫
⎝
⎛-ε
=++n n y y h s .ε为精度要求，
若s<0.75,步长折半；若s>1.5步长加倍。

6.亚当斯方法(Adams) 1).显式Adams 方法：
记：),(k k k y x f f =；
二阶显式Adams 方法：][211-+-+=k k k k f f 3h
y y ;
三阶显式Adams 方法：][2211--++-+=k k k k k f 5f 16f 231h
y y ;
四阶显式Adams 方法：]5559379[24
1231k k k k k k f f f f h
y y +-+-+=---+.
2). 隐式Adams 方法
二阶隐式Adams 方法：][211++-+=k k k k f f 3h
y y ;
三阶隐式Adams 方法：]8[2111-++-++=k k k k k f f f 51h
y y ;
四阶隐式Adams 方法：)9195(24
1121+--+++-+=k k k k k k f f f f h
y y
3）.Adams 预报-校正系统：
先用显式格式作为预测值，再用隐式格式来校正。

预测值： )9375955(24
3211---+-+-+=k k k k k k f f f f h
y y
校正值：
()),9195(24
11121++--+++-+
=k k k k k k k y x f f f f h
y y 4).改进的Adams 预报-校正系统：
预测：)9375955(243211---+-+-+=k k k k k k f f f f h
y p
改进：)(k k 1k 1k p c 270251
p m -+=++
校正值：
())519,9(24
21111--++++-++
=k k k k k k k f f f m x f h
y c 改进：)(1k 1k 1k 1k p c 270
19
c y ++++--= 7.收敛与稳定性
对于固定的kh x x k +=0，如果数值解y k 当,0h →（同时+∞→n ）时趋向于准确解y (x k )，则称该数值方法方法是收敛的。

如果一种数值方法，在节点值y k 上大小为δ的扰动，于以后各节点值y m (m >k )上产生的偏差均不超过δ，则称数值方法是稳定的.
一般，隐式格式较显式格式有较好的稳定性。

8.刚性方程组：
考虑n 阶常微分方程组：)(x Ay y φ+='，⋯⋯⋯（*）
若矩阵A 的特征值λ1,λ2,⋯λn 的实部Re(λi )<0，i=1,2,⋯,n.，则称)
Re()
Re(i n
i 1i n i 1min max S λλ=
≤≤≤≤为*式
的刚性比，当刚性比很大时称*式为刚性方程组。

刚性方程组需采用稳定性较好的方法。

二. MATLAB 的有关命令
[t ，x]=solver （’f’,ts,x0,options ）
t:输出的自变量值, x: 输出的函数值；
solver :包括ode45、 ode23 、ode113、ode15、sode23s.
非刚性系统：
ode23：用2阶(及3阶)龙格-库塔算法。

ode45：用4阶(及5阶)龙格-库塔-算法。

ode113(Adams-Bashforth-Moulton PECE)多步方法。

刚性系统：
ode15s(Gear 方法)多步方法
ode23s(二阶modified Rosenbroack formula)单步。

ode23t(trapezoidal rule)solve DAEs. ode23tb(TR-BDF2)低精度算法。

f ：由待解方程写成的m-文件名
ts=[t0，tf]，t0、tf 为自变量的初值和终值 x0：函数的初值
options=odeset(’reltol’,rt,’abstol’,at ), rt ，at ：分别为设定的相对误差和绝对误差. (缺省时设定相对误差10-3, 绝对误差10-6) 注意:
1、在解n 个未知函数的方程组时，x0和x 均为n 维向量，m-文件中的待解方程组应以x 的分量形式写成.
2、使用Matlab 软件求数值解时，高阶微分方程必须等价地变换成一阶微分方程组.
四.重点算法Matlab 程序
Euler 格式、
function[x,y]=naeuler(dyfun,xspan,y0,h) x=xspan(1):h:xspan(2); y(1)=y0;
for n=1:length(x)-1
y(n+1)=y(n)+h*feval(dyfun,x(n),y(n)); end
x=x';y=y';
隐式Euler 格式
function[x,y]=naeulerb(dyfun,xspan,y0,h) x=xspan(1):h:xspan(2); y(1)=y0;
for n=1:length(x)-1
y(n+1)=iter(dyfun,x(n+1),y(n),h);
end
x=x';y=y';
function y=iter(dyfun,x,y,h)
y0=y;e=1e-4;K=1e+4;
y=y+h*feval(dyfun,x,y);
y1=y+2*e;k=1;
while abs(y-y1)>e
y1=y;
y=y0+h*feval(dyfun,x,y);
k=k+1;if k>K,error('迭代发散');end
end
改进Euler格式
function[x,y]=naeulerg(dyfun,xspan,y0,h)
x=xspan(1):h:xspan(2);
y(1)=y0;
for n=1:length(x)-1
K1=feval(dyfun,x(n),y(n));
y(n+1)=y(n)+h*K1;
K2=feval(dyfun,x(n+1),y(n+1));
y(n+1)=y(n)+h*(K1+K2)/2;
end
x=x';y=y';
4阶经典Runge-Kutta格式
function[x,y]=nak4(dyfun,xspan,y0,h)
x=xspan(1):h:xspan(2);
y(1)=y0;
for n=1:length(x)-1
K1=feval(dyfun,x(n),y(n));
K2=feval(dyfun,x(n)+h/2,y(n)+h/2*K1);
K3=feval(dyfun,x(n)+h/2,y(n)+h/2*K2);
K4=feval(dyfun,x(n+1),y(n)+h*K3);
y(n+1)=y(n)+h*(K1+2*K2+2*K3+K4)/6; end
x=x';y=y';
变步长4阶经典Runge-Kutta格式
function[x,y]=nark4v(dyfun,xspan,y0,e,h)
if nargin<5,h=(xspan(2)-xspan(1))/10;end
n=1;x(n)=xspan(1);y(n)=y0;
[y1,y2]=comput(dyfun,x(n),y(n),h);
while x(n)<xspan(2)-eps;
if abs(y2-y1)/10>e
while abs(y2-y1)/10>e
h=h/2;
[y1,y2]=comput(dyfun,x(n),y(n),h);
end else
while abs(y2-y1)/10<=e h=2*h;
[y1,y2]=comput(dyfun,x(n),y(n),h); end
h=h/2;h=min(h,xspan(2)-x(n)); [y1,y2]=comput(dyfun,x(n),y(n),h); end
n=n+1;x(n)=x(n-1)+h;y(n)=y2;
[y1,y2]=comput(dyfun,x(n),y(n),h); end
x=x';y=y';
function [y1,y2]=comput(dyfun,x,y,h); y1=rk4(dyfun,x,y,h); y21=rk4(dyfun,x,y,h/2);
y2=rk4(dyfun,x+h/2,y21,h/2); function y=rk4(dyfun,x,y,h); K1=feval(dyfun,x,y);
K2=feval(dyfun,x+h/2,y+h/2*K1); K3=feval(dyfun,x+h/2,y+h/2*K2); K4=feval(dyfun,x+h,y+h*K3); y=y+h*(K1+2*K2+2*K3+K4)/6;
五.试验例题
例1 不同方法的精度比较
用多种方法解下述初值问题，并与其准确解x e y x +=-进行比较。

⎩
⎨⎧=≤<++-='1)0(6
.00,1y x x y y
解：取步长h=0.1, x k =kh (k=0,1,⋯,6)。

各方法计算结果及绝对误差见表(1)、(2)、(3)。

表（1）
x n 欧拉公式 改进的欧拉公式 四阶标准龙格—库塔公式
y n 误差 y n 误差 y n 误差
0.0 1.000000 0 1 0 1 0 0.1 1.000000 1.0050000 1.00483750 0.2 1.010000 1.019025 1.01873090 0.3 1.029000 1.041218 1.04081842 0.4 1.056100 1.070802 1.07032029 0.5 1.090490 1.107076 1.10653093 0.6 1.131441 1.149404 1.14881193
表（2） 四阶阿当姆斯公式
n x n
显示公式隐式公式①
y n误差y n误差
3 0.3 取自准确解— 1.04081801
4 0.4 1.07032292 1.07031966
5 0.5 1.10653548 1.10653041
6 0.6 1.14881481 1.14881101
备注y1, y2, y3 取自准确解y1, y2 取自准确解
①本题关于y为线性，代入隐式公式，可解出y n+1，因此可直接求隐式公式之解。

表（3）四阶阿当姆斯预测—校正公式
n x n
显示公式隐式公式①
Y n误差y n误差
4 0.4 1.07031992 1.07032014
5 0.5 1.10653027 1.10653077
6 0.6 1.14881103 1.14881175
备注y1, y2, y3由四阶标准龙格—库塔公式提供
对比以上各表数据知，在相同步长下求解同一问题时，方法阶数越高，解的精度也越高，一阶的欧拉公式精度最低，四阶标准龙格—库塔公式的精度大大高于改进的欧拉公式。

四阶阿当姆斯方法、显式的精度略低于隐式。

同是阿当姆斯预测—校正系统，带误差修正的精度又高于不带误差修正的。

从计算量看，四阶龙格—库塔法计算量最大，每前进一步要计算4次函数值f，而带误差修正的阿当姆斯预测—校正法与它精度差不多的，每前进一步只计算两次函数值f，所以后者是可取的。

但它是四步法，必须用其它方法提供出始值，程序略复杂些。

例2 (步长的计算结果的影响)用欧拉公式求下述初值问题在x = 1处的近似解，并与准确解y (1)=1比较。

解分别取步长h=10-1、h=10-2、h=10-3、h=10-4计算，结果见下表。

由表中数据可见，
当时结果完全失真，而取比计算量增加了十倍，但解
的精度却基本一样，可见取太浪费计算量了。

步长h y (1) 的计算值
上述结果差异很大的原因在于欧拉公式的绝对稳定区间为（－2, 0）步长h 应满足
)0,2(-∈∂∂h y
f
，对本题，x x y y x f 2)(
1000),(2+--=
1000-=∂∂y
f
，故应取h 满
足
即 。

可见取 时欧拉公式是数值不稳
定的，导致结果失真，而取
和
都满足稳定性要求，可用于求解。

此例说明，求解微分方程数值解一定选取注意步长，过大则导致解的失真，过小又会使计算量大增。

究竟取多大步长才合适，不仅取决于所采用的数值方法，还决定于待解微分方程本身的特性。

例3：(Matlab 不同命令之差异)取h=0.2分别用不同的Matlab 命令解方程
⎩
⎨
⎧=='11y 2t/y
-y y )( (0<t<4). 其准确解为:t y 21+= 解：
>> [t,y]=ode45(odefun,[0,4],1);n=length(t);e=sqrt(sum((sqrt(1+2*t)-y).^2/n));[n,e]
ans =45.00000000000000 0.00020257808813
>> [t,y]=ode23(odefun,[0,4],1);n=length(t);e=sqrt(sum((sqrt(1+2*t)-y).^2/n));[n,e] ans =13.00000000000000 0.19048360113013
>> [t,y]=ode113(odefun,[0,4],1);n=length(t);e=sqrt(sum((sqrt(1+2*t)-y).^2/n));[n,e] ans =18.00000000000000 0.00973277413962
>> [t,y]=ode23t(odefun,[0,4],1);n=length(t);e=sqrt(sum((sqrt(1+2*t)-y).^2/n));[n,e] ans =18.00000000000000 0.03919668151270
>> [t,y]=ode23s(odefun,[0,4],1);n=length(t);e=sqrt(sum((sqrt(1+2*t)-y).^2/n));[n,e] ans =81.00000000000000 2.54374270001824
>> [t,y]=ode23tb(odefun,[0,4],1);n=length(t);e=sqrt(sum((sqrt(1+2*t)-y).^2/n));[n,e] ans =15.00000000000000 0.24308870757312
>> [t,y]=ode15s(odefun,[0,4],1);n=length(t);e=sqrt(sum((sqrt(1+2*t)-y).^2/n));[n,e] ans =22.00000000000000 0.45509043060378
可见ode45精度较高，计算量较大，ode23计算量小，但误差大，ode113适中，用刚性方程组揭发解非刚性问题不适合，特别ode23s 计算量大且但误差大。

例4：解刚性方程组： ⎩
⎨⎧=='=-='10z -100z,z 2
1y 99.99z,-01y 0y )()(.
其准确解为:x x x
e z e e
y 10010001.0,---=+=就h=0.01、h=0.02使用改进欧拉格式及方程组，并在一张
图上绘出其图形。

function[x,y]=naeuler2s(dyfun,xspan,y0,h) x=xspan(1):h:xspan(2);
y=zeros(length(y0),length(x)); y(:,1)=y0(:);
for n=1:length(x)-1
K1=feval(dyfun,x(n),y(:,n)); y(:,n+1)=y(:,n)+h*K1;
K2=feval(dyfun,x(n+1),y(:,n+1)); y(:,n+1)=y(:,n)+h*(K1+K2)/2; end
x=x';y=y'; %用欧拉格式解常微分方程组的m 文件
function f=dyfun(x,y)
f(1)=-0.01*y(1)-99.99*y(2); f(2)=-100*y(2);
f=f(:); %函数m 文件
clear;
[x,y]=naeuler2s(@dyfun,[0 500],[2,1],0.01); plot(x,y),axis([-50 500 -0.5 2])
gtext('y(x) h=0.01'), gtext('z(x) h=0.01') %步长为0.01时解方程并画图 clear;
[x,y]=naeuler2s(@dyfun,[0 500],[2,1],0.02); hold on
plot(x,y),gtext('y(x) h=0.02'), gtext('z(x) h=0.02') %步长为0.02时解方程并画图。

clear;
x=0:0.01:500;
y=exp(-0.01*x)+exp(-100*x);z=exp(-100*x); hold on
plot(x,y,'r',x,z,'r'),gtext('y(x) 真根'), gtext('z(x) 真根') %画出真根图形
运行结果真根与h=0,01时结果重合，h=0.02时结果是错误的。

故刚性方程组应采用稳定性较好方法。

例6：数值实验（摘自蔡）观察欧拉显示方法的收敛性。

内容：取h=0.1, 0.01,⋯用欧拉显式方法求解：用欧拉法
解方程⎩⎨⎧=='1
1y xy y 1/3
)( (0<x<1).
计算到(2)y ，并与精确解3/2
2()(2)/3y x x ⎡⎤=+⎣⎦
比较。

MATLAB 程序
t=1; y=1;h=0.1 for k=1:10;
y=y+h*t*y^(1/3) t=t+h end
运行结果
h y t 0.1 0.01 0.001 0.0001
y=((t 2+2)/3)3/2
(
准确解)
1 1 1 1 1 1
1.1 1.1 1.1175 1.1068 1.1068
1.10681660630838 1.2
1.2136
1.2393
1.2277
1.2279
1.22787959356620
1.3 1.3416 1.3762 1.3639 1.3641 1.36413599028836 1.4 1.4849 1.5294 1.5162 1.5165 1.51656453868604 1.5 1.6447 1.6998 1.6857 1.6861 1.68617060118373 1.6 1.8217 1.8883 1.8734 1.8739 1.87398185690257 1.7
2.017 2.0962 2.0804 2.0810 2.08104468957300 1.8 2.2319 2.3243 2.3076 2.3083 2.30842116960842 1.9 2.4671 2.5739 2.5563 2.5571 2.55718653993016 2.0
2.7239
2.8177
2.8274 2.8283
2.82842712474619
欧拉显示方法是收敛的。

例7：实验目的：比较不同算法用于“非光滑”解时的结果。

内容：用欧拉显式方法和2阶龙格−库塔算法求解：sin ,(0)1,y kx y '== 计算(2)y π并在一张图上画出解曲线。

k=6时的程序如下： clear
t=0;y=1;h=pi/60; for k=0.1:0.1:2;
y=y+h*abs(sin(6*t)) t=t+h
plot(t,y,'o'), hold on, end
t0=0;tf=1/3*pi;x0=[1];k=pi/60
[t,x]=ode23('ex545f',[t0:k:tf],1)
plot(t,x), hold on,
for i=0:pi/60:pi/3
if i<=pi/6;
y1=7/6-1/6*cos(6*i)
else
y1=9/6+1/6*cos(6*i)
end
plot(i,y1,'*'),hold on,
end
例8 实验目的：观察欧拉显式方法的数值不稳定性。

内容：取h= 0.05, 0.04, 0.03, ⋯, 用欧拉显式方法求解：
'=-=并画出曲线，同样用欧拉隐式方法重复求解上述问题。

y y y
50,(0)100,
解：欧拉显式方法程序略。

'是y的线性函数故用欧拉隐式方法可直接解出y。

程序如下
因50y
y-
=
clear;
h=0.01;
t=0;y=100;
for k=0:h:1;
y=y/(1+50*h);
t=t+h;
plot(t,y,'b *'), hold on,
end
gtext('隐t=0.01')
运行结果：
图中标示为“隐00x ”为隐式欧拉格式结果，其余是欧拉格式的运行结果，可见隐式欧拉格式较欧拉格式稳定。

例9：解2阶微分方程⎪⎪
⎩⎪⎪⎨
⎧==π-=+....,)(0
y 00y 2t 1y y 2设初始时间t0=0，终止时间tf=3π， 解：将方程改写为：⎪⎩⎪⎨⎧π
-+-==2t 1x x x x 2
122
1..
，
写为矩阵表达式：⎪⎪⎭
⎫
⎝⎛=⎪⎪⎭⎫ ⎝⎛⎪⎪⎭⎫ ⎝⎛π-⎪⎪⎭⎫ ⎝⎛+⎪⎪⎭⎫ ⎝⎛⎪⎪⎭⎫ ⎝
⎛-=⎪⎪⎭⎫ ⎝⎛000x 0x 2t 110x x 0110x x 2122121)()(;..。

本题解析解：2
2
2
t cont 12
1y ππ-
-+
=))((
MATLAB 程序
function xdot=ex545f(t,x) u=1-(t.^2)/(pi^2);
xdot=[0 1;-1 0]*x+[0 1]'*u; %建立ex545f.m 函数程序。

clf,t0=0;tf=3*pi;x0t=[0;0]; %给出初始值。

[t,x]=ode23('ex545f',[t0,tf],x0t) %调用MATLAB 中的数值积分2阶3阶龙—库公式。

y=x(:,1); %y 是x 的第2列。

for I=1:length(t);
y2(I)=(1+2/(pi^2))*(1-cos(t(I)))-t(I)^2/(pi^2) end %在数值积分的序列点上计算解析解的值。

u=1-(t.^2)/(pi^2);
clf,plot(t,y,'-',t,u,'+',t,y2,'o') % 绘出数值解、解析解的图形。

legend('数值解','输入量','解析解')
运行结果解析解与数值解看不出差别，是因为本数值积分是按默认精度10-3,如果用长格式输出y 、y2可看出它们的微小误差。

例10 解微分方程组.
⎪⎪
⎩⎪
⎪⎨⎧
===-=-==1
)0(,1)0(,0)0(51.0'''321213312321y y y y y y y y y y y y
解 1、建立m-文件rigid.m 如下：
function dy=rigid(t,y) dy=zeros(3,1); dy(1)=y(2)*y(3); dy(2)=-y(1)*y(3);
dy(3)=-0.51*y(1)*y(2);
2、取t0=0，tf=12，输入命令：
>> [T,Y]=ode45('rigid',[0 12],[0 1 1]); plot(T,Y(:,1),'-',T,Y(:,2),'*',T,Y(:,3),'+') 3、结果如图.
【例11】采用最简洁格式的ODE 文件和解算指令，研究围绕地球旋转的卫星轨道。

（1）问题的形成 轨道上运动的卫星，在Newton 第二定律
22dt
r
d m ma F ==
，和万有引力定律
程序见数学建模高教出版社数学建模与
数学实验To Matlab （ff5）
r r mM G F E 3-=作用下，有r r M G dt r d a E
322-==。

即3r x GM a E x -=，3r
y
GM a E y -=，而
22y x r +=。

这里)/(10672.62211kg m N G ⋅⨯=-是引力常数，)(1097.524kg M E ⨯=是地球
的质量。

又假定卫星以初速度)/(4000)0(s m v y =在)(102.4)0(7
m x ⨯-=处进入轨道。

（2）构成一阶微分方程组
令[][
][]T
T
y
x T
y x y x
v v y x
y y y y Y ''===43
21
，则
⎥
⎥⎥⎥
⎥⎥⎥⎦⎤
⎢⎢⎢
⎢⎢⎢
⎢⎣
⎡+⋅-+⋅-=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡''''='2
/322
2
2/322
14
34
3
21
)()()(y x y GM y x y GM y y y y y y t Y E
E
(5.14.2.1-5)
初始条件为
[]
T
y v x Y )0(00)0()0(=
(5.14.2.1-6)
（1）根据式(5.14.2.1-5)编写最简洁格式的ODE 文件 [DYDt2.m]
function Yd=DYDt2(t,Y,flag,G,ME) % flag 按ODE 文件格式规定，必须是第三输入宗量。

对它的赋值由ode45指令自动产生。

第4、5宗量是被传递的参数 switch flag
case '' %按规定：这里必须是空串。

在此为"真"时，完成以下导数计算。

X=Y(1:2);V=Y(3:4);r=sqrt(sum(X.^2));Yd=[V;-G*ME*X/r^3]; otherwise
error(['Unknown flag ''' flag '''.']); end
（2）对微分方程进行解算
G=6.672e-11;ME=5.97e24;vy0=4000;x0=-4.2e 7;t0=0;tf=60*60*24*9;
tspan=[t0,tf]; %指定解算微分方程时间区间 Y0=[x0;0;0;vy0]; %按式(5.14.2.1-6)给定初值向量
[t,YY]=ode45('DYDt2',tspan,Y0,[],G ,ME); %第4宗量取缺省设置,第5、6宗量是被传递参数
X=YY(:,1); %输出Y 第一列是位移数据x(t) Y=YY(:,2); %输出Y 第二列是位移数据y(t) plot(X,Y,'b','Linewidth',2); hold on
axis('image') %保证x 、y
轴等长刻度，且坐
标框恰包容图形
[XE,YE,ZE] = sphere(10); %产生单位球面数据
RE=0.64e7; %地球半径
XE=RE*XE;YE=RE*YE;ZE=0*ZE; %坐标纸上的地球平面数据
mesh(XE,YE,ZE),hold off %绘地球示意图
【*例12】带事件设置的ODE文件及主程序编写演示。

本例将以较高精度计算卫星经过近地点和远地点的时间，并在图上标志。

（1）ODE文件的编写
[DYDt3.m]
function varargout=DYDt3(t,Y,flag,G,ME,tspan,Y0)
% DYDt3.m 供主程序调用的ODE函数文件
% 本文件自带三个子函数：f,fi,fev。

% t, Y 分别是自变量和一阶函数向量，是最基本的输入宗量。

% flag 第三输入宗量，它专供解算指令（如ode45）作调用通知。

% 在运行中，解算指令会根据需要向flag发不同的字符串。

% varargout 是"变"输出宗量。

它由变维的元胞数组构成。

每个元胞中可以存放指令
% 所产生的任意形式的数据。

switch flag
case '' % 必须用空串符。

情况为"真"时，计算导数 dY/dt = f(t,Y)。

varargout{1} = f(t,Y,G,ME); %输出为一个元胞，容纳f子函数的一个输出Yd
case 'init' % 必须用'init'，情况为"真"时，传送计算区间、初值、设置参数。

[varargout{1:3}] = fi(tspan,Y0); %输出为三个元胞，容纳fi子函数的三个输出量
case 'events' % 必须用'events'，情况为"真"时，设置事件性质。

[varargout{1:3}] = fev(t,Y,Y0); %输出三个元胞，容纳fev子函数的三个输出量
otherwise
error(['Unknown flag ''' flag '''.']);
end
% ------------------------------------------------------------------
function Yd = f(t,Y,G,ME)
%计算导数子函数，被"父"函数DYDt3调用。

X=Y(1:2);V=Y(3:4);r=sqrt(sum(X.^2));Yd=[V; -G*ME*X/r^3];
% ------------------------------------------------------------------
function [ts,y0,options] = fi(tspan,Y0)
%设置时间区间、初值、算法参数子函数，被"父"函数DYDt3调用。

ts=tspan;y0 = Y0;
options = odeset('Events','on','Reltol',1e-5,'Abstol',1e-4);
%开动ode45的"事件判断"功能，设置相对误差和绝对误差。

% ------------------------------------------------------------------
function [value,isterminal,direction] = fev(t,Y,Y0)
%事件判断子函数，被"父"函数DYDt3调用。

dDSQdt = 2 * ((Y(1:2)-Y0(1:2))' * Y(3:4));
%dDSQdt是离初始点的二乘距离u的时间导数du/dt，而u=(x-x0)^2+(y-y0)^2 。

value = [dDSQdt; dDSQdt]; %定义两个穿越0的事件
direction = [1; -1]; %第一事件：以渐增方式穿越0。

第二事件：以渐减方式穿越0
isterminal = [1; 0]; %第一事件发生后，终止计算；而第二事件发生后，继续计算。

（2）运行以下主程序
G=6.672e-11;ME=5.97e24;vy0=4000; x0=-4.2e7;t0=0;tf=60*60*24*9; tspan=[t0,tf];Y0=[x0;0;0;vy0];
[t,YY,Te,Ye,Ie]=ode45('DYDt3',[],[],[],G,ME,tspan,Y0); %
<3>
X=YY(:,1);Y=YY(:,2);
plot(X,Y,'b','Linewidth',2);hold on text(0,6e7,'轨道','Color','b')
%产生蓝色文字注释
axis('image'); %保证x 、y 轴等长刻度，且坐标框恰包容图形在三个事件发生点上画标记 plot(Ye(1,1),0.4e7+Ye(1,2),'r^','MarkerSize',10) plot(Ye(2,1),0.4e7+Ye(2,2),'bv','MarkerSize',10) plot(Ye(3,1),-0.4e7+Ye(3,2),'b^
','MarkerSize',10)
%把轨道的半周期和全周期标在图上
text(0.8*Ye(3,1),-2e7+Ye(3,2),['t3=' sprintf('%6.0f',Te(3))]) text(0.8*Ye(2,1),1.5e7+Ye(2,2),['t2=' sprintf('%6.0f',Te(2))]) %在x-y 坐标上画地球 [XE,YE,ZE]=
sphere(10);RE=0.64e7;XE=RE*XE;Y E=RE*YE;ZE=0*ZE; mesh(XE,YE,ZE)
text(1e7,1e7,'地球','Color','r'),
hold off
%产生红色文字注释
例13 范例：状态转移方程组模型
计算机网络可靠性分析问题随着计算机通信网络系统特别是Internet 网络日益广泛的应用，提高系统的可靠性意义重大．研究和分析具有实用性的高可靠计算机通信网络系统，是国际上非常活跃的一个研究方向．
（1） 问题及假设
我们将研究的计算机通信网络系统称为无冗余的防火墙协议系统．计算机随时可能发生三个事件一无故障、间歇故障和永久故障．因此，计算机一般处于三种工作状态：
无故障工作、带故障工作和不工作．这三种状态之间的转移过程如图3．9所示．要求建立该系统的状态转移模型，并进行可靠性分析． (2) 分析与模型
该问题属于状态转移问题，利用马尔科夫状态转移原理，用)(1t P ,)(2t P 和)(3t P 分别表示系统处于无故障工作、带故障工作和不工作三种状态的概率，则有以下状态转移方程组：
⎪⎪⎪⎪
⎩⎪
⎪
⎪
⎪⎨⎧+++=++-=++-=+++-=)()()()2/()
()()]([)()2/()
()()()()()(]2/)2/[()
(21321221211t P t P dt t dP t P t P dt t dP t P t P t P t P dt t dP t p t p t p t t p t t p λλλλλλγλγλλγλλλ
初始条件]0,0,1[)]0(),0(),0([321=P P P ,参数取值
1.0~01.0:,10~10:,10~10:3445γλλ----t p 这是一个带参数的微分方程组模型．
(3) 模型求解
求解这个带参数的微分方程组模型有两种方法，一是用特征根法求解析解，另一个是用数值解法求数值解．分别求解如下：
假定模型中的参数取下限值，即
01.0,10,104
5===--γλλt p ①.解析解法
利用常系数线性微分方程组的特征根法． MATLAB 程序；
lp=10^(-5);lt=10^(-4);gm=0.01;
A=[-(lp+lt),gm,0;lt/2,-(gm+lp+lt),0;lp+lt/2,lp+lt,0]; [l,v]=eig(A); c=inv(1)*[1;0;0] 输出结果：
特征值（0，-0.0102,-0.0001） 与之对应的特征向量(0 0 1);
（0．7053 -0．7089 0．0035）； （-0．7053 -0．0035 0．7089）
根据特征值和特征向量写出其解析解,利用初始条件]0,0,1[)]0(),0(),0([321=P P P 得到特解
t t e e t P 0001.00102.0199503724.0004937.0)(--+= t t e e t P 000.00102.020049378.00049623.0)(--+-= t t e e t P 0001.00102.0300011612.10000245.01)(---+= ②．数值解法
MATLAB 程序如下
首先编写m 文件（eqs0．m ）
function xdot=eqs0(t,p,flag,lp,lt,gm)
a=[-(lp=lt),gm,0;lt/2,-(gm+lp+lt),0;lp+lt/2,lp+lt,0];
p=[p(1);p(2);p(3)] xdot=a*p;
在工作空间执行以下程序：
ts=[0
1000];p0=[1;0;0];lp=10^(-5);lt=10^(-4);gm=0.01;
[t,p]=ode23(’eqs0’,ts,p0,[],lp,lt,gm) plot(t,1-p(:,));
xlabe l(‘时间t(小时)’);ylabel(‘可靠度R(t)’);
050010001500200025003000-2.5
-2-1.5
-1-0.5
0.5
11.5
2title(‘参数取值lp=0.00001;lt=0.00001;gm=0.01’);grid 输出结果如图3.10
计算表明：当t=1000小时情况下，].,.,.[)](),(),([058400047093690t p t p t p 321= 显示系统工作的概率为94.18%，从图3.10中还可以看出系统可靠性变化更加细致的情况 , 何时
可靠度达到99%，何时达到98%等等。

③.敏感性分析
当参数取值在以下范围 1.0~01.0:,10~10:,10~10:3445γλλ----t p 内变化时，系统的可靠性情况如何呢？留给读者自己分析
例14 ⎪⎩
⎪⎨⎧===---0)0(';2)0(0
)1(1000222x x x dt
dx
x dt x d 解: 令 y 1=x ，y 2=y 1’ 则微分方程变为一阶微分方程组： ⎪⎩
⎪⎨⎧
==--==0
)0(,2)0()1(1000''21122
1221y y y y y y y y
1、建立m-文件vdp1000.m 如下： function dy=vdp1000(t,y) dy=zeros(2,1);
dy(1)=y(2);
dy(2)=1000*(1-y(1)^2)*y(2)-y(1);
2、取t0=0，tf=3000，输入命令： [T,Y]=ode15s('vdp1000',[0 3000],[2 0]); plot(T,Y(:,1),'-')
3、结果如图
程序见数学建模高教出版社数学建模与
数学实验To Matlab （ff4）
导弹追踪问题
设位于坐标原点的甲舰向位于x 轴上点A(1, 0)处的乙舰发射导弹，导弹头始终对准乙舰.如果乙舰以最大的速度v0(是常数)沿平行于y 轴的直线行驶，导弹的速度是5v0，求导弹运行的曲线方程.又乙舰行驶多远时，导弹将它击中？ 解法一（解析法）
假设导弹在t 时刻的位置为P(x(t), y(t))，乙舰位于),1(0t v Q
由于导弹头始终对准乙舰，故此时直线PQ 就是导弹的轨迹曲线弧OP 在点P 处的切线，
即有 x y
t v y --=1'0
即 y y x t v +-=')1(0 (1)
又根据题意，弧OP 的长度为AQ 的5倍， 即
t v dx y x
00
25'1=
+⎰
(2)
由(1),(2)消去t 整理得模型: 21
(1)"1' (3)5
x y y -=+ 初值条件为: 0)0(=y 0)0('=y
解即为导弹的运行轨迹: 24
5
)1(125)1(8556
54+-+--=x x y
当1=x 时245=y ，即当乙舰航行到点)24
5
,1(处时被导弹击中.
被击中时间为:0
0245v v y t ==. 若v 0=1, 则在t=0.21处被击中. 轨迹图见程序chase1 解法二(数值解)
令y1=y,y2=y1’，将方程（3）化为一阶微分方程组。

21(1)''15x y y -=+ ⇒ 12
2
21'1'1/(1)5y y y y x =⎧⎪
⎨=+-⎪⎩
1.建立m-文件eq1.m
function dy=eq1(x,y) dy=zeros(2,1); dy(1)=y(2);
dy(2)=1/5*sqrt(1+y(1)^2)/(1-x);
2. 取x0=0，xf=0.9999，建立主程序ff6.m 如下： x0=0，xf=0.9999
[x,y]=ode15s('eq1',[x0 xf],[0 0]); plot(x,y(:,1),’b.') hold on y=0:0.01:2; plot(1,y,’b*')
结论: 导弹大致在（1，0.2）处击中乙舰 程序数学建模高教出版社数学建模与数学实验To Matlab(ff6)
六.习题
1. 实验目的：利用实际例子筛选有效算法。

内容：将不同算法用于下面的方程，选出你最满意的方法。

也可以自己设计新方法，看是否比我们所讲方法好。

1).22()/,(1)1y y xy x y '=+=, 计算到x=2.
2). 21
4,(2)12
y x y x y '=+-=-, 计算到x=4. 3). 1(sin ),(2/)/2y y y y x y ππ-'=-=, 计算到x=2. 4).sin ,(1)1y y x y '==, 计算到x=2.
提示：先确定“满意”的尺度和标准，可以比较算法的效率、精确度等。

2. 实验目的：检查各种算法的长期行为。

内容：给定方程组()(),()(),(0)0,(0).
x t ay t y t bx t x y b '=⎧⎪'=-⎪
⎨=⎪⎪=⎩的解是x −y 平面上一个椭圆，利用你已经知道的算法，取足够
小的步长，计算上述方程的轨道，看那种算法能够保持椭圆轨道不变。

(计算的时间步要足够多)。

3.实验目的：初步认识刚性微分方程。

内容：任选一种你知道的显式方法，取不同步长，求解
2000999.751000.25
(0)0,(0)2dy u v dx dv u v dt
u v ⎧=-++⎪⎪⎪=-⎨⎪==-⎪⎪⎩
4. 实验课题（二）变步长龙格—库塔法
用变步长四阶标准龙格—库塔法解下列初值问题，并与准确解进行比较。

取h = 0.1，允许误差取ε=10-4及ε=10-6分别计算。

（1）⎪⎩⎪⎨⎧=≤<+='1)0(20,219.0y x x y y （2）⎪⎩⎪⎨⎧=≤<+-='1
)0(20,21y x x xy y 准确解：（1）45.0)21(-+=x y ； （2）212)1(-+=x y 。

5. 实验一 求解下列微分方程（组）
1．简单微分方程
;0)1(2)1(,0)0(,42)1)(2;
1)1(,2)0(,2)1='-=-=''+==-=-''y y y y y x y y x y y 2. 特殊微分方程
l ）Duffing’s 方程．03=-+u u εμ．可以变形为，3,x x y y x ε+-==，给出其相平面
图，参数ε的取值可以考虑-1/4,0.1和1/4等．
2）单摆运动方程0sin 8=+u μ，用图形表示其解．
3．微分方程组
1）1） 线性微分方程组X ’＝AX ，其中矩阵A 的定义如下:
⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡-⎥⎦⎤⎢⎣⎡-1214340106312 2)非线性微分方程组
||4,4);
,)2xy x y y xy y x x b x y x y x a +-=-==--=
6. 将下列二阶方程化为一阶方程组
（1）
（2）
7. 1926年荷兰数学家范德堡(van der Pol)构造了研究无线电热离子振荡器性态的模
型：
其中
是正的物理常数，这是电路理论中用到的非线性微分方程，如果引入新变
量： )3
(,3
21y y y y y y -+'==μ 可以将范德堡方程化为一阶常微分方程组的形式，试利用新变量推导一阶常微分方程组。

8. 分别用欧拉法和改进的欧拉法解下面初值问题⎩⎨⎧-='=<≤=-''1
)1()1(3.11,023y y y y 取h = 0.1计算，并与准确解2
1-=
x y 相比较。

9实验课题（四）应用题
一个初始重量为1350公斤(包括燃料1080公斤)的小火箭，点燃后垂直向上运动，火箭中的燃料以每秒18公斤的常速燃烧，提供3150公斤的动力。

假定大气阻力与速度平方成正比D =kv 2，且设空气动力阻力系数 k =0.039公斤⋅秒/米2。

试求燃料烧尽之前的位移、速度、加速度的变化规律。

方案I 用变步长四阶标准龙格—库塔法、步长h = 0.1。

方案II 用哈明方法，步长h = 0.1。

提示：设火箭位移为y 。

作用在火箭上的力有：推力F = 3150; 时
刻t 时火箭的重量引起的重力W =1350-18t ；大气阻力D=0.039(y ')2。

根据牛顿第二定律得运动微分方程 2)(y w kg g w gF y '--='' 初始条件y (0)= y '(0)=0。

求解范围；t ∈[0,T]，其中T = 1080/18=60（秒）。