刚性常微分问题的数值解法及编程

创新实验论文
题目：刚性常微分问题的数值解法
课程名称: _________ 创新实验______________ 学院：_____________ 理学院______________ 专业：数学与应用数学_________ 年级：__________ 应数131 _________ 学号：1307010239 234 236
姓名：袁蕊张蕾刘霖
指导教师：___________ 罗贤兵_____________ 2015 年07月14日
目录
第一章绪论 (3)
选题背景 (3)
刚性问题的算法 (3)
引言 (4)
第二章刚性问题 (5)
第三章预备知识 (8)
第四章计算实验 (15)
附页 (20)
第一章绪论
自然界和工程技术的很多现象，其数学模型是常微分方程(组)的初值问题，普通的常微分方程的数值解法已经比较成熟，理论比较完整，也有许多方法可供选择。

但，有一类常微分方程组，求解值时遇到相当大的困难，这类常微分方程组解的分量有的变化很快，有的变化缓慢，常常出现这种现象：变化快的分量很快趋于它的稳定值，而变化慢的分量缓慢趋于它的稳定值。

从数值解的观点看来，当变化快时应该用小步长积分，当变化快的分量已趋于稳定，或者说已经没有变化快的分量时就应该用较大的步长积分，但是理论和实际都说明，很多方法特别是显示方法的步长任不能放大，否者便出现数值不稳定现象，即误差急剧增加，已致掩盖了真解，使求解过程无法继续进行。

常微分方程组的这种性质叫做刚性，我们考虑一阶常微分方程初值问题的数值解法。

誉f(x,y) a*b (1.1)
.y(a)二y。

常微分方程的解能用初等函数、特殊函数或它们的级数与积分形式表达的非常之少，用解析办法只能求解线性常系数等特征类型的常微分方程。

在实际问题中归结出来的求解微分方程的方法只要依靠数值解法。

所谓数值解法，就是通过某种离散化办法，将微分方程转化为差分方程来求解。

求方程( 2-1 )的数值解，
即对一系列离散节点 a =X o <X1 < < x n <X n 1 <
建立求y(X n)的近似值y n的递推格式，由此求得解y(x)在各节点的近似值，
n=0, 1,2，…。

相邻两个节点的间距h二人q - x n称为步长，这时节点x^ x0 nh。

因
此，这样得到的数值解法也称为差分方法。

初值问题(1.1 )的数值解法，可区分为两大类：
(1)单步法：此类方法在计算Xm上的近似值y n 1时只用到了前一点禺上的信
息。

女口Euler法、Runge-Kutta法和Taylor级数法这就是这类方法的典型代表。

(2 )多步法：此类方法在计算时，除了需要x n点的信息外，还需要
x n4,Xnj,…前面若干个点上的信息。

线性多步法是这类方法的典型代表。

本文讨论的是几种隐式方法(向后欧拉，梯形公式，改进欧拉)和隐式
Run ge-Kutta
选取的步长h必须很小, 线性常微分方程初值问题第二章刚性问题
满足h v 1/100,才能保证绝对稳定性要求。

对于非
若初值问题是稳定的，即
V= f(x, y(x)),a 兰x 兰b,
、y(x°) = y°,y° E R m,
dy/dx <0.用欧拉法进行数值求解时，h应满足
cf卄
1+ —— h <1。

若
M
在方程组的情况下,
cf
cy
,h应满足h=2/M。

例如一阶常系数线性方程组
Uy
dx y
y(a)二y。

这里的A=a j ss，y=力山…y s T.记A的特征值为「，匕,…丄，对稳定的初值
问题应满足Re - <0.用欧拉法数值求解时，为了保证计算的稳定性，h的选取2
应满足h <------
max打
1&总maxRe 為
由下面的例子可以知道，当比值很大时，h很小，这时计算不熟
mn〔Re却很多，耗时很长，给实际计算带来了很大的困难。

