[y（0）=（1,0,-1「
（1.2）
C21
19-20"
19
-2120
A的特征值分别为
<40
—40—40』
i）。
方程（2-11）
的解为
也Ay
11
y_!
（x）=-exp（一2x）+3（cos40x+sin40x）exp（-40x）
（1.3）
11
*y2（x）=—exp（-2x）-一（cos40x+sin40x）exp（」40x）
22
y3（x）=-（cos40x—sin40x）exp（-40x）
我们对解式（1.3）编程作图，可以看出这组解在开始时刻变化激烈，随后逐渐进入稳态，对应于入2,打的分量在解中的作用随时间x的推移越来越显得无足
轻重。
解式（1.3）程序和图像曲线（图1）如下面所示解程序：
h1=ezplot（'（1/2）*exp（-2*x）+（1/2）*（cos（40*x）+sin（40*x））*exp（-40*x）'）;axis（[01-11.5]）;
holdon;
h2=ezplot（'（1/2）*exp（-2*x）-（1/2）*（cos（40*x）+sin（40*x））*exp（-40*x）'）;axis（[01-11.5]）;
set（h2,'Color','r'）
holdon;
h3=ezplot（'-（cos（40*x）-sin（40*x））*exp（-40*x）'）;
axis（[01-11.5]）;
set（h3,'Color','k'）
holdonlegend（'y1','y2','y3'）
图像：
图一
由于在开始的一段时间量x,解曲线变化激烈，对方程进行数值求解时，自然要求数值有较高的精度，而对较大的时间量x,解曲线变化缓慢，因此，对数
值方法的精度不必苛刻的要求，但就数值方法的稳定性而言，它并不随时间量x
的大小而改变。
例如对初值问题（1.2）用欧拉折现法，步长必须满足hv—^=辺壯0.035,这样小的步长对于较大的求解区间是难以接受的。
我们
maxl^40
看到，步长主要受特征值九2=k3=40/2的限制，如前所述，正是这两个特征
定义1.若线性系统式（2-11）中A的特征值r满足条件
（1）Rei：
：
0,i=0,1,,m
maxRei
（2）R=J_|i〉〉1
称式（1.2）为刚性方程，比值R为刚性比。
第三章预备知识隐式方法
1.1欧拉公式（显示）
考虑初值问题
化。
设h詁-x/i=1,2…），将式（1.4）离散化的办法有Taylor展开法、数值
微分法及数值积分法如果在点Xn将y（Xn1）=y（Xn，h）做Taylor展开，得
y（Xn1），并注意到y（Xn）二f（Xn,y（Xj），便得
在（1.8）若用yn近似替代y（Xn）、yn1近似替代yX彳），同样得到递推公式（1.7）.
如果在［Xn,Xn1］上对y"=f（X,y（X））积分，得
y（Xn1）-y（Xn）二f（X,y（X））dX.（1.9）
那么对上式右端积分用左矩形求积公式，并用yn近似替代y（Xn）、yn1近似替代
y（Xn1），也可得递推公式（1.7）.
我们知道，在xOy平面上，微分方程dy=f（x,y）的解称为积分曲线，积分
dX
曲线上一点（x,y）的切线斜率等于函数f（x,y）D的值。
如果D中每一点，都画上一条以f（x,y）在这点的值为斜率并指向x增加方向的有向线段，即在D上作出了一个由方程dy=:
f（x,y）确定的方向场，那么方程的解f=y（x）.从几何dx
上看，就是位于此方向场中的曲线，它在所经过的每一点都与方向场的该点的方向一致。
从初始点P0（xo,yo）出发，过这点的积分曲线为y=y（x），斜率为f（X0,y°）.设在x=x°附近y（x）可用过Po点的切线近似表示，切线方程为y=y°•f（Xo,y°）（x-Xo）.当x=Xi时，y（xj的近似值为y°•f（心yo）（%-x°），并记为yi，这就是得到x=%时计算y（xj的近似公式
y二yif（洛，％）&-为）.
当x=X2时，y（X2）的近似值为yi•f（xi,yJ（X2-xj，并记为y•于是就得到当
X=X2的近似公式
丫2*1f（为，％）卜2-Xi）.
重复上面方法，一般可得当X=Xn彳的计算y（XnJ的近似公式
yni二ynf（Xn,yn）（Xni—Xn）.
如果h=Xj-Xi/i=i,2…），则上面公式就是（i.7）.将P°,P,…,Pn连续起来，就得到一条折线，所以Euler方法又称为折线法。
由公式（i.4）看出，已知yo便可算出力.已知yi，便可算出七，如此继续下去，这种只用前一步的值yk便可计算出y「的递推公式称为单步法。
i.2向后欧拉公式（隐式）
yn替代
若在（1.6）中，其右边的积分由数值积分的右矩形公式近似，并用
y（Xn）替代y（Xn1），则可得到
（i.io）
乂=y（Xo）
■:
yn1=ynhf（Xni,yni）
并称（i.io）为后退Euler公式。
i.3梯形公式（隐式）
yn替代
若在公式（i.6）中，其右边积分用数值梯形积分公式近似，并用
y（Xn），yni替代f（Xni），则可得到梯形方法公式
h
yn1-yn2[f（xn,yn^f（xn1,yn-1）]（1.
梯形方法同头退Euler方法一样都是隐式单步法。
对于隐式方法，通常采用迭代法。
对后退Euler方法，先Euler方法计算y.i，并将它作为初值，即
yn0）i.=y0hf（xn,yn），再把它代入（1.7）的右端，便得到后退Euler方法的迭
代公式为
<（0）
（1.12）
yn41=yn+hf（xn,yn）
（k+）丄-七/（k）
yn*=yn+hf（Xn卅ynd
同样地，仍用Euler方法提供初始值，梯形方法的迭代公式为
yn）i*nhf（Xn,yn）
（1.13）
k=0,1,
h
yn「=yn+?
[f（Xn,yn）+f（禺杯y^）]
1.4改进欧拉（隐式）
yn1=ynhf（Xn,yn）
yn+=yn+§（f（Xn』n）+f（Xn*,yn卅）
1.5Runge-Kutta法
龙格-库塔方法的基本思路是想办法计算f（x,y）在某点上的函数值，然后对这些函数值做数值线性组合，构造出一个近似的计算公式；再把近似的计算公式和解的泰勒展开式相比较，便得前面的若干项相吻合，从而达到较高的精度，一
般的显示R-K方法的形式如下：
r
yn41=yn灼iki
*=hf（Xn,Yn）（1.14）
i-4
ki=hf（Xn+dh,yn+迟0jkj）,（2MiMr）
Ljm
当式（1.14）的布局截断误差达到O（hp1）时成公式（1.14）为p阶r段R-K方法。
类似于显示方法我们可以得到Runge-Kattu隐式方法
r
yn1二ynh''iki
（i=1,2,r）
i#
常用隐式R-K方法:
1级2阶中点公式：
yn^=yn+hfxn+^h,——）
l22丿
2级2阶梯形公式：
yn1=ynhfXn,ynfX.1,『n1
hyn^=yn+?
（Ki+K2）
2级4阶R-K公式：
<6412
K2=fxn+A^h,yn+3_2<3知+丄如
6124
在Matlab中，专门提供勒求解微分方程的ode函数，如ode45,ode23,ode113,ode15s,ode23s等等。
Ode求解器分为变步长和固定布场两种求解模式，其中，ode45是采用R-K法的变步长求解器，和它一样的还有ode23。
目前，ode45是
使用最多的求解函数，主要用于求解非刚性常微分方程。
（如果求解时长时间没
有相应，贝U方程是刚性的，可以换用ode23求解），ode函数的调用方式大都类
似，以ode45为例，常用的语法格式有：
[T,Y]=ode45（odefun,tspan,yO）;
函数
问题类型
精确度
说明
ode45
非刚性
中等
米用算法为4-5阶Runge-Kutta法，大多数情况下首选的函数
ode23
非刚性
低
基于Bogacki-Shampine2-3阶Runge-Kutta公式，在精确要求不咼的场合，以及对于轻度刚性方程，ode23的效率可能好于ode45
ode113
非刚性
低到高
基于变阶次Adams-Bashforth-MoutlonPECE算法。
在对误差要求严格的场合或者输入参数doefun代表的函数本身计算量很大情况比ode45效率咼。
ode113可以看成一个多步解算器，因为它会利用前几次时间节点上的解计器当前时间节点的解。
因此它不适应与非连续系统。
ode15s
刚性
低到中
基于数值差分公式（后向差分公式，BDFs也叫Gear方
法），因此效率不是很高。
同ode113一样，ode15s也是一个多步计算器。
当ode45求解失败，或者非常慢，并且怀疑问题是刚性的，或者求解微分代数问题时可以考虑用ode15s
ode23s
刚性
低
基于修正的二阶Rosenbrock公式。
由于是单步解算器，当精度要求不咼时，匕效率可能会咼于ode15s。
他可以
解决一些ode15s求解起来效率不太咼的刚性冋题。
ode23t
适度刚
性
低
ode23可以用来求解微分代数方程。
ode23tb
刚性
低
当方程是刚性的，并且求解要求精度不高时可以使用。
第四章实验计算
例如：
刚性常微分方程
y（x）=-15y（x）
：
y（0）=1
这个初值问题的解析解为yx二e」5x，那么它在0.1处的解为?
解：
理论解:
1.向后欧拉（隐式）源程序
functiony=Euler_2（fun,x0,y0,xN,N）
x=zeros（1,N+1）;y=zeros（1,N+1）;x
（1）=x0;y
（1）=y0;
h=（xN-x0）/N;
forn=1:
N
x（n+1）=x（n）+h;
z0=y（n）+h*feval（fun,x（n）,y（n））;fork=1:
3
z1=y（n）+h*feval（fun,x（n+1）,z0）;ifabs（z1-z0）<1e-3break;
endz0=z1;
end
y（n+1）=z1;
end
T=[x',y']
functionz=f（x,y）
z=-15*y;
迭代100次
j电rrrn<»mrtpupwvtruw*■v
»Euler_2（#0,130,L100）,
Columns池throith卵
D.Z50O5159X29SE：
0.鬲血閘彌轩闾fl.：
5O1$23W6O5S59
Colukns99throujhlGl
1232525K00721430.2290miS54B»270.22570560；634C3
■■_■
迭代1OOOO次
fx»Euler.aC,0,1,04h10000）
Columns9997through9998
0.2232891920L19670.223255703657172
Columns9999thrmigh10000
0.223222*********0.22318S742OU328
Column10001
0.223155268724772
源程序作图比较：
（红色线是迭代的结果，蓝色线是准确值）functiony=Euler_2（fun,xO,yO,xN,N）
x=zeros（1,N+1）;y=zeros（1,N+1）;
x
（1）=xO;y
（1）=yO;h=（xN-xO）/N;
forn=1:
N
x（n+1）=x（n）+h;
zO=y（n）+h*feval（fun,x（n）,y（n））;
fork=1:
3
z1=y（n）+h*feval（fun,x（n+1）,z0）;
ifabs（z1-z0）<1e-3
break;
end
z0=z1;
end
y（n+1）=z1;
end
T=[x,y];
plot（x,y,'r'）
CommandWindow
5uler_S（'0,1,（J.13100）;
axiS{[0D.10.2（LG）:
holdw;
色Let「（一15#ie■'j[OjC.1]J;
axis（CO0.10.20,引）
由图可看出当取端点较小迭代次数较大时，误差较小。
2.梯形公式（隐式）源代码
functiony=Euler_3（fun,xO,yO,xN,N）
x=zeros（1,N+1）;y=zeros（1,N+1）;
x
（1）=x0;y
（1）=y0;h=（xN-xO）/N;
forn=1:
N
x（n+1）=x（n）+h;
z0=y（n）+h*feval（fun,x（n）,y（n））;
fork=1:
3
z1=y（n）+h/2*（feval（fun,x（n）,y（n））+feval（fun,x（n+1）,z0））
J
ifabs（z1-z0）<1e-3
break;
end
z0=z1;
end
y（n+1）=z1;
end
T=[x',y']
functionz=f（x,y）
z=-15*y;
迭代次数100次
Euler3（JzJ鼻0〉0.1j
100）;
0.0950000000000001
0.240521461194936
0.0960000000000001
0.236940697931545
0,0970000000000001
0.23341324329109
0.0980000000000001
0.229933303631593
0.0390000000000001
0.226515097136278
0.1
0.223142853627662,
迭代次数10000次时
Euler3Cz,0,1,0,1,10000）
C^luuis5HOthjoufhSJS6
AlM晰加354
i.zzn&isTiim
ColwuisSWTthiMchIKK）]
氐洱32融衲TIM
fl,2232M5925：
1⑸軒】W49犬1
■
叽229】63血皿啪
«.223IN16UOM；S
3.改进欧拉公式（隐式）源代码
functiony=Euler_4（fun,xO,yO,xN,N）
x=zeros（1,N+1）;y=zeros（1,N+1）;
x
（1）=x0;y
（1）=y0;h=（xN-xO）/N;
forn=1:
N
x（n+1）=x（n）+h;
z0=y（n）+h*feval（fun,x（n）,y（n））;
y（n+1）=y（n）+h/2*（feval（fun,x（n）,y（n））+feval（fun,x（n+1）,z0））;end
T=[x',y']
functionz=f（x,y）
z=-15*y;
迭代次数100次时
Euler4Cz\0;l30.b100）
迭代次数10000次时
Euler_4CzJ』0?
1』0.1,10000}
C^lunfit9990thicMh9的百
d盟翼9軸就：
1林：
九2口4捕10：
：
如3枱042X31$帥射曲3也M口仙0：
第14篦40{餡斯4苏1H別齋
C»lwuu999Tfhitujh10^0]
4h223»40：
HT14SI也就扯酋0.£»]»?
]I04S3T1«.?
231634^433103»3l30ieH03«：
4
R-K隐式4阶方法
1.源程序
functionR=rk4（f,a,b,ya,yb,N）;
T=a:
h:
b;
Y
（1）=ya;
forj=1:
N
K1=feval（f,T（j）+h*（3+3A（1/2））/6,Y（j）+k1*1/4+k2*（3+2*3A（1/2）））/12;
K2=feval（f,t（j）+h*（3-3A（1/3））/6,y（j）+k1*（3-2*3A（1/2））/12+k2*1/4）;
Y（j+1）=Y（j）+（k1+k2）/2;
end
R=[T'Y']
functionz=f（x,y）
z=-15*y;
>>rk4（'f',0,0.1,1000）
0*240521461184936
0.236940697931545
0.23341324329109
0,229938303631593
0.226515097136278
0.223142853627662.
0*0950000000000001
0.0960000000000001
0,0970000000000001
0,0980000000000001
.0990********
0.1
总结：
经多次试验发现向后欧拉，梯形公式，改进欧拉公式都是取端点越小迭代次数越大时，越接近于准确值。
当放大端点值时也应增加迭代次数。
在取端点值一样且迭代次数一样时，梯形公式和改进欧拉公式比向后欧拉更接近于准确值。
但在计算端点较大时，他们不太实适用，因为耗时长。
Runge-kutta解此刚性常微分方程组比较使适用。
所以隐式R-K方法是解刚性常微分一般方法。
【1】黄云清
【2】张德峰
【3】袁兆鼎
参考文献
舒适等著数值计算方法2009
著Matlab数值分析与应用2007
等著刚性常微分方程初值问题的数值解法1987
论文
成绩
评
语
年
月
日
升级会员 