第九章 常微分方程数值解
/* 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
)) (,( ) ( 1 1 0 1 x y x f h y x y ? ?
) 1,...,0 ( ),( 1 1 1 ? ? ? ? ? ? ? n i y x f h y y i i i i
由于未知数 yi+1 同时出现在等式的两边，不能直接得到，故
称为 隐式 /* implicit */ 欧拉公式，而前者称为 显式 /* explicit */
欧拉公式。
一般先用显式计算一个初值，再 迭代 求解。
? 隐式 欧拉法的局部截断误差,
11 )( ?? ?? iii yxyR )()( 322 hOxy ih ?????
即隐式欧拉公式具有 1 阶精度。
Hey! Isn’t the leading term of the local
truncation error of 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 x2 x1
))(,(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 i i i i y x f h y y ? ? ?
Step 2,再将 代入 隐式 梯形公式的右边作 校正,得到 1 ? i y
)],( ),( [ 2 1 1 1 ? ? ? ? ? ? i i i i i i y x f y x f h y y
注,此法亦称为 预测 -校正法 /* 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 ??????
将改进欧拉法推广为,
),(
),(
] [
1 2
1
2 2 1 1 1
phK y ph x f K
y x f K
K K h y y
i i
i i
i i
? ? ?
?
? ? ? ? ? ?
),(),(),(
),(),(
),()(
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
21 ??? ??p
注意到,就是改进的欧拉法。
Q,为获得更高的精度，应该如何进一步推广？
其中 ?i ( i = 1,…,m ),?i ( i
= 2,…,m ) 和 ?ij ( i = 2,…,
m; j = 1,…,i?1 ) 均为待定
系数，确定这些系数的
步骤与前面相似。
§ 2 Runge-Kutta Method
),..,(
...,.,
),(
),(
),(
],.,[
1 1 2 2 1 1
2 32 1 31 3 3
1 21 2 2
1
2 2 1 1 1
? ?
?
? ? ? ? ? ?
? ? ? ?
? ? ?
?
? ? ? ? ?
m m m m m m i m
i i
i i
i i
m m i i
hK hK hK y h x f K
hK hK y h x f K
hK y h x f K
y x f K
K K K h y y
? ? ? ?
? ? ?
? ?
? ? ?
? 最常用为四级 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年给出了计算量与可达到的最高精
度阶数的关系,
7 5 3
可达到的最高精度
6 4 2 每步须算 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
2 1 0 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
)) (,( ) ( 1 1 0 1 x y x f h y x y ? ?
) 1,...,0 ( ),( 1 1 1 ? ? ? ? ? ? ? n i y x f h y y i i i i
由于未知数 yi+1 同时出现在等式的两边，不能直接得到，故
称为 隐式 /* implicit */ 欧拉公式，而前者称为 显式 /* explicit */
欧拉公式。
一般先用显式计算一个初值，再 迭代 求解。
? 隐式 欧拉法的局部截断误差,
11 )( ?? ?? iii yxyR )()( 322 hOxy ih ?????
即隐式欧拉公式具有 1 阶精度。
Hey! Isn’t the leading term of the local
truncation error of 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 x2 x1
))(,(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 i i i i y x f h y y ? ? ?
Step 2,再将 代入 隐式 梯形公式的右边作 校正,得到 1 ? i y
)],( ),( [ 2 1 1 1 ? ? ? ? ? ? i i i i i i y x f y x f h y y
注,此法亦称为 预测 -校正法 /* 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 ??????
将改进欧拉法推广为,
),(
),(
] [
1 2
1
2 2 1 1 1
phK y ph x f K
y x f K
K K h y y
i i
i i
i i
? ? ?
?
? ? ? ? ? ?
),(),(),(
),(),(
),()(
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
21 ??? ??p
注意到,就是改进的欧拉法。
Q,为获得更高的精度，应该如何进一步推广？
其中 ?i ( i = 1,…,m ),?i ( i
= 2,…,m ) 和 ?ij ( i = 2,…,
m; j = 1,…,i?1 ) 均为待定
系数，确定这些系数的
步骤与前面相似。
§ 2 Runge-Kutta Method
),..,(
...,.,
),(
),(
),(
],.,[
1 1 2 2 1 1
2 32 1 31 3 3
1 21 2 2
1
2 2 1 1 1
? ?
?
? ? ? ? ? ?
? ? ? ?
? ? ?
?
? ? ? ? ?
m m m m m m i m
i i
i i
i i
m m i i
hK hK hK y h x f K
hK hK y h x f K
hK y h x f K
y x f K
K K K h y y
? ? ? ?
? ? ?
? ?
? ? ?
? 最常用为四级 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年给出了计算量与可达到的最高精
度阶数的关系,
7 5 3
可达到的最高精度
6 4 2 每步须算 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
2 1 0 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