例，某一物理现象可归结为一个线性方程组
其中x为时间变量，而A =
■ i
= -1, ' 2 ~ -40 1 ■ i ,，3 二-40 1 - <dx
[y(0)=(1,0,-1「
(1.2 )
C2119 -20"
19-21 20A的特征值分别为
<40—40 —
40』
i)。

方程(2-11 ) 的解为
也Ay
1 1
y_! (x )= -exp( 一2x) +3(cos40x +sin 40x)exp( -40x)
1 1
(1.3 ) * y2(x) = — exp(-2x) -一(cos40x +sin 40x)exp(」40x)
2 2
y3(x)= -(cos40x—sin 40x)exp(-40x)
我们对解式（1.3 ）编程作图，可以看出这组解在开始时刻变化激烈，随后逐渐进入稳态，对应于入2,打的分量在解中的作用随时间x的推移越来越显得无足
轻重。

解式(1.3 )程序和图像曲线(图1)如下面所示解程序：
h1=ezplot('(1/2)*exp(-2*x)+(1/2)*(cos(40*x)+si n(40*x))*exp(-40*x)'); axis([0 1 -1
1.5]);
hold on;
h2=ezplot('(1/2)*exp(-2*x)-(1/2)*(cos(40*x)+si n(40*x))*exp(-40*x)'); axis([0 1 -1
1.5]);
set(h2,'Color','r')
hold on;
h3=ezplot('-(cos(40*x)-si n(40*x))*exp(-40*x)');
axis([0 1 -1 1.5]);
set(h3,'Color','k')
hold on lege nd('y1','y2','y3')
图像：
图一
由于在开始的一段时间量x ,解曲线变化激烈，对方程进行数值求解时，自 然要
求数值有较高的精度，而对较大的时间量 x,解曲线变化缓慢，因此，对数
值方法的精度不必苛刻的要求，但就数值方法的稳定性而言，它并不随时间量 x 的大小而改变。

例如对初值问题（1.2 ）用欧拉折现法，步长必须满足 hv —^ =辺壯
0.035,这样小的步长对于较大的求解区间是难以接受的。

我们
maxl^ 40
看到，步长主要受特征值 九2 = k 3 =40/2的限制，如前所述，正是这两个特征
定义1.若线性系统式（2-11 ）中A 的特征值r 满足条件
（1） Re i ：： 0,i =0,1, ,m
max Re i
（2） R
=J_| i 〉〉1
称式（1.2 ）为刚性方程，比值R 为刚性比。

值，在微分方程中随时间量 x 的增大而显得作用减小，这种矛盾完全是比值
-{cos(40 x)-sm(40 x)) exp<40 x)
m ir n
第三章预备知识隐式方法
在（1.8 ）若用y n 近似替代y （X n ）、y n 1近似替代yX 彳），同样得到递推公式（1.7 ）.
如果在［X n ,X n 1］上对y " = f (X, y( X))积分，得
y （X n 1） - y （X n ）二 f （X, y （X ））dX. （ 1.9 ）
那么对上式右端积分用左矩形求积公式，并用
y n 近似替代y （X n ）、y n 1近似替代
y （X n 1），也可得递推公式（1.7 ）. 我们知道，在xOy 平面上，微分方程dy = f （x,y ）的解称为积分曲线，积分
dX
曲线上一点（x,y ）的切线斜率等于函数f （x,y ） D 的值。

如果D 中每一点，都 画上一条以f （x,y ）在这点的值为斜率并指向x 增加方向的有向线段，即在 D 上作出了一个由方程dy =: f (x,
1.1欧拉公式（显示）
考虑初值问题
i y )
y(x o ) =y 。

(1.4)
(1.5)
为了求得它在等距离散点 <X 1 < X 2 <…< X n <…上的数值解，首先将（1.1 ）离散
化。

设h 詁-x/i =1,2…），将式（1.4 ）离散化的办法有Taylor 展开法、数值
微分法及数值积分法如果在点X n 将y （X n 1） = y （X n ，h ）做Taylor 展开，得
y( X n 1)= y(X n ) hy (X n )
h 2 2! y ( n), \ ■ (X n ,X n .1) (1.6)
那么当h 充分小时, h 2 略去误差项-y （ n ），用yn 近似值代替皿）、
y n -1近似代替 y （X n 1），并注意到 y （X n ）二 f （X n , y （Xj ），便得
X = y(x o ),
Vn+ = y^hf (X n ,y n ),
n = 0,1； (1.7)
b a
其中X n =xo • nh,h =
.用（1.4）求解（1.1 ）的方法称为Euler 方法。

