常微分方程数值解法

用分段的折线逼近函数,此为 “折线法”而非“切线法”， 除第一个点是曲线上的切线，
其它都不是。
2、Euler方法的误差估计
1)局部截断误差。 在一步中产生的误差而非累积误差：
~
T x y y

n1
n1
n1
其中
~
y
是当
y
n

y(
x
)
n
（精确解！）时
n1
由Euler法求出的值，即y 无误差！ n
y y x , y x y 则得： h f
f ,
n1
n2
n
n
n1
n1
同样与Euler法结合，形成迭代算法，对n 0，1，2，
y y x , y 0 hf
n1
n
n
n
y y x , y x y

k 1
推出总体误差与步长的关系。
由微分方程解的存在唯一性,自然假定 （ f x，y）
充分光滑，或满足 Lipschitz条件：
f
x,ny源自xn

f
x
,
n
y
n

L
yxn
y n
第 n 步 的 总 体 截 断 误 差 记 为
en y
xn
y n
则 对 n 1 步：
e x y x y y y T y y ~ ~
用yn1, yn代替y(xn1), y(xn ), 对右端积分采用 取左端点的矩形公式
则有
xn1 xn
f
(x,
y)dx

hf
(xn ,
yn )
yn1 yn hf (xn , yn ) (n 0,1, )
C 在xn 附近 y(x) 的 Taylor 展开：
y(xn
h)
k 1
n1
hf
n
k
n1 n1
k 0,1, 2,
向后Euler 方法收敛条件与截断误差
y y x y

0

n1
hf
n
,
n
n
y y x , y

k 1
n1
hf
n
k
n1 n1
k 0,1, 2,
y y x y x y k1 n1
Euler 方法以Oh 速率收敛：h 0, eN 0.
4、微分方程数值解的稳定性
稳定性分析，对计算误差： y y*
n 1
n 1
n 1
其中 y* 是 y 的近似计算值，误差积累会淹没真值？
n 1
n 1
定义：一种数值方法求解 y' = y, 其中 是复常数。对给
定的步长 h 0，若计算误差 在计算 y k 1，2， 时
dy
yy 0
x x0
y x x y dx x0,y0
(x )f( , )
0
0
0
0
等步长 为h，则
x1
x0

h，可由切线算出 y 1
：
y 1

y 0
hf（x0，y
）
0
逐步计算出y

y（x）在
x ，点值 n 1
：
y n 1

y n

hf（xn，y
）
n
n 0，1，2，
hf
n
,
nn

y(
x
)
0

y 0

y(
x
)
0

y 0
1，几何意义。
n 0,1,
由
x
,
0
y
0
出发取曲线 y y x 的切线（存在！），则 斜率
dy
f
dx
x
,
0
y
0
x
,
0
y
0
由于 f
x
,
0
y
0
及
x
,
0
y
0
已知，必有切线方程。
由点斜式写出切线方程：
将
x y( ) n 1
在 xn 点 Taylor 展开：
x y( ) n 1

yxn

h
y(xn)
hf
(xn ,
y( xn))

h2 y//
2
x x
n
n 1
~
y
x yn hf
,
n
yn
n1
~
y

y
(
x
)
n

hf
x
,
n
y

x
n

n1
k
h f
n1
k
,
f
n1 n1
k 1
,
n1 n1
y y k
k 1
hL
n1
n1
收敛条件 0 hL 1
T h 局部截断误差
O 2 ，整体截断误差，
n1
当 0<hL<1 时：O h
向后Euler 法的稳定性
对 y/ y 用向后Euler法：y y h y
n
n
n1
y y x y hf ( , ) n 0, 1, 2,
n1
n
n
n
Taylor展开法不仅可得到求数值解的公式，且容易估计
截断误差。
§2 尤拉（Eular)方法
一、 Euler 公式 计算公式：
y y x y dy

dx

f
(x, y)
a

x

b



n1
n1
h f n2
n
f
n
, k k 0,1,2,
n 1
n 1
梯形公式的收敛性
由迭代算法，对n 0，1，2，
y y x , y

0

n1
hf
n
nn
y y x ,y x y

k 1

n1
n

h 2

f
f
nn
,
n1
y y
n1
n1
n1
n1
n1
~

n1
n1
n1
n1
n1
y y 以 下 估 计 ~ n1 n1
其 中~y

y

x
n


h
f

x
,
n
y

x
n


n1
总体方法误差(2)
~
y
y x
y
hf
n1
n
x
,
n
y

x
n

y hf n

0，
则：
eN

TN
e 1 hL N 1

T T
N 1 hL

N 1
1hL
2
T N 2

1hL
N 1
T1
由局部截断误差 T n O h2 ，则
eN T N T 1 hL n1 1hL N 1 T 1

y
(
x
n
)

h
y
/
(
x
)
n

h2
2
y x //
( ) n

y(xn) hf
(
x
,
n
y
(
x
))
n

