力学中的计算方法(常微分方程)

定义 若某算法的局部截断误差为 O(hp+1)，则称该算法有 p Ri 的主项 阶精度。 /* leading term */
欧拉法的局部截断误差：

h2 2
泰勒展开
Ri y( xi 1 ) yi 1 [ y( xi ) hy( xi ) h2 y( xi ) O( h3 )] [ yi hf ( xi , yi )]
f ( xi h, yi hK 3 )
§2 Runge-Kutta Method
注：
龙格-库塔法的主要运算在于计算 Ki 的值，即计算 f 的
值。Butcher 于1965年给出了计算量与可达到的最高精 度阶数的关系：
每步须算Ki 的个数
2
3
4
5
6
7
O ( h6 )
n8
O( hn2 )
4 3 4 5 可达到的最高精度 O( h2 ) O( h ) O( h ) O( h ) O( h )
由于龙格-库塔法的导出基于泰勒展开，故精度主要受
解函数的光滑性影响。对于光滑性不太好的解，最好 采用低阶算法而将步长h 取小。
变步长的龙格 - 库塔法
舍入误差↑
h →h
截断误差↓
与„xn, xn+1‟有关 以四阶龙格 - 库塔法为例： (h) 5 y ∵ Rn(xn+1)～O（h ），从xn出发计算xn+1的近似值 n 1
第三章 常微分方程数值解
/* Numerical Methods for Ordinary Differential Equations */
考虑一阶常微分方程的初值问题 /* Initial-Value Problem */:
dy f ( x, y) dx y ( a ) y0 x [a , b ]
3 要求 Ri y( xi 1 ) yi 1 O( h ) ，则必须有：
1 1 2 1 , 2 p 2
这里有 3 个未知 数， 2 个方程。
存在无穷多个解。所有满足上式的格式统称为2阶龙格 - 库 塔格式。注意到，p 1, 1 2 1 就是改进的欧拉法。
假设 yi 1 y( xi 1 ), yi y( xi ) ，则可以导出 Ri y( xi 1 ) yi 1 O(h3 ) 即中点公式具有 2 阶精度。
§1 Euler’s Method
方 法 显式欧拉 隐式欧拉 梯形公式 中点公式

简单 稳定性最好 精度提高 精度提高, 显式
h x0 记为 y( x1 ) y( x0 ) hy( x0 ) y0 h f ( x0 , y0 ) y1
x1
定义 在假设 yi = y(xi)，即第 i 步计算是精确的前提下，考 虑的截断误差 Ri = y(xi+1) yi+1 称为局部截断误差 /* local
truncation error */。
... ... Km f ( xi m h, y m1 hK1 m2 hK2 ... m m1 hKm1 )
其中i ( i = 1, …, m )，i ( i = 2, …, m ) 和 ij ( i = 2, …, m; j = 1, …, i1 ) 均为待定 系数，确定这些系数的 步骤与前面相似。
2
y( xi ) O( h3 )
欧拉法具有 1 阶精度。
§1 Euler’s Method
欧拉公式的改进：

隐式欧拉法 /* implicit Euler method */
y( x1 )
向后差商近似导数
y ( x1 ) y0 h f ( x1 , y ( x1 ))
y( x1 ) y( x0 ) h
§1 欧拉方法 /* Euler’s Method */
欧拉公式： 亦称为欧拉折线法 /* Euler’s polygonal arc method*/ y( x1 ) y( x0 ) 向前差商近似导数 y( x0 )
yi 1 yi h f ( xi , yi ) ( i 0, ... , n 1)
( h) 5 y( x n 1 ) yn Ch 1
将h折半，从xn出发计算两步得xn+1的近似值 y
y( x
(h ) 2 n1 n1 (h) 2 y ( x n 1 ) y n 1 ( h) y ( x n1 ) yn 1
) y
h 5 2C ( ) 2
2
Q: 为获得更高的精度，应该如何进一步推广？
§2 Runge-Kutta Method
yi1 yi h[ 1 K1 2 K2 ... m Km] K1 f ( xi , yi ) K2 f ( xi 2 h, yi 21 hK1 ) K3 f ( xi 3 h, yi 31 hK1 32 hK2 )
只要 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存 在唯一解。 要计算出解函数 y(x) 在一系列节点 a = x0< x1<…< xn= b 处的近似值 yi y( xi ) ( i 1, ... , n) 节点间距 hi xi 1 xi (i 0, ... , n 1) 为步长，通常采用等距节点， 即取 hi = h (常数)。
Step 2: 将 K2 代入第1式，得到
yi 1 y i h 1 y( x i ) 2[ y( x i ) phy( x i ) O ( h 2 )] y i (1 2 )h y( x i ) 2 ph2 y( x i ) O ( h 3 )
h yi [ f ( x i , yi ) f ( x i 1 , yi 1 )] 2
y i 1
h yi f ( xi , yi ) f xi 1 , yi h f ( xi , yi ) 2
( i 0, ... , n 1)
注：此法亦称为预测-校正法 /* predictor-corrector method */。 可以证明该算法具有 2 阶精度，同时可以看到它是个单 步递推格式，比隐式公式的迭代求解过程简单。后面将 看到，它的稳定性高于显式欧拉法。
隐式欧拉法的局部截断误差：
Ri y( xi 1 ) yi 1 h2 y( xi ) O( h3 )
2
即隐式欧拉公式具有 1 阶精度。
梯形公式 /* trapezoid formula */
§1 Euler’s Method
— 显、隐式两种算法的平均
h yi 1 yi [ f ( xi , yi ) f ( xi 1 , yi 1 )] ( i 0, ... , n 1) 2
§1 Euler’s Method
改进欧拉法 /* modified Euler’s method */
Step 1: 先用显式欧拉公式作预测，算出 y i 1 y i h f ( x i , y i ) Step 2: 再将 yi 1 代入隐式梯形公式的右边作校正，得到
y i 1
Step 1: 将 K2 在 ( xi , yi ) 点作 Taylor 展开
K 2 f ( x i ph, yi phK1 ) f ( x i , yi ) phf x ( x i , yi ) phK1 f y ( x i , yi ) O( h2 )
y( xi ) phy( xi ) O( h2 )
y i 1 K1 K2

