微分方程数值解法

定义 若某算法的局部截断误差为 En1 O(h p1) ，则称该算法
有p 阶精度。
欧拉法的局部截断误差,由Taylor展开:
En1

y(xn1)
yn1
[y(xn )
hy(xn )
h2 2
y(xn ) O(h3)][yn
hf
( xn ,
yn )]

h2 2
2. Taylor展开法
将y( xn h)在点xn作Taylor展开
y( xn
h)
y( xn ) hy'( xn )
h2 2!
y''( )
忽略高阶项,取近似值可得到Euler公式 yn1 yn h f ( xn , yn ) (n 0,1, 2, ... )
3. 数值积分法
stable */。
一般分析时为简单起见，只考虑试验方程 /* test equation */
y y
当步长取为 h 时，将某算法应用于上式，并假设只在初值 产生误差 0 y0 y0 ，则若此误差以后逐步衰减，就称该
算法相对于 h h 绝对稳定，h 的全体构成绝对稳定区域。
yn1 yn h f ( xn1, yn1 ) (n 0,1, 2, ... )
由于未知数 yi+1 同时出现在等式的两边，不能直接得到，故 称为隐式 /* implicit */ 欧拉公式，而前者称为显式 /* explicit */ 欧拉公式。
一般先用显式计算一个初值，再迭代求解。
整体截断误差除与 xn步计算有关外,还与 xn1, , x1的计算
有关 分析计算中的某一步,显式单步法的一般形式可写为:
yn1 yn hQ(xn , yn , h)
其中 Q(xn , yn , h) 称为增量函数。如对于Euler公式其增量函 数 Q(xn , yn , h) f (xn , yn ) y'(xn )
梯形公式 /* trapezoid formula */ — 显、隐式两种算法的平均
yn1

yn

h 2
[f
(xn
,
y
n
)

f
(xn1
,
y
n1
)]
(n 0,1, 2 L )
中点欧拉公式 /* midpoint formula */
中心差商近似导数
y( x1 )
y( x2 ) y( x0 ) 2h
lim(1 h)1/ h e
h0
y0e xn y( xn )
✓
关于整体截断误差与局部截断误差的关系,有如下定理
定理:对IVP(1)式的单步法 yn1 yn hQ( xn , yn , h,)
若局部截断误差为O(hp1) ( p 1) ,且函数 Q( xn , yn , h) 对y
我们称算法A 比算法B 稳定，就是指 A 的绝对稳定区域比 B 的大。
例：考察显式欧拉法 yn1 yn h yn (1 h)n1 y0
0 y0 y0
yn1 (1 h)n1 y0
n1 yn1 yn1 (1 h)n10
比较近似数值解：
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
由此可见,Euler公式的近似值接近方程的精确值.
二、构造初值问题数值方法的基本途径
以Euler法为例说明构造IVP问题数值方法的三种基本途径
将 方程y' f ( x, y)在 区间[ xn , xn1]上 积分
xn1 y'dx xn1 f ( x, y)dx (n 0,1,L )
xn
xn
用yn1 , yn代替y( xn1 ), y( xn ), 对右端积分采用 取左端点的矩形公式
则有
xn1 xn
f ( x, y)dx h f ( xn , yn )
ìy' = -y + x +1 í îy(0) =1
取步长h 0.1,计算到x 0.5
x³0
解 : f ( x, y) y x 1,由Euler公式
yn1 yn h( yn xn 1)
代入 h 0.1, 有 yn1 0.9 yn 0.1( xn 1), 依次计算得
例：就初值问题
y y

y(0)
y0
考察欧拉显式格式的收敛性。
解：该问题的精确解为 y( x) y0e x 欧拉公式为 yn1 yn h yn (1 h) yn
yn (1 h)n y0
对任意固定的 x = xn = nh ，有
yn y0 (1 h)xn / h y0[(1 h)1/ h ] xn
件，即存在与 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 xi1 xi (i 0, ... , n 1) 为步长，通常采用等距节点， 即取 hi = h (常数)。
yn1 yn h f ( xn , yn ) (n 0,1, 2, ... )
三、Euler法的改进及梯形公式
隐式欧拉法 /* implicit Euler method */
向后差商近似导数
y( x1 )
y( x1 ) h
y( x0 )
x0
x1
y( x1 ) y0 h f ( x1, y( x1 ))
例：用尤拉公式和改进的尤拉公式解初值问题

y'Βιβλιοθήκη y2x y
y(0) 1.
(0 x 1);
准确解:
y 1 2x
解：取步长h 0.1，
尤拉公式为：
yn1

yn h( yn

2 xn ). yn

y
p


yn
h( yn

2 xn ); yn
改进的尤拉公式为： yc

yn
h( y p

2 xn1 ); yp

yn1

1 2 ( yp

yc ).
计算结果略。
习题1: 分别用上述两种数值方法求出数值解，并与精确解比较。
四、单步法的误差分析和稳定性
1. 整体截断误差和局部截断误差 整体截断误差:数值解 yn和精确解 y( xn ) 之差
en y( xn ) yn
满足Lipschitz条件,即存在L>0,使得
Q(x, y, h) Q(x, y, h) L y y
对一切 y和y 成立，则该方法收敛，且有整体截断误 差
en O(h p )
注：由该定理可知整体截断误差总比局部截断误差低一阶
例：考察改进Euler法的收敛性。对改进的Euler法
Q( x, y, h) 1 f ( x, y) f ( x h, y hf ( x, y))
限定h即可知Q满足Lipschitz条件,故改进的Euler法收敛.
3. 稳定性
y( x) 30 y( x)
例：考察初值问题

y(0) 1
在区间[0, 0.5]上的解。
分别用欧拉显、隐式格式和改进的欧拉格式计算数值解,h=0.1。
节点 xi 0.0 0.1
0.2 0.3 0.4
求解(1)最基本的方法是单步法
单步法:从初值 y0 开始,依次求出 y1, y2 ,, 后一步的值 yn1 只依赖前一步的 yn
典型的单步法是Euler(欧拉)方法,其计算格式是:
yn1 yn h f ( xn , yn ) (n 0,1, 2, ... )
例:求解常微分方程初值问题
2 于是有
Q( x, y, h) Q( x, y, h) 1 { f ( x, y) f ( x, y) 2 f ( xh, y hf ( x, y)) f ( x h, y hf ( x, y)) }
设L为f关于y的Lipschitz常数,则由上式可得
Q(x, y, h) Q(x, y, h) L(1 h/ 2) y y
1. 数值微分法,用差商代替微商（Euler折线）
向前差商近似导数
y( x0 )
y( x1 ) h
y( x0 )
y( x1 ) y( x0 ) hy( x0 ) y0 h f ( x0 , y0 ) 记为 y1 x0
x1
yn1 yn h f ( xn , yn ) (n 0,1, 2, ... )
第一章 常微分方程初值问题 数值解法
第一节 基本原理，Euler法
一、初值问题的数值解
考虑一阶常微分方程的初值问题 /* Initial-Value Problem */:
dy f ( x, y) x [a, b] dx
(1)
y(a) y0
只要 f (x, y) 在[a, b] R1 上连续，且关于 y 满足 Lipschitz 条
y( x2 ) y( x0 ) 2h f ( x1 , y( x1 ))
x0