常微分方程数值解法ppt

(
h p1 p 1)!
y(
p1) (
)
即(7.2-1)为p阶公式，上述公式称为Taylor公式。
注：利用Tayler公式构造，不实用，高阶导数f (i)不易计算。
16
7.2.2 RungeKutta方法
1. 基本思想
因为
y(xi1) y(xi )
xi 1 xi
f (x, y(x))dx
= y(xi) + hf (，y())
2
7.0 基本概念
1. 一阶常微分方程的初值问题
y(x) f (x, y(x)), x (a, b)
y(a)
y0
(7.0-1)
注：若f 在D = {a x b , | y |<+}内连续，且满足Lip条件：
L 0，使
|f (x，y1) – f (x，y2)| L|y1 – y2|
则(7.0-1)的连续可微解y(x)在[a，b]上唯一存在。
18
3. 常数的确定
确定的原则是使精度尽可能高。
以二阶为例：
yi1 yi h(c1K1 c2 K2 )
K1
f (xi , yi )
K2 f ( xi a2h, yi b21hK1)
(希望y(xi+1) – yi+1 = O(hp)的阶数p尽可能高)
首先：
y( xi1 )
y(xi )
| (x, y, h) (x, ~y, h) | L | y ~y |
且其局部截断误差为hp+1阶，则其(整体)截断误差为hp阶，
即该数值解公式为p阶方式。
注：① 局部截断误差较易估计, 定理7.1-1表明：
若ei = O(hp+1) 则i = O(hp).
② Euler局部截断误差为 ei1
h2 2
f(xi，yi) = yi + (2 + 0.1i)yi2
于是Euler计算公式为
yi+1 = yi + 0.1[yi + (2 + 0.1i)yi2]，i = 0，1，2，3，4 注：Euler方法精度较低
ex91.m
13
例2：用改进Euler方法求解初值问题：
y'
(
x)
Байду номын сангаас
1 x
( y( x)
y2 (x)),
y( )
O(h2 ) ,
所以一阶精度。
向后Euler法也是一阶精度。
③ 梯形公式为二阶精度。
12
例1：用Euler方法求解初值问题：
y' (x) y(x) (1 x) y2 (x), 1 x 1.5
y(1) 1
取步长h
=
0.1，并与准确解
y(x)
1 x
比较
解：因为xi = 1 + 0.1i，而f(x，y) = y + (1+x)y2，故
y( xi1 )
y(xi )
h 2
[
f
(
xi
,
y(xi ))
f
( xi1,
y( xi1 ))]
yi1
yi
h[ 2
f
(xi ,
yi )
f
( xi 1 ,
yi 1 )]
(7.1-4)
称(7.1-4)为梯形公式隐式公式。
显化：预估值：
yi
1
yi
hf
(xi ,
yi )
校正值：
yi
1
yi
h[ 2
f
(xi , yi)
hy' (xi )
h2 2!
y" (xi ) ...
hp p!
y( p) (xi )
h p1 y ( p1) ( )
( p 1)!
当h充分小时，略去Taylor公式余项，并以yi、yi+1分别代替
y(xi)、y(xi+1)，得到差分方程：
h2 yi1 yi hf (xi , yi ) 2! f '(xi , yi )
c1 c2 1
c2a2
1 2
b21 1 a2
特例：c1 = 0 c2 = 1，a2 = 1/2，b21 = 1/2，得：
yi1
yi
hK 2
K1
f (xi , yi )
K2
f (xi
h 2
,
yi
h 2
K1 )
称为中点公式。
(7.2-7)
22
4. 最常用的R-K公式 —— 标准4阶R-K公式
4
③ 二阶常微分方程y''(x) = f (x，y(x)，y'(x))可设为一阶常微 分方程组的初值问题： 引进新的未知函数z(x) = y'(x)，则
y '(x) z(x)
z
'(
x)
f
(x,
y(x),
z(x))
其初始条件为：
y(a) y0
z(a)
y'0
称为一阶微分方程组的初值问题，方法类似。 ④ 边界问题，常用差分方法解。
( p 1)!
当h充分小时，略去Taylor公式余项，并以yi、yi+1分别代替
y(xi)、y(xi+1)，得到差分方程：
yi1
yi
hf
(
xi
,
yi
)
h2 2!
f ' (xi , yi ) ...
hp p!
f ( p1) ( xi , yi )
(7.2-1)
其局部截断误差为：
y( xi1)
~yi1
可望得到较高精度的数值解，从而避免求f 的高阶导数。
17
2. RK公式
p
yi1 yi h c j K j
K1
j 1
f (xi , yi )
j 1
Kj
f (xi
a jh, yi
h bjsKs ),
s1
2 j p
(7.2-4)
其中Kj为y = y(x)在xi + ajh (0 aj 1)处的斜率预测值。 aj，bjs，cj为特定常数。
hp p!
f
( p1) (xi ,
yi )
其局部截断误差为：
y( xi1)
~yi1
(
h p1 p 1)!
y ( p1) ( )
(7.2-1)
15
y( xi1) y( xi h)
y(xi ) hy' (xi )
h2 2!
y" (xi ) ...
hp p!
y( p) (xi )
h p1 y ( p1) ( )
(1 1
yi ) yi 0.1i
(1 yi1 ) yi1 1 0.1(i 1)
i = 0，1，2，3，4
注：改进Euler方法精度比Euler方法精度高
14
7.2 RungeKutta方法
7.2.1 构造高阶单步法的直接方法
由Taylor公式：y(xi1) y(xi h)
y(xi )
(7.0-2)
3
2. 初值问题的数值解
称(7.0-1)的解y(x)在节点xi处的近似值 yi y(xi) a < x1 <x2 < ... < xn = b.
为其数值解，其求解方法称为数值方法。 注： ① 考虑等距节点: xi = a + ih，h = (b – a)/n. ② 从初始条件y(a) = y0出发，依次逐个计算y1，y2，…，yn的 值，称为步进法。 两种：单步法、多步法。
y( xi
)
y( xi
).
得差分方程：
yi1 yi h
f (xi , yi )
yi+1 = yi + hf (xi，yi). 即为Euler公式。
若记
y(xi1 ) y(xi ) h
y(xi1 )
f (xi1, y(xi1 ))
yi+1 = yi + hf (xi+1，yi+1).
称为向后Euler法。
数值计算方法
第七章 常微分方程数值解法
1
在工程和科学计算中，所建立的各种常微分方程 的初值或边值问题，除很少几类的特殊方程能给出解 析解，绝大多数的方程是很难甚至不可能给出解析解 的，其主要原因在于积分工具的局限性。
y(x) x2 y2
因此，人们转向用数值方法去解常微分方程，并 获得相当大的成功，讨论和研究常微分方程的数值解 法是有重要意义的。
1
x
1.5
y(1) 0.5
取步长h = 0.1，并与准确解
y(x) x 1 x
比较
解：xi
=
1
+
0.1i，f
(xi ,
yi )
1 xi
( yi
yi2 )
(1 yi ) yi 1 0.1i
于是改进Euler法的计算公式为
yi
1
yi
0.1(1 yi ) yi 1 0.1i
yi
1
yi
0.1 2
= y(xi) + hK
其中K = f (，y())称为y(x)在[xi，xi+1]上的平均斜率。
若取 K1 = f (xi，y(xi)) ——Euler公式
取 K2 = f (xi+1，y(xi+1)) —— 向后Euler公式 一阶精度
取
1 2
(
K1
K2
)
—— 梯形公式
二阶精度
猜想：若能多预测几个点的斜率，再取其加权平均作为K，
f
( xi1,
yi1 )]
称(7.1-5)为改进的Euler公式(显示格式)
(7.1-5)
9
4. 几何意义 Euler法折线法 改进Euler法平均斜率折线法
10
7.1.2 截断误差与代数精度
定义7.1-1
① 称 i = y(xi) – yi 为数值解yi的(整体)截断误差。
② 若yk = y(xk)，k = 0,1,2,…,i – 1. 由求解公式得数值解 ~yi y(xi )，称 ei y(xi ) ~yi为yi的局部截 断误差。
y(xi+1) y(xi) + hf (xi+1，y(xi+1))
yi+1 = yi + hf (xi，yi) 向后Euler 公式
(7.1-3)
8
② 用梯形公式：
xi1 xi
f
(x,
y(x))dx
h[ 2
f
(xi ,
y(xi ))
f
( xi 1,
h3 y(xi1))] 12 f
( ,
y( ))
y(xi1) y(xi )
xi1 f (x, y(x))dx
xi
取不同的数值积分可得不同的求解公式，如：
① 用矩形公式：
xi 1 xi
f
(x,
y(x))dx
hf
(xi ,
y(xi ))
hf
( xi 1 ,
y(xi1 ))
y(xi+1) y(xi) + hf (xi，y(xi))
yi+1 = yi + hf (xi，yi) Euler 公式
x
y
)f y
(x0, y0 )
)
K2 = f (xi，yi) + a2hf x'(xi，yi) + b21hK1 f y'(xi，yi) + O(h2).
代入(7.2-5)得：
yi+1 = yi + hc1 f (xi，yi) + hc2 f (xi，yi) + h2c2[a2 f x'(xi，yi) + b21K1 f y'(xi，yi)] + O(h3)
hy' (xi )
1 h2 2!
y" (xi )
O(h3 )
(7.2-5)
19
另一方面： 将K2在(xi，yi)处展开。
yi1 yi h(c1K1 c2 K2 )
K1
f (xi , yi )
K2 f ( xi a2h, yi b21hK1)
(f
( x0
x, y0
y)
f
(x0, y0 ) (x
注：局部截断误差ei是指单步计算产生的误差，而(整体)截
断误差i则考虑到每步误差对下一步的影响。
11
定义7.1-2 若求解公式的（整体）截断误差为O(hp)，则称该 方法是p阶方法，或是p阶精度。<p阶公式>
定理7.1-1 设数值解公式：yi+1 = yi + h(xi，yi，h)中的函数
(x，y，h)关于y满足Lipschitz条件：
(7.1-2)
注：① Euler法为显式，向后Euler法为隐式——须解出yi+1.
② 可用迭代法yi+1 (k+1) = yi + hf (xi+1，yi+1(k))
k = 0,1,2,… 解得yi+1 ，其中yi+1(0) = yi + hf (xi，yi).
7
3. 对(7.0-1)两边取积分得
c1 c2 1
c2a2
1 2
b21 1 a2
特例：a2 = 1 c1 = c2 = 1/2，b21 = 1，得2阶R-K公式：
yi
1
K1
yi
h 2
(K1
f (xi , yi )
K2 )
改进欧拉公式
K2
f (xi
h, yi
hK1)
21
希望：ei+1 = y(xi+1) – yi+1 = O(h3). 则应：
5
7.1 初值问题数值解法的构造及其精度
7.1.1 构造方法
对于(7.0-1)可借助Taylor展开（导数法）、差商法、积分法 实现离散化来构造求积公式：
1. 设y C[a，b]将y(xi+1) = y(xi+h)在xi处展开
y(xi1 )
y(xi ) hy(xi )
h2 2
y( )
y(xi ) hf
= yi + h(c1 + c2) f (xi，yi) + c2a2h2[f x'(xi，yi) + (b21/a2) f (xi，yi) f y'(xi，yi)]+O(h3)
（希望）
yi h(c1 c2 ) y' (xi ) c2a2h2 y" (xi ) O(h3 )
20
希望：ei+1 = y(xi+1) – yi+1 = O(h3). 则应：
(xi , y(xi ))
h2 2
y( )
y(xi+1) yi+hf (xi，yi) 其中yi y(xi).
称yi+1 = yi + hf (xi，yi). i = 0,1,2, …,n – 1
为Euler求解公式，（Euler法）
[xi，xi+1]
(7.1-1)
6
2. 用差商来表示：
y( xi1) h