1 1 yi h K 1 K 2 2 2 f ( xi , yi ) f ( x i h, y i hK 1 )
步长一定是一个h 吗？
§2 Runge-Kutta Method
将改进欧拉法推广为：
yi 1 K1 K2 yi h [1 K 1 2 K 2 ] f ( x i , yi ) f ( xi ph, yi phK 1 )
注：的确有局部截断误差 Ri y( xi 1 ) yi 1 O( h3 ) ， 需要22 个初值 y0和 y1来启动递推 即梯形公式具有 阶精度，比欧拉方法有了进步。 过程，这样的算法称为双步法 /* double-step 但注意到该公式是隐式公式，计算时不得不用到 method */，而前面的三种算法都是单步法 迭代法，其迭代收敛性与欧拉公式相似。 /* single-step method */。


§2 Runge-Kutta Method
Step 3: 将 yi+1 与 y( xi+1 ) 在 xi 点的泰勒展开作比较
yi 1 yi (1 2 )h y( xi ) 2 ph2 y( xi ) O( h3 )
h2 y( xi 1 ) y( xi ) hy( xi ) y( xi ) O( h3 ) 2
x0
x1
y i 1 y i h f ( x i 1 , yi 1 ) ( i 0, ... , n 1)
由于未知数 yi+1 同时出现在等式的两边，不能直接得到，故 称为隐式 /* implicit */ 欧拉公式，而前者称为显式 /* explicit */ 欧拉公式。 一般先用显式计算一个初值，再迭代求解。

精度低 精度低, 计算量大 计算量大 多一个初值, 可能影响精度
Can’t you givethink me a formula Do OK, you Well, callyet mewithout greedy… let’s with all the advantages any it possible? make it of the disadvantages? possible.
d f ( x, y) dx 首先希望能确定系数 1、2、p，使得到的算法格式有 2阶 dy 精度，即在 yi y( xi ) 的前提假设下，使得 f x ( x, y) f y ( x, y) dx Ri y( xi 1 ) yi 1 O( h3 ) f x ( x, y) f y ( x, y) f ( x, y) y( x )
中点欧拉公式 /* midpoint formula */
中心差商近似Hale Waihona Puke Baidu数
y( x1 )
y( x2 ) y( x0 ) 2h f ( x1 , y( x1 ))
y( x 2 ) y( x0 ) 2h
x0 x1 x2
yi 1 yi 1 2h f ( xi , yi ) i 1, ... , n 1
§2 龙格 - 库塔法 /* Runge-Kutta Method */
建立高精度的单步递推格式。
单步递推法的基本思想是从 ( xi , yi ) 点出发，以某一斜 率沿直线达到 ( xi+1 , yi+1 ) 点。欧拉法及其各种变形所 能达到的最高精度为2阶。
考察改进的欧拉法，可以将其改写为： 斜率 一定取K1 K2 的平均值吗？
最常用为四级4阶经典龙格-库塔法 /* Classical Runge-Kutta Method */ ：
y i 1 K1 K2 K3 K4 yi h ( K1 2K 2 2K 3 K 4 ) 6 f ( xi , yi )
h f ( xi h , y K1 ) i 2 2 h f ( xi h , y K2 ) i 2 2