刚性微分方程组隐式龙格库塔方法

毕业设计
题目：刚性系统的隐式RK方法
学院：数理学院
专业名称：信息与计算科学
学号： ************
学生姓名：**
指导教师：***
2016年05月15日
摘要
本文主要介绍单步隐式Runge – Kutta方法,简要的介绍了Gauss型隐式Runge – Kutta方法、Radau型隐式Runge – Kutta方法和Lobatto型隐式Runge – Kutta方法。

并利用这些基本的隐式Runge – Kutta方法来对刚性方程组进行数值求解，并将隐式Runge – Kutta方法与显式经典Runge – Kutta方法求解的结果进行对比，说明两种数值解法的优缺点。

关键词：刚性系统隐式Runge – Kutta方法单步方法Newton迭代法
Abstract
This paper mainly introduces the Implicit Runge-Kutta Methods and a simple description of Gauss implicit Runge-Kutta method ,Radau implicit Runge-Kutta method and Lobatto implicit Runge-Kutta method . These basic Implicit Runge-Kutta methods are used to solve the stiff equations. These implicit Runge-Kutta methods iare compared with the classical explicit Runge-Kutta method. This paper explain the advantages and disadvantages of the two kind of numerical methods.
Keywords: Stiff system Implicit Runge-Kutta method One step method Newton iterative method
目录
1.绪论 (1)
1.1刚性方程 (1)
1.2隐式RK方法的研究意义 (2)
1.3 RK方法的研究现状 (3)
2.单步RK方法的收敛性和稳定性 (5)
2.1单步RK方法的收敛性 (5)
2.2单步RK方法的稳定性 (6)
3.三类隐式RK方法 (8)
3.1引言 (8)
3.2 Gauss型隐式RK方法 (9)
3.3 Radau型隐式RK方法 (10)
3.4 Lobatto型隐式RK方法 (11)
4隐式RK方法的实现 (13)
4.1非线性系统的改进 (13)
4.2简化的Newton迭代法 (13)
5数值实验与结果分析 (15)
参考文献 (18)
附录 (21)
1.绪论1.1刚性方程
对于一般的线性常系数系统
y′=Ay+φ(t) A为m×m的矩阵，特征值为λi(i=1,2,⋯,m)。

定义1[23]若一个系统满足
（1）Reλi<0, i=1,2,⋯,m
（2）max
i |Reλi|min
i
|Reλi|
⁄=R≫1
其中R为刚性比，则这个系统称为刚性系统。

定义2[27]若线性系统
y′=Ay x∈[0,T]
或非线性系统
y′=f(x,y) x∈[0,T]
的矩阵A或Jacobi矩阵ðfðy
⁄的特征值λi满足
max
1≤i≤m
|Reλi|≫1
则其是刚性的。

