专题九隐式龙格-库塔法.pdf

点
c1
,L
cs
就是拉道求积节点，它分别为
Ps*
(
t
)
+
P* s −1
(t
)
的零点或
Ps*
(
t
)
−
Ps
* −1
(
t
)
的
零点．这类方法计算节点及系数 c,b, A 的步骤是：
(1)求
Ps*
(t
)
+
P* s −1
(t
)
=
0
或
Ps*
(t
)
−
P* s −1
(t
)
=
0
的实根，得到
c1 ,L cs
，前者含
c1 = 0 ，后者含 cs = 1．
对(9.1)中方程从 tn 到 tn+1 积分可得
∫ ( ) y (tn+1 ) − y (tn ) =
f tn+1
tn
t, y (t ) dt
若用 s 个节点的求积公式计算右端积分，即
∫ ( ) ∑ ( ) ( ) ( ) f tn+1 t, y t tn
s
dt ≈ h bi f tn + cih, y tn + cih
由此得到
∑s
i =1
bicik −1
=
1 k
,k
= 1, 2,Lq
(9.9)
当 f (tn + hu) = (hu)q 时，求积公式(9.8)不精确成立，此时
( ) ( ) y tn+1 − yn+1 = O hq+1
若求积节点 ci (i = 1, 2,L, s) 已经给定，对(4.7)从 tn 到 tn,i = tn + cih 求积，则得
如果当 j ≥ i 时 aij = 0 ，公式(4.2b)中 ki 可由 k1,L, ki−1 得到，此时公式(9.2a)和
(9.2b)为显式方法，若当 j > i 时，aij = 0 而 aii = 0 ，公式称为半隐式 RK 方法，此
时(9.2b)右端 f 含 ki ，因此要求得 ki ，需要解 m 阶非线性方程组（假定 f 关于 y 为
∫ ∑ =
h
⎡ ⎢⎣
1 0
f
(tn
+ hu) du −
s
bi
i =1
f
(tn
+ cih)⎤⎥⎦
∫ ∑ 1 0
f
(tn
+ hu) du
≈
s i =1
bi
f
(tn
+ cih)
(9.8)
具有 q −1次代数精确度，即(9.8)对 f (tn + hu) = (hu)k−1 ， k = 1, 2,L, q 精确成立，
∫ ∫ y (tn + cih) =
( ) f tn +cih t dt = h ci
tn
0
f (tn + hu) du
用 c1, c2 ,L, cs 为节点的求积公式得到
∫ ∑ ( ) ( ) f ci 0
s
tn + hu du ≈ aij f
j =1
tn + c jh ,i = 1, 2,L, s
( ) 这样得到的隐式 RK 方法局部截断误差 ln+1 = ( ) y tn+1 − yn+1 = O h2s+1 ，故方法
的阶 p = 2s ，对 s 级 2s 阶的隐式 RK 方法(4.2)其节点及系数 c,b, A 可由以下 3 步 得到．
(1)求节点 c1, c2,L, cs ，由于[0,1] 上的高斯求积公式节点为区间[0,1] 的 s 阶勒 让德多项式 Ps* (t ) 的 s 个零点，而[−1,1] 上的 s 阶勒让德多项式为
j =1
为了建立求积系数 bi 与节点系数 ci 之间的关系，可通过对特殊初值问题积分得
到，考虑
y ' = f (t), y (tn ) = 0
(9.7)
由 tn 到 tn+1 积分得
( ) ∫ ( ) y tn+1
=
f tn+1
tn
t dt
令 t = tn + hu ，并用 s 个节点的求积公式得到
数 aij (i, j = 1, 2Ls) ．上面讨论表明利用数值求积公式建立(9.1)的隐式 RK 方法，
其系数应满足关系(9.9)及(9.11)，因为 s ≥ 1 ，故(9.9)(9.11)对 k = 1总是成立的，于 是有
s
s
∑ ∑ bi = 1, aij = ci ,i = 1, 2,L, s
s
∑ yn+1 = yn + h biki i =1
∑ ki =
⎛ f ⎜ tn + cih, yn + h
s
aij
k
j
⎞ ⎟
，
i
=
1,
2,L
s
⎝
j =1
⎠
(9.2a) (9.2b)
tn = t0 + nh (n = 0,1,L) ， h > 0 称为积分步长， yn 为(9.1)的精确解 y (tn ) 的数
, j = 1, 2,Ls
对 l = 1, 2,L, s 成立．
（9.16）
这样得到的公式很多，下面只给出 s = 2, p = 3 的两个公式．当 s = 2 时，由
P2* (t ) +
P1*
(t)
=
6t 2
−
4t
=
0
求得
c1
=
0, c2
=
2 3
，再由（9.9）求得
b1
=
1 4
, b2
=
3 4
，
最后由（9.11）求得 a11
i =1
(9.6)
这里 ci 为求积节点，bi 为求积系数，若令 yn = y (tn ) ，Yi ≈ y (tn + ) cih ，yn+1 ≈ y (tn+1 ) ，
则由(9.6)可得
s
∑ y (tn+1 ) − y (tn ) = h bi f (tn + cih,Yi ) i =1
若对(9.1)从 tn 到 tn + cih 分，则得
由以上 3 步求得的 s 级 2s 阶隐式 RK 方法如下： s = 1, p = 2 的隐式 RK 方法公式
11 22
1
yn+1
=
yn
+
hf
⎛ ⎜⎝
tn
+
1 2
h,
1 2
(
yn
+
yn+1
)
⎞ ⎟⎠
s = 2, p = 4 的隐式 RK 方法是
（9.13）
1− 3 26 1+ 3 26
1 4 1+ 3 46
专题九 隐式龙格-库塔法
9.1 龙格-库塔法的一般结构 显式龙格-库塔方法的绝对稳定域是有限区域，它不适用于求解刚性方程，
因此要研究隐式 RK 方法．对于初值问题
y ' = f (t, y), y (t0 ) = y0
(9.1)
y ∈ Rm , f : R × Rm → Rm , 隐式 RK 方法用公式可表示为
(2)由（9.9）求解得 b1,L,bs ．
(3)由（9.11）求解得 aij (i, j = 1, 2,L, s) ．
还可由下式（9.16）对 l = 1, 2,L, s 求得另一组 aij (i, j = 1, 2,L, s) ．
∑ ( ) s
i =1
bi
ci
l
a −1 ij
=
1 l
b
j
1− cjl
Yi = yn + h aij f tn + c jh,Yj ， i = 1, 2,Ls j =1
(9.4b)
计算时，先由(9.4b)求出Y1,LYs ，再由(9.4a)求出 yn+1 ．若假定在 tn 点(9.1)的解 y (tn )
与数值解 yn 相等，即 yn = y (tn ) ．由(9.4)求出 tn+1 = tn + h 的数值解 yn+1 与精确解
( ) Ps ( x) =
1 ds 2s s! dxs
⎡ ⎢⎣
x2 −1
s⎤ ⎥⎦
令 x = 2t −1，则得 Ps* (t ) = Ps (2t −1) ，当 s = 1, 2,3 时有
P1* (t ) = 2t −1 P2* (t ) = 6t2 − 6t +1 P3* (t ) = 20t3 − 30t2 +12t −1
ห้องสมุดไป่ตู้
(9.17b)
如 果 公 式 节 点 由 P2* (t ) − P1* (t ) = 6t2 − 8t + 2 = 2(3t −1)(t −1) = 0 求 得
c1
=
1 3
, c2
= 1 ，再由（9.9）及（9.11）求得以下公式
1 3
−
15 1212
1
3 4
1 4
3/ 4 1/ 4
(9.18)
称为 Radau IIA 公式． （III）基于洛巴托求积公式的隐式 RK 方法
p = 2s ．对显式方法，只要确定公式的阶 p 就可利用泰勒展开比较系数得到公式
的系数 b, c 及 A ．对隐式公式，理论上讲也可类似求出公式系数，但由于泰勒展
开及其计算的复杂，这样推导是不可取的．因此一般不采用泰勒展开的方法． 9.2 基于数值求积公式的隐式 RK 方法
为了得到隐式 RK 方法(9.2)中的系数 b, c, A ，仍从 m = 1的初值问题(9.1)出发，
Ps* (t ) = 0 的 s 个实根 c1,L, cs 即为公式节点．
(2)对 k = 1, 2,L, s ，由(4.9)得到关于 b1,L,bs 的线性方程组．求此方程组的解
则得 b1,L, bs ．
(3)分别对 i = 1, 2,L, s ，求(4.11)的解得 aij ( j = 1, 2,L, s) ．
1/ 2
1− 3 46
1 4
1/ 2
s = 3, p = 6 的隐式 RK 方法是
（9.14）
5 − 15 10
1 2
5 + 15 10
5 36 1+ 3 46
1/ 2
10 − 3 15 45
2 9
10 + 3 15 45
25 − 6 15 180
10 − 3 15 72
5 / 36
5 /18
1/18
5 /18
（II）基于拉道求积公式的隐式 RK 方法
（9.15）
若求积公式节点中包含求积区间[0,1] 的一个端点 c1 = 0 或 cs = 1，这种高斯型
求积公式就是拉道求积公式，它具有 2s − 2 次代数精确度．故 q = 2s −1 ，
( ) ln+1 = ( ) y tn+1 − yn+1 = O h2s ， p = 2s −1，可得 s 级 2s −1阶隐式 RK 方法，其节
非线性）．一般情形公式(9.2b)为隐式 RK 方法，在(9.2b)中， k1,L, ks 要通过求解
s × m 阶的方程组才能得到，计算比较复杂． 公式(9.2)还可改写为更便于推导的形式
s
∑ yn+1 = yn + h bi f (tn + cih,Yi ) i =1
(9.4a)
∑ ( ) s
值解，公式(9.2a)和(9.2b)称为 s 级的 RK 方法，其中系数 bi , ci 及 aij (i, j = 1, 2L s) 都
( ) 是实数， ci 称为节点， bi 称为权系数， A =
aij
∈ Rs×s 为系数矩阵，公式(9.2a)
s×s
和(9.2b)通常用以下的格式表示
或
(9.3)
(9.10)
此公式至少具有 s −1次代数精确度，即对 f (tn + hu ) = (hu)k−1 , k = 1, 2,L, s 精确成
立，由此得到
∑s
j =1
aijc j k −1
=
1 k
cik , i
= 1, 2,L, s
(9.11)
对 k = 1, 2,L, s 成立，它表明当求积节点 ci (i = 1, 2,L, s) 确定后，由(9.11)可求出系
y (tn+1 ) 之差 ln+1 ，当 h → 0 时，如果有整数 p > 0 ，使
( ) ( ) ln+1 = y tn+1 − yn+1 = O h p+1
(9.5)
则称方法(9.2)是 p 阶相容的， ln+1 为公式(9.2)的局部截断误差，对 s 级的显式 RK
方法的阶 p ≤ s ，当 s ≤ 4 时有 p = s ．而对隐式 RK 方法的阶 p 最高可达到
i =1
j =1
(9.12)
这就是隐式 RK 方法的相容性条件． 根据以上讨论可利用高斯型求积公式建立隐式 RK 方法，下面给出 3 种隐式
RK 方法． （I）基于高斯求积公式的隐式 RK 方法
若用 s 点的高斯求积公式．其代数精确度为 2s −1，即 q −1 = 2s −1，故 q = 2s ，
=
a12
=
0
及 a21
=
a22
=
1 3
，于是得公式
0
0
0
211 333
1/ 4 3/ 4
(9.17a)
这公式称为 Radau IA 公式，它不是 A − 稳定的，一般不用．另一个 A − 稳定的
Radau IA 公式，其 aij 由（9.16）求出的，公式为
0
1 4
−
1 4
2 15 3 4 12
1/ 4 3/ 4
当区间为[0,1] ，其两端均为求积节点，对 s ≥ 2 的洛巴托求积公式，代数精
于是有
∫ ∫ y (tn+1 ) =
tn+1 tn
f
(t ) dt
=
h
1 0
f
(tn
+ hu) du
s
∑ ≈ h bi f (tn + cih) = yn+1 i=0
若求积公式
( ) ∫ ( ) ∑ ( ) y tn+1
− yn+1 =
f tn+1
tn
t
s
dt − h bi f
i =1
tn + cih
∫ ( ) y (tn + cih) − y (tn ) =
f tn +cih
tn
t, y (t ) dt,i = 1,Ls
若右端积分用求积系数为 aij (i, j = 1, 2L s) 的求积公式，则得
∑ ( ) s
Yi = yn + h aij f tn + cih,Yj , i = 1,L s