第二讲 常微分方程数值求解

第二讲

常微分方程数值求解

—— Euler

——

主要内容

 Euler 法与改进的 Euler 法

 算法分析：误差与收敛性

 Runge-Kutta 法

 一阶常微分方程组与高阶方程

 应用举例

 Matlab 相关函数

初值问题

 考虑一维经典初值问题

问题：如何计算 y(b) ^{的近似值？}

0

( , ) ( )

dy f x y dx

y a y

 



 

 x [ , ]a b

Euler

两边求积分可得

左点矩形公式 Euler 公式

( , )

b b

a a

dy dx f x y dx dx 

 

( ) ( ) ^b ( , )

y b  y a 



( ) ( ) ^b ( , )

y b  y a 



( ) ( ) ( , ( ))

y a b a f a y a

  

0 ( ) ( , 0)

y b a f a y

  

Euler

分割区间，在每个小区间上使用 Euler 公式，逐步递推

( ) f x

a x x b _x

y

o

x0  x₁ x₂

 

 

 

Euler

 Euler 法

1 0 1 0 0 0 1

( ) ( ) ( , ) 

y x  y  x  x f x y y

2 1 2 1 1 1 2

( ) ( ) ( , ) 

y x  y  x  x f x y y

1 1 1 1

( n) n ( n n ) ( n , n )  n

y x



1 ( 1 ) ( , )

k k k k k k

y _  y  x _  x f x y

0 1 2, , , ..., 1

k  n 

几何含义

举例

例：用 Euler 法解初值问题

解： Matlab 代码见 Euler.m, Example1.m

真解

2

0 1 0 1

[ , ] ( )

dy x

dx y y x

y

  

 

 



2 1

y  x 

梯形法

计算定积分时采用梯形公式，即

梯形公式

缺点：公式右端含有 y_k+1 ，求解较困难

梯形法

1

( 1) ( ) ^k ( , )

k

x

k k x

y x _  y x 



 

1

1 1

( _k ) x^k 2 x^k ( _k , ( _k )) ( _k , ( _k )) y x ^  f x y x f x _ y x _

  

 

1

1 1 1

2 ( , ) ( , )

k k

k k k k k k

x x

y _ y ^  f x y f x _ y _

  

0 1 2, , , ..., 1

k  n 

改进的 Euler 法

 先用 Euler 法计算出 y_k+1 的近似值，然后代入后端项中

改进的 Euler 法

 第一步称为预估，第二步称为校正

 通常也写为：

其中：

 

1

1 1 1

2 ( , ) ( ,  )

k k

k k k k k k

x x

y _ y ^  f x y f x _ y _

  

1 ( 1 ) ( , )

_{k k k k k k}

y _  y  x _  x f x y

0 1 2, , , ..., 1

k  n 

1 ( _k, _k ), K  f x y

 

1 1 2

1

k k 2 k

y _  y  h K  K

2 ( _{k k} , _{k k} 1)

K  f x  h y  h K

1

k k k

h  x _  x

举例

例：用 Euler 法和改进的 Euler 法解初值问题

解： Matlab 代码见

EulerM.m, Example2.m 真解

[0,1]

12 (0) 2

x y y

dy y

dx x

 

 

 

12 2 2^x 6

y  x  e 

主要内容

 Euler 法与改进的 Euler 法

 算法分析：误差与收敛性

 Runge-Kutta 法

 一阶常微分方程组与高阶方程

 应用举例

 Matlab 相关函数

单步法

 单步法与多步法

 如果在计算 y_k+1 时，只需用到 y_k 的值，则称为单步法

 如果在计算 y_k+1 时，需用到 y_k , y_k-1 , . . . , y_k-r+1 的值

（即需要使用前 r 的点的值），则称为多步法

 解初值问题的等步长单步法的一般形式

 如果  中含有 y_k+1 ，则称算法为隐式的，否则为显式

的 显式单步法的一般形式

称  为增量函数

 

1  , , 1,

k k k k k k

y _  y  h x y y _ h

 

1  , ,

k k k k k

y _  y  h x y h

局部截断误差

定义：设 y(x) 是初值问题的精确解，则称

为显式单步法的局部截断误差。

 R_k+1 是局部的，因为这里假定 y_k 是精确的，即 y_k=y(x_k)

例：采用等步长的 Euler 法和梯形公式 Euler 法的局部截断误 差

Euler 法的局部截断误差：

梯形公式 Euler 法的局部截断误差：

 

1 ( 1) ( )  , ( ),

k k k k k k k

R _  y x _  y x  h x y x h

2 3

1 2 ''( ) ( )

k k