定理1（解的存在性与唯一性）
（1）对于所有(x,y)∈D，函数f(x,y)是连续的；
（2）对于任何(x,y),(x,y∗)∈D，存在常数L，是函数满足
‖f(x,y)−f(x,y∗)‖≤L‖y−y∗‖则初值问题
{ẏ=f(x,y)(a≤x≤b) y(a)=η
有唯一解。

其中y=(y1,y2,⋯,y m)T，D={(x,y)|a≤x≤b,−∞<a<b<∞;−∞<y i<∞,i=1,2,⋯,m }。

其中L被称为Lipscℎitz常数
定义3 如果一个常微分系统的Lipscℎitz常数L很大(大于20)，则它是刚性的。

1.2隐式RK方法的研究意义
在常微分方程及常微分方程组的数值解法中，Runge – Kutta方法是目前应用最为广泛的数值解法之一，同时又具有误差小，精度高的特点。

尽管显式Runge – Kutta方法能够非常准确、快速的给出大部分常微分方程组的数值解。

但是在化学、自动控制电力系统等领域中，会出现一些病态的常微分方程组，也就是刚性方程组。

刚性方程组对于数值解法的稳定性要求苛刻，比如方程组
{y1=−0.01y1−99.99y2, y1(0)=2
y2=−100y2, y2(0)=1将其表示为矩阵形式：
[y1
y2]=[
−0.01−99.99
0−100
]⌈
y1
y2⌉
令
A=[−0.01−99.99
0−100
]
发现A特征值为：λ1=−100，λ2=−0.01，刚性比s=|λ1|
|λ2|
=10000≫1。

方程组的解为：
{y1(x)=e−100x+e−0.01x y2(x)=e−100x
解由快瞬态和慢瞬态两部分构成。

由于慢瞬态的部分，y1(x)衰减变得十分缓慢。

当自变量变到x=391时，函数值还未下降到初值的1%，求解区间至少取为(0,391)。

另一方面，由于快瞬态的部分，y2(x)衰减的非常快，因此步长要取得非常小。

从绝对稳定性的方面来看，如果用四阶显式经典RK方法求解，绝对稳定区间要求λℎϵ(−2.78,0)，则要求ℎ< 0.0278。

这样，在(0,391)上就要计算14 065步，计算量巨大，因此计算区间(0,391)内的解时，舍入误差积累会特别严重。

例如取求解区间为[0,1]，用不同步长ℎ来计算y1(1)和y2(1)的值。

利用四阶显式经典RK方法求解如下：
ℎy1(1)y2(1)
0.04 2.9802322e+17 2.9802322e+17
0.029.9004983e-01 1.3929556e-24
快瞬态慢瞬态
0.01 9.9004983e -01 2.5300364e -43 0.001 9.9004983e -01 3.7204130e -44 0.0001 9.9004983e -01 3.7200760e -44 真值
9.9004983e -01
3.7200760e -44
很显然，保留八位有效数字的情况下，要保持良好的精度，步长要取得非常小，这就增加了计算量。

而随着步长的加大，误差也会越来越大。

当步长加大到绝对稳定与之外时（即ℎ>0.0278），计算结果就完全不可信了。

对于刚性方程组，显式方法已经远远不能胜任了，一般采用绝对稳定性更好的方法（如隐式Runge – Kutta 方法）进行求解，本文采用单步隐式Runge – Kutta 方法，而对于隐式方法中的级值得求解，本文采用Newton 迭代法。

1.3 RK 方法的研究现状
研究基于标准模型方程的Runge – Kutta 方法的常见形式为：
{
y i+1=y i +ℎ∑αj k j r
j=1
k 1=f (x i ,y i ) k j =f(x i +λi ℎ,y i +ℎ∑μjl k l j−1
l=1
) (j =2,3,…r)
（显式）
{
y i+1
=y i +ℎ∑αj k j
r
j=1
k j =f(x i +λi ℎ,y i +ℎ∑μjl k l r
l=1
) (j =1,2,3,…r)
（隐式）
因为Runge – Kutta 方法是比较成熟的常微分方程数值解法。

所以如今主要是对于经典的Runge – Kutta 方法进行完善和扩充，在一定的条件下，提高级数以提高精度。

或者是将Runge – Kutta 方法与某些领域结合使用。

在1994年，费景高[1]给出了一种显式Runge – Kutta 并行方法，从而实现Runge – Kutta 方法在多处理机上的应用。

1997年，Enenkel 和Jackson [2]实现了Runge – Kutta 方法的对角隐式
并行改进。

1999年，廖文远和李庆扬[5]给出了一类求解刚性常微分方程的半隐式多步Runge – Kutta方法。

2000年，张诚坚和余洪兵[3]针对非线性延迟系统构造了一类并行预校算法，给出其算法的局部误差估计，数值实验表明该算法是有效的，且具一定的可比性。

2003年，李爱雄[4]通过对传统单支方法的计算格式进行改造，得到了解DDEs的两类单支并行算法单支并行预校算法和单支并行算法，并对方法的收敛性和稳定性做出了分析。

2008年，庞丽君和朱永忠[6]给出了一类随机微分方程Runge – Kutta方法的指数稳定性。

2.单步RK 方法的收敛性和稳定性
2.1单步RK 方法的收敛性
对于常微分初值问题
{ẏ=f (x,y ) (a ≤x ≤b )
y (a )=η
(1)
的单步显式公式
{
y i+1=y i +ℎφ(x i ,y i ,ℎ) (i =0,1,⋯,n −1)
y 0=η
(2)
局部截断误差可以表示为
R i+1=y (x i+1)−[y (x i )+ℎφ(x i ,y i ,ℎ)] (i =0,1,⋯,n −1) (3) 定理2[16]：
设y(x)为(1)的解，{y i }i=0n 为(2)的解。

如果： （1）存在常数c ，使得
|R i+1|≤cℎp+1 (i =0,1,2,⋯,n −1)
（2）存在a >0，使得
max (x,y )∈D σ0≤a≤ℎ
|ðφ(x,y,ℎ)
ðy
|≤L 其中D σ={(x,y )|a ≤x ≤b,y (x )−σ≤y ≤y (x )+σ} 记c 0=c
L [e L (b−a )−1]，则当ℎ≤min {a,√σ
c 0
p
}时，有
E(ℎ)≤c 0ℎp
证:由(3)得
y (x i+1)=y (x i )+ℎφ(x i ,y(x i ),ℎ)+R i+1 (i =0,1,⋯,n −1) (4) 将(4)与(2)相减
y (x i+1)−y (x i )=y (x i )−y i +ℎ[φ(x i ,y (x i ),ℎ)−φ(x i ,y i ,ℎ)]+
R i+1 (i =0,1,⋯,n −1)
由y (x 0)−y 0=0知道，在i =0时，|y (x i )−y i |≤c 0ℎp 成立。

现在假设0≤i ≤k −1时也是成立的。

在ℎ≤√σ
c 0p
时有：
|y (x i )−y i |≤c 0(√σ
c 0
p )p
=σ (i =0,1,⋯,n −1)
进而
|φ(x i ,y (x i ),ℎ)−φ(x i ,y i ,ℎ)| =|ðφ(x i ,ηi ,ℎ)ðy (y (x i −y i ))|
≤L |y (x i −y i )| 0≤i ≤k −1
其中ηi 是y (x i )和y i 之间的数。

于是定理结合条件与(4)式，可以得到
|y (x i+1)−y i+1|
≤|y (x i )−y i |+ℎ |φ(x i ,y (x i ),ℎ)−φ(x i ,y i ,ℎ)|+| R i+1| ≤|y (x i )−y i |+Lℎ|y (x i )−y i |+cℎp+1
=(1+Lℎ)|y (x i )−y i |+cℎp+1 0≤i ≤k −1
从而
|y(x k −y k )|≤(1+Lℎ
)k |
y (x 0)−y 0|+(1+Lℎ)k −1(1+Lℎ)−1
cℎp+1
因为y (x 0)−y 0=0及(1+Lℎ)k ≤e Lkℎ≤e L(b−a)，得到
|y (x k )−y k |≤c
L
[e L (b−a )−1]ℎp =c 0ℎp
所以当i =k 时|y (x i )−y i |≤c 0ℎp 也是成立的。

（证毕）
对于单步显式格式而言，当f(x,y)和
ðf(x,y)ðy
在D σ内连续时，那么定理1的条件
（2）总是满足的。

所以满足上述条件的单步显式公式，只要也满足相容性条件
φ(x,y (x ),0)=f(x,y(x))
那么它一定是收敛的。

2.2单步RK 方法的稳定性
定义4 对于初值问题(1),若{y i }i=0n 是由(2)得到，{z i }i=0n
是下面扰动问题的解：
{z i+1=z i +ℎ[φ(x i ,z i ,ℎ)+δi+1] (i =0,1,2,⋯,n −1)z 0=η+δ0
|δi|≤ε时，如果存在正常数C，ε0及ℎ0，使所有的ε∈(0,ε0]，ℎ∈(0,ℎ0]，当max
0≤i≤n
有
|y i−z i|≤Cε
max
0≤i≤n
就称(2)是稳定的或者零稳定的。

定理3[16]在定理1的条件下，单步显式公式(2)是稳定的。

3.三类隐式RK 方法
3.1引言
对于初值问题(1)，隐式Runge – Kutta 方法的一般形式为
{
y i+1
=y i +ℎ∑b j k j
r
j=1
(5)
k j =f (x i +c i ℎ,y i +ℎ∑a jl k l r
l=1
) (j =1,2,3,…r ) (6) 其中，x i =x 0+iℎ，n =0,1,2,⋯，为数轴上的离散点列；ℎ为步长，y i 为解y(x i )的近似值；c 1，c 2，…，c r 称为隐式Runge – Kutta 方法的步长；b 1，b 2，…，b r 为权系数。

令A =(a jl )，j =1,2,3,…r 。

称A 为系数矩阵，因此上式可以简化表示为
使用Butcℎer 可以看到，如果A 是一个主对角元素均为零的下三角矩阵，那么上表就可以表示一个显式Runge – Kutta 方法。

如果A 是一个主对角线非零的下三角矩阵，相应的Runge – Kutta 方法就是半隐式方法，(5)式等号左右就含有级值k 1，k 2，⋯，k r ，计算k i 时要解一个包含r 个未知量k i 的非线性方程组。

如果A 是满足A ij (i ≤j)不全为零，则相应的方法就是隐式方法。

若 系统是m 维的，对于隐式Runge -Kutta 方法实现的每一个递推步，均需要求解一个m ×r 维的非线性方程组，一般用牛顿迭代法求解。

例如
这是Kutta得到的3级3阶显式
与
分别是3级4阶的半隐式Runge – Kutta公式和隐式Runge – Kutta公式。

要求出具体的Runge – Kutta公式，一般有两种办法。

一种是将(6)在点(x i,y i)处进行泰勒展开，并且代入到(5)中，再与y(x i+ℎ)在x i处的泰勒展开式进行对比，从而确定参数c，b，A。

另一种方法就是将微分方程转化为等价的积分方程，用数值积分得到表达式，再进行对比得到参数。

基于后一种方法，可以构造得到以下三种不同的隐式Runge-Kutta方法。

3.2 Gauss型隐式RK方法[17]
设c1，c2，…，c s为P s(2c−1)=0的根，这里P s(t)是[0,1]上的s次Legendre多项式，0<c i<1，i=1,⋯,s。

求s级的2s阶的Gauss型Runge – Kutta公式的参数c，b，A需要以下步骤:
（1）求出s次的Legendre多项式P s(2c−1)=0的s个零点c1，c2，…，c s；
（2）由线性方程组
∑b j c j k−1
s
j=1
=1k⁄ k=1,⋯,s
确定系数b j，j=1,⋯,s；
（3）计算系数矩阵
A=CW−1
在这个基础上，就给出了s=1,2,⋯,5，p=2s的一系列Gauss型隐式Runge – Kutta 方法。

定理4[9][10]如果一个数值方法应用于方程y′=λy,其中λ为复常数，则对于所有λ，Reλ<0和对于固定步长ℎ>0，当n→∞时，得到的数值解趋于零，则称这种方法是A稳定的。

定理5[9][10] s级Gauss型隐式Runge – Kutta方法是2s阶的，其稳定函数是(s,s)Padé近似且是A稳定的。

以下列出s=
3.3 Radau型隐式RK方法[17]
这种方法的参数c，b，A需要下面几个步骤求得：
（1）求多项式P s(t)−P s−1(t)的零点c1，c2，…，c s，并指定c1>0，c s=1；（2）计算系数b k，
b k=∫
P s(t)−P s−1(t)
(t−c k)[P s′(c k)−P s−1′(c k)] 1
dt k=1,⋯,s （3）计算系数矩阵A，
A=CW−1
例如s=3，p=5的RadauⅡA型Runge – Kutta方法之一为：
定理6[31]s级RadauⅠA型Runge – Kutta方法和s级RadauⅡA型Runge – Kutta方法是2s−1阶的，稳定函数是(s−1,s)次对角Padé近似。

这两种都是A稳定的。

3.4 Lobatto型隐式RK方法[17]
这种方法的参数c，b，A需要下面几个步骤求得：
（1）求多项式P s(t)−P s−2(t)的零点c1，c2，…，c s，并指定c1=0，c s=1；（2）计算系数b k，
b k=∫
P s(t)−P s−2(t)
(t−c k)[P s′(c k)−P s−2′(c k)] 1
dt k=1,⋯,s
（3）计算系数矩阵A，
A=D−1(W T)−1(N−C)T D
列出4级6阶RadauⅢA型隐式Runge – Kutta方法
定理7[31]s级LobattoⅢA，ⅢB，ⅢC型隐式Runge – Kutta方法是2s−2阶
的，LobattoⅢA和ⅢB型Runge – Kutta方法的方法的稳定函数是(s−1,s−1)对角Padé近似，LobattoⅢC型隐式Runge – Kutta方法是(s−2,s)次对角Padé近似。

所以这些方法都是A稳定的。

4隐式RK 方法的实现
4.1非线性系统的改进
对于s 级隐式隐式Runge – Kutta 方法
{
Y i
=y n +ℎ∑a ij f(x n +c j ℎ,Y j
)s
j=1 i =1,⋯,s (7)
y n+1=y n +ℎ∑b j f(x n +c j ℎ,Y j )s
j=1
(8)
令
z i =Y i −y (9)
于是(7)变为
z i =∑a ij f(x n +c j ℎ,y n +z j )s
j=1
(10)
当得到的(10)解z 1，⋯，z 1时，由显式(8)可以很快得到解y n+1。

若隐式Runge – Kutta 方法的系数矩阵式非奇异的，那么(10)可以写为
[z 1⋮z s ]=A [ℎf(x n +c 1ℎ,y n +z 1)
⋮ℎf(x n +c s ℎ,y n +z s )] (11)
所以(8)可以看成与
y n+1=y n +∑d i z i s
i=1
等价的，其中
(d 1,⋯,d s )=(b 1,…,b s )A −1 (12)
比如s =3，p =5的Radau ⅡA 型Runge – Kutta 方法中，d =(0,0,1)。

4.2简化的Newton 迭代法
对于一般的非线性微分方程，系统(10)必须用迭代的方法求解。

本文采用Newton 迭代法。

Newton 迭代法应用于系统(10)，每次迭代都需要一个系数矩阵
[ I −ℎa 11ðf ðy (x n +c 1ℎ,y n +z 1)⋯ℎa 1s ðf
ðy (x n +c s ℎ,y n +z s )⋮−ℎa s1ðf ðy (x n +c 1ℎ,y n +z 1)⋯ℎa ss ðf ðy (x n +c s ℎ,y n +z s )]
的线性方程组。

定义5 若s 阶矩阵B =(b ij )，m 阶矩阵A =(a ij )，且有
B ⊗A =[b 11A
⋯b 1s A ⋮
⋱⋮b s1A
⋯
b ss A
] 则称运算⊗是B 与A 的Kronecker 积。

为了简化(10)的计算，我们用近似值
J ≈
ðf
ðy (x n ,y n
) 代替Jacobi 矩阵ðf ðy (x n +c i ℎ,y n +z i )⁄的值，于是方程(10)的简化Newton 迭代法变为
(I −ℎA ⊗J )△Z (k )=−Z (k )+ℎ(A ⊗I)F(Z (k))
Z (k+1)=Z (k)+△Z (k) (13)
这里Z
(k+1)
=
(z 1(k ),⋯,z s (k)
)T 是解的第k 次近似，△Z
(k )
=(△
(Z 1(k )
),⋯,△
(Z s (k)))T
是增量，F(Z
(k)
)是F(Z
(k )
)=(f(x n +c 1ℎ,y n +z 1(k)
),⋯,f(x n +c s ℎ,y n +z s
(k)
))
T
的缩写。

每一次的迭代需要s 次由函数f 的求值和一个N ∙S 维线性方程组的求解。

因为矩阵(I −ℎA ⊗J )对于所有的迭代是相同的，所以其LU 分解在一个积分步内只要进行一次，进而减少了计算时间。

5数值实验与结果分析
问题：
{
y 1′=10−7×y 1+sin (x ) y 1(0)=1y 2
′=10
−3
×y 2+cos (x ) y 2(0)=1
写成矩阵形式就是：
(y 1
′y 2
′)=(
10−7
00
10−3)(y 1y 2
)+(sin (x )cos (x ))
很明显该方程组的刚性比R =10−3
10−7=104≫1,因此是个刚性方程组。

取步长为0.1，在区间[0,1]内分别用四级四阶经典显式Runge – Kutta 方法
{
y i+1=y i +1
6(k 1+2k 2+2k 3+k 4)
k 1=f (x i ,y i ) k 2=f (x i +12ℎ,y i +1
2ℎk 1) k 3=f (x i +23ℎ,y i +23
ℎk 2) k 4=f (x i +ℎ,y i +ℎk 3)
和二级四阶隐式Runge – Kutta 方法
{
y i+1=y i +1
2(k 1+k 2)
k 1=f(x i +3−√36ℎ,y i +14k 1ℎ+3−2√3
12k 2ℎ)k 2=f(x i +3+√36ℎ,y i +14k 2ℎ+3+2√312k 1ℎ)
进行求解。

得到结果如下 表一：
真值结果
x i y 1 y 2
0.0000000e+00 1.0000000e+00 1.0000000e+00 1.0000000e -01 1.0049958e+00 1.0999384e+00 2.0000000e -01 1.0199334e+00 1.1988893e+00 3.0000000e -01 1.0446635e+00 1.2958649e+00
4.0000000e-01 1.0789390e+00 1.3898974e+00
5.0000000e-01 1.1224175e+00 1.4800481e+00
6.0000000e-01 1.1746644e+00 1.5654174e+00
7.0000000e-01 1.2351579e+00 1.6451531e+00
8.0000000e-01 1.3032934e+00 1.7184598e+00
9.0000000e-01 1.3783901e+00 1.7846058e+00
1.0000000e+00 1.4596978e+00 1.8429313e+00
表二：
隐式结果
x i y1y2
0.0000000e+00 1.0000000e+00 1.0000000e+00
1.0000000e-01 1.0049958e+00 1.0999384e+00
2.0000000e-01 1.0199334e+00 1.1988893e+00
3.0000000e-01 1.0446635e+00 1.2958649e+00
4.0000000e-01 1.0789390e+00 1.3898974e+00
5.0000000e-01 1.1224175e+00 1.4800481e+00
6.0000000e-01 1.1746644e+00 1.5654173e+00
7.0000000e-01 1.2351579e+00 1.6451531e+00
8.0000000e-01 1.3032934e+00 1.7184598e+00
9.0000000e-01 1.3783901e+00 1.7846058e+00
1.0000000e+00 1.4596978e+00 1.8429313e+00
表三：
显式结果
x i y1y2
0.0000000e+00 1.0000000e+00 1.0000000e+00
1.0000000e-01 1.0049958e+00 1.0999336e+00
2.0000000e-01 1.0199334e+00 1.1988802e+00
3.0000000e-01 1.0446635e+00 1.2958521e+00
4.0000000e-01 1.0789390e+00 1.3898814e+00
5.0000000e-01 1.1224175e+00 1.4800297e+00
6.0000000e-01 1.1746645e+00 1.5653972e+00
7.0000000e-01 1.2351579e+00 1.6451319e+00
8.0000000e-01 1.3032934e+00 1.7184382e+00
9.0000000e-01 1.3783901e+00 1.7845846e+00
1.0000000e+00 1.4596978e+00 1.8429111e+00
很明显的看到，在保留八位小数的情况下，二级四阶隐式Runge – Kutta方法与真值无异，能够精确到小数点后七位以上。

而经典Runge – Kutta只能精确到小数点后四位。

因此对于刚性方程组，相同步长下，隐式Runge – Kutta方法的精度比显式Runge – Kutta方法高得多。

并且由于步长小，在实际的求解区间里面，涉及的递推次数少，从而函数的计算量就小。

参考文献
[1]费景高. 并行显式Runge—Kutta 公式的实现[J]. 计算机工程与设计, 1994 (5):
34-40.
[2] Enenkel R F, Jackson K R. DIMSEMs-diagonally implicit single-eigenvalue methods
for the numerical solution of stiff ODEs on parallel computers[J]. Advances in Computational Mathematics, 1997, 7(1-2): 97-133.
[3] 张诚坚, 余红兵. 非线性延迟系统的一类并行预校算法[J]. 华中理工大学学报,
2000, 28(11): 104-106.
[4] 李爱雄, 刘伟丰, 张诚坚, 等. 求解DDEs 的多导龙格库塔方法的渐近稳定性
Asymptotic stability of multi-derivative Runge-Kutta method for delay differential equations[J]. 华中科技大学学报(自然科学版), 2002, 6: 108-110.
[5] 廖文远, 李庆扬. 一类求解刚性常微分方程的半隐式多步RK 方法[J]. 清华大
学学报: 自然科学版, 1999, 39(6): 1-4.
[6]庞立君, 朱永忠. 类随机微分方程Runge. Kutta 方法的指数稳定性[J]. 河海大
学学报(自然科学版), 2008, 36(3).
[7]Curtiss C F, Hirschfelder J O. Integration of stiff equations[J]. Proc. Nat. Acad. Sci,
1952, 38(235): 1.
[8]Gear C W. 常微分方程初值问题的数值解法[J]. 1978.
[9] Butcher J C. Implicit runge-kutta processes[J]. Mathematics of Computation, 1964,
18(85): 50-64.
[10] Ehle B L. High order A-stable methods for the numerical solution of systems of
DE's[J]. BIT Numerical Mathematics, 1968, 8(4): 276-278.
[11] Burrage K, Butcher J C, Chipman F H. An implementation of singly-implicit
Runge-Kutta methods[J]. BIT Numerical Mathematics, 1980, 20(3): 326-340. [12] Shampine L F. Implementation of implicit formulas for the solution of ODEs[J].
SIAM Journal on Scientific and Statistical Computing, 1980, 1(1): 103-118. [13] Lindberg B. On smoothing and extrapolation for the trapezoidal rule[J]. BIT
Numerical Mathematics, 1971, 11(1): 29-52.
[14] Lindberg B. IMPEX 2-a procedure for solution of systems of stiff differential
equations[R]. CM-P0*******, 1973.
[15] Bader G, Deuflhard P. A semi-implicit mid-point rule for stiff systems of ordinary
differential equations[J]. Numerische Mathematik, 1983, 41(3): 373-398.
[16] 孙志忠, 袁慰平, 闻震初. 数值分析[M]. 东南大学出版社, 2002.
[17] 刘德贵, , 费景高. 刚性大系统数字仿真方法[M]. 河南科学技术出版社, 1996.
[18] Cockburn B, Shu C W. The Runge–Kutta discontinuous Galerkin method for
conservation laws V: multidimensional systems[J]. Journal of Computational Physics, 1998, 141(2): 199-224.
[19] Monovasilis T, Kalogiratou Z, Simos T E. Exponentially fitted symplectic runge–
Kutta–Nyström methods[J]. Appl. Math. Inf. Sci, 2013, 7(1): 81-85.
[20] Hochbruck M, Pazur T. Implicit Runge--Kutta Methods and Discontinuous
Galerkin Discretizations for Linear Maxwell's Equations[J]. SIAM Journal on Numerical Analysis, 2015, 53(1): 485-507.
[21] Papadopoulos D F, Simos T E. A modified Runge-Kutta-Nyström method by using
phase lag properties for the numerical solution of orbital problems[J]. Applied Mathematics & Information Sciences, 2013, 7(2): 433-437.
[22] Zhong X, Shu C W. A simple weighted essentially nonoscillatory limiter for Runge–
Kutta discontinuous Galerkin methods[J]. Journal of Computational Physics, 2013, 232(1): 397-415.
[23] Lambert J D, Lambert J D. Computational methods in ordinary differential
equations[M]. London: Wiley, 1973..
[24] Zhu J, Zhong X, Shu C W, et al. Runge–Kutta discontinuous Galerkin method using
a new type of WENO limiters on unstructured meshes[J]. Journal of Computational
Physics, 2013, 248: 200-220.
[25] Jayakumar T, Maheshkumar D, Kanagarajan K. Numerical solution of fuzzy
differential equations by Runge-Kutta method of order five[J]. International Journal of Applied Mathematical Science, 2012, 6: 2989-3002.
[26] Dimarco G, Pareschi L. Asymptotic preserving implicit-explicit Runge--Kutta
methods for nonlinear kinetic equations[J]. SIAM Journal on Numerical Analysis, 2013, 51(2): 1064-1087.
[27] Miranker W L, Miranker A. Numerical Methods for Stiff Equations: And Singular
Perturbation Problems[M]. Springer Science & Business Media, 2001.
[28] Hadi M, Anderson M, Husein A. Using 4th order Runge-Kutta method for solving a
twisted Skyrme string equation[C]//THE 4TH INTERNATIONAL CONFERENCE ON THEORETICAL AND APPLIED PHYSICS (ICTAP) 2014. AIP Publishing, 2016, 1719: 030005.
[29] Jimenez J C, Sotolongo A, Sanchez-Bornot J M. Locally linearized runge kutta
method of dormand and prince[J]. Applied Mathematics and Computation, 2014, 247: 589-606.
[30] Jimenez J C, Sotolongo A, Sanchez-Bornot J M. Locally linearized runge kutta
method of dormand and prince[J]. Applied Mathematics and Computation, 2014, 247: 589-606.
[31] Wanner G, Hairer E. Solving ordinary differential equations II[M]. Springer-Verlag,
Berlin, 1991.
附录
绪论程序：
a=0;b=1; %求解区间
f1=@(t,x,y)(-0.01*x-99.99*y);
f2=@(t,x,y)(-100*y);
h=0.0001; %步长
T=zeros(1,((b-a)/h)+1);
X=zeros(1,((b-a)/h)+1);
Y=zeros(1,((b-a)/h)+1);
T=a:h:b;
X(1)=2; %初值
Y(1)=1;
for j=1:(b-a)/h
k11=feval(f1,T(j),X(j),Y(j));
k12=feval(f2,T(j),X(j),Y(j));
k21=feval(f1,T(j)+h/2,X(j)+h/2*k11,Y(j)+h/2*k12);
k22=feval(f2,T(j)+h/2,X(j)+h/2*k11,Y(j)+h/2*k12);
k31=feval(f1,T(j)+h/2,X(j)+h/2*k21,Y(j)+h/2*k22);
k32=feval(f2,T(j)+h/2,X(j)+h/2*k21,Y(j)+h/2*k22);
k41=feval(f1,T(j)+h,X(j)+h*k31,Y(j)+h*k32);
k42=feval(f2,T(j)+h,X(j)+h*k31,Y(j)+h*k32);
X(j+1)=X(j)+h*(k11+2*k21+2*k31+k41)/6;
Y(j+1)=Y(j)+h*(k12+2*k22+2*k32+k42)/6;
R=[T' X' Y'];
end
save 0.0001结果.txt R –ascii
真值程序：
y1=exp(-100)+exp(-0.01);
y2=exp(-100);
R=[1 y1 y2];
save 真值结果.txt R –ascii
第五章程序：
syms x;
a=0;b=1;
h=0.1;
f1=dsolve('Dy=0.0000001*y+sin(x)','y(0)=1','x');
f2=dsolve('Dy=0.001*y+cos(x)','y(0)=1','x');
X=a:h:b;
Y1=double(subs(f1,{x},{X}));
Y2=double(subs(f2,{x},{X}));
R1=[X' Y1' Y2'];
xlswrite('结果',R1,'真值')
format long
f1=@(x,y1)(0.0000001*y1+sin(x));
f2=@(x,y2)(0.001*y2+cos(x));
a=0;b=1;
h=0.1;
X=zeros(1,(b-a)/h+1);
Y1=zeros(1,(b-a)/h+1);
Y2=zeros(1,(b-a)/h+1);
X=a:h:b;
Y1(1)=1;
Y2(1)=1;
syms k1 k2;
for j=1:(b-a)/h
x10=[1 1];
x20=[1 1];
tol=1e-10;
x1=x10(1);x2=x20(1);
y1=x10(2);y2=x20(2);
K1=feval(@(k1,k2)(k1-feval(f1,X(j)+h*((3-sqrt(3))/6),Y1(j)+h*(1/4)*k1+h*((3-2*sqrt( 3))/12)*k2)),x1,y1);
K2=feval(@(k1,k2)(k2-feval(f1,X(j)+h*((3+sqrt(3))/6),Y1(j)+h*((3+2*sqrt(3))/12)*k1+ h*(1/4)*k2)),x1,y1); % 迭代函数
G1=feval(@(k1,k2)(k1-feval(f2,X(j)+h*((3-sqrt(3))/6),Y2(j)+h*(1/4)*k1+h*((3-2*sqrt( 3))/12)*k2)),x2,y2);
G2=feval(@(k1,k2)(k2-feval(f2,X(j)+h*((3+sqrt(3))/6),Y2(j)+h*((3+2*sqrt(3))/12)*k1+ h*(1/4)*k2)),x2,y2); % 迭代函数
F1=[K1 K2];F2=[G1 G2];
J1=double(subs(jacobian([k1-feval(f1,X(j)+h*((3-sqrt(3))/6),Y1(j)+h*(1/4)*k1+h*((3-2 *sqrt(3))/12)*k2);k2-feval(f1,X(j)+h*((3+sqrt(3))/6),Y1(j)+h*((3+2*sqrt(3))/12)*k1+h* (1/4)*k2)],[k1 k2]),{k1,k2},{x1,y1}));
J2=double(subs(jacobian([k1-feval(f2,X(j)+h*((3-sqrt(3))/6),Y2(j)+h*(1/4)*k1+h*((3-2 *sqrt(3))/12)*k2);k2-feval(f2,X(j)+h*((3+sqrt(3))/6),Y2(j)+h*((3+2*sqrt(3))/12)*k1+h* (1/4)*k2)],[k1 k2]),{k1,k2},{x2,y2}));
x11=x10-F1/J1;
x21=x20-F2/J2;
while (norm(x11-x10)>tol)
x10=x11;
x1=x10(1);
y1=x10(2);
K1=feval(@(k1,k2)(k1-feval(f1,X(j)+h*((3-sqrt(3))/6),Y1(j)+h*(1/4)*k1+h*((3-2*sqrt(
3))/12)*k2)),x1,y1);
K2=feval(@(k1,k2)(k2-feval(f1,X(j)+h*((3+sqrt(3))/6),Y1(j)+h*((3+2*sqrt(3))/12)*k1+ h*(1/4)*k2)),x1,y1);
F1=[K1 K2];
J1=double(subs(jacobian([k1-feval(f1,X(j)+h*((3-sqrt(3))/6),Y1(j)+h*(1/4)*k1+h*((3-2 *sqrt(3))/12)*k2);k2-feval(f1,X(j)+h*((3+sqrt(3))/6),Y1(j)+h*((3+2*sqrt(3))/12)*k1+h* (1/4)*k2)],[k1 k2]),{k1,k2},{x1,y1}));
x11=x10-F1/J1;
end
while (norm(x21-x20)>tol)
x20=x21;
x2=x20(1);
y2=x20(2);
G1=feval(@(k1,k2)(k1-feval(f2,X(j)+h*((3-sqrt(3))/6),Y2(j)+h*(1/4)*k1+h*((3-2*sqrt( 3))/12)*k2)),x2,y2);
G2=feval(@(k1,k2)(k2-feval(f2,X(j)+h*((3+sqrt(3))/6),Y2(j)+h*((3+2*sqrt(3))/12)*k1+ h*(1/4)*k2)),x2,y2);
F2=[G1 G2];
J2=double(subs(jacobian([k1-feval(f2,X(j)+h*((3-sqrt(3))/6),Y2(j)+h*(1/4)*k1+h*((3-2 *sqrt(3))/12)*k2);k2-feval(f2,X(j)+h*((3+sqrt(3))/6),Y2(j)+h*((3+2*sqrt(3))/12)*k1+h* (1/4)*k2)],[k1 k2]),{k1,k2},{x2,y2}));
x21=x20-F2/J2;
end
Y1(j+1)=Y1(j)+h*(x11(1)+x11(2))/2;
Y2(j+1)=Y2(j)+h*(x21(1)+x21(2))/2;
end
R2=[X' Y1' Y2'];
xlswrite('结果',R2,'隐式结果')
f1=@(x,y1)(0.0000001*y1+sin(x));
f2=@(x,y2)(0.001*y2+cos(x));
a=0;b=1;
h=0.1;
X=zeros(1,(b-a)/h+1);
Y1=zeros(1,(b-a)/h+1);
Y2=zeros(1,(b-a)/h+1);
X=a:h:b;
Y1(1)=1;
Y2(1)=1;
for j=1:(b-a)/h
k11=feval(f1,X(j),Y1(j));
k12=feval(f2,X(j),Y2(j));
k21=feval(f1,X(j)+h/2,Y1(j)+h/2*k11);
k22=feval(f2,X(j)+h/2,Y2(j)+h/2*k11);
k31=feval(f1,X(j)+h/2,Y1(j)+h/2*k21);
k32=feval(f2,X(j)+h/2,Y2(j)+h/2*k21);
k41=feval(f1,X(j)+h,Y1(j)+h*k31);
k42=feval(f2,X(j)+h,Y2(j)+h*k31);
Y1(j+1)=Y1(j)+h*(k11+2*k21+2*k31+k41)/6;
Y2(j+1)=Y2(j)+h*(k12+2*k22+2*k32+k42)/6; end
R3=[X' Y1' Y2'];
xlswrite('结果',R3,'显式结果')。