슬라이드 1

Size: px

Start display at page:

Download "슬라이드 1"

棉伸门
5 years ago
Views:

1 08-09 年度第一学期计算方法 (B) 童伟华管理科研楼 05 室 E-mal: 中国科学技术大学数学科学学院 ttp://mat.ustc.edu.cn/

2 第八章常微分方程数值解

3 微分方程数值解在科学研究或工程领域中, 有许多数学模型都是通过微分方程来描述的, 求解微分方程是非常重要的关键的问题微分方程按自变量的个数可分为 : 常微分方程偏微分方程微分方程按定解条件可分为 : 初值问题边值问题 3

4 微分方程数值解研究微分方程解析解的学科 : 数学物理方程, 但是绝大部分的微分方程是没有解析解的数值求解微分方程没有统一的算法, 针对不同类型的微分方程, 需要设计特定的算法目前, 研究数值求解微分方程的方法是热门的课题, 正在迅速发展之中, 常见的方法 : 有限差分法有限元方法有限体积法边界元方法谱方法 4

5 微分方程数值解常微分方程的初值问题 : 为了使解存在唯一, 一般需要对函数件 ( 初值问题解的存在唯一性 ) 若函数即 dy = f d y( a) = (, y) y 0 上存在且唯一, [ a, b] a b, < y < [ ab, ] f( y, ) f( y, ) Ly y f( y, ) f( y, ) 加限制条在条带上连续, 且满足 Lpsctz 条件,, 则初值问题的解在区间 5

6 微分方程数值解常微分方程的解是一个函数, 但是, 计算机没有办法直接对函数进行运算常微分方程的数值解采用数值离散的方法, 即在一系列离散点列上, 求未知函数在这些点上函数值的近似基本步骤如下 : () 对区间进行分割 : ( 即对 I : a = 0 < < < m = b 函数定义域进行离散 ), 目标是求解 { y= y ( )} 的值 () 对微分方程进行离散, 建立关于 { y } 的方程, 一般要求满足 : 解的存在唯一性稳定性收敛性相容性 (3) 解关于 { } 方程, 求出 { } 的值 y y 6

7 微分方程数值解主要问题 : 如何对定义域进行离散? 如何对微分方程进行离散? 收敛性问题, 即步长充分小时, 所得到的数值解能否逼近问题的真解误差估计稳定性问题, 即舍入误差在以后各步的计算中, 是否会无限制扩大计算效率并行计算 7

8 Euler 公式对定义域向前差商公式作等距剖分, 即 [ ab, ] b a I : = a +, =, = 0,,, m m y ( + ) y ( ) = y'( ) + y''( ξ) y ( + ) y ( ) = f(, y ( )) + y''( ξ) y ( + ) = y ( ) + f (, y ( )) + y''( ξ) y = y + f (, y ) + 8

9 Euler 公式定义 : 在假设第步计算是精确的前提下, 考虑截断误 p 差 T+ = y ( + ) y+, 称 T + 为局部截断误差若 T = O ( + + ), 则称方法是 p阶相容的, 简称相容向前差商公式的局部截断误差 : T = y ( ) y = O = y ( ) + f (, y ( )) + y''( ξ) y+ = y( ) + f (, y( )) + y ''( ξ) y f (, y)) = y ( ) ''( ξ ) 9

10 Euler 公式整体截断误差和收敛性 : 考虑局部截断误差的积累和传播 e = y ( ) y = y( ) + f (, y( )) + y ''( ξ) y f (, y)) y ( ) y + f(, y ( )) f(, y) + T + e + L y( ) y + T + ( + L) e + T, T = ma T ( ) ( + L) ( + L) e + T + T = ( + L) e + (( + L) + ) T j ( ) ( ) ( + L) ( + L) e + T + ( + L) + T j ( ) = ( + L) e + ( + L) + ( + L) + T 3 0

11 Euler 公式整体截断误差和收敛性 : ( ) + L e + + L L + T + ( ) 0 ( ) ( ) ( + L) = ( + ) + ( + L) + + L e0 T + + ( + L) ( + L) e0 + T L + T ( + L) e0 + L 是阶方法 ( + ) L T ( b a) L T e e0 + e e0 + L L T ( b a) L 当 e0 = 0, 0时, T= O ( ) e e, 即 Euler 公式 L 是收敛的

12 Euler 公式稳定性 : 误差在以后各步的计算中不会无限制扩大下面考虑简单情况 : 仅初值有误差, 而其他计算步骤无误差设是初值有误差后的计算值, 则 { z } y+ = y + f (, y) z = z + f (, z ) + e y z e + f(, y ) f(, z ) e + L y z = e ( + L) e ( + L) e e C e + ( + ) L 可以看出, 向前差商公式关于初值是稳定的当初始误差充分小, 以后各步的误差也充分小

13 Euler 公式向后差商公式 y ( + ) y ( ) = y'( + ) + y''( ξ) y ( + ) y ( ) = f( +, y ( + )) + y''( ξ) y ( + ) = y ( ) + f ( +, y ( + )) + y''( ξ) y = y + f (, y ) 隐式格式, 需要迭代求解 3

14 Euler 公式 Pcard 迭代格式 : (0) y+ = y + f (, y), k = 0,,, ( k+ ) ( k) y+ = y + f ( +, y+ ), 记 ( y) = y + f (, y), 则当充分小时, φ + ' φ ( y) = f y(, y) L < 从而迭代收敛 4

15 Euler 公式中心差商公式 y ( + ) y ( ) = y'( ) + y''( ξ) 3! y ( + ) y ( ) = f(, y ( )) + y''( ξ) 3! y = y + f (, y ) + 多步格式, 二阶格式, 数值不稳定 5

16 Euler 公式基于数值积分的近似公式若取若取若取 dy + + = f (, y), [ a, b] y( + ) y( ) = f (, y) d y '( t) d d = y'( t) y'( ) = f(, y ( )) + +, 则有 y '( t) d ( ) y '( ) = f (, y( )) y = y + f (, y ) + + y'( t) y'( ) = f(, y ( )) + + +, 则有 y '( t) d ( ) y '( ) = f (, y( )) y = y + f (, y ) y'( td ) ( f(, y ( )) + f( +, y ( + )), 则有 y+ = y + ( f(, y) + f( +, y+ )) y+ = y + f (, y) 预估 - 校正公式 y+ = y + ( f(, y) + f( +, y+ )) 6

17 Euler 公式向前差商的 Euler 算法 7

18 Runge-Kutta 方法 Taylor 级数法 ( 高阶单步方法 ): 取取 k k+ ( k) ( k+ ) y ( + ) = y ( ) + y '( ) + + y ( ) + y ( ξ) k! ( k + )! y'( ) = f( y, ) y''( ) = f( y, ) + fy( y, ) y' y'''( ) = y+ = y + ( f, y) + ( f(, y) + fy(, y) f(, y) ) +! k = k =, 向前差商的 Euler 公式, 可得 y+ = y + f(, y) + f(, y) + fy(, y) f(, y) ( ) 8

19 Runge-Kutta 方法 Taylor 级数法 : 要大量计算复合函数的全导数 Runge-Kutta 方法 : 一个点上的导数值可以用它邻近的一些点上的函数值来近似表示 ( 数值微分的思想 )+ Taylor 级数法基本思想 : d y'( ) y ( ) f(, y( )) ck k k k j j k ( k ) + + = j j j k! j= j! d j= f(, y ) j k = f( + a, y+ bk) a j j jl l l= j j = = l= b jl 9

20 Runge-Kutta 方法基本做法 : 利用二元函数的 Taylor 展开形式 f( +, y+ l) = l f(, y) k 0 k! + = y 比较方程两边幂次相同的项系数, 得到关于的方程组取, 作 Taylor 展开 k = cf(, y ( )) + cf( + ay, ( ) + bf (, y ( ))) cf(, y ( )) + cf(, y ( )) + caf (, y ( )) + cbf y f y O (, ( )) y(, ( )) + ( ) ( ) 与 f(, y) + f(, y) + fy(, y) f(, y) 比较得 : c+ c = ca = / cb = / 4 个未知数,3 个方程, 有无穷多组解 k a, b, c j jl j 0

21 Runge-Kutta 方法修正的 Euler 法 ( 中点法 ): 取 c = 0, c =, a = b = / y+ y k = + k = f(, y) k = f ( + /, y + k/) c = /, c = /, a = b = y = + y + k ( + k)/ k = f(, y) k = f ( +, y + k) c = / 4, c = 3 / 4, a = b = / 3 y+ y k ( 3 k)/4 = + + k = f(, y) k = f ( + /3, y + k/3) 二阶 Runge-Kutta 公式 : 取 Henu 公式 : 取

22 Runge-Kutta 方法三阶 Kutta 公式 : y = y + k ( + 4 k + + k3) / 6 k = f(, y) k = f ( + /, y + k/) k3 = f ( +, y k+ k) 三阶 Henu 公式 : y = y + k ( k3) / 4 k = f(, y) k = f ( + /3, y + k/3) k3 = f ( + /3, y + k /3)

23 Runge-Kutta 方法四阶 Runge-Kutta 公式 : y = y + k ( + k + k k4)/6 k = f(, y) k = f ( + /, y + k/) k3 = f ( + /, y + k /) k4 = f ( +, y + k3) 四阶 Kutta 公式 : y = y + k ( + 3k + 3 k k4) / 8 k = f(, y) k = f ( + /3, y + k/3) k3 = f ( + /3, y k/3 + k) k4 = f ( +, y + k k + k3) 3

24 Runge-Kutta 方法步长的自适应 : 设置一个误差容限,() 假如超出误差容限, 则否决这一步并减小步长 ;() 假如符合误差容限, 则接受这一步并选取适合下一步的步长改变步长的策略 : 步长加倍或减半根据阶的信息选取适当的步长 e c p+ p T y + e 0.8* + = < T e y 嵌入 Runge-Kutta 对 : 一个阶和另一个必要的计算 p p + 阶, 共享 4

25 Runge-Kutta 方法嵌入 Runge-Kutta 4/5 对公式 k = f(, y ) k = f ( +, y + k) k3 = f ( +, y + k+ k) k4 = f ( +, y + k k + k3) k5 = f ( +, y + k 8 k + k3 k4) k6 = f ( +, y k+ k k3+ k4 k5) y+ = y + k+ k3+ k4 k z+ = y + k+ k3+ k4 k5 + k

26 Runge-Kutta 方法步长控制所需的误差估计为 : 8 97 e = z y = k k k + k + k 步长自适应策略 : e () 若误差测试 < T成功, 则用 z 取代, 程序进 y + y + 入下一步 ; T y p+ () 否则, 用 + = 0.8* 再尝试一遍 ( 如果重复 e 失败, 则步长减半直至成功 ) 在 Matlab 中, 就使用这种方法, 譬如 ode3,ode45 等命令 ( 分别使用嵌入 Bogack-Sampne /3 对嵌入 Dormand-Prnce 4/5 对 ) 6

27 Runge-Kutta 方法 Runge-Kutta 算法 7

28 线性多步法线性步方法的一般形式 : k k α jym+ j= β jfm+ j, j= 0 j= 0 需要给定个初始值 k y, y,, yk k 0 才能启动迭代基本思想 : 微分方程化为积分方程, 用数值积分近似 dy = f y d ya ( ) = y0 (, ), ( ) ( * = ) + * '( ) a b y y y tdt 8

29 线性多步法 * 取分别为, 有,, + p y = y + y d= y + a f y + T ( ) ( ) + '( ) ( ) (, ( )) p + p p j j j + j= 0 显式公式 : 用数值积分节点,,, q构造插值多项 + 式近似 y'( ), 在区间 [ p, + ] 上计算数值积分 y '( ) d p 隐式公式 : 用数值积分节点 +,,, + q构造插值多项 + 式近似 y'( ), 在区间 [ p, + ] 上计算数值积分 y '( ) d p p 控制积分区间, q控制插值节点 q 9

30 线性多步法 Adams 格式 : 取例 : 构造 p p = 0 = 0, q= 的显式格式 [ ] + ( + ) = ( ) + 0( ) '( ) + ( ) '( ) + ( ) + + a0 = l0( ) d = d = + + a = l( ) d = d = + = y( ) + a y'( ) + a y'( ) + T + = R( ) d = y y l y l y R d T y = y (3) + y ( η ) 5! f y f y 3 [ 3 (, ) (, )] 3 (3) ( )( ) d = y ( ξ ) 二阶显式 Adams 公式 30

31 线性多步法三阶显式 Adams 公式 : [ 3 (, ) 6 (, ) + 5 (, )] 3 4 (4) f y f y f y y+ = y +, T+ = y ( ξ ) 8 四阶显式 Adams 公式 : y + = y (5) T + = y ( ξ ) [ + ] f y f y f y f y 55 (, ) 59 (, ) 37 (, ) 9 ( 3, 3), 三阶隐式 Adams 公式 : 四阶隐式 Adams 公式 : 4 [ 5 (, ) + 8 (, ) (, )] 4 (4) f + y+ f y f y y+ = y +, T+ = y ( ξ ) 4 y + = y (5) T + = y ( ξ ) [ + + ] f y f y f y f y 9 ( +, + ) 9(, ) 5 (, ) (, ), 4 3

32 线性多步法例 : 构造 p =, q= 的显式格式 [ ] + ( + ) ( ) 0( ) '( ) ( ) '( ) ( ) '( ) ( ) a0 l0( ) d d y = y + l y + l y + l y + R d a = y( ) + a y'( ) + a y'( ) + a y'( ) + T ( )( ) 7 = = = ( )( ) 3 = ( )( ) l ( ) d = d = ( )( ) ( )( ) = = = ( )( ) a l( ) d d (4) y ( η ) = = 3! 7 f(, y) f(, y ) + f(, y ) y+ = y T+ R( ) d ( )( )( ) d [ ] 3

33 线性多步法例 : 构造 p =, q= 的隐式格式 [ ] + ( + ) ( ) 0( ) '( + ) ( ) '( ) ( ) '( ) ( ) a0 l0( ) d d y = y + l y + l y + l y + R d a = y( ) + a y'( ) + a y'( ) + a y'( ) + T ( )( ) 3 = = = ( )( ) 4 = + + ( + )( ) l ( ) d = d = 0 ( )( ) ( )( ) 9 = = = ( )( ) a l( ) d d + (4) y ( η ) = = 3! [ 3 f( +, y+ ) + 9 f(, y ) ] y+ = y T+ R( ) d ( + )( )( ) d 33

34 线性多步法为避免迭代, 可用预估 - 校正公式 [ 7 f(, y ) f(, y ) + f(, y )] y = y f(, y ) + 9 f(, y ) + + y+ = y

35 常微分方程组的数值解法设 m个一阶方程组成的常微分方程初值问题 dy = f( y,,, ym) d dym = fm( y,,, ym), a b d y( a) = η ym( a) = ηm 写成向量形式 y( ) f(, y( ),, ym( )) dy (, ( )) y( ) f( y, ( ),, y( )), ( ), (, ( )) m = F Y d Y F Y = = Y( a) = η y ( ) f (, y ( ),, y ( )) m m m 35

36 常微分方程组的数值解法各种方法都可以直接运用过来例 : 两个一阶微分方程组成的微分方程组 dy = f( yz,, ) d dz = gyz (,, ), a b d ya ( ) = y0 za ( ) = z0 Euler 公式 y z n+ n+ = = y z n n + f ( + g( n n,, y y n n, z, z n n ) ) 36

37 常微分方程组的数值解法 Runge-Kutta 公式 yn+ yn = + K + K + K 3+ K 4 zn+ zn 6 ( ) f( n, yn, zn) K = g ( n, yn, z) n () () f( n +, yn + K, zn + K ) K = () () g ( n +, yn + K, zn + K ) () () f( n +, yn + K, zn + K ) K 3 = () () g ( n +, yn + K, zn + K ) K f ( +, y + K, z + K ) = (,, ) () () n n 3 n 3 4 () () g n + yn + K3 zn + K3 37

38 常微分方程组的数值解法高阶微分方程 : 引入辅助变量, 可化为一阶微分方程组设 m阶常微分方程 m d y = f yy m d ya ( ) = η y'( a) η = y ( m ) ( a) 引入辅助变量 y = dy d dy d y = m y = y m ( m ) (,, ',, y ) = η m, a b dy = y d dym = f( y,,, ym), a b d y( a) = η ym( a) = ηm 38

39 常微分方程的稳定性定义 : 若一个离散变量的 k 步方法, 在 Lpsctz 条件下, 存在常数 C 和, 使得以中任意值为步长 0 > 0 (0, 0 ], 通过任意两组初值 { } k u m m = 0和 { } k v m m = 0得到两组离散值 { u } N m = 和 { } N T t v m m= k( N 0 = ), 都成立则称该方法是稳定的 ma u v C ma u v, m m m m k m N 0 m k 数值方法的稳定性 : 初值产生误差的传播 ( 或者离散解关于初值的连续依赖性 ) 数值方法的稳定性不仅于数值方法有关, 而且与微分方程本身有关如果微分方程本身是不稳定的, 那就没理由要求数值方法稳定因此, 数值方法的稳定性概念是建立在微分方程稳定的基础上的 39 m k

40 常微分方程的稳定性为了比较不同方法的性质, 取典型的微分方程 dy = λy, Re( λ) < 0 d ya ( ) = y0 将一般的差分方程 ( 数值格式 ): k j= 0 a j y n+ j = j= 0 n+ j 应用到典型的微分方程, 可得 : k 对于给定的初始误差 e 的形式 k b j f ( j n+ j j n+ j j= 0 j= 0 k k, y n+ j ay = λ by, Re( λ) < 0 j n+ j j n+ j j= 0 j= 0 k ) 0, e,, e k ae = λ be, Re( λ) < 0, 误差方程具有一样 40

41 常微分方程的稳定性定义 : 差分方程称为绝对稳定的, 若差分方程作用到典型的微分方程 dy d = λy, Re( λ) < 0, 对任意的初值, 总存在左半复平面上的一个区域, 当在这个区域时, 差分方程的解趋于 0 ( ), 这个区域称为稳定区域例 : 向前 Euler 公式的稳定性 y = y + λy e = e + λe n+ n n n+ n n e e n+ = + λ < n Img Re 4

42 常微分方程的稳定性例 : 向后 Euler 公式的稳定性 y = y + λ y y = y λ e+ = e λ e e λ n+ = < n 定义 : 当绝对稳定区域是左半平面时, 则称该数值方法是无条件绝对稳定的一般来说, 隐式格式的绝对稳定性比同阶的显式法的好 Img 0 Re 4

43 常微分方程的稳定性例 :3 阶 Runge-Kutta 公式的稳定性 y K K K n+ 3 = y = λy n = λy = λy n n n + [ K 6 ( + λ) + 4K [ + λ + ( λ) ] + K 3 yn+ = yn + λ+ ( λ) + ( λ) 6 e ( ) ( ) e 6 n+ 3 = + λ+ λ + λ < n 3 ] 可以证明 :Runge-Kutta 方法和隐式 Adams 方法都是绝对稳定的 k = - k = k = 4 k = 3 Img Re

常微分方程的数值解法 - Numerical solution of ordinary differential equation

常微分方程的数值解法 - Numerical solution of ordinary differential equation 常微分方程的数值解法 Numerical solution of ordinary differential equation 张晓平 2018 年 12 月 17 日武汉大学数学与统计学院 Table of contents 1. 一般概念 2. 欧拉方法 3. 龙格 - 库塔方法 (Runge-Kutta method) 1 一般概念一般概念 1. 常微分方程的求解问题在实践中经常遇到, 但我们只知道一些特殊类型的常微分方程的解析解