计算方法(12)第八章 常微分方程(1)

y ( xn )
1.000000 1.004837 1.018731 1.040818 1.070320 1.106531
y ( xn ) yn
1.6 10 4
2.9 10 4
4.0 10 4
4.8 10 4
5.5 104
四、数值方法的误差估计、收敛和稳定性 1. 整体截断误差和局部截断误差 整体截断误差: 数值解 yn 和精确解y( xn )之差 n y( xn ) yn 整体截断误差除与 xn 步计算有关外,还与 xn1 ,, x1 的计算有关. 分析计算中的某一步,显式单步法的一般形式可写为:
取步长 h 0.1, 计算到 x 0.5
解 : f ( x, y ) y x 1,由Euler 公式
yn1 yn h( yn xn 1)
代入h 0.1, 有yn1 0.9 yn 0.1( xn 1), 依次算得果如下:
n0 1 2 3 4 5 xn 0 0.1 0.2 0.3 0.4 0.5 yn 1.0 1.0 1.01 1.029 1.0561 1.09049
只依靠前一步的 yn ，是一种逐点求解的离散化方法。
典型的单步法是Euler(欧拉)方法,其计算格式是:
yn1 yn h f ( xn , yn ) (n 0,1,2, ... )
例8-1:求解常微分方程初值问题 y' y x 1 x0 y ( 0) 1
yn1 yn hQ( xn , yn , h)
其中 Q( xn , yn , h) 称为增量函数。如对于Euler公式其 增量函数 Q( x , y , h) f ( x , y ) y' ( x )
n n n n n
定义 9-1 9-1）式的解，用单 8-1 设 y ( x) 是初值问题（ 8-1
2 ' x ' y
h [ yn f ( xn , yn )h [ f ( xn , yn ) f ( xn , yn ) f ( xn , yn )] O( h 3 )] 2
Rn y( xn 1 ) yn 1 [ y( xn ) hy( xn )
h2 2
公式是 p 阶的， p 是正整数。
) ，则称此数值
欧拉法的局部截断误差,由Taylor展开: Rn y( xn1 ) yn1
[ y( xn ) hy( xn )
h2 2
y( xn ) O(h 3 )] [ yn hf ( xn , yn )]
欧拉法具有 1 阶精度。
(0) y Step 2: 再将 n1 代入隐式梯形公式的右边作校正，
得到
h (0) yn 1 yn ( f ( xn , yn ) f ( xn 1 , yn 1 ) 2
(0) yn 1 yn hf ( xn , yn ) h (0) y y ( f ( x , y ) f ( x , y n 1 n n n n 1 n 1 ) 2
欧拉折线法
数值分析
三、Euler公式的改进及梯形公式 1、隐式欧拉法 向后差商近似导数
y( x1 ) y( x1 ) y( x0 ) h
x0 x1
y ( x1 ) y0 h f ( x1 , y ( x1 ))
yn1 yn h f ( xn1 , yn1 ) (n 0,1,2, ... )
此法亦称为预测-校正法 。一方面它有较高精度，同时 可以看到它是个单步递推格式，比隐式公式的迭代 求解过程简单。后面将看到，它的稳定性高于显式 欧拉法。
y n 1 h yn f (xn , y n ) f xn 1 , y n h f (xn , y n ) 2
用 y1 y0 f ( x0 , y0 )( x1 x0 ) y0 hf ( x0 , y0 ) 近似代替 y ( x1 )
1. 数值微分法,用差商代替微商 得到Euler公式
亦称为欧拉折线法
yn1 yn h f ( xn , yn ) (n 0,1,2, ... )
h2 y( xn ) O( h3 ) 2
类似可以证明欧拉预估-校正公式具有2阶精度 h yn 1 yn ( f ( xn , yn ) f ( xn 1 , yn hf ( xn , yn )) 2 上式在( xn , yn )处的Taylor展开式为
yn 1 h yn { f ( xn , yn ) [ f ( xn , yn ) hf x' ( xn , yn ) 2 hf y' ( xn , yn ) f ( xn , yn )]} O( h3 )
步法计算到第 n 步没有误差，即 yn y ( xn ) ，则
为单步法在点 xn 1 处的局部截断误差。
Rn y( xn1 ) yn1 y( xn1 ) y( xn ) hQ( xn , y( xn ), h)
p 1
O (h 定义8-2若数值公式的局部截断误差 E R n n1
3、 中点欧拉公式
中心差商近似导数
y( x 2 ) y( x0 ) y( x1 ) 2h
y( x2 ) y( x0 ) 2h f ( x1 , y( x1 ))
yn1 yn1 2hf(xn ,yn ) n 1, 2,
x0 x1 x2
2、梯形公式
h y n 1 y n [f(xn , y n ) f(xn 1 , y n 1 )] (n 0,1, 2 ) 2
可得
h yn1 yn h (1 )( xn yn ) 1 2
yn 1 0.905 yn 0.095 xn 0.1
以 h 0.1 代入得
xn
0 0.1 0.2 0.3 0.4 0.5
yn
1.000000 1.004762 1.018594 1.040633 1.070096 1.106278
2 h ' ' 3 yn f ( xn , yn )h [ f x ( xn , yn ) f y ( xn , yn ) f ( xn , yn )] O(h ) 2 Rn y( xn 1 ) yn 1
[ y( xn ) hy( xn )
h2 2
y( xn ) O( h 3 )]
(0) yn 1
(n 0, 1, 2
)
例8-2 用欧拉预估-校正公式 解初值问题 y ' y x 1 y(0) 1 取 h 0.1 ，计算到 x 0.5 。 解：利用
h yn 1 yn ( f ( xn , yn ) f ( xn 1 , yn hf ( xn , yn )) 2
由于未知数 yn+1 同时出现在等式的两边，不能直接得 到，故称为隐式 欧拉公式，而前者称为显式 欧拉公 式。 一般先用显式计算一个初值，再迭代求解。
2、梯形公式
h y n 1 y n [f(xn , y n ) f(xn 1 , y n 1 )] (n 0,1, 2 ) 2 ------- 显、隐式两种算法的平均
直接解微分方程可得精确解 : y f ( x ) x e x , 故x5 0.5, y(0.5) 1.106531
由此可见,Euler公式的近似值接近方程的精确值.
二、构造初值问题数值方法的基本途径 以Euler法为例说明构造IVP问题数值方法的三种基 本途径 1. 数值微分法,用差商代替微商
向前差商近似导数
记为
y( x1 ) y( x0 ) h y( x1 ) y( x0 ) hy( x0 ) y0 h f ( x0 , y0 ) y( x0 )
x0
பைடு நூலகம்
x1
y1
x0 x1
过点 ( x0 , y0 ) ，以 f ( x0 , y0 ) 为切线斜率的
切线方程为
y y0 f ( x0 , y0 )( x x0 )
(0) yn 1 yn hf ( xn , yn ) ( k 1) h (k ) y y ( f ( x , y ) f ( x , y n n n n 1 n 1 )) n1 2
(k 0,1, )
4、欧拉预估-校正公式 Step 1: 先用显式欧拉公式作预测，算出 (0) yn 1 yn hf ( xn , yn )
要计算出解函数 y(x) 在一系列节点 a = x0< x1<…< xn= b 处的近似值 yi y( xi ) (i 1, ... , n) 节点间距 hi xi 1 xi (i 0, ... , n 1) 为步长，通常采用等距节点， 即取 hi = h (常数)。
求解(8-1) (8-2)最基本的方法是单步法 单步法:从初值 y0 开始,依次求出 y1 , y2 , ,后一步的值 yn1
y yn f ( xn , yn )( x xn )
将x =xn +1处该直线上的函数值作为y( xn +1 )的近似值yn +1， 则
yn1 yn f ( xn , yn )( xn1 xn ) yn h f ( xn , yn ) (n 0,1,2, ... )
f ( x , y )dx h f ( xn , yn )
yn1 yn h f ( xn , yn ) (n 0,1,2, ... )
数值分析
4. 几何方法
设已求得解y =y( x )在xn处的近似值yn，则过点( xn , yn ) 以f ( xn , yn )为斜率的直线方程为
y( xn ) O( h 3 )]