Transcript 改进的Euler 法收敛

第二讲
常微分方程（组）
数值求解
1
主要内容
 Euler 法与改进的 Euler 法
 算法分析：误差与收敛性
 Runge-Kutta 法
 一阶常微分方程组与高阶方程
 应用举例
 Matlab 相关函数
2
初值问题
 考虑一维经典初值问题
 dy
 f ( x, y)

 dx

 y( a )  y0
x  [a, b]
问题：如何计算 y(b) 的近似值？
3
Euler 公式
 Euler 公式
两边求积分可得

b
a
dy
dx 
dx

b
a
f ( x , y) dx
b
y( b)  y( a)   f ( x, y) dx
a
b
y( b)  y( a)   f ( x, y) dx
a
 y( a)  ( b  a) f ( a, y( a))
 y0  (b  a) f (a, y0 )
左点矩形公式
Euler 公式
4
Euler 法
 Euler 法
分割区间，在每个小区间上使用 Euler 公式，逐步递推
x0 
x1
x2
x i 1
y
f ( x)
o a
xi 1 xi
xn1
xi
b
x
 xn
5
Euler 法
 Euler 法
y( x1 )  y0  ( x1  x0 ) f ( x0 , y0 )
y1
y( x2 )  y1  ( x2  x1 ) f ( x1 , y1 )
y2
y( xn )  yn1  ( xn  xn1 ) f ( xn1, yn1 )
yn
yk 1  yk  ( xk 1  xk ) f ( xk , yk )
k  0, 1, 2, ..., n  1
6
几何含义
7
举例
例：用 Euler 法求解蹦极模型
cd 2
 dv
 g
v

m
 dt

 v ( 0)  0
t  [0, T ]
参数取值：g = 9.81, m = 68.1, cd = 0.25
解：见 bungee.m
8
梯形法
 梯形法
计算定积分时采用梯形公式，即
y( xk 1 )  y( xk )  
xk 1
xk
f ( x, y) dx
梯形公式
xk 1  xk
 y( x k ) 
f ( xk , y( xk ))  f ( xk 1 , y( xk 1 )) 

2
yk 1
xk 1  xk
 yk 
f ( xk , yk )  f ( xk 1 , yk 1 ) 

2
k  0, 1, 2, ..., n  1
梯形法
缺点：公式右端含有 yk+1，求解较困难
9
改进的Euler法
 改进的 Euler 法
 先用 Euler 法计算出 yk+1 的近似值，然后代入后端项中
yk 1  yk  ( xk 1  xk ) f ( xk , yk )
yk 1
xk 1  xk
 yk 
f ( xk , yk )  f ( xk 1 , yk 1 ) 

2
k  0, 1, 2, ..., n  1
 第一步称为预估，第二步称为校正
改进的 Euler 法
 通常也写为：
1
yk 1  yk  hk  K1  K 2 
hk  xk 1  xk
2
其中： K1  f ( xk , yk ), K2  f ( xk  hk , yk  hk K1 )
10
举例
例：用 Euler 法和改进的 Euler 法解初值问题
2x
 dy
 y

y
 dx
 y(0)  1

x  [0, 1]
解：见 Euler.m, EulerM.m,
Example2.m
11
主要内容
 Euler 法与改进的 Euler 法
 算法分析：误差与收敛性
 Runge-Kutta 法
 一阶常微分方程组与高阶方程
 应用举例
 Matlab 相关函数
12
单步法
 单步法与多步法
 如果在计算 yk+1 时，只需用到 yk 的值，则称为单步法
 如果在计算 yk+1 时，需用到 yk , yk-1 , . . . , yk-r+1 的值（即
需要使用前 r 的点的值），则称为多步法
 解初值问题的等步长单步法的一般形式
yk1  yk  h    xk , yk , yk1, hk 
 如果  中含有 yk+1 ，则称算法为隐式的，否则为显式的
 显式单步法的一般形式
yk1  yk  h    xk , yk , hk 
称  为增量函数
13
局部截断误差
定义：设 y(x) 是初值问题的精确解，则称
Rk 1  y( xk 1 )  y( xk )  hk   xk , y( xk ), hk 
为显式单步法的局部截断误差。
 Rk+1 是局部的，因为这里假定 yk 是精确的，即 yk=y(xk)
例：采用等步长的 Euler 法和梯形公式Euler法的局部截断误差
Euler 法的局部截断误差： Rk 1
h2

y ''( xk )  O ( h3 )
2
梯形公式 Euler法的局部截断误差：
Rk 1
h3

