第五章 常微分方程数值解
§ 5.1 问题的提出是其解析解。
，
、边值问题是其解析解。
、初始问题常微分方程定解问题
btaRtx
xbxxax
btatxf
dt
dx
ttRtx
xtx
tttxf
dt
dx
n
ba
n
,)(
2)-(5
)(,)(
),(
2
,)(
1)-(5
)(
),,(
1
0
00
0
。形式，称此为标准模式均可化为一般
),(
0),,,,( )(
txf
dt
dx
xxxtF n
)()1(
,
)()1(
12
2
1
2
2
1
21
2
2
2
tfxxx
dt
dx
x
dt
dx
dt
dx
xxx
tfx
dt
dx
x
dt
xd
，则引入例如
).(~)(
.,,,,,,
)(
1010
txtx
xxxttt
tx
nn
得近似值逼近的方法得到然后在利用数的近似值求在给定时刻
。数值解法就是得解析解对于非线性问题很难求
然后按精度要求截断。
，即或采用台劳展开的方法然后数值积分。
，
积分方程：）微分方程（初始问题算法的构造是
4)-(5 )()()(
3)-(5,,2,1,0)),(()()(
1-5 *
1
1
1
1
nnnnn
nnn
t
t
nn
txhtxtx
tth
ndtttxftxtx
n
n
则为多步算法。则称之为单步算法。否即算法是只用到前面一步信息定义：若计算
),,(
,
1
1
nnnnnn
nn
htxhxx
xx
.,,,
1)(
1
,0
1:1
)()(
1
)1()2(
1
)2(
1
)( ),(
1
)(
1)1(,0)0(
10),()1(
*
10
10
1
2
11
2
11
2
1
2
n
nn
iiiiiiii
iiiiiii
xxx
xx
h
x
ni
tfxxx
h
xxxx
h
xxx
h
txxx
h
tx
xx
ttfxxxx
由此解出例子：有限元法等。差分法的法很多，有差分方法、边值问题的算法构造方
本章只介绍初始问题的解法，原因是初始问题系统科学中研究的最重要的问题。
§ 5.2 初始问题的单步法
!!
2
)(
)),(()()(
~
)(
~
)(
),(,
,2,1,0
1
2
1
1
111
1
1
nn
nnnnn
nnn
nnnnnn
nnnn
h
x
T
ttxfhtxtx
txtxT
txfftth
n
fhxx
局部截断误差：
其中，
、显式欧拉法一、欧拉法
的李布希兹常数为其中，
总体误差：
fL
e r rLhTe r r
xtxLhxtx
xtxLhxtxxtx
txfttxfhxtx
xtx
xtxtxtx
xtxe r r
nnnn
nnnnn
nnnnnnn
nnnnnnn
nn
nnnn
nnn
)1(
)()1()(
~
)()()(
~
)),()),((()(
)(
~
)(
~
)(
~
)(
)(
11
11
11
11
1111
111
!! )(
)1(
2
1)1(
2
)1())1((
2
)1(
2
,
2
),,)(m a x (),:1,m a x (
1
)1(2
1
22
1
0
1
0
22
22
1
22
1
102
hoe r r
he
L
M
hL
hL
h
M
e r r
e r rhLhLh
M
e r rhLh
M
e r r
h
M
T
ttttxMnkhh
n
hLn
n
n
n
n
k
k
nn
n
nk
及则设以上说明显式欧拉法是一阶收敛的算法
2、隐式欧拉法
。二阶收敛的方法，即
）、改进的欧拉法（隐式一阶收敛的方法，即
)(
,2,1,0
)],(),([
3
)(
,2,1,0
),(
2
1
111
1
111
hoe r r
n
txftxfhxx
hoe r r
n
txfhxx
n
nnnnnnn
n
nnnnn
二、龙格 -库塔法（ RK)
我们的目的是构造更高阶的单步法！
阶方法。即此单步法为步法的整体误差则此单截断误差结论：若单步法的局部
phoe r r
hoT
p
n
p
n
).(
),(
1
1
1
( *) )()(
.:0),),((
)),(()),((
:0,10,,
)),(()()(
2
0
1
2
0
1
1
1
1
)o( hkctxtx
mjrrxfk
)o( hrrxfcdtttxf
miahatrtth
dtttxftxtx
m
m
j
jjnn
jjj
m
m
j
jjj
t
t
iininn
t
t
nn
n
n
n
n
令利用拉格朗日求积公式分析：由
( ** ),1
),)((
)(
)()()(
)()()(
),),((
1
1
1
1
2
1
11
111
00
mj
hatkbhtxfk
bkbhtx
hokaahtx
hkaarxrx
ttxfktr
jn
j
i
ijinj
ji
j
i
ijin
j
i
iiin
jjjjj
nnn
待定
，则取的显式单步法：
阶可构造出）可知：适当选取参数）和式（由式（ 1,,***?mbac jijj
m
j
jjnn kcxx
0
1
,1
),(
),(
1
1
0
mj
hatkbhxfk
txfk
jn
j
i
ijinj
nn
称此结构为 m+1级龙格 -库塔法。它最多是 m+1阶的。
一些常用的显式 RK方法
1、二级二阶 RK法改进的欧拉法（ p277),中点公式（ p278),休恩方法（ p278)等。
2、三级三阶 RK法休恩三阶方法（ P279),库塔三阶方法（ P279)等
3、四级四阶 RK法经典四级四阶 RK法（ P280),基尔四级四阶 RK法 (P280)
三、单步法的数值稳定性为复常数试验模型：稳定性时，都针对同一数值依赖性，在讨论方法的小很困难。为摆脱这种估计其大太依赖于由于时，将产生误差计算则用单步法
，即实际得到近似值产生误差若计算
,
),,(),,?(),,(
)],,?(),,([?
),,(
,?
,?)(
111
1
1
xx
txfhtxhtx
htxhtxhxx
xx
x
htxhxx
xx
xx
nnnn
nnnnnn
nnn
n
nnnn
nnn
nnnn
这样做的根据是,1）对实验模型数值不稳定的方法，不可用；
2）一般的初始问题在其解的存在区域内，可局部线性化转化为试验方程。
几个常用算法的数值稳定性
.)
2
0
))(
02
11,
)1(
)1(
1
1
1
限制性（还要满足算法数值稳定限制，而且不仅要满足算法收敛（步长
），（间：为实数时，绝对稳定区当绝对稳定区域式为用于试验模型的计算公显式定区域、显式欧拉法的绝对稳
h
hoh
h
h
h
xhxhxx
E u l e r
nn
nnnn
是绝对稳定的。
都稳定的算法，只要显然，凡是稳定的。称此算法是面的整个左半平面，则复平的绝对稳定区域包含了定义：一个算法如果它限制即可！（只需考虑满足算法收敛步长
），（间：为实数时，绝对稳定区当绝对稳定区域式为用于试验模型的计算公隐式定区域、隐式欧拉法的绝对稳
0h0,)R e (-A
))(
0
11,
1
1
1
1
2
1
111
A
h
hoh
h
h
h
x
h
xxhxx
E ul e r
nn
nnnnn
限制即可！（只需考虑满足算法收敛步长
），（间：为实数时，绝对稳定区当稳定的。是绝对稳定区域
）（
式为用于试验模型的计算公改进的稳定区域、改进的欧拉法的绝对
))(
0
,1
2
2
,
2
2
2
2
2
3
2
1
111
hoh
h
A
h
h
h
h
x
h
h
xxx
h
xx
E u l e r
nn
nnnnnn
的限制。还要考虑数值稳定性（
限制，而且（不仅考虑满足算法收敛步长
），（间：为实数时，绝对稳定区当
（绝对稳定区域
（
（
公式为法用于试验模型的计算经典四级四阶法的绝对稳定区域、经典四级四阶
)
78.2
0
))(
078.2
1)
!4
)(
!3
)(
!2
)(
1,
)
!4
)(
!3
)(
!2
)(
1
)
!4
)(
!3
)(
!2
)(
1
4
2
432
432
1
432
1
h
hoh
h
hhh
h
hhh
h
x
hhh
hx
RK
RK
nn
nn
实例分析
.1)0(,20,20 xtxdtdx
*分别用显式 Euler法、改进的 Euler法、经典 4级 4阶 RK法计算 *
考虑数值稳定性对步长 h的限制：
20,20 t
.13 9.00,078.244
.1.00,02
.1.00,02
hhhRK
hhhE ul e r
hhhE ul e r
即满足法的步长阶级即满足法的步长改进的即满足法的步长显式
。分别取为步长 15.0,1.0,05.0h
四、变步长和误差控制算法从前面的分析可知，求解常微分方程初始问题步长的选择是关键。然而，事先确定一个既保证收敛又数值稳定的步长很困难，另外方程的解可能在求解区间的某些部分变化平缓，可用大步长，而在其他部分变化剧烈，要用小步长。所以讨论变步长的算法是必要的。
对于一个选定的算法，我们可以事先分析确定保证数值稳定的步长 h的一个可选范围 [h_min,h_max]。接下来，我们讨论变步长的误差控制算法。
)(
,.(
111
11
呢使如何确定待定）处的近似值
。下面求时的近似值阶方法已经得到设用某种
nnn
nnnnnn
nn
xtxe r r
hxhhtt
xtp
)).(,
),(
)(
),(
00
nn txt
tx
xtx
txfx
过（
的解析解为设问题
).,
),(
~
)(
~
),
~
(
~
nn
nn
xt
tx
xtx
txfx
过（
的解析解为设问题
“一步误差”。
阶方法的为称
）
阶单步法方法采用
pd
hohxc
xtxd
hoT
hxthxx
p
n
pp
n
nnn
p
n
nnnn
1
21
111
1
1
1
)()(
(
~
)(
),,(
为正的常数）。则可以证明：只要 cce r rhd n
n
n (,
1
1
11
1
*
21
11
2
11
21
11
111
1
)()(
)()(
~
)()()(
~
.,
1,
的估 计1
nn
n
pp
nnn
p
nn
pp
nnn
nnnn
n
n
xxd
hohtcxx
hoxtx
hohtcxtx
xxhtt
pphh
、d
则处的近似值阶单步法分别计算阶和用设当前步长重新执行前面的计算。是合适的，否则说明显然，若
，
为步长，此时不一定合适，可以考虑当前的
、一步误差控制算法
,,1
))((
2
/1
*
1
*
11
1*
1
1
1
hhh
d
h
h
d
h
d
dhtcd
hhhh
n
p
n
pnn
p
n
p
nn
nn
1,
,,
1
),(
111
11
1
1
nnn
nn
nnn
nnnn
nn
nn
xxd
xpxp
htt
e n d ;tbhhtbif
hh
b,xt
计算阶法算：阶法算：
）（
计算终点算法：已知
）返回（
停止。
）（
）（
1
3;
)(/
1;
)(2
1
111
1
e l s e
btif
e nd
hh
dabsh
e l s e;h;hn nxx
ε/hdi f abs
n
nn
p
n
nnnn
n
function [t,x]=weifenfchqiujie(f_name,t0,x0,h_min,h_max,epsl,t1)
%解 dx/dt=f(x,t),x(t0)=x0,终点 t1.用变步长的算法。
t(1)=t0;
x(1)=x0;%初始值
h=(h_max+h_min)/2;%h选在绝对稳定区间内
n=1;
while 1
if t(n)-t1==0,break,end;
if t1-t(n)<h
h=t1-t(n);
end
k1=feval(f_name,x(n),t(n));
k2=feval(f_name,x(n)+h/4*k1,t(n)+h/4);
k3=feval(f_name,x(n)+3/32*h*k1+9/32*h*k2,t(n)+3/8*h);
k4=feval(f_name,x(n)+1932/2197*h*k1-
7200/2197*h*k2+7296/2197*h*k3,t(n)+12/13*h);
k5=feval(f_name,x(n)+439/216*h*k1-8*h*k2+3680/513*h*k3-
845/4104*h*k4,t(n)+h);
k6=feval(f_name,x(n)-8/27*h*k1+2*h*k2-
3544/2565*h*k3+1859/4104*h*k4-11/40*h*k5,t(n)+1/2*h);
x(n+1)=x(n)+h*(16/135*k1+6656/12825*k3+28561/56430*k4-
9/50*k5+2/55*k6); %5阶 RK-Fehlberg法
x1(n+1)=x(n)+h*(25/216*k1+1408/2565*k3+2197/4104*k4-
1/5*k5);%4阶 RK-Fehlberg法
%d=h*(1/360*k1-128/4275*k3-
2197/75240*k4+1/50*k5+2/55*k6);%一步误差
d=x1(n+1)-x(n+1);
if (abs(d)/h)>epsl
seita=0.84*sqrt(sqrt(epsl*h/abs(d)));%调整步长参数
h=seita*h;
else
t(n+1)=t(n)+h;
n=n+1;
end
end
实例 1
.1)0(,10,50 xtxdtdx
function f=ffff(x,t)
f=-50*x;
>> [t,x]=weifenfchqiujie('ffff',0,1,0.1,0.5,0.00001,1);
>> plot(t,x)
实例 2
.0)0(,20,1 21 2 xtttxdtdx
function f=ffff2(x,t)
f=1-2*t*x/(1+t^2);
>> [t,x]=weifenfchqiujie('ffff2',0,0,0.01,1,0.00001,2);
>> plot(t,x)
五、单步法的小结
1、单步法计算格式简单，是自开始的算法。
2、步长 h的选择原则上受到两个条件的限制,一是要保证算法收敛，二是要保证算法绝对稳定。
2、单步法的显格式绝对稳定区间相对比较窄，而隐格式绝对稳定区间相对比较宽，有的可是 A-稳定的。
3、利用隐格式绝对稳定区间相对比较宽的特点,可有效得到刚性问题的解。
§ 5.3 初始问题的线性多步法称为隐式。称为显式，如果步线性多步法。或写为则称为如步法。如果计算公式形这样的计算公式称为多用到前若干步的近似值如果计算
00
,,,,
00
10
11
1
0
1
0
1
11
1)-3-(5
n
k
j
jnj
k
j
jnjn
k
j
jnj
k
j
jnj
npnpnn
fhfhxx
k
fhx
xxxx
线性多步法的好处就是：少量的函数值计算可得到高阶收敛的方法！！！
一、一些常用的多步法
1,Adams方法
720
251
)937
5955(
24
44
8
3
)5
1623(
12
33
12
5
)3(
2
22
32
11
2
11
11
1
nn
nnnn
n
nnnn
nnnn
p
ff
ff
h
xx
f
ff
h
xx
ff
h
xx
cp
公式k
Adams显式方法绝对稳定区间
（ -1,0）
（ -6/11,0）
（ -3/10,0）
1 60
3
-
)191 06
2 646 46
2 51(
7 20
54
7 20
19
)5
199(
2443
24
1
)8
5(
1232
12
1
)(
2
21
32
1
11
21
11
1
11
11
1
nn
nn
nnn
nn
nnnn
nn
nnn
nnnn
p
ff
ff
f
h
xx
ff
ff
h
xx
ff
f
h
xx
ff
h
xx
cp 公式k
Adams隐式方法绝对稳定区间
（ -6,0）
（ -3,0）
（ -90/49,0）
)0,(
2,Gear方法（隐式）
1
32
11
12
11
111
1
25
12
25
3
25
36
25
48
44
11
6
11
2
11
9
11
18
33
3
2
3
1
3
4
22
n
nn
nnn
nn
nnn
nnnn
p
hf
xx
xxx
hfx
xxx
hfxxx
cp
25
16
公式k
绝对稳定区域
)90( 0A
)'2788( 0A
)'1473( 0A
二、实际计算要考虑的问题
1、多步法不能自开始，需要事先设置初始表头。例如 Adams
3步 3阶显式方法：
)51623(12 211 nnnnn fffhxx
.,,210 xxx需要表头表头的计算推荐用高阶的显式单步法，比如经典 4级 4阶
RK法。
2、隐式多步法不仅需要表头，而且每步要解非线性方程（组）。
当我们讨论的问题不是刚性问题时，可以采用预估 -校正策略：
111 3
2
3
1
3
4
nnnn hfxxx
例如 2步 2阶 Gear方法方法阶显式级。例如用同阶显式方法作预估形成表头法计算阶级用经典
A d a m s
xxxRK
22)2(
.,,44)1( 101
e nd
txffffxx
e nd
txhfxxx
mkf or
ff
h
xx
nif or
txfftxff
ii
m
ii
i
k
iii
k
i
ii
);,(;;
),(
3
2
3
1
3
4
:1
)3(
2
:1
),();,(
11110
)(
11
1
)1(
11
)(
1
01
)0(
1
111000
校正预估
m为自定的校正次数
§ 5.4 微分方程组和 刚性问题的解法
)'1,3()0(
2
2
1 0 0 19 9 9
9 9 91 0 0 1
),(
,)',(
.50
,1)0(,21 0 0 19 9 9
,3)0(,29 9 91 0 0 1
21
221
2
121
1
x
xtxf
dt
dx
xxx
t
xxx
dt
dx
xxx
dt
dx
则令例如此时，前面所讨论的方法均可用于解方程组。
本节着重讲解刚性方程组及其解法！！！
什么是刚性方程组？
的特征值。矩阵是其中刚性比成立：称为刚性方程组，若对问题：
x
f
J ac ob it
ttts
nkt
Ttt
Rxxtx
Ttttxf
dt
dx
k
k
nk
k
nk
k
n
)(
.1))(R e (m i n/))(R e (m a x)()2(
:1,0))(R e ()1(
],[
.)(
),,(
11
0
00
0
一般来说,s(t)为 O(10)时，称为一般刚性问题，当 s(t）超过
)10( 6o 称为严重刚性问题。
才能保证数值稳定。取为在求解区间内，步长方法计算上述问题，则阶级采用
0,0 0 1 4
2 0 0 0/78.2
44
h
RK
。是比较严重的刚性问题所以刚性比，矩阵的特征值上例,10,22 0 0 0 321 sJ a c o b i
实际上，计算到 t=0.01时快变部分已经不起作用，接下来我们应该取大的步长，然而由于稳定性要求仍然要用小步长。这就是求解刚性问题所面临的矛盾！
解决矛盾的途径在于使用的方法最好是 A-稳定的，即 h不受数值稳定性限制。但遗憾的是：早已证明，显式的方法不可能是 A-稳定的，隐式的多步法也只有不超过 2阶的方法才可能是
A-稳定的，而高阶的隐式 单步法计算量又太大 。
1)(
1)(
22 0 0 0
2
22 0 0 0
1
tt
tt
eetx
eetx
其解析解为：
。都含快变部分 tetxtx 2 0 0 021 )(),(?
看来,A-稳定的要求对于实际问题是过于严格了。因此提出了 稳定的概念。)(?A
一个解微分方程的数值算法称为 稳定的，如果它的 绝对稳定区域包含一个无限的三角形区域 。
)(?A
]2,0[
综合上述，解刚性方程组的算法一定是隐式的。
解刚性方程组的算法分为两大类：高阶的隐式单步法和隐式多步法。
1、常用的 A-稳定的隐式 RK法
)
2
,
2
(
)1(
1
h
tk
h
xfk
hkxx
RK
nn
nn
法（中点法）一级二阶
)),(
2
(
),(
)(
2
)2(
212
1
211
htkk
h
xfk
txfk
kk
h
xx
nn
nn
nn
RK 法二级二阶
)
6
3
2
),
3
32
(
4
(
)
6
3
2
),
3
32
(
4
(
)(
2
)3(
212
211
211
htkk
h
xfk
htkk
h
xfk
kk
h
xx
nn
nn
nn
）（）（
）（）（
法二级四阶
1
1
1
1
RK
2、常用的线性多步法就是 2,3,4步的 Gear方法，它是解刚性方程组的有效方法。
例 1：用 3步 Gear方法解本节开始的问题。
[t,x]=gearfa('fff',0,[3,1]',5,50,0.0001);
h=0.1
t=0:0.003;
例 2：用经典 4级 4阶 RK法计算本节开始的问题。
前面已经讨论 h<0.0014才能保证数值稳定。下面比较
h=0.001,0.00138,0.0014的计算结果。
h=0.001
[t,x]=rk44('fff',0,[3,1]',5,0.001);
t=0:0.003;
[t,x]=rk44('fff',0,[3,1]',5,0.00138);
t=0:0.003;
[t,x]=rk44('fff',0,[3,1]',5,0.0014);
t=0:0.003;
解析解的曲线：
1)(
1)(
22 0 0 0
2
22 0 0 0
1
tt
tt
eetx
eetx
t=0:0.001:5; x1=exp(-2000*t)+exp(-2*t)+1;
x2=-exp(-2000*t)+exp(-2*t)+1; plot(t,x1,t,x2);
t=0:0.0001:0.003;
§ 5.1 问题的提出是其解析解。
，
、边值问题是其解析解。
、初始问题常微分方程定解问题
btaRtx
xbxxax
btatxf
dt
dx
ttRtx
xtx
tttxf
dt
dx
n
ba
n
,)(
2)-(5
)(,)(
),(
2
,)(
1)-(5
)(
),,(
1
0
00
0
。形式，称此为标准模式均可化为一般
),(
0),,,,( )(
txf
dt
dx
xxxtF n
)()1(
,
)()1(
12
2
1
2
2
1
21
2
2
2
tfxxx
dt
dx
x
dt
dx
dt
dx
xxx
tfx
dt
dx
x
dt
xd
，则引入例如
).(~)(
.,,,,,,
)(
1010
txtx
xxxttt
tx
nn
得近似值逼近的方法得到然后在利用数的近似值求在给定时刻
。数值解法就是得解析解对于非线性问题很难求
然后按精度要求截断。
，即或采用台劳展开的方法然后数值积分。
，
积分方程：）微分方程（初始问题算法的构造是
4)-(5 )()()(
3)-(5,,2,1,0)),(()()(
1-5 *
1
1
1
1
nnnnn
nnn
t
t
nn
txhtxtx
tth
ndtttxftxtx
n
n
则为多步算法。则称之为单步算法。否即算法是只用到前面一步信息定义：若计算
),,(
,
1
1
nnnnnn
nn
htxhxx
xx
.,,,
1)(
1
,0
1:1
)()(
1
)1()2(
1
)2(
1
)( ),(
1
)(
1)1(,0)0(
10),()1(
*
10
10
1
2
11
2
11
2
1
2
n
nn
iiiiiiii
iiiiiii
xxx
xx
h
x
ni
tfxxx
h
xxxx
h
xxx
h
txxx
h
tx
xx
ttfxxxx
由此解出例子：有限元法等。差分法的法很多，有差分方法、边值问题的算法构造方
本章只介绍初始问题的解法，原因是初始问题系统科学中研究的最重要的问题。
§ 5.2 初始问题的单步法
!!
2
)(
)),(()()(
~
)(
~
)(
),(,
,2,1,0
1
2
1
1
111
1
1
nn
nnnnn
nnn
nnnnnn
nnnn
h
x
T
ttxfhtxtx
txtxT
txfftth
n
fhxx
局部截断误差：
其中，
、显式欧拉法一、欧拉法
的李布希兹常数为其中，
总体误差：
fL
e r rLhTe r r
xtxLhxtx
xtxLhxtxxtx
txfttxfhxtx
xtx
xtxtxtx
xtxe r r
nnnn
nnnnn
nnnnnnn
nnnnnnn
nn
nnnn
nnn
)1(
)()1()(
~
)()()(
~
)),()),((()(
)(
~
)(
~
)(
~
)(
)(
11
11
11
11
1111
111
!! )(
)1(
2
1)1(
2
)1())1((
2
)1(
2
,
2
),,)(m a x (),:1,m a x (
1
)1(2
1
22
1
0
1
0
22
22
1
22
1
102
hoe r r
he
L
M
hL
hL
h
M
e r r
e r rhLhLh
M
e r rhLh
M
e r r
h
M
T
ttttxMnkhh
n
hLn
n
n
n
n
k
k
nn
n
nk
及则设以上说明显式欧拉法是一阶收敛的算法
2、隐式欧拉法
。二阶收敛的方法，即
）、改进的欧拉法（隐式一阶收敛的方法，即
)(
,2,1,0
)],(),([
3
)(
,2,1,0
),(
2
1
111
1
111
hoe r r
n
txftxfhxx
hoe r r
n
txfhxx
n
nnnnnnn
n
nnnnn
二、龙格 -库塔法（ RK)
我们的目的是构造更高阶的单步法！
阶方法。即此单步法为步法的整体误差则此单截断误差结论：若单步法的局部
phoe r r
hoT
p
n
p
n
).(
),(
1
1
1
( *) )()(
.:0),),((
)),(()),((
:0,10,,
)),(()()(
2
0
1
2
0
1
1
1
1
)o( hkctxtx
mjrrxfk
)o( hrrxfcdtttxf
miahatrtth
dtttxftxtx
m
m
j
jjnn
jjj
m
m
j
jjj
t
t
iininn
t
t
nn
n
n
n
n
令利用拉格朗日求积公式分析：由
( ** ),1
),)((
)(
)()()(
)()()(
),),((
1
1
1
1
2
1
11
111
00
mj
hatkbhtxfk
bkbhtx
hokaahtx
hkaarxrx
ttxfktr
jn
j
i
ijinj
ji
j
i
ijin
j
i
iiin
jjjjj
nnn
待定
，则取的显式单步法：
阶可构造出）可知：适当选取参数）和式（由式（ 1,,***?mbac jijj
m
j
jjnn kcxx
0
1
,1
),(
),(
1
1
0
mj
hatkbhxfk
txfk
jn
j
i
ijinj
nn
称此结构为 m+1级龙格 -库塔法。它最多是 m+1阶的。
一些常用的显式 RK方法
1、二级二阶 RK法改进的欧拉法（ p277),中点公式（ p278),休恩方法（ p278)等。
2、三级三阶 RK法休恩三阶方法（ P279),库塔三阶方法（ P279)等
3、四级四阶 RK法经典四级四阶 RK法（ P280),基尔四级四阶 RK法 (P280)
三、单步法的数值稳定性为复常数试验模型：稳定性时，都针对同一数值依赖性，在讨论方法的小很困难。为摆脱这种估计其大太依赖于由于时，将产生误差计算则用单步法
，即实际得到近似值产生误差若计算
,
),,(),,?(),,(
)],,?(),,([?
),,(
,?
,?)(
111
1
1
xx
txfhtxhtx
htxhtxhxx
xx
x
htxhxx
xx
xx
nnnn
nnnnnn
nnn
n
nnnn
nnn
nnnn
这样做的根据是,1）对实验模型数值不稳定的方法，不可用；
2）一般的初始问题在其解的存在区域内，可局部线性化转化为试验方程。
几个常用算法的数值稳定性
.)
2
0
))(
02
11,
)1(
)1(
1
1
1
限制性（还要满足算法数值稳定限制，而且不仅要满足算法收敛（步长
），（间：为实数时，绝对稳定区当绝对稳定区域式为用于试验模型的计算公显式定区域、显式欧拉法的绝对稳
h
hoh
h
h
h
xhxhxx
E u l e r
nn
nnnn
是绝对稳定的。
都稳定的算法，只要显然，凡是稳定的。称此算法是面的整个左半平面，则复平的绝对稳定区域包含了定义：一个算法如果它限制即可！（只需考虑满足算法收敛步长
），（间：为实数时，绝对稳定区当绝对稳定区域式为用于试验模型的计算公隐式定区域、隐式欧拉法的绝对稳
0h0,)R e (-A
))(
0
11,
1
1
1
1
2
1
111
A
h
hoh
h
h
h
x
h
xxhxx
E ul e r
nn
nnnnn
限制即可！（只需考虑满足算法收敛步长
），（间：为实数时，绝对稳定区当稳定的。是绝对稳定区域
）（
式为用于试验模型的计算公改进的稳定区域、改进的欧拉法的绝对
))(
0
,1
2
2
,
2
2
2
2
2
3
2
1
111
hoh
h
A
h
h
h
h
x
h
h
xxx
h
xx
E u l e r
nn
nnnnnn
的限制。还要考虑数值稳定性（
限制，而且（不仅考虑满足算法收敛步长
），（间：为实数时，绝对稳定区当
（绝对稳定区域
（
（
公式为法用于试验模型的计算经典四级四阶法的绝对稳定区域、经典四级四阶
)
78.2
0
))(
078.2
1)
!4
)(
!3
)(
!2
)(
1,
)
!4
)(
!3
)(
!2
)(
1
)
!4
)(
!3
)(
!2
)(
1
4
2
432
432
1
432
1
h
hoh
h
hhh
h
hhh
h
x
hhh
hx
RK
RK
nn
nn
实例分析
.1)0(,20,20 xtxdtdx
*分别用显式 Euler法、改进的 Euler法、经典 4级 4阶 RK法计算 *
考虑数值稳定性对步长 h的限制：
20,20 t
.13 9.00,078.244
.1.00,02
.1.00,02
hhhRK
hhhE ul e r
hhhE ul e r
即满足法的步长阶级即满足法的步长改进的即满足法的步长显式
。分别取为步长 15.0,1.0,05.0h
四、变步长和误差控制算法从前面的分析可知，求解常微分方程初始问题步长的选择是关键。然而，事先确定一个既保证收敛又数值稳定的步长很困难，另外方程的解可能在求解区间的某些部分变化平缓，可用大步长，而在其他部分变化剧烈，要用小步长。所以讨论变步长的算法是必要的。
对于一个选定的算法，我们可以事先分析确定保证数值稳定的步长 h的一个可选范围 [h_min,h_max]。接下来，我们讨论变步长的误差控制算法。
)(
,.(
111
11
呢使如何确定待定）处的近似值
。下面求时的近似值阶方法已经得到设用某种
nnn
nnnnnn
nn
xtxe r r
hxhhtt
xtp
)).(,
),(
)(
),(
00
nn txt
tx
xtx
txfx
过（
的解析解为设问题
).,
),(
~
)(
~
),
~
(
~
nn
nn
xt
tx
xtx
txfx
过（
的解析解为设问题
“一步误差”。
阶方法的为称
）
阶单步法方法采用
pd
hohxc
xtxd
hoT
hxthxx
p
n
pp
n
nnn
p
n
nnnn
1
21
111
1
1
1
)()(
(
~
)(
),,(
为正的常数）。则可以证明：只要 cce r rhd n
n
n (,
1
1
11
1
*
21
11
2
11
21
11
111
1
)()(
)()(
~
)()()(
~
.,
1,
的估 计1
nn
n
pp
nnn
p
nn
pp
nnn
nnnn
n
n
xxd
hohtcxx
hoxtx
hohtcxtx
xxhtt
pphh
、d
则处的近似值阶单步法分别计算阶和用设当前步长重新执行前面的计算。是合适的，否则说明显然，若
，
为步长，此时不一定合适，可以考虑当前的
、一步误差控制算法
,,1
))((
2
/1
*
1
*
11
1*
1
1
1
hhh
d
h
h
d
h
d
dhtcd
hhhh
n
p
n
pnn
p
n
p
nn
nn
1,
,,
1
),(
111
11
1
1
nnn
nn
nnn
nnnn
nn
nn
xxd
xpxp
htt
e n d ;tbhhtbif
hh
b,xt
计算阶法算：阶法算：
）（
计算终点算法：已知
）返回（
停止。
）（
）（
1
3;
)(/
1;
)(2
1
111
1
e l s e
btif
e nd
hh
dabsh
e l s e;h;hn nxx
ε/hdi f abs
n
nn
p
n
nnnn
n
function [t,x]=weifenfchqiujie(f_name,t0,x0,h_min,h_max,epsl,t1)
%解 dx/dt=f(x,t),x(t0)=x0,终点 t1.用变步长的算法。
t(1)=t0;
x(1)=x0;%初始值
h=(h_max+h_min)/2;%h选在绝对稳定区间内
n=1;
while 1
if t(n)-t1==0,break,end;
if t1-t(n)<h
h=t1-t(n);
end
k1=feval(f_name,x(n),t(n));
k2=feval(f_name,x(n)+h/4*k1,t(n)+h/4);
k3=feval(f_name,x(n)+3/32*h*k1+9/32*h*k2,t(n)+3/8*h);
k4=feval(f_name,x(n)+1932/2197*h*k1-
7200/2197*h*k2+7296/2197*h*k3,t(n)+12/13*h);
k5=feval(f_name,x(n)+439/216*h*k1-8*h*k2+3680/513*h*k3-
845/4104*h*k4,t(n)+h);
k6=feval(f_name,x(n)-8/27*h*k1+2*h*k2-
3544/2565*h*k3+1859/4104*h*k4-11/40*h*k5,t(n)+1/2*h);
x(n+1)=x(n)+h*(16/135*k1+6656/12825*k3+28561/56430*k4-
9/50*k5+2/55*k6); %5阶 RK-Fehlberg法
x1(n+1)=x(n)+h*(25/216*k1+1408/2565*k3+2197/4104*k4-
1/5*k5);%4阶 RK-Fehlberg法
%d=h*(1/360*k1-128/4275*k3-
2197/75240*k4+1/50*k5+2/55*k6);%一步误差
d=x1(n+1)-x(n+1);
if (abs(d)/h)>epsl
seita=0.84*sqrt(sqrt(epsl*h/abs(d)));%调整步长参数
h=seita*h;
else
t(n+1)=t(n)+h;
n=n+1;
end
end
实例 1
.1)0(,10,50 xtxdtdx
function f=ffff(x,t)
f=-50*x;
>> [t,x]=weifenfchqiujie('ffff',0,1,0.1,0.5,0.00001,1);
>> plot(t,x)
实例 2
.0)0(,20,1 21 2 xtttxdtdx
function f=ffff2(x,t)
f=1-2*t*x/(1+t^2);
>> [t,x]=weifenfchqiujie('ffff2',0,0,0.01,1,0.00001,2);
>> plot(t,x)
五、单步法的小结
1、单步法计算格式简单，是自开始的算法。
2、步长 h的选择原则上受到两个条件的限制,一是要保证算法收敛，二是要保证算法绝对稳定。
2、单步法的显格式绝对稳定区间相对比较窄，而隐格式绝对稳定区间相对比较宽，有的可是 A-稳定的。
3、利用隐格式绝对稳定区间相对比较宽的特点,可有效得到刚性问题的解。
§ 5.3 初始问题的线性多步法称为隐式。称为显式，如果步线性多步法。或写为则称为如步法。如果计算公式形这样的计算公式称为多用到前若干步的近似值如果计算
00
,,,,
00
10
11
1
0
1
0
1
11
1)-3-(5
n
k
j
jnj
k
j
jnjn
k
j
jnj
k
j
jnj
npnpnn
fhfhxx
k
fhx
xxxx
线性多步法的好处就是：少量的函数值计算可得到高阶收敛的方法！！！
一、一些常用的多步法
1,Adams方法
720
251
)937
5955(
24
44
8
3
)5
1623(
12
33
12
5
)3(
2
22
32
11
2
11
11
1
nn
nnnn
n
nnnn
nnnn
p
ff
ff
h
xx
f
ff
h
xx
ff
h
xx
cp
公式k
Adams显式方法绝对稳定区间
（ -1,0）
（ -6/11,0）
（ -3/10,0）
1 60
3
-
)191 06
2 646 46
2 51(
7 20
54
7 20
19
)5
199(
2443
24
1
)8
5(
1232
12
1
)(
2
21
32
1
11
21
11
1
11
11
1
nn
nn
nnn
nn
nnnn
nn
nnn
nnnn
p
ff
ff
f
h
xx
ff
ff
h
xx
ff
f
h
xx
ff
h
xx
cp 公式k
Adams隐式方法绝对稳定区间
（ -6,0）
（ -3,0）
（ -90/49,0）
)0,(
2,Gear方法（隐式）
1
32
11
12
11
111
1
25
12
25
3
25
36
25
48
44
11
6
11
2
11
9
11
18
33
3
2
3
1
3
4
22
n
nn
nnn
nn
nnn
nnnn
p
hf
xx
xxx
hfx
xxx
hfxxx
cp
25
16
公式k
绝对稳定区域
)90( 0A
)'2788( 0A
)'1473( 0A
二、实际计算要考虑的问题
1、多步法不能自开始，需要事先设置初始表头。例如 Adams
3步 3阶显式方法：
)51623(12 211 nnnnn fffhxx
.,,210 xxx需要表头表头的计算推荐用高阶的显式单步法，比如经典 4级 4阶
RK法。
2、隐式多步法不仅需要表头，而且每步要解非线性方程（组）。
当我们讨论的问题不是刚性问题时，可以采用预估 -校正策略：
111 3
2
3
1
3
4
nnnn hfxxx
例如 2步 2阶 Gear方法方法阶显式级。例如用同阶显式方法作预估形成表头法计算阶级用经典
A d a m s
xxxRK
22)2(
.,,44)1( 101
e nd
txffffxx
e nd
txhfxxx
mkf or
ff
h
xx
nif or
txfftxff
ii
m
ii
i
k
iii
k
i
ii
);,(;;
),(
3
2
3
1
3
4
:1
)3(
2
:1
),();,(
11110
)(
11
1
)1(
11
)(
1
01
)0(
1
111000
校正预估
m为自定的校正次数
§ 5.4 微分方程组和 刚性问题的解法
)'1,3()0(
2
2
1 0 0 19 9 9
9 9 91 0 0 1
),(
,)',(
.50
,1)0(,21 0 0 19 9 9
,3)0(,29 9 91 0 0 1
21
221
2
121
1
x
xtxf
dt
dx
xxx
t
xxx
dt
dx
xxx
dt
dx
则令例如此时，前面所讨论的方法均可用于解方程组。
本节着重讲解刚性方程组及其解法！！！
什么是刚性方程组？
的特征值。矩阵是其中刚性比成立：称为刚性方程组，若对问题：
x
f
J ac ob it
ttts
nkt
Ttt
Rxxtx
Ttttxf
dt
dx
k
k
nk
k
nk
k
n
)(
.1))(R e (m i n/))(R e (m a x)()2(
:1,0))(R e ()1(
],[
.)(
),,(
11
0
00
0
一般来说,s(t)为 O(10)时，称为一般刚性问题，当 s(t）超过
)10( 6o 称为严重刚性问题。
才能保证数值稳定。取为在求解区间内，步长方法计算上述问题，则阶级采用
0,0 0 1 4
2 0 0 0/78.2
44
h
RK
。是比较严重的刚性问题所以刚性比，矩阵的特征值上例,10,22 0 0 0 321 sJ a c o b i
实际上，计算到 t=0.01时快变部分已经不起作用，接下来我们应该取大的步长，然而由于稳定性要求仍然要用小步长。这就是求解刚性问题所面临的矛盾！
解决矛盾的途径在于使用的方法最好是 A-稳定的，即 h不受数值稳定性限制。但遗憾的是：早已证明，显式的方法不可能是 A-稳定的，隐式的多步法也只有不超过 2阶的方法才可能是
A-稳定的，而高阶的隐式 单步法计算量又太大 。
1)(
1)(
22 0 0 0
2
22 0 0 0
1
tt
tt
eetx
eetx
其解析解为：
。都含快变部分 tetxtx 2 0 0 021 )(),(?
看来,A-稳定的要求对于实际问题是过于严格了。因此提出了 稳定的概念。)(?A
一个解微分方程的数值算法称为 稳定的，如果它的 绝对稳定区域包含一个无限的三角形区域 。
)(?A
]2,0[
综合上述，解刚性方程组的算法一定是隐式的。
解刚性方程组的算法分为两大类：高阶的隐式单步法和隐式多步法。
1、常用的 A-稳定的隐式 RK法
)
2
,
2
(
)1(
1
h
tk
h
xfk
hkxx
RK
nn
nn
法（中点法）一级二阶
)),(
2
(
),(
)(
2
)2(
212
1
211
htkk
h
xfk
txfk
kk
h
xx
nn
nn
nn
RK 法二级二阶
)
6
3
2
),
3
32
(
4
(
)
6
3
2
),
3
32
(
4
(
)(
2
)3(
212
211
211
htkk
h
xfk
htkk
h
xfk
kk
h
xx
nn
nn
nn
）（）（
）（）（
法二级四阶
1
1
1
1
RK
2、常用的线性多步法就是 2,3,4步的 Gear方法，它是解刚性方程组的有效方法。
例 1：用 3步 Gear方法解本节开始的问题。
[t,x]=gearfa('fff',0,[3,1]',5,50,0.0001);
h=0.1
t=0:0.003;
例 2：用经典 4级 4阶 RK法计算本节开始的问题。
前面已经讨论 h<0.0014才能保证数值稳定。下面比较
h=0.001,0.00138,0.0014的计算结果。
h=0.001
[t,x]=rk44('fff',0,[3,1]',5,0.001);
t=0:0.003;
[t,x]=rk44('fff',0,[3,1]',5,0.00138);
t=0:0.003;
[t,x]=rk44('fff',0,[3,1]',5,0.0014);
t=0:0.003;
解析解的曲线：
1)(
1)(
22 0 0 0
2
22 0 0 0
1
tt
tt
eetx
eetx
t=0:0.001:5; x1=exp(-2000*t)+exp(-2*t)+1;
x2=-exp(-2000*t)+exp(-2*t)+1; plot(t,x1,t,x2);
t=0:0.0001:0.003;