N 如果利用差商代似替代微商，那么可得
y(X n 1) - y(X n )
h
y (X n )二 f (X n , y (X n )). (1.8 )
y)确定的方向场，那么方程的解f=y (x).从几何dx
上看，就是位于此方向场中的曲线，它在所经过的每一点都与方向场的该点的方向一致。

从初始点P0(x o, y o)出发，过这点的积分曲线为y=y(x)，斜率为f(X0,y°).设在x =x°附近y(x)可用过P o点的切线近似表示，切线方程为y = y°• f (X o, y°)(x-X o).当x=X i时，y(xj 的近似值为y°• f (心y o)(% - x°)，并记为y i，这就是得到x =%时计算y(xj的近似公式
y 二y i f (洛，％)&-为).
当x =X2时，y(X2)的近似值为y i • f (x i, yJ(X2 - xj，并记为y•于是就得到当
X =X2的近似公式
丫2 *1 f (为，％)卜2 -X i).
重复上面方法，一般可得当X =X n彳的计算y(X n J的近似公式
y n i 二y n f (X n,y n)(X n i —X n).
如果h=X j -Xi/i =i,2…)，则上面公式就是(i.7 ).将P°,P,…,P N连续起来，就得到一条折线，所以Euler方法又称为折线法。

由公式(i.4 )看出，已知y o便可算出力.已知y i，便可算出七，如此继续下去，这种只用前一步的值y k便可计算出y「的递推公式称为单步法。

i.2向后欧拉公式(隐式)
若在(1.6 )中，其右边的积分由数值积分的右矩形公式近似，并用y(X n)替代y(X n 1)，则可得到
乂 = y(X o)
■:
y n 1 =y n hf(X n i, y n i)
并称(i.io)为后退Euler公式。

i.3梯形公式(隐式)
若在公式(i.6 )中，其右边积分用数值梯形积分公式近似，并用
y(X n)，y n i替代f (X n i)，则可得到梯形方法公式
y n替代(i.io)
y n替代
h y n1 - y n 2[f(x n ,y n^ f ( x n 1 , y n -1)] (1.I II
梯形方法同头退Euler 方法一样都是隐式单步法。

对于隐式方法，通常采用 迭代法。

对后退Euler 方法，先Euler 方法计算y . i ，并将它作为初值，即
y n 0)
i .=y 0 h f (x n , y n )，再把它代入(1.7 )的右端，便得到后退 Euler 方法的迭 代公式为
< (0) y n41 = y n + hf (x n , y n )
(k+) 丄-七/ (k) y n* =yn +hf (X n 卅 ynd
同样地，仍用Euler 方法提供初始值，梯形方法的迭代公式为
y n )i *n hf(X n ,y n )
h y n 「=yn +? [ f (X n , y n ) + f (禺杯 y^)]
1.4改进欧拉(隐式)
y n 1 =y n hf (X n ,y n )
< h ----
yn+ = y n +§(f (X n 』n ) + f (X n*,y n 卅)
1.5Runge-Kutta 法 龙格-库塔方法的基本思路是想办法计算f(x,y)在某点上的函数值，然后对 这些函数值做数值线性组合，构造出一个近似的计算公式；再把近似的计算公式 和解的泰勒展开式相比较，便得前面的若干项相吻合，从而达到较高的精度， 一
般的显示R-K 方法的形式如下：
r
y n41 = y n 灼 i k i
* =hf (X n , Y n )
(1.14)
i -4 k i =hf (X n +d h, y n +迟 0j k j ),(2Mi Mr) L jm
当式(1.14 )的布局截断误差达到O(h p1)时成公式(1.14 )为p 阶r 段R-K 方 法。

类似于显示方法我们可以得到 Runge-Kattu 隐式方法
r
y n 1 二 y n h'
'i k
i i# 常用隐式R-K 方法:
I r
匕=f X n +qh, y n +Gh 送丙 k j
< j 』
(1.12)
(1.13) k = 0,1,
(i= 1,2, r )
1级2阶中点公式：y n^ = y n+hf x n+^h,——)
l 2 2丿
2级2阶梯形公式：y n 1 =y n h f X n,y n f X. 1,『n 1
h yn^ =y n
+?(K i +K2 )
2级4阶R-K公式：<Ki = f冷+出3人小+丄知+ 3*厶3如
< 6 4 12
K2 = f x n+ A^h,y n +3_2<3知+丄如
6 12 4
在Matlab中，专门提供勒求解微分方程的ode函数，如ode45,ode23, ode113, ode15s, ode23s等等。

Ode求解器分为变步长和固定布场两种求解模式，其中，ode45是采用R-K法的变步长求解器，和它一样的还有ode23。

目前，ode45是
使用最多的求解函数，主要用于求解非刚性常微分方程。

（如果求解时长时间没
有相应，贝U方程是刚性的，可以换用ode23求解），ode函数的调用方式大都类
第四章实验计算
例如：刚性常微分方程
y (x) = -15y(x)
：y(0) = 1
这个初值问题的解析解为yx二e」5x，那么它在0.1处的解为? 解：
理论解:
Command Window
5uler_S (' 0, 1, (J. 13 100); axi S{[0 D. 1 0. 2 (LG): hold w;
色Let「(一 1 5#IE ■' j [Oj C. 1 ] J ; axis (CO 0. 1 0. 2 0,引)
由图可看出当取端点较小迭代次数较大时，误差较小。

2. 梯形公式(隐式)源代码
fun ctio n y=Euler_3(fu n,xO,yO,xN,N)
x=zeros(1,N+1);y=zeros(1,N+1);
x(1)=x0;y(1)=y0; h=(xN-xO)/N;
for n=1:N
x( n+1)=x( n)+h;
z0=y( n)+h*feval(fu n,x( n),y( n));
for k=1:3
z1=y( n)+h/2*(feval(fu n, x( n),y( n))+feval(fu n,x( n+1),z0))
J
if abs(z1-z0)<1e-3
break;
end
z0=z1;
end
y(n+1)=z1;
end
T=[x',y']
迭代次数10000次时
Euler_4 C z J 』0? 1』0. 1, 10000}
C^lunfit 9990 thicMh 9的百
d 盟翼9軸就：1林：
九2口4捕10：：如3枱 042X31$帥射曲3 也M 口仙0：第14篦4 0{餡斯4苏1H 別齋
C»lwuu 999T fhitujh 10^0] 4H 223»40：HT14SI 也就扯酋0. £»]»?]I04S3T1 «.?231634^433103 »3l30ieH03«：4
R-K 隐式4阶方法
1.源程序
function R=rk4(f,a,b,ya,yb,N);
T=a:h:b;
Y(1)=ya;
for j=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 ']
fun cti on z=f(x,y)
z=-15*y;
>>rk4( ' f ' ,0,0.1,1000)
0* 0950000000000001
0. 0960000000000001
0, 0970000000000001
0,0980000000000001
0.0990000000000001
0.1 总结：经多次试验发现向后欧拉，梯形公式，改进欧拉公式都是取端点越小 迭代次数越大时，越接近于准确值。

当放大端点值时也应增加迭代次数。

在取端 点值一样且迭代次数一样时，梯形公式和改进欧拉公式比向后欧拉更接近于准确 值。

但在计算端点较大时，他们不太实适用，因为耗时长。

Runge-kutta 解此刚 性常微分方程组比较使适用。

所以隐式 R-K 方法是解刚性常微分一般方法。

Coltuns $2 thr ougli &S
山255393B9S574285 C«l\uu^ 轴 throuch IQI
0L £?9938303S3I593
G 251^917199390^8 0.24：S46L482O8443
0. 2265IS097L36278 0.223142853627662
0*240521461184936 0.236940697931545 0.23341324329109 0,229938303631593 0.226515097136278 0. 223142853627662.
参考文献 舒适等著 数值计算方法 2009 著Matlab 数值分析与应用 2007 等著 刚性常微分方程初值问题的数值解法 1987
【1】黄云清
【2】张德峰。