h2
2
y x //
( ) n
dy


dx

f (x, y)


y
(
x
)
0

y 0
a

x

b



y x x 取 h 的线性部分，且 y( ) 得 y( ) 的近似值：
但每次迭代均多算一次函数值—提高精度的计算代价。
梯形公式的稳定性
y y x , y x y n1
h f n2
f
n
n
,
n1
n 1
对y/ y
用梯形公式
：y n1

y n

h2（ yn

y） n1
h
n 1
n
2
n
n1
1

x
,
n
y
n
n1
y xn yn h f
x
,
n
yxn
f
x
,
n
y
n
由Lipschitz条件
1 hL
y xn
y n

(1 hL)
en
e T e

n1
n1 1 hL n 对一切n 成立。
对取定 N ， 由
e0

y
x0

y 0
y x y 代替，则：
n1
n1
y y
n1
n
h
f
xn
,
y n
y y hf
n1
n
xn
,
y n
n 0,1, 2,
B. 数值积分
用数值积分方法离散化 :
xn1 dydx xn1 f (x, y)dx (n 0,1, )
xn dx
xn
dx
x
y(1) 3
只要函数f (x, y)适当光滑连续，且关于y满足李普 希兹(Lipschitz)条件，即存在常数L，使得
f (x, y) f (x, y) L y y 由常微分方程理论知，初值问题的解必存在且唯一。
微分方程的数值解：设方程问题的解y(x)的存在区间 是[a,b]，令a= x0< x1<…< xn =b,其中hk=xk+1-xk , 如是等距 节点h=(b-a)/n , h称为步长。 y(x)的解析表达式不容易得到或根本无法得到，我们用数 值方法求得y(x)在每个节点xk上y(xk)的近似值，用yk表示，
即yk≈y(xk)，这样y0 , y1 ,...,yn称为微分方程的数值解。
微分方程离散化常用方法
A 用差商代替微商
dy y
x x dx x x xn, yn
y
n1

n f (xn , y(xn ))
n1
n
用h xn1 xn ,
y n

yxn,
N1

1hL k
T h O 2 N K
N 1

1hL k
k 0
k 0
1hL N 1

1 hL 1
O h2 O h


lim h0
1hL
N
lim h0
1hL
xN x0 h
e L x N x0 与步长 h 无关常数
2424242424adams19如果令由方程组可解得相应的线性多步公式为称其为四阶隐式公式其局部截断误差为利用数值积分方法求线性多步公式基本思想是首先将初值问题化成等价的积分形式用过节点次插值多项式代替求积分例如时过节点的三次插值多项式为其中555937对上式用代替这就是四阶显式公式
第九章 常微分方程数值解法
k
n1

k 0,1,2,
y y x y x y k1 n1
k h f n1 2
k
,
f
n1 n1
k 1
,
n1 n1
y y hL
k
k 1

2
n1
n1
收敛条件 0 hL 1 2
梯形公式比 Euler法的局部与总体误差均高一阶，eN 0 h2 ，
T x y h y 则
y
n1
~
2

n1
n1 2
//
x x
n
n1
令 M 2 max y// (x) , y(x) 充分光滑，则： a xb
T M h h n 1
2
O 2 22
3、 总体方法误差(1)
递推方法：从任意两相邻步的总体误差关系

n1
n
1 h

1
2
h
2


1
Re (h)
1 4
h2

2
2

1
Re (h)

1 4
h2

2

当Re (h) 0时
n1
1，梯形公式是A - 稳定的。

n
四、改进的尤拉公式
梯形公式虽然提高了精度，但使算法复杂。而在实际计算中 只迭代一次，这样建立的预测—校正系统称作改进的尤拉公式。
n 1
n
n 1
误差公式： h
n1
n
n1
则

1
n1
1 h

n
1
1h 1h
1 2
1
1

1
1
1h h2 2 12Re(h)h2 2 2
只要 Re (h) 0 则
n1 1
y, xn 1

x y

n1
h
yxn
y y x y
为隐式算法：


n1
hf
n
,
n1 n1
y x0

y 0
为避免解非线性方程，与Euler 法结合：
y y x y

0

n1
hf
n
,
n
n
y y x , y

n
nk
不产生增大的误差，即 ，称对 h 与 这种方法是
n+k
n
绝对稳定的。
对 ，h 的允许范围内是绝对稳定的，
则称绝对稳定区域。
Euler法的绝对稳定区域
y/ y 的 Euler 算法：
Im(λh)
y y h y
n 1
n
n
计算值：y* y* h y*
n 1
预测
yn1 yn hf ( xn , yn );
校正
yn 1

yn

h 2
[
f
( xn
,
yn
)

f ( xn1, y n1)],
这一计算公式也可表为
h
yn 1

yn

[f 2
( xn , yn )
f
( xn
h, yn
hf
(xn , yn ))],
n
因此向后 Euler 法是A 稳定的。但收敛