数值分析复习---第八章 常微分方程数值解法

§8.2 欧拉(Euler)法
欧拉方法的收敛性
将 y(xn+1)在 xn 点 Taylor展开, ∃ξn ∈ (xn , xn+1),
y ( xn +1 )
=
y ( xn
+ h)
=
y(xn ) + hy′(xn ) +
h2 2
y′′(ξn )
=
y(xn ) + hf (xn , y(xn )) +
9/29
郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis
§8.2 欧拉(Euler)法
欧拉方法的收敛性
定义 若给定方法的局部截断误差满足 | Tn+1 |= O(h p+1),
则称该方法是 P 阶的，或称为具有 P 阶精度。
10/29
郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis
解的存在唯一性（“常微分方程”理论）：只要 f (x, y)
在[a, b] × R1 上连续，且关于 y 满足 Lipschitz 条件，即存在与 x, y 无关的常数 L 使
| f (x, y1)− f (x, y2)| ≤ L| y1 − y2 |
对任意定义在 [a, b] 上的 y1(x) 和 y2(x) 都成立，则上述IVP存 在唯一解。
§8.3 改进欧拉(Euler)方法
先用欧拉公式(8.2)求出一个初步的近似值 y n +1 ,称为预测值, 它的精度不高, 再用梯形公式对它校正 一次,即迭代一次,求得yn+1,称为校正值, 这种预测-校 正方法称为改进的欧拉公式：
欧拉方法的收敛性
对yn+1, yn+1应用欧拉公式(8.2)得
| yn+1 − yn+1 |=| y(xn ) + h f (xn , y(xn )) − ( yn + h f (xn , yn )) |
≤| y(xn ) − yn | +h | f (xn , y(xn )) − f (xn , yn ) |
∵
y n +1
=
yn
+
h[ 2
f
(xn ,
yn ) +
f (x, y) = − y ∴
f (xn+1 , yn+1 )]
yn+1
=
yn
+
h 2
[−
yn
− yn+1]
整理成显式
y n +1
=
⎜⎛ ⎝
2 2
− +
h h
⎟⎞ ⎠
y
n
反复迭代,得到
y n+1
=
⎜⎛ ⎝
2 2
− +
h h
⎟⎞ ⎠
y
n
=
⎜⎛ 2 − h ⎟⎞ 2 ⎝2+h⎠
向前欧拉公式和向后欧拉公式：
yn+1 = yn + h f (xn , yn )， yn+1 = yn + h f (xn+1, yn+1)。
6/29
郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis
§8.2 欧拉(Euler)法
若对积分用梯形公式，则得
[ ] y(xn+1)
h2 2
y′′(ξn )
(1)
假定已知准确值 y(xn ),利用欧拉公式(8.2)，定义
yn+1 y(xn ) + h f (xn , y(xn ))
(2)
8/29
郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis
§8.2 欧拉(Euler)法
局部截断误差
xn+1点向后差商近似导数
y′(xn+1)
≈
y(xn+1) − h
y(xn )
y(xn+1) ≈ y(xn ) + hy′(xn+1)
↑ y(xn ) ≈ yn代入
y( xn+1) ≈ yn+1 = yn + h f ( xn+1, yn+1)
14/29
隐式或后退欧拉公式
郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis
≤ T + (1+ hL)[T + (1+ hL) | en−1 |]
≤ T + (1+ hL)T + (1+ hL)2 | en−1 |≤ "
≤ T + (1+ hL)T + (1+ hL)2T +" + (1+ hL)nT
+ (1+ hL)n+1 | e0 |
≤ (1+ hL)n+1 −1T = (1+ hL)n+1 −1 O(h2 )
§8.2 欧拉(Euler)法
整体截断误差
记en+1 = y(xn+1) − yn+1.
因为计算 yn+1 时, 用到的y1, y2,",yn是 y(x1), y(x2 ),",
y(xn ) 的近似值,每步产生的误差会累积到计算y(xn+1)的误差
中,因此 en+1 与 y1, y2,",yn 都有关,称为整体截断误差.
2/29
郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis
§8.2 欧拉(Euler)法
通常取 xi+1 − xi = hi = h (常数),则Euler法的计算格式
⎧ ⎨ ⎩
yi+1 = yi + y0 = y(x0 )
hf
(
xi
,
yi
)
i=0,1,…,n
§8.2 欧拉(Euler)法
若用向后差商近似导数，即
y′(xn+1) = f (xn+1, y(xn+1))
y′( xn+1 ) ≈
y( xn+1 ) − h
y( xn )
y(xn+1) − h
y( xn )
≈
f (xn+1, y(xn+1))
y( xn+1 ) ≈ y( xn ) + hf ( xn+1 , y( xn+1 ))
17/29
郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis
例8.2.3
对初值问题
⎧ ⎨ ⎩
y′ + y (0 )
y= =1
0
证明用梯形公式求得的近似解为
yn
=
⎜⎛ 2 − ⎝2+
h ⎟⎞ n h⎠
并证明当步长h→0时,yn收敛于精确解 e − x
证明: 解初值问题的梯形公式为
§8.2 欧拉(Euler)法 向后欧拉公式 比较欧拉显式公式和隐式公式及其局部截断误差
显
式 yn+1 = yn + h f (xn , yn ) n = 0, 1,...
公 式
Tn+1
=
y(xn+1) −
yn+1
=
h2 2
y′′(xn ) + O(h3)
隐 式
yn+1 = yn + h f ( xn+1, yn+1)
y n −1
=
⎜⎛ ⎝
2 2
− +
h h
⎟⎞ ⎠
3
yn−2
= ...
=
⎜⎛ ⎝
2 2
− +
h h
⎟⎞ ⎠
n +1
y
0
∵
y0 = 1
∴yn
=
⎛ ⎜⎝
2 2
− +
h h
⎞n ⎟⎠
→
e−x
18/29
郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis
欧拉显 式公式
欧拉隐 式公式
# 梯形公式 /* trapezoid formula */ — 显、隐式两种算法的平均
yn+1
=
yn
+
h[ f 2
(xn ,
yn ) +
f
(xn+1, yn+1)]
注：的确有局部截断误差 Rn+1 = y(xn+1) − yn+1 = O(h3 ) ，
即梯形公式具有2 阶精度，比欧拉方法有了进步。但 注意到该公式是隐式公式，计算时不得不用到迭代 法，其迭代收敛性与欧拉公式相似。
公 式
Tn+1
=
y(xn+1) −
yn+1
=
−
h2 2
y′′(xn ) + O(h3)
16/29
郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis
若将这两种方法进行算术平均，即可消除误差 的主要部分/*leading term*/
而获得更高的精度,称为梯形法
( xn ,
y( xn ))
y(xn+1) ≈ y(xn) + hy′(xn) = y(xn) + hf (xn, y(xn))
⎧ ⎨ ⎩
yn+1 y0 =
= yn + y(a)
hf
(
xn
,
yn
)
差分方程初值问题 向前Euler方法
4/29
郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis
梯形 公式
欧拉法小结
公式
局部截断误差
( ) h2 y(2)
2
xn
− h2 2
y(2)
( xn )
( ) h3 y(3)
3
xn
精显稳 度 隐 定 步数
性 1阶 显 差 单步
1阶 隐 好 单步
2阶 隐 好 单步
19/29
郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis
§8.2 欧拉(Euler)法 向后欧拉公式
由于未知数 yn+1 同时出现在等式的两边，故称为隐式 /* implicit */ 欧拉公式，而前者称为显式 /* explicit */ 欧拉公 式。隐式公式不能直接求解，一般需要用Euler显式公式 得到初值，然后用Euler隐式公式迭代求解。因此隐式公 式较显式公式计算复杂，但稳定性好（后面分析）。
隐式欧拉公式中的未知数 yn+1 可通过以下迭代法 求解：
⎧ ⎨ ⎩
y（0) n +1
=
y = (k +1) n +1
yn + h yn + h
f (xn , yn )
f
( xn+1,
y(k) n +1
)
15/29
郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis
−
y(
xn
)
≈
h 2
f (xn, y(xn))+
f (xn+1, y(xn+1))
[ ] ⎪⎧
⎨
yn+1
≈
yn
+
h 2
f (xn, yn ) + f (xn+1, yn+1)
⎪⎩ y0 = y(a)
梯形欧拉公式
7/29
郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis
| en+1 |=| y(xn+1) − yn+1 |≤| y(xn+1) − yn+1 | + | yn+1 − yn+1 |
≤| Tn+1 | + | yn+1 − yn+1 | .
(3)
11/29
郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis
§8.2 欧拉(Euler)法
( 8.2 )
3/29
郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis
§8.2 欧拉(Euler)法
(1) 用差商近似导数
y′(xn) = f (xn, y(xn))
y′(xn) ≈
y(xn+1) − h
y( xn )
y(xn+1) − h
y( xn )
≈
Βιβλιοθήκη Baidu
f
Tn+1 = y(xn+1) − yn+1
= y(xn+1) − ( y(xn ) + hf (xn , y(xn )))
=
h2 2
y′′(ξn ).
称为局部截断误差
令 M = max | y′′(x) |,则 a≤ x≤b
| Tn+1
|=
h2 2
|
y′′(ξn ) |≤
h2 2
M
= O(h2 ).
(1+ hL) −1
hL
≤ O(h).
由此知，当 h → 0, 有 | en+1 |→ 0, 欧拉公式是一阶收敛的。
13/29
郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis
§8.2 欧拉(Euler)法 向后欧拉公式
#隐式欧拉法或向后欧拉法
#/* implicit Euler method or backward Euler method*/
≤| y(xn ) − yn | +hL | y(xn ) − yn | (李普希兹条件)
≤ (1+ hL) | en | .
(4)
12/29
郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis
记
T
=
max k
| Tk
|, 将(4)代入(3)得
| en+1 |≤ T + (1+ hL) | en |
郑州大学研究生课程 （2011-2012学年第一学期）
数值分析 Numerical Analysis
习题课 第八章 常微分方程数值解法
一、要点回顾
待求解的问题：一阶常微分方程的初值问题 /* Initial-Value Problem */:
⎪⎧ dy = f ( x, y) x ∈[a, b] ⎨ dx ⎪⎩ y(a) = y0
xn
xn
∫ y( xn+1 ) − y( xn ) =
xn+1 f ( x, y( x))dx
xn
分别用左矩形和右矩形公式，即
f (x, y) ≈ f (xn , y(xn )), f (x, y) ≈ f (xn+1, y(xn+1)), 代替上式右端的积分,并注意 y(xn ) ≈ yn,y(xn+1) ≈ yn+1分别得到
⎧ ⎨ ⎩
yn+1 y0 =
= yn + y(a)
hf
(
xn+1,
yn+1
)
向后Euler方法
5/29
郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis
§8.2 欧拉(Euler)法
∫ ∫ （2）用数值积分方法
xn+1 y′( x)dx = xn+1 f ( x, y( x))dx