常微分方程初值问题的的数值解法

本章讨论常微分方程初值问题的数值解法
2
考虑一阶常微分方程的初值问题
⎧ dy ⎪ = f ( x, y ) ⎨ dx ⎪ ⎩ y (a ) = y0
x ∈ [a, b]
只要 f (x, y) 在[a, b] × R1 上连续，且关于 y 满足 Lipschitz 条 件，即存在与 x, y 无关的常数 L 使对任意x∈[a, b] ，和y1， y2 ∈ R1 都有 | f ( x, y1) − f ( x, y2 ) | ≤ L| y1 − y2 | 在唯一解。 成立， 则上述问题存
⎧ ⎪ ⎨ ⎪ ⎩ y n +1 = yn + hf ( xn , yn ), h yn +1 = yn + [ f ( xn , yn ) + f ( xn +1 , y n +1 )] 2
改进的Euler方法：y0=1,
y1=y0+hf (x0, y0) =1.1, y1=1+01./2 ×[(1−2 ×0/1)+(1.1−2 ×0.1/1.1)] =1.095909, …… y11=…… y11=1.737869.
1 yn +1 = yn + h[ f ( xn , yn ) + f ( xn +1 , yn +1 )] 2
12
称之为梯形公式。这是一个隐式的计算公式，欲求的yn+1需 解一个方程。
3．截断误差
定义 在假设 yn = y(xn)，即第 n 步计算是精确的前提下，考 虑的截断误差 εn+1 = y(xn+1) − yn+1 称为局部截断误差
⎧ y n +1 = y n + k1 ⎨ ⎩k1 = hf ( xn ,y n )
n = 0,1,2, ... n = 0,1,2, ...
考察改进的欧拉法，可以将其改写为：
⎧ y n +1 = yn + hf ( xn , yn ) ⎪ ⎨ h ⎪ yn +1 = yn + [ f ( xn , yn ) + f ( xn +1 , y n +1 )] 2 ⎩ 1 ⎧ ⎪ y n +1 = y n + 2 ( k 1 + k 2 ) ⎪ = hf ( x n , y n ) ⎨ k1 ⎪k = hf ( x n + h , y n + k 1 ) 2 ⎪ ⎩
第八章 常微分方程初值问题的数值解法
实际中，很多问题的数学模型都是微分方程。利用微 分方程理论，我们可以研究它们的一些性质。但是，只有 极少数特殊的方程有解析解。对于绝大部分的微分方程是 没有办法求出它的解析解的。 常微分方程作为微分方程的基本类型之一，在自然界 与工程界有很广泛的应用。很多问题的数学表述都可以归 结为常微分方程的定解问题。很多偏微分方程问题，也可 以化为常微分方程问题来近似求解。
求解方法即采用离散化方法，（通常用数值积分、微 分、泰勒展开等）将上述初值问题化成关于离散变量 的相应问题。把这个相应问题的解yn作为y(xn)的近似 值。这样求得的yn就是上述初值问题在节点xn上的数 值解。 不同的离散化导致不同的求解方法。
6
1. 欧拉(Euler)方法(折线法)
几何 意义
a
x1 x2 x3 b
3
方程的几何意义:寻找积分曲线
y0
a
b
4
求解思路：
要计算出解函数 y(x)比较困难，改求曲线上的离散点 a = x0< x1<…< xn= b处的近似函数值
y i ≈ y ( xi ) ( i = 1, ... , n )
y(x)
a
x1
x2
x3
b
5
通常采用等距节点，步长
h = xi+1 − xi (i = 0, ... , n −1)
20
龙格 - 库塔法
高精度的单步递推格式
单步递推法的基本思想是从 ( xi , yi ) 点出发，以某一斜 率沿直线达到 ( xi+1 , yi+1 ) 点。欧拉法及其各种变形所 能达到的最高精度为2阶。
21
考察欧拉法，可以将其改写为：
y n +1 = y n + hf ( xn ,y n )
定义 若某算法的局部截断误差为 O(hp+1)，则称该算法有p 阶精度或称p阶方法。
13
) 欧拉法的局部截断误差：
将 y ( x n +1 )在 x n点展开成泰勒级数： y ( x n +1 ) = y ( x n ) + h y ′ ( x n ) +
h2 2
y ′′(ξ )
欧拉公式： y n +1 = y n + hf ( x n , y n )
n = 0,1,2, ...
8
例 用欧拉法求初值问题
0 .9 ⎧ y ⎪ y' = − 1+ 2x ⎨ ⎪ y( x0 ) = 1 x0 = 0 ⎩
当h = 0.02时在区间[0, 0.10]上的数值解。 方程真解：y( x ) = (1 + 2 x ) −0.45
y n +1 = y n + hf ( xn ,y n )
f ( xn , y n ) 1 h ⎞ ⎛ f ⎜ x n + h, y n + k1 ⎟ 2 2 ⎠ ⎝ 1 h ⎞ ⎛ f ⎜ xn + h, y n + k 2 ⎟ 2 2 ⎠ ⎝ f (x n + h, y n + hk 3 )
2 ⎧k1 = xn − yn ⎪ 2 = + − ( yn + 0.05k1 ) ( 0 . 05 ) k x ⎪ 2 n ⎨ 2 = + − ( yn + 0.05k 2 ) ( 0 . 05 ) k x ⎪ 3 n ⎪k = ( x + 0.1) 2 − ( y + 0.1k ) n n 3 ⎩ 4
T1 = y( xn+1 ) − yn+1 =
h2 2
y′′(ξ )
若y ( x)在[ a, b]充分光滑，令 M = max y′′( x) ，
a < x <b
h2 则 T1 ≤ M ⋅ = O (h 2 ) 2
欧拉法具有 1 阶精度。
14
) 梯形法的局部截断误差：
T2 = y( xn+1 ) − yn+1 = O(h3 )
改进欧拉法yn 1 1.095909 1.184096 1.260201 1.343360 1.416102 1.482956 1.552515 1.616476 1.678168 1.737869 1.795822 1.852242 1.907323 1.961253 2.014207
准确解 1 1.095445 1.183216 1.264911 1.341641 1.414214 1.483240 1.549193 1.612452 1.673320 1.732051 1.788854 1.843909 1.897367 1.949359 2.000000
问题: 为获得更高的精度，应该如何进一步推广？
23
将改进欧拉法推广为：
r ⎧ ⎪ y n +1 = y n + ∑ ω i k i i =1 ⎪ ⎪ ⎨ k 1 = hf ( x n , y n ) ⎪ i −1 ⎪ k i = hf ( x n + α i h , y n + ∑ β ij k j ), ( i = 2 ,3, " , r ) ⎪ j =1 ⎩
y ( xn +1 ) ≈ y ( xn ) + hf ( xn , y ( xn ))
yn +1 = yn + hf ( xn , yn ) ——Euler公式
2. 用梯形公式
∫
则有
xn+1
xn
h f ( x, y ( x))dx ≈ { f ( xn , y ( xn )) + f ( xn +1 , y ( xn +1 ))} 2
22
类似的有三阶Runge-Kutta公式：
⎧ ⎪ y n +1 ⎪ ⎪ k1 = ⎨ ⎪k = ⎪ 2 ⎪k = ⎩ 3 h = yn + (k1 + 4 k 2 + k 3 ) 6 f (xn , yn ) 1 1 f (xn + h, yn + hk 1 ) 2 2 f ( x n + h , y n − hk 1 + 2 hk 2 )
3
)
截断误差阶为 O(h5)。
25
例题: 取步长h = 0.1， 用四阶龙格―库塔法解初值问题的数 值解 y′ = x2 – y (0≤x≤1) y(0) = 1
解 ：初值问题的精确解为 y=x2−2x+2−e− x
用四阶龙格―库塔公式
⎧ k1 = ⎪ ⎪k 2 = ⎪ ⎨ ⎪k = ⎪ 3 ⎪ ⎩k 4 =
0 0.0005 0.0005 0.0014 0.0018 0.0021
10
2．梯形公式
一阶方程的初值问题与积分方程
y ( xn +1 ) = y ( xn ) + ∫
是等价的
xn+1
xn
f ( x, y ( x))dx
借助于数值积分，求y(xn+1)的值
11
1. 用左矩形公式
∫
xn+1
xn
f ( x, y ( x))dx ≈ hf ( xn , y ( xn ))
梯形法具有 2 阶精度。