计算物理-6

y ( x) = y ( xn ) + y′( xn )( x − xn ) + O(( x − xn ) 2 ) y′′( xn ) y ( x) = y ( xn ) + y′( xn )( x − xn ) + ( x − xn ) 2 + O(( x − xn )3 ) 2
ε n +1 = O(h 2 )
常微分方程初值问题 ⎧ y′ = f ( x, y ), x0 < x < x f ⎨ ⎩ y ( x0 ) = y0 解函数y(x), x0≤x≤xf. 解析法：常系数线性微分方程可用特征根法， 一些特殊的微分方程可以用分离变量法等解. 绝大多数常微分方程都难以求得精确解.
⎧ y′ = x + y , 0 < x < 1 ⎨ ⎩ y (0) = 1
隐式Euler格式
y n +1 = y n + h( y n +1 − 2 x n +1 / y n +1 )
改进Euler格式
改进Euler法精度明显高
~ y n +1 = y n + h( y n − 2 x n / y n ) ⎧ ⎪ ⎨ y = y + h ( y − 2x / y + ~ ~ − y 2 x / n +1 n n n n n +1 n +1 y n +1 ) ⎪ 2 ⎩
§6.2.2 Runge-Kutta法
原理：Lagrange微分中值定理
y ( xn +1 ) − y ( xn ) ⇒ y ( xn +1 ) = y ( xn ) + hK * K = y′(ξ ) = h
*
平均斜率
问题已转化为如何对K*进行近似计算。
*
K
xn
xn+1
龙格-库塔(Runge-Kutta)法的基本思想 Euler公式可改写成：
yi = y ( xi )
或用向前差 商近似导数
y ( x n +1 ) − y ( x n ) y ′( x n ) ≈ h
x0 x1
y ( xn +1 ) ≈ y ( xn ) + hy ′( xn ) = yn + h f ( xn , yn )
y i +1 = y i + h f ( x i , y i ) ( i = 0, ... , n − 1)
f ( x + Δx ) − 2 f ( x ) + f ( x − Δ x ) 2 ′ ′ = f ( x ) + O (( Δ x ) ) 2 Δx
这说明二阶中心差商的精度也为二阶
6.2 常微分方程的数值解法
§6.2.1 Euler法及其改进 §6.2.2 Runge-Kutta法 §6.2.3 收敛性与稳定性（略） §6.2.4 微分方程组与高阶微分方程
二阶差商多取中心式，即
Δ2 y f ( x + Δx) − 2 f ( x) + f ( x − Δx) = 2 Δx (Δx) 2
以上是一元函数的差分与差商。多元函数f(x,y,…)的
差商与导数之间的误差表明差商逼近导数的程度，称为逼近误差。 由函数的Taylor展开，可以得到逼近误差相对于自变量差分（增量） 的量级，称为用差商代替导数的精度，简称为差商的精度。
改进欧拉法
Step 1: 先用显式欧拉公式作预测，算出 y i +1 = y i + h f ( x i , y i ) Step 2: 再将 yi +1
代入隐式梯形公式的右边作校正，得到
y i +1
h = yi + [ f ( x i , yi ) + f ( x i +1 , yi +1 )] 2
⎧ y i + 1 = y i + hK 1 ⎨ ⎩K 1 = f (xi, yi )
则yi+1的表达式与y(xi+1)的Taylor展开式的前两项 2 O ( h ) 。 完全相同,即局部截断误差为 改进的Euler公式 又可改写成：
h ⎧ ⎪ y i +1 = y i + 2 ( K 1 + K 2 ) ⎪ ⎨K 1 = f (xi , yi ) ⎪ K = f ( x , y + hK ) i+1 i 1 ⎪ 2 ⎩
Δy = f ( x + 1 1 Δx ) − f ( x − Δ x ) 2 2
（1-3）
（1-4）
Δx >0
上面谈的是一阶导数，对应的称为一阶差分。对一阶差分再作一 阶差分，所得到的称为二阶差分，记为 Δ2 y 。
以向前差分为例，有
Δ2 y = Δ(Δy ) = Δ[ f ( x + Δx) − f ( x)] = Δf ( x + Δx) − Δf ( x) = [ f ( x + 2Δx) − f ( x + Δx)] − [ f ( x + Δx) − f ( x)] = f ( x + 2Δx) − 2 f ( x + Δx) + f ( x)
2. 欧拉方法的导出 把区间［a,b］ 分为n个小区间
N等分
步长为 hi = ( xi +1 - xi )
节点 xi = a + ihi ,一般取 hi = h ( = (b − a ) / n )即等距
要计算出解函数 y(x) 在一系列节点
a = x0 < x1 < L < xn = b
处的近似值
一阶向后差商也具有一阶精度。
将
f ( x + Δx )
与
f ( x − Δx )
的Taylor展开式相减可得:
f ( x + Δx ) − f ( x − Δx ) = f ′( x) + O ((Δx) 2 ) 2 Δx
可见一阶中心差商具有二阶精度。
将 f ( x + Δx ) 与
f ( x − Δx ) 的Taylor展开式相加可得
2
§6.2.1 Euler法及其改进
1. 数值解
取步长h，记xn=x0+nh，n=1,2, …,N, 求得各节点 xn处函数y(xn)的近似值yn ，称y0，y1, …，yn, … yN为常微分方程的数值解。
yn y (x n )
y0 x1
x0
xn
各种Euler法格式
Euler格式 y n +1 = y n + hf ( x n , y n )，n=0,1,2, … 隐式Euler格式 y n +1 = y n + hf ( x n +1 , y n +1 ) ，n=0,1,2, … 二步Euler格式 y n +1 = y n −1 + 2hf ( x n , y n ) ，n=0,1,2, … 原理(1)微分方程⇒数值微分 原理(2)微分方程⇒积分方程⇒数值积分
y i +1
h = yi + [ f ( x i , yi ) + f ( x i +1 , yi + h f ( x i , yi ))] 2
( i = 0, ... , n − 1)
注：此法亦称为预测-校正法 。可以证明该算法具有 2 阶精 度，同时可以看到它是个单步递推格式，比隐式公式的 迭代求解过程简单。后面将看到，它的稳定性高于显式 欧拉法。
(Δx) 2 (Δx)3 (Δx) 4 IV f ( x − Δx) = f ( x) − Δx ⋅ f ′( x) + f ′′( x) − f ′′′( x) + f ( x) + O((Δx) 5 ), 2! 3! 4! f ( x) − f ( x − Δx) ∴ = f ′( x) + O(Δx) Δx
一阶中心差商为
Δy = Δx f (x + 1 1 Δx) − f ( x − Δx) 2 2 Δx
或
Δy f ( x + Δ x ) − f ( x − Δ x ) = 2 Δx Δx
一阶中心差商为
Δy = Δx f (x + 1 1 Δx) − f ( x − Δx) 2 2 Δx
或
Δy f ( x + Δ x ) − f ( x − Δ x ) = Δx 2 Δx
依此类推，任何阶差分都可由其低一阶的差分再作一阶差分得到。 例如n 阶前差分为
函数的差分与自变量的差分之比，即为函数对自变量的差商。 一阶向前差商为
Δy f ( x + Δx) − f ( x) = Δx Δx
（1-7）
一阶向后差商为
Δy f ( x) − f ( x − Δx) = Δx Δx
（1-8）
f ( x, y ) = f ( x n , y n ) + ( x − x n ) f x ( x n , y n ) + ( y − y n ) f y ( x n , y n )
+ O(( x − x n ) 2 + ( y − y n ) 2 )
3 Euler法的改进
原理：微分方程⇒积分 方程⇒数值积分⇒梯形 公式 梯形格式：
dy dx
是函数对自变量的导数，又称微商；
Δy 、 Δx 分别称为函数 及自变量的差分，
Δy 为函数对自变量的差商。 Δx
一阶导数，对应的称为一阶差分: 向前差分 向后差分 中心差分
Δy = f ( x + Δx ) − f ( x )
（1-2）
Δy = f ( x ) − f ( x − Δx )