R _  h y x  O h

3 4

1 12 '''( ) ( )

k k

R _   h y x  O h

相容性与阶

定义：^{若增量函数}  ^{(x, y, h)} 在 h=0 处连续，且满足 则称显式单步法是相容的。

 等步长显式单步法的相容性与阶

定义：设 p 是使得下式成立的最大正整数

则称显式单步法是 p 阶的。

例：等步长的 Euler 法是 1 阶的，

等步长的梯形公式 Euler 法是 2 阶的。





 x y, ,  f x y( , )

( ^p 1) Rk  O h ^

收敛性

定义：若设 y(x) 是解析解， y_k 是 x_k 处的数值解，若 则称算法是 收敛的。

定理：^{若增量函数}  ^{(x, y, h)} 关于 x, y, h 均满足 Lipschitz 条件，则显式单步法的收敛性与相容性等价。

 如果不是等步长算法，则 h 是指小区间长度的最大值

0 0

lim _k ( _k )

h y y x

   (k  1 2, , ..., )n

收敛性

定理：^{设显式单步法}

是 p 阶的，且增量函数  关于 y 满足 Lipschitz 条件，又设 初值是精确的，即 y₀=y(x₀) ，则算法整体误差为

即算法是收敛的，且收敛阶是 p 。

推论：设 f(x, y) 关于 y 满足 Lipschitz 条件，则 (1) Euler 法收敛

(2) 改进的 Euler 法收敛

 

1  , ,

k k k k k k

y _  y  h x y h

( ) ( ^p)

k k

y  y x  O h

主要内容

 Euler 法与改进的 Euler 法

 算法分析：误差与收敛性

 Runge-Kutta 法

 一阶常微分方程组与高阶方程

 应用举例

 Matlab 相关函数

Runge-Kutta

 算法的精度取决于数值积分的精度

 一般来说，积分点越多，精度可能越高（但不宜太多）

 由于 无法获得，因此需要用近似方法计 算

 采用高精度的数值积分方法就可能构造出高精度的算法

1

1 d

( ) ( ) ^k ( , )

k

x

k k x

y x _  y x 



1 2, , ...,

k  n

 

1

1

( , ) d , ( )

k k

x r

i k i k i

x i

f x y x h  f x  h y x  h











( _{k i} )

y x   h

R-K

 显式 R-K 法的一般格式

 这里 是参数

 参数选取准则：使得公式具有尽可能高的精度

 工具： Taylor 展开

1, 2, ...,

k  n

1

1 r

k k i i

i

y _ y h  K



 



 

1 _k , _k

K  f x y

1

1

, ⁱ

i k i k ij j

j

K f x  h y h ^  K



 

    





, ,

i i ij

  

r = 1

课堂板书

利用 Taylor 展开，可得局部截断误差：

利当 时，达到最高阶数 =1 这事实上就是 Euler 法

 

1 1 1 1 ,

k k k k k

y _  y  h K  y  hf x y

 

1 ( 1) ( ) 1 ,

k k k k k

R _  y x _  y x  hf x y

2 3

1

'( ) 1 ''( ) ( ) '( )

k 2 k k

y x h y x h O h  y x h

   

2 3

1

(1 ) '( ) 1 ''( ) ( )

k 2 k

y x h y x h O h



   

1 1

 

r = 2

课堂板书

利用 Taylor 展开，可得局部截断误差：

令

 

1 1 1 2 2

k k

y _  y  h  K   K

 

1 1 1 2 '( )

k k

R _    y x h

' 2

2 2

1 ( , ( ))

2   f x y x^{x k k} h

 

  

 

1 _k, _k ,

K  f x y K2  f x





' 2 3

2 21

1 ( , ( )) ( , ( )) ( ) 2   f x y x^{y k k} f x y x^{k k} h O h

 

   

1 2

1    0, 1 _{2 2}

2    0, 1 _{2 21} 2    0

r = 2

解不惟一，设

公式一：令 a = 1/2

这事实上就是改进的 Euler 法

1 1 a,

   ₂ 1 2a ,

  ₂₁ 1

  2a

2 a 0

  

 

1 1 2

k k 2

y _  y  h K  K

 

1 _k, _k ,

K  f x y K2  f x





r = 2

公式二：令 a = 1

中点公式 Euler 法

注： r = 2 时，显式 R-K 法的阶至多只能达到二

阶

1 2

k k

y _  y  hK

 

1 _k, _k ,

K  f x y ₂ 1 1 ₁

2 , 2

k k

K  f x  h y  hK 

r = 3

课堂板书

一个常用的三阶 R-K 法

 

1 1 4 2 3

k k 6

y _  y  h K  K  K

 

1 _k, _k ,

K  f x y

2 1

1 1

2 , 2

k k

K  f x  h y  hK 

 

3 _k , _k 1 2 2

K  f x  h y  hK  hK

四阶经典 R-K 法

设 f(x, y) 关于 y 满足 Lipschitz 条件，则该算法是收敛 的，且具有 4 阶收敛阶。

1, 2, ...,

k  n

 

1 1 2 2 2 3 4

k k 6

y _  y  h K  K  K  K

 

1 _k , _k

K  f x y

2 1

1 ,

2 2

k k

K  f x  h y  h K 

3 2

1 ,

2 2

k k

K  f x  h y  h K 

 

4 _k , _k 3

K  f x  h y  hK

举例

例：用经典 R-K 方法解初值问题

解： Example3.m ， RK4.m

[0,1]

12 (0) 2

x y y

dy y

dx x

 

 

 

主要内容

 Euler 法与改进的 Euler 法

 算法分析：误差与收敛性

 Runge-Kutta 法

 一阶常微分方程组与高阶方程

 应用举例

 Matlab 相关函数

一阶常微分方程组

其中

0 0

( , )

( )

dY F x Y dx

Y x Y

 



 



x0  a x  b

1 2

m

y Y y

y

 

 

 

  

 

 



0 1

0 2 0

0 m

y Y y

y

 

 

 

  

 

 



1 2

m

f F f

f

 

 

 

  

 

 



四阶 R-K 法

1, 2, ...,

k  n

 

1 _k , _k

K  F x Y

2 1

1 ,

2 2

k k

K  F x  h Y  h K 

3 2

1 ,

2 2

k k

K  F x  h Y  h K 

 

4 _k , _k 3

K  F x  h Y  hK

 

1 1 2 2 2 3 4

k k 6

Y _  Y  h K  K  K  K

举例

例： m=2 时的四阶 R-K 方法

四阶 R-K 方法

1, 2, ...,

k  n

 

1 1 2 2 2 3 4

k k 6

y _  y  h K  K  K  K

0 0 0 0

,

' ( , , )

' ( , , )

( ) ( )

y f x y z z g x y z

y x y z x z

 

 

  



 

1 1 2 2 2 3 4

k k 6

z _  z  h L  L  L  L

举例

 

1 _k, _k, _k

K  f x y z

2 1 1

1 , ,

2 2 2

k k k

h h

K  f x  h y  K z  L 

 

3 2 2

1 , ,

2 2 2

k k k

h h

K  f x  h y  K z  L 

 

4 _k , _k 3, _k 3

K  f x  h y  hK z  hL

 

1 _k, _k, _k

L  g x y z

2 1 1

1 , ,

2 2 2

k k k

h h

L  g x  h y  K z  L 

3 2 2

1 , ,

2 2 2

k k k

h h

L  g x  h y  K z  L 

 

4 _k , _k 3, _k 3

L  g x  h y  hK z  hL

高阶方程

 处理方法：化为一阶方程组

引入变量：

( 1)

( 1) ( 1)

0 0 0 0 0 0

( , , ', ..., )

( ) , '( ) ', ... , ( )

m

m m

m m

d y f x y y y dx

y x y y x y y x y



 

 



   



( 1) 1 , 2 ', ..., _m ^m z  y z  y z  y ^

1 2

2 3

1 2

( 1)

' '

' ( , , , ..., )

( ) , ( ) ', ... , ( )

m m

m

z z z z

z f x z z z

z x y z x y z x y ^

 

 



 

   





举例： Van der Pol

例：求解 Van der Pol 方程

Balthasar van der Pol 是荷兰的一位著名 电子工程师。在 20 世纪 20 年代至 30 年代 之间，他开创了现代实验动力学。 1927 年，为了描述在电子电路中三极管的震荡 效应，第一次推导出了著名的 van der Pol 方程

Balthasar van der Pol 1889 -- 1959

van der Pol 震荡系统作为一种经典的自激

震荡系统，已经作为一种重要的数学模型并 且在更加复杂的动力系统的建模中广泛使 用。

(1 2) 0

y    y y  y

 

举例： Van der Pol

例：求解 Van der Pol 方程

转化为一阶方程组：

初值： ，求解区间 [0,100]

以 u=5 为例，编程求解： vandelpol.m, vandelpol_main.m

2 3 2

1 0

3

d y dy dy

dt   ^ ^ dt ^  dt ^  y

1 2

3

2 1 2 2

' ' 1

3 z z

z z  z z

 

      

  



(0) 1, '(0) 1

y  y 