第九章 常微分方程数值解
/* Numerical Methods for Ordinary Differential Equations */
考虑 一阶 常微分方程的 初值问题 /* Initial-Value Problem */:
0)(
],[),(
yay
baxyxf
dx
dy
只要 f (x,y) 在 [a,b]? R1 上连续，且关于 y 满足 Lipschitz 条件,即存在与 x,y 无关的常数 L 使对任意定义在 [a,b] 上的 y1(x) 和 y2(x) 都成立，则上述 IVP存在唯一解 。
|||),(),(| 2121 yyLyxfyxf
要计算出解函数 y(x) 在一系列节点 a = x0< x1<…< xn= b
处的近似值 ),...,1()( nixyy ii
节点间距 为步长，通常采用 等距节点,
即取 hi = h (常数 )。
)1,.,,,0(1 nixxh iii
§ 1 欧拉方法 /* Euler’s Method */
欧拉公式：
x0 x1
向前差商近似导数
h
xyxyxy )()()( 01
0
),()()()( 000001 yxfhyxyhxyxy 1y
记为
)1,...,0(),(1 niyxfhyy iiii
定义 在假设 yi = y(xi),即第 i 步计算是精确的前提下,考虑的截断误差 Ri = y(xi+1)? yi+1 称为 局部截断误差 /* local
truncation error */。
定义 若某算法的局部截断误差为 O(hp+1),则称该算法有 p
阶精度 。
欧拉法的局部截断误差：
)],([)]()()()([)( 3211 2 iiiihiiiii yxhfyhOxyxyhxyyxyR
)()( 322 hOxy ih 欧拉法具有 1 阶精度。
Ri 的 主项
/* leading term */
亦称为 欧拉折线法
/* Euler’s polygonal arc method*/
§ 1 Euler’s Method
欧拉公式的改进：
隐式欧拉法 /* implicit Euler method */
向后差商近似导数
h
xyxyxy )()()( 01
1
x0 x1
))(,()( 1101 xyxfhyxy
)1,...,0(),( 111 niyxfhyy iiii
由于未知数 yi+1同时出现在等式的两边，不能直接得到，故称为 隐式 /* implicit */ 欧拉公式，而前者称为 显式 /* explicit */
欧拉公式。
一般先用显式计算一个初值，再 迭代 求解。
隐式 欧拉法的局部截断误差：
11 )( iii yxyR )()( 322 hOxy ih
即隐式欧拉公式具有 1 阶精度。
Hey! Isn’t the leading term of the local
truncation errorof Euler’s method?
Seems that we can make a good
use of it …
)(22 ih xy
§ 1 Euler’s Method
梯形公式 /* trapezoid formula */ — 显、隐式两种算法的 平均
)1,.,,,0()],(),([2 111 niyxfyxfhyy iiiiii
注,的确有局部截断误差,
即梯形公式具有 2 阶精度，比欧拉方法有了进步。
但注意到该公式是 隐式 公式，计算时不得不用到迭代法，其迭代收敛性与欧拉公式相似。
)()( 311 hOyxyR iii
中点欧拉公式 /* midpoint formula */
中心差商近似导数
h
xyxyxy
2
)()()( 02
1
x0 x2x1
))(,(2)()( 1102 xyxfhxyxy
1,.,,,1),(211 niyxfhyy iiii
假设,则可以导出即中点公式具有 2 阶精度。
)(),( 11 iiii xyyxyy )()( 311 hOyxyR iii
需要 2个初值 y0和 y1来启动递推过程，这样的算法称为 双步法 /* double-step
method */，而前面的三种算法都是 单步法
/* single-step method */。
方 法
§ 1 Euler’s Method
显式欧拉隐式欧拉梯形公式中点公式简单 精度低稳定性最好 精度低,计算量大精度提高 计算量大精度提高,显式 多一个初值,可能影响精度
Can’t you give me a formula
with all the advantages yet without any
of the disadvantages?
Do you think
it possible?Well,call me greedy…OK,let’s make it
possible.
改进欧拉法 /* modified Euler’s method */
Step 1,先用 显式 欧拉公式作 预测,算出 ),(1 iiii yxfhyy
Step 2,再将 代入 隐式 梯形公式的右边作 校正,得到1?iy
)],(),([2 111 iiiiii yxfyxfhyy
注,此法亦称为 预测 -校正法 /* predictor-corrector method */。
可以证明该算法具有 2 阶精度，同时可以看到它是个 单步 递推格式，比隐式公式的迭代求解过程 简单 。后面将看到，它的 稳定性高 于显式欧拉法。
)1,...,0(),(,),(2 11 niyxfhyxfyxfhyy iiiiiiii
§ 1 Euler’s Method
§ 2 龙格 - 库塔法 /* Runge-Kutta Method */
建立高精度的单步递推格式。
单步递推法的 基本思想 是从 ( xi,yi ) 点出发，以 某一斜率 沿直线达到 ( xi+1,yi+1 ) 点。欧拉法及其各种变形所能达到的最高精度为 2阶 。
考察改进的欧拉法，可以将其改写为：
),(
),(
2
1
2
1
12
1
211
hKyhxfK
yxfK
KKhyy
ii
ii
ii
斜率一定取 K1 K2
的 平均值 吗？
步长一定是一个 h 吗？
§ 2 Runge-Kutta Method
首先希望能确定系数?1,?2,p，使得到的算法格式有 2阶精度，即在 的前提假设下，使得)(
ii xyy?
)()( 311 hOyxyR iii
Step 1,将 K2 在 ( xi,yi ) 点作 Taylor 展开
)(),(),(),(
),(
2
1
12
hOyxfp hKyxp hfyxf
p hKyphxfK
iiyiixii
ii
)()()( 2hOxyphxy ii
将改进欧拉法推广为：
),(
),(
][
12
1
22111
phKyphxfK
yxfK
KKhyy
ii
ii
ii
),(),(),(
),(),(
),()(
yxfyxfyxf
dx
dy
yxfyxf
yxf
dx
d
xy
yx
yx
Step 2,将 K2 代入第 1式，得到
)()()()(
)]()()([)(
32
221
2
211
hOxyphxyhy
hOxyphxyxyhyy
iii
iiiii
§ 2 Runge-Kutta Method
Step 3,将 yi+1 与 y( xi+1 ) 在 xi 点的 泰勒 展开作比较
)()()()( 322211 hOxyphxyhyy iiii
)()(2)()()( 3
2
1 hOxy
hxyhxyxy
iiii
要求,则必须有,)()( 3
11 hOyxyR iii
2
1,1
221 p
这里有 个未知数,个方程。
3
2
存在 无穷多个解 。所有满足上式的格式统称为 2阶龙格 - 库塔格式 。
2
1,1
21p
注意到,就是改进的欧拉法。
Q,为获得更高的精度，应该如何进一步推广？
其中?i ( i = 1,…,m ),?i ( i
= 2,…,m ) 和?ij ( i = 2,…,
m; j = 1,…,i?1 ) 均为待定系数，确定这些系数的步骤与前面相似。
§ 2 Runge-Kutta Method
)...,(
......
),(
),(
),(
]...[
112211
23213133
12122
1
22111
mm mmmmim
ii
ii
ii
mmii
hKhKhKyhxfK
hKhKyhxfK
hKyhxfK
yxfK
KKKhyy
最常用为四级 4阶 经典龙格 -库塔法 /* Classical Runge-Kutta
Method */,
),(
),(
),(
),(
)22(
34
2223
1222
1
432161
hKyhxfK
KyxfK
KyxfK
yxfK
KKKKyy
ii
h
i
h
i
h
i
h
i
ii
h
ii
§ 2 Runge-Kutta Method
注：
龙格 -库塔法 的主要运算在于计算 Ki 的值，即计算 f 的值。 Butcher 于 1965年给出了计算量与可达到的最高精度阶数的关系：
753
可达到的最高精度
642每步须算 Ki 的个数
)( 2hO )( 3hO )( 4hO )( 5hO )( 6hO)( 4hO )( 2?nhO
8?n
由于龙格 -库塔法的导出基于泰勒展开，故精度主要受解函数的光滑性影响。对于光滑性不太好的解，最好采用 低阶算法 而将步长 h取小 。
HW,
p.202 #1,2
§ 2 Runge-Kutta Method
Lab 15,Runge-Kutta (Order Four)
Use Runge-Kutta method of order four to approximate the solution of the
initial-value problem
,,and (1)
You are supposed to write a function
void RK4 ( double (*f)( ),double a,double b,double y0,int n,FILE *outfile)
to approximate the solution of Problem (1) with y’ = f,x in [a,b],and the
initial value of y being y0,Output the approximating values of y on the
n+1 equally spaced grid points from a to b to outfile.
Input
There is no input file,Instead,you must hand in your function in a *.h file.
The rule of naming the *.h file is the same as that of naming the *.c or
*.cpp files.
Output (?represents a space)
For each test case,you are supposed to print n+1 lines,and each line is in
the following format,fprintf( outfile,,%8.4f?%12.8e\n”,x,y );
),( yxfy ],[ bax? 0)( yay?
§ 2 Runge-Kutta Method
Sample Judge Program
#include <stdio.h>
#include <math.h>
#include "98115001_15.h"
double f0(double x,double y)
{ return (y?x*x+1.0); }
void main( )
{
FILE *outfile = fopen("out.txt","w");
int n;
double a,b,y0;
a = 0.0; b = 2.0; y0 = 0.50; n = 10;
RK4(f0,a,b,y0,n,outfile);
fprintf(outfile,"\n");
fclose(outfile);
}
Sample Output
(?represents a space)
0.0000?5.00000000e?001
0.2000?8.29293333e?001
0.4000?1.21407621e+000
0.6000?1.64892202e+000
0.8000?2.12720268e+000
1.0000?2.64082269e+000
1.2000?3.17989417e+000
1.4000?3.73234007e+000
1.6000?4.28340950e+000
1.8000?4.81508569e+000
2.0000?5.30536300e+000
§ 3 收敛性与稳定性 /* Convergency and Stability */
收敛性 /* Convergency */
定义 若某算法对于任意固定的 x = xi = x0 + i h,当 h?0
( 同时 i) 时有 yi? y( xi ),则称该算法是 收敛 的 。
例,就初值问题 考察欧拉显式格式的收敛性。
0)0( yy
yy?
解,该问题的精确解为 xeyxy?
0)(?
欧拉公式为
iiii yhyhyy )1(1 0)1( yhy ii
对任意固定的 x = xi = i h,有
ii xhhxi hyhyy ])1[()1( /10/0
eh hh /10 )1(l im
)(0 ix xyey i
§ 3 Convergency and Stability
稳定性 /* Stability */
例,考察初值问题 在区间 [0,0.5]上的解。
分别用欧拉显、隐式格式和改进的欧拉格式计算数值解。
1)0(
)(30)(
y
xyxy
0.0
0.1
0.2
0.3
0.4
0.5
精确解改进欧拉法欧拉隐式欧拉显式节点 xi xey 30
1.0000
2.0000
4.0000
8.0000
1.6000?101
3.2000?101
1.0000
2.5000?10?1
6.2500?10?2
1.5625?10?2
3.9063?10?3
9.7656?10?4
1.0000
2.5000
6.2500
1.5626?101
3.9063?101
9.7656?101
1.0000
4.9787?10?2
2.4788?10?3
1.2341?10?4
6.1442?10?6
3.0590?10?7
What is wrong!
An Engineer complains,"Math
theorems are so unstable that
a small perturbation on the conditions
will cause a crash on the conclusions!"
§ 3 Convergency and Stability
定义 若某算法在计算过程中任一步产生的误差在以后的计算中都 逐步衰减,则称该算法是 绝对稳定的 /*absolutely
stable */。
一般分析时为简单起见，只考虑 试验方程 /* test equation */
yy
常数，可以是复数当步长取为 h 时，将某算法应用于上式，并假设只在初值产生误差,则若此误差以后逐步衰减，就称该算法相对于 绝对稳定,的全体构成 绝对稳定区域 。
我们称 算法 A 比算法 B 稳定,就是指 A 的绝对稳定区域比 B
的 大 。
000 yy
h? h? h
§ 3 Convergency and Stability
例,考察显式欧拉法
011 )1( yhyhyy iiii
000 yy 011 )1( yhy ii
01111 )1( iiii hyy
由此可见，要保证初始误差?0 以后逐步衰减，
必须满足：
hh
1|1| h
0-1-2 Re
Img
例,考察隐式欧拉法 11 iii yhyy?
ii yhy
1
1
1 0
1
1 1
1
i
i h
可见绝对稳定区域为,1|1| h
210 Re
Img
注,一般来说，隐式欧拉法的绝对稳定性比同阶的显式法的好。
§ 3 Convergency and Stability
例,隐式龙格 -库塔法
),...,1(
)...,(
]...[
11
111
mj
hKhKyhxfK
KKhyy
mmjjijij
mmii
而 显式 1~ 4 阶方法的绝对稳定区域为
)2,2( 11
11
KhyhxfK
hKyy
ii
ii其中 2阶方法 的绝对稳定区域为
0 Re
Img
k=1
k=2
k=3
k=4
-1-2-3
-
-
-
1
2
3
Re
Img
无条件稳定
HW,
p.202 #6
/* Numerical Methods for Ordinary Differential Equations */
考虑 一阶 常微分方程的 初值问题 /* Initial-Value Problem */:
0)(
],[),(
yay
baxyxf
dx
dy
只要 f (x,y) 在 [a,b]? R1 上连续，且关于 y 满足 Lipschitz 条件,即存在与 x,y 无关的常数 L 使对任意定义在 [a,b] 上的 y1(x) 和 y2(x) 都成立，则上述 IVP存在唯一解 。
|||),(),(| 2121 yyLyxfyxf
要计算出解函数 y(x) 在一系列节点 a = x0< x1<…< xn= b
处的近似值 ),...,1()( nixyy ii
节点间距 为步长，通常采用 等距节点,
即取 hi = h (常数 )。
)1,.,,,0(1 nixxh iii
§ 1 欧拉方法 /* Euler’s Method */
欧拉公式：
x0 x1
向前差商近似导数
h
xyxyxy )()()( 01
0
),()()()( 000001 yxfhyxyhxyxy 1y
记为
)1,...,0(),(1 niyxfhyy iiii
定义 在假设 yi = y(xi),即第 i 步计算是精确的前提下,考虑的截断误差 Ri = y(xi+1)? yi+1 称为 局部截断误差 /* local
truncation error */。
定义 若某算法的局部截断误差为 O(hp+1),则称该算法有 p
阶精度 。
欧拉法的局部截断误差：
)],([)]()()()([)( 3211 2 iiiihiiiii yxhfyhOxyxyhxyyxyR
)()( 322 hOxy ih 欧拉法具有 1 阶精度。
Ri 的 主项
/* leading term */
亦称为 欧拉折线法
/* Euler’s polygonal arc method*/
§ 1 Euler’s Method
欧拉公式的改进：
隐式欧拉法 /* implicit Euler method */
向后差商近似导数
h
xyxyxy )()()( 01
1
x0 x1
))(,()( 1101 xyxfhyxy
)1,...,0(),( 111 niyxfhyy iiii
由于未知数 yi+1同时出现在等式的两边，不能直接得到，故称为 隐式 /* implicit */ 欧拉公式，而前者称为 显式 /* explicit */
欧拉公式。
一般先用显式计算一个初值，再 迭代 求解。
隐式 欧拉法的局部截断误差：
11 )( iii yxyR )()( 322 hOxy ih
即隐式欧拉公式具有 1 阶精度。
Hey! Isn’t the leading term of the local
truncation errorof Euler’s method?
Seems that we can make a good
use of it …
)(22 ih xy
§ 1 Euler’s Method
梯形公式 /* trapezoid formula */ — 显、隐式两种算法的 平均
)1,.,,,0()],(),([2 111 niyxfyxfhyy iiiiii
注,的确有局部截断误差,
即梯形公式具有 2 阶精度，比欧拉方法有了进步。
但注意到该公式是 隐式 公式，计算时不得不用到迭代法，其迭代收敛性与欧拉公式相似。
)()( 311 hOyxyR iii
中点欧拉公式 /* midpoint formula */
中心差商近似导数
h
xyxyxy
2
)()()( 02
1
x0 x2x1
))(,(2)()( 1102 xyxfhxyxy
1,.,,,1),(211 niyxfhyy iiii
假设,则可以导出即中点公式具有 2 阶精度。
)(),( 11 iiii xyyxyy )()( 311 hOyxyR iii
需要 2个初值 y0和 y1来启动递推过程，这样的算法称为 双步法 /* double-step
method */，而前面的三种算法都是 单步法
/* single-step method */。
方 法
§ 1 Euler’s Method
显式欧拉隐式欧拉梯形公式中点公式简单 精度低稳定性最好 精度低,计算量大精度提高 计算量大精度提高,显式 多一个初值,可能影响精度
Can’t you give me a formula
with all the advantages yet without any
of the disadvantages?
Do you think
it possible?Well,call me greedy…OK,let’s make it
possible.
改进欧拉法 /* modified Euler’s method */
Step 1,先用 显式 欧拉公式作 预测,算出 ),(1 iiii yxfhyy
Step 2,再将 代入 隐式 梯形公式的右边作 校正,得到1?iy
)],(),([2 111 iiiiii yxfyxfhyy
注,此法亦称为 预测 -校正法 /* predictor-corrector method */。
可以证明该算法具有 2 阶精度，同时可以看到它是个 单步 递推格式，比隐式公式的迭代求解过程 简单 。后面将看到，它的 稳定性高 于显式欧拉法。
)1,...,0(),(,),(2 11 niyxfhyxfyxfhyy iiiiiiii
§ 1 Euler’s Method
§ 2 龙格 - 库塔法 /* Runge-Kutta Method */
建立高精度的单步递推格式。
单步递推法的 基本思想 是从 ( xi,yi ) 点出发，以 某一斜率 沿直线达到 ( xi+1,yi+1 ) 点。欧拉法及其各种变形所能达到的最高精度为 2阶 。
考察改进的欧拉法，可以将其改写为：
),(
),(
2
1
2
1
12
1
211
hKyhxfK
yxfK
KKhyy
ii
ii
ii
斜率一定取 K1 K2
的 平均值 吗？
步长一定是一个 h 吗？
§ 2 Runge-Kutta Method
首先希望能确定系数?1,?2,p，使得到的算法格式有 2阶精度，即在 的前提假设下，使得)(
ii xyy?
)()( 311 hOyxyR iii
Step 1,将 K2 在 ( xi,yi ) 点作 Taylor 展开
)(),(),(),(
),(
2
1
12
hOyxfp hKyxp hfyxf
p hKyphxfK
iiyiixii
ii
)()()( 2hOxyphxy ii
将改进欧拉法推广为：
),(
),(
][
12
1
22111
phKyphxfK
yxfK
KKhyy
ii
ii
ii
),(),(),(
),(),(
),()(
yxfyxfyxf
dx
dy
yxfyxf
yxf
dx
d
xy
yx
yx
Step 2,将 K2 代入第 1式，得到
)()()()(
)]()()([)(
32
221
2
211
hOxyphxyhy
hOxyphxyxyhyy
iii
iiiii
§ 2 Runge-Kutta Method
Step 3,将 yi+1 与 y( xi+1 ) 在 xi 点的 泰勒 展开作比较
)()()()( 322211 hOxyphxyhyy iiii
)()(2)()()( 3
2
1 hOxy
hxyhxyxy
iiii
要求,则必须有,)()( 3
11 hOyxyR iii
2
1,1
221 p
这里有 个未知数,个方程。
3
2
存在 无穷多个解 。所有满足上式的格式统称为 2阶龙格 - 库塔格式 。
2
1,1
21p
注意到,就是改进的欧拉法。
Q,为获得更高的精度，应该如何进一步推广？
其中?i ( i = 1,…,m ),?i ( i
= 2,…,m ) 和?ij ( i = 2,…,
m; j = 1,…,i?1 ) 均为待定系数，确定这些系数的步骤与前面相似。
§ 2 Runge-Kutta Method
)...,(
......
),(
),(
),(
]...[
112211
23213133
12122
1
22111
mm mmmmim
ii
ii
ii
mmii
hKhKhKyhxfK
hKhKyhxfK
hKyhxfK
yxfK
KKKhyy
最常用为四级 4阶 经典龙格 -库塔法 /* Classical Runge-Kutta
Method */,
),(
),(
),(
),(
)22(
34
2223
1222
1
432161
hKyhxfK
KyxfK
KyxfK
yxfK
KKKKyy
ii
h
i
h
i
h
i
h
i
ii
h
ii
§ 2 Runge-Kutta Method
注：
龙格 -库塔法 的主要运算在于计算 Ki 的值，即计算 f 的值。 Butcher 于 1965年给出了计算量与可达到的最高精度阶数的关系：
753
可达到的最高精度
642每步须算 Ki 的个数
)( 2hO )( 3hO )( 4hO )( 5hO )( 6hO)( 4hO )( 2?nhO
8?n
由于龙格 -库塔法的导出基于泰勒展开，故精度主要受解函数的光滑性影响。对于光滑性不太好的解，最好采用 低阶算法 而将步长 h取小 。
HW,
p.202 #1,2
§ 2 Runge-Kutta Method
Lab 15,Runge-Kutta (Order Four)
Use Runge-Kutta method of order four to approximate the solution of the
initial-value problem
,,and (1)
You are supposed to write a function
void RK4 ( double (*f)( ),double a,double b,double y0,int n,FILE *outfile)
to approximate the solution of Problem (1) with y’ = f,x in [a,b],and the
initial value of y being y0,Output the approximating values of y on the
n+1 equally spaced grid points from a to b to outfile.
Input
There is no input file,Instead,you must hand in your function in a *.h file.
The rule of naming the *.h file is the same as that of naming the *.c or
*.cpp files.
Output (?represents a space)
For each test case,you are supposed to print n+1 lines,and each line is in
the following format,fprintf( outfile,,%8.4f?%12.8e\n”,x,y );
),( yxfy ],[ bax? 0)( yay?
§ 2 Runge-Kutta Method
Sample Judge Program
#include <stdio.h>
#include <math.h>
#include "98115001_15.h"
double f0(double x,double y)
{ return (y?x*x+1.0); }
void main( )
{
FILE *outfile = fopen("out.txt","w");
int n;
double a,b,y0;
a = 0.0; b = 2.0; y0 = 0.50; n = 10;
RK4(f0,a,b,y0,n,outfile);
fprintf(outfile,"\n");
fclose(outfile);
}
Sample Output
(?represents a space)
0.0000?5.00000000e?001
0.2000?8.29293333e?001
0.4000?1.21407621e+000
0.6000?1.64892202e+000
0.8000?2.12720268e+000
1.0000?2.64082269e+000
1.2000?3.17989417e+000
1.4000?3.73234007e+000
1.6000?4.28340950e+000
1.8000?4.81508569e+000
2.0000?5.30536300e+000
§ 3 收敛性与稳定性 /* Convergency and Stability */
收敛性 /* Convergency */
定义 若某算法对于任意固定的 x = xi = x0 + i h,当 h?0
( 同时 i) 时有 yi? y( xi ),则称该算法是 收敛 的 。
例,就初值问题 考察欧拉显式格式的收敛性。
0)0( yy
yy?
解,该问题的精确解为 xeyxy?
0)(?
欧拉公式为
iiii yhyhyy )1(1 0)1( yhy ii
对任意固定的 x = xi = i h,有
ii xhhxi hyhyy ])1[()1( /10/0
eh hh /10 )1(l im
)(0 ix xyey i
§ 3 Convergency and Stability
稳定性 /* Stability */
例,考察初值问题 在区间 [0,0.5]上的解。
分别用欧拉显、隐式格式和改进的欧拉格式计算数值解。
1)0(
)(30)(
y
xyxy
0.0
0.1
0.2
0.3
0.4
0.5
精确解改进欧拉法欧拉隐式欧拉显式节点 xi xey 30
1.0000
2.0000
4.0000
8.0000
1.6000?101
3.2000?101
1.0000
2.5000?10?1
6.2500?10?2
1.5625?10?2
3.9063?10?3
9.7656?10?4
1.0000
2.5000
6.2500
1.5626?101
3.9063?101
9.7656?101
1.0000
4.9787?10?2
2.4788?10?3
1.2341?10?4
6.1442?10?6
3.0590?10?7
What is wrong!
An Engineer complains,"Math
theorems are so unstable that
a small perturbation on the conditions
will cause a crash on the conclusions!"
§ 3 Convergency and Stability
定义 若某算法在计算过程中任一步产生的误差在以后的计算中都 逐步衰减,则称该算法是 绝对稳定的 /*absolutely
stable */。
一般分析时为简单起见，只考虑 试验方程 /* test equation */
yy
常数，可以是复数当步长取为 h 时，将某算法应用于上式，并假设只在初值产生误差,则若此误差以后逐步衰减，就称该算法相对于 绝对稳定,的全体构成 绝对稳定区域 。
我们称 算法 A 比算法 B 稳定,就是指 A 的绝对稳定区域比 B
的 大 。
000 yy
h? h? h
§ 3 Convergency and Stability
例,考察显式欧拉法
011 )1( yhyhyy iiii
000 yy 011 )1( yhy ii
01111 )1( iiii hyy
由此可见，要保证初始误差?0 以后逐步衰减，
必须满足：
hh
1|1| h
0-1-2 Re
Img
例,考察隐式欧拉法 11 iii yhyy?
ii yhy
1
1
1 0
1
1 1
1
i
i h
可见绝对稳定区域为,1|1| h
210 Re
Img
注,一般来说，隐式欧拉法的绝对稳定性比同阶的显式法的好。
§ 3 Convergency and Stability
例,隐式龙格 -库塔法
),...,1(
)...,(
]...[
11
111
mj
hKhKyhxfK
KKhyy
mmjjijij
mmii
而 显式 1~ 4 阶方法的绝对稳定区域为
)2,2( 11
11
KhyhxfK
hKyy
ii
ii其中 2阶方法 的绝对稳定区域为
0 Re
Img
k=1
k=2
k=3
k=4
-1-2-3
-
-
-
1
2
3
Re
Img
无条件稳定
HW,
p.202 #6