y '''( xk )  O ( h4 )
12
14
相容性与阶
 等步长显式单步法的相容性与阶
定义：若增量函数 (x, y, h) 在 h=0 处连续，且满足
  x , y, 0   f ( x , y )
则称显式单步法是相容的。
定义：设 p 是使得下式成立的最大正整数
Rk  O ( h p 1 )
则称显式单步法是 p 阶的。
例：等步长的 Euler 法是 1 阶的，
等步长的梯形公式 Euler 法是 2 阶的。
15
收敛性
定义：若设 y(x) 是解析解，yk 是 xk 处的数值解，若
lim yk  y( xk )  0
h0
( k  1, 2, ..., n)
则称算法是 收敛的。
 如果不是等步长算法，则 h 是指小区间长度的最大值
定理：若增量函数 (x, y, h) 关于 x, y, h 均满足 Lipschitz 条
件，则显式单步法的收敛性与相容性等价。
16
收敛性
定理：设显式单步法
yk 1  yk  hk   xk , yk , hk 
是 p 阶的，且增量函数  关于 y 满足 Lipschitz 条件，又设初
值是精确的，即 y0=y(x0)，则算法整体误差为
yk  y( xk )  O( h p )
即算法是收敛的，且收敛阶是 p。
推论：设 f(x, y) 关于 y 满足 Lipschitz 条件，则
(1) Euler 法收敛
(2) 改进的 Euler 法收敛
17
主要内容
 Euler 法与改进的 Euler 法
 算法分析：误差与收敛性
 Runge-Kutta 法
 一阶常微分方程组与高阶方程
 应用举例
 Matlab 相关函数
18
Runge-Kutta法
y( xk 1 )  y( xk )  
xk 1
xk
f ( x, y) dx
k  1, 2, ..., n
 算法的精度取决于数值积分的精度

xk  1
xk
r
f ( x , y ) dx  h i f  xk  i h, y( xk  i h) 
i 1
 采用高精度的数值积分方法就可能构造出高精度的算法
 一般来说，积分点越多，精度可能越高（但不宜太多）
 由于 y( xk  i h) 无法获得，因此需要用近似方法计算
19
R-K法
 显式 R-K 法的一般格式
r
yk 1  yk  h  i K i
k  1, 2, ..., n
i 1
K1  f  xk , yk 
i 1


K i  f  xk  i h, yk  h  ij K j 
j 1


 这里
i  2, ..., r
 i , i , ij 是参数
 参数选取准则：使得公式具有尽可能高的精度
 采用的工具：Taylor 展开
20
R-K法
 一阶 R-K 法
 二阶 R-K 法
 三阶 R-K 法
课堂板书
21
经典R-K法
 经典 R-K 方法
yk 1
h
 yk   K1  2 K 2  2 K 3  K 4 
6
k  1, 2, ..., n
K1  f  xk , yk 
1
h


K 2  f  xk  h, yk  K1 
2
2


1
h


K 3  f  xk  h, yk  K 2 
2
2


K4  f  xk  h, yk  hK3 
设 f(x, y) 关于 y 满足 Lipschitz 条件，则该算法是收敛的，
且具有 4 阶收敛阶。
22
举例
例：用经典 R-K 方法解初值问题
2x
 dy
 y

y
 dx
 y(0)  1

x  [0, 1]
解：见 Example3.m，RK4.m
23
主要内容
 Euler 法与改进的 Euler 法
 算法分析：误差与收敛性
 Runge-Kutta 法
 一阶常微分方程组与高阶方程
 应用举例
 Matlab 相关函数
24
一阶方程组
 一阶常微分方程组
 dY
 F ( x,Y )

 dx

Y ( x0 )  Y0
 y1 
y 
其中 Y   2 
 
 
 ym 
x0  a  x  b
 f1 
f 
2 

F
 
 
 fm 
 y10 
 0
y2 

Y0 
 
 0
 ym 
25
四阶 R-K 法
 四阶 R-K 法
Yk 1
h
 Yk   K1  2 K 2  2 K 3  K 4 
6
k  1, 2, ..., n
K1  F  xk ,Yk 
1
h


K 2  F  xk  h, Yk  K1 
2
2


1
h


K 3  F  xk  h, Yk  K 2 
2
2


K4  F  xk  h,Yk  hK3 
26
举例
例：m=2 时的四阶 R-K 方法
 y '  f ( x , y, z )

 z '  g ( x , y, z )
 y( x )  y , z ( x )  z

0
0
0
0
四阶 R-K 方法
h
yk 1  yk   K1  2 K 2  2 K 3  K 4 
6
h
zk 1  zk   L1  2 L2  2 L3  L4 
6
k  1, 2, ..., n
27
举例
K1  f  xk , yk , zk 
L1  g  xk , yk , zk 
1
h
h 

K 2  f  xk  h, yk  K1 , zk  L1 
2
2
2 

1
h
h 

L2  g  xk  h, yk  K1 , zk  L1 
2
2
2 

1
h
h 

K 3  f  xk  h, yk  K 2 , zk  L2 
2
2
2 

1
h
h 

L3  g  xk  h, yk  K 2 , zk  L2 
2
2
2 

K4  f  xk  h, yk  hK3 , zk  hL3 
L4  g  xk  h, yk  hK3 , zk  hL3 
28
高阶方程
 处理方法：化为一阶方程组
 dm y
( m 1)

f
(
x
,
y
,
y
',
...,
y
)

m
 dx
 y( x )  y , y '( x )  y ', ... , y( m1) ( x )  y( m1)

0
0
0
0
0
0
引入变量： z  y, z  y ', ..., z  y( m1)
1
2
m
 z1 '  z2

z2 '  z3



 z '  f ( x , z , z , ..., z )
1
2
m
 m
( m 1)

z
(
x
)

y
,
z
(
x
)

y
',
...
,
z
(
x
)

y
 1 0
0
2
0
0
m
0
0
29
举例：Van der Pol
例：求解 Van der Pol 方程
Balthasar van der Pol 是荷兰的一位著名电
子工程师。在20世纪20年代至30年代之间，
他开创了现代实验动力学。1927年，为了
描述在电子电路中三极管的震荡效应，第
一次推导出了著名的van der Pol方程
y   (1  y 2 ) y  y  0
Balthasar van der Pol
1889 -- 1959
van der Pol 震荡系统作为一种经典的自激震
荡系统，已经作为一种重要的数学模型并且
在更加复杂的动力系统的建模中广泛使用。
30
举例：Van der Pol
例：求解 Van der Pol 方程
3

d y
1  dy 
dy 
   
 y0

2
 3  dt 

dt
dt


2
初值：y(0)  1, y '(0)  1 ，求解区间 [0,100]
转化为一阶方程组：
 z1 '  z2


1 3

z
'


z


z

z
1
2
3 2
 2



以 u=5为例，编程求解：vandelpol.m, vandelpol_main.m
31