数值模拟导论——第十三讲 多步法收敛性 Jacob White 合作伙伴Deepak Ramaswamy, MichalRewienski, and Karen Veroy 概要 多步法的小时间步问题 局部切断误差 选择系数 不收敛法 稳定并连续则收敛 下一章讨论大时间步问题 两个时间刻度实例的绝对稳定性 震荡器 () () () , d xt f xt ut dt = 多步法 基本方程式 通用符号 非线性差分方程 k步法： () () 00 ?? , kk lj lj jjlj ax t f x utβ ?? ? == =? ∑∑ 离散点的解 时间离散化 多步法系数 () ( ) 11 1 ?? ? , ll l l xx tfxut ?? ? ?=? () ( ) 1 ?? ? , ll l l xx tf xut ? ?=? 多步法 基本方程式 通用算法 多步法方程： () () 00 ?? , kk lj lj jjlj ax t f x utβ ?? ? == =? ∑∑ 前欧拉近似法： ( ) ( ) ( ) ( )( ) 111 , ll ll xt xt tf xt ut ??? ≈? 前欧拉离散方程： 多步法系数： 01 01 1, 1, 1, 0, 1k α αββ= ==?== 后欧拉法离散方程： 多步法系数： 01 01 1, 1, 1, 1, 0k α αββ= ==?== 捕捉法离散方程： () () () () ( ) 11 1 ?? ? ? ,, , 2 ll l l ll t x x fxut fx ut ?? ? ? ?= 多步法系数：01 0 1 11 1, 1, 1, , 22 k αα β β= ==?= = 多步法 基本方程式 定义及考察 多步法方程： () () 00 ?? , kk lj lj jjlj ax t f x utβ ?? ? == =? ∑∑ 1）若 0β ≠ 则多步法是绝对的 2）k步法在 ''xs f s和 之前用到k 3）需要范数化，都有 0 1α＝ 4）k步法有2k＋1个自由系数 需要高的精度时如何选择合适的系数？ () () 00 ?? kk lj lj jj d xt Axt ax t x dt β ? ? == =?=? ∑∑ 多步法离散化 () () 1 00 ?? kk lj lj jj Eyt xt ay t E AEyβ ? ?? == =? =? ∑∑ 令 多步法 简化问题分析 标量ODE： () () () 0 , 0 C d vt vt v v dt λλ= =∈ 为什么会出现这样的简单测试问题？ z非线性分析有很多微妙的现象没有揭示 z对于多步法标量等同于矢量 1 00 kk lj lj jj n ay t y λ β λ ? ? == ?? ?? ?=? ?? ?? ∑∑ null离散方程 多步法 简单问题分析 标量ODE：() () () 0 , 0 C d vt vt v v dt λλ= =∈ 多步法标量公式： 00 ?? kk l j l j jj av t vβλ ? ? == =? ∑∑ 对于所有 Cλ∈ 都必须考虑λ的值 多步法 收敛分析 收敛定义 定义：对于给定的任何初始条件，用多步法求解[0,T]上的初值问题是收敛的。 ( ) 0, ? 0 max 0 l T l t tvvt ?? ∈ ?? ? ?? ? →??→当时 多步法 收敛分析 P阶收敛 定义：对于任一给定的λ和任何初始条件，用多步法求解[0,T]上的初 值问题是P阶收敛的。 对于所有 0 tt? ?小于给定的的 () () 0, ? max p l T l t vvlt Ct ?? ∈ ?? ? ?? ? ?≤? 前欧拉法和后欧拉法都是1阶收敛的，梯形法则是2阶收敛的 多步法 收敛分析 实例反应方程 前欧拉法和后欧拉法，误差∞Δt 后欧拉法，误差∞（Δt）2 多步法 收敛分析 收敛的两个条件 1）局部条件：“一步法”误差较小（连续性） 可以用泰勒奇数证明 2）全局条件：单步误差不会快速增加（稳定性） 在这种情况下，所有一步法（K＝1）都是稳定的 多步法（K>1）需要详细分析 多步法 收敛分析 全局误差等式 00 ?? 0 kk lj lj jj av t vβλ ?? == ??= ∑∑多步法公式： 精确解基本上和多步法一致： () () 00 ? kk l jlj j lj jj d av t t vt e dt β ?? == ??= ∑∑ 局部切断误差（LTE） 全局误差： ( ) ? ll l E vt v≡? 差分方程将局部切断误差和全局误差联系起来 ()( ) ( ) 1 00 11 ll lkl kk atEatE atEeλβ λβ λβ ?? ? ?+?? ++?? =null 1 ?? ? 0 ll l vvtvλ + ? ?? = 前欧拉法 收敛分析 前欧拉法的连续性 前欧拉法的定义 分解精确解v（t）并展开得 () () ()() () 2 2 2 1 2 l e dv l t t d v vl t vlt t dt dt τ?? +?? ??? = nullnullnullnullnullnullnull 其中 l e 为局部切断误差，其范围： () 2 l eCt≤? ，当C＝0.5时 d vv dt λ= ( ) ,1ltl tτ∈ ?+?? ? ? ? [] ( ) 2 0, 2 max T dv dt τ τ ∈ 前欧拉法 收敛分析 全局误差方程 前欧拉法定义 1 ? l v + ? l v= + ? l tvλ? 利用局部切断误差定义 ( ) ( ) ( )() 1 l vl t vlt t lt eλ+?= ??? ?+ 上式左右相减得全局误差方程 1 () lll E ItEeλ + = +? + 利用 l e的数量值及其范围得 ( ) () 2 1 1 lll l EItEe tECtλλ + ≤+? + ≤+? + ? 前欧拉法 收敛分析 差分方程一个非常有利的跳跃 一个λ跳跃得差分方程的解 若() 10 1, 0, 0 ll uubuεε + ≤+ + = > 那么 l l e ub ε ε ≤ 为了便于证明，将 l u 写成幂级数及求和形式 () () () 1 0 11 1 11 l j l l j ubb ε ε ε ? = ?+ ≤+ = ?+ ∑ 一步法 收敛分析 差分方程一个非常有利的跳跃 最后，由于 () () 11 l l ee ε ε εε+ ≤?+ ≤ () () () 11 1 1 11 ll l l e ubbb ε εε ε ?+ + ? ≤=≤ ?+ 一步法 收敛分析 回顾前欧拉法的收敛分析 应用λ并消去Δt得 null () 2 2 1() lt ll l b e EtECtCt t λ ε λ λ ? + ?? ?? ≤+? +?≤ ? ? ?? nullnullnull 最后，由于 ,lt T?≤ [] 0, max T l lT C E et λ λ ∈ ≤ ?? 前欧拉法 收敛分析 考查前欧拉法分析 [] 0, max T l lT C Ee t λ λ ∈ ≤ ?? z前欧拉法时一阶收敛 z跳跃随时间间隔程指数增长 zC与精确解得二次倒数相关 z跳跃随时间增长程指数增加 前欧拉法 收敛分析 不稳定反应的精确解和前欧拉法图 前欧拉法误差随时间延长而增大 前欧拉法 收敛分析 求解反应方程的前欧拉法误差 正如限制所约束的那样，误差随时间增长程指数增加 前欧拉法 收敛分析 电路问题中的精确解和前欧拉法解 图 前欧拉法误差并不总是随时间延长而增大 前欧拉法 收敛分析 求解电路问题方程的前欧拉法的误 差 误差并不是总是随时间程指数增长！ 跳跃是保守的 多步法 使得局部切断误差足够小 精确约束 局部切断误差： () () 00 kk l jlj j lj jj d avt t vt e dt β ?? == ??= ∑∑ () () 1p p d vt t vt pt dt ? =? = 若 注：不能根据 () () d vt vt dt λ= 式推导 () () () () 1 00 kj kj kk pp jj d vt vt dt akj ttpkj teβ ? ? ? == ???? ?? = ∑∑ nullnullnullnullnullnullnullnullnullnullnullnullnullnull 局部切断误差 多步法 使得局部切断误差足够小 精确约束 () () () () () 1 00 1 00 kk pp jj kk pp p k jj akjt t pkjt takj pkj e β β ? == ? == ???? ?? = ?? ????= ?? ?? ∑∑ ∑∑ 对于所有 0 ,p p≤ 若() () 1 00 0 kk pp jj ak j pk jβ ? == ?? ? ??= ?? ?? ∑∑ () () () 0 1 00 kk p l jlj j lj jj d avt t vt e C t dt β + ?? == ?? = = ? ∑∑ 那么 对于( ) p vt t=若 () ( ) ( ) 1 00 0 kk pp p jj takj pkjβ ? == ?? ? ?? ? = ??∑∑ 那么 0 k e = 对于所有光滑v(t)在t处都有一个局部精确泰勒级数 多步法 使得局部切断误差足够小 精度约束k=2实例 精度拘束： () () 1 00 0 kk pp jj ak j pk jβ ? == ?? ? ??= ?? ?? ∑∑ 对于k＝2，产生系数的5x6维系统方程式 0 1 3 4 p p p p = = = = 0 1 2 0 1 2 1 1 1 0 0 0 0 2 1 0 -1 -1 -1 0 4 1 0 -4 -2 0 0 8 1 0 -12 -3 0 0 16 1 0 -32 -4 0 0 a a a β β β ?? ???? ?? ???? ?? ?? = ?? ???? ?? ?? ?? ???? ?? ?? 注：总存在0 i a = ∑ 多步法 使得局部切断误差足够小 精度约束k=2实例 0 1 3 4 p p p p = = = = 0 1 2 0 1 2 1 1 1 0 0 0 0 2 1 0 -1 -1 -1 0 4 1 0 -4 -2 0 0 8 1 0 -12 -3 0 0 16 1 0 -32 -4 0 0 a a a β β β ?? ? ??? ?? ? ??? ?? ? ?? ? ???= ?? ? ?? ? ??? ?? ? ?? ? ??? ?? ?? 对于k＝2精度约束 前欧拉法 01 2 0 1 2 1, 1, 0, 0, 1, 0aa aβ ββ==?= === 前欧拉法满足p=0和p=1但是不满足p=2 2 ()Ct??局部切断误差＝ 梯形捕捉法 01 2 0 1 2 1, 1, 0, 0, 0, 0aa aβ ββ==?= = = = 梯形捕捉法满足p＝0，p=1，p=2但是不满足p＝3 3 ()Ct??局部切断误差＝ 后欧拉法 01 2 0 1 2 1, 1, 0, 0, 0, 0aa aβ ββ==?= = = = 后欧拉法满足p=0和p=1但是不满足p=2 2 ()Ct??局部切断误差＝ 多步法 使得局部切断误差较小 精度约束k=2实例,生成方法 首先介绍一个规范化的问题，例如 0 1a = 1 2 0 1 2 1 1 0 0 0 1 1 0 -1 -1 -1 2 1 0 -4 -2 0 4 1 0 -12 -3 0 8 1 0 -32 -4 0 16 a a β β β ??? ? ??? ??? ??? ? ??? ??? ??? =? ??? ? ??? ??? ??? ? ? ????? 用两步法求解最小的局部切断误差 01 2 0 1 2 14 1 1, 0, 1, , , 333 aaa βββ== =?= = = 满足所有五个精度约束局部切断误差＝ 5 ()Ct? 用两步直接法求解最小的局部切断误差 01 2 0 1 2 1, 4, 5, 0, 4, 2aaa β ββ= ==?=== 只能满足四个精度局部切断误差＝ 4 ()Ct? 多步法 使得局部切断误差足够小 前欧拉法，后欧拉法，梯形捕捉法 和“最”直接法的局部切断误差图 前欧拉法 使得局部切断误差足够小 前欧拉法，后欧拉法，梯形捕捉法 和“最”直接法的全局误差图 前欧拉法 使得局部切断误差足够小 前欧拉法，后欧拉法，梯形捕捉法 和“最”直接法的全局误差图 前欧拉法 多步法的稳定性 差分方程 为什么“最”2步直接法不能收敛？ 多步法差分方程 全局误差 我们使得局部切断误差较小，全局误差怎么变得这样大？ 求解差分方程补充 考虑一般的k阶差分方程 1 01 ll lkl k ax ax ax u ?? + ++ =null 其中必须已知k个初始条件 01 01 , , , k k xxx x x x ?? = ==null 明显地，当方程由以下公式更新时 () 10 11 1 0 1 k k xaxaxu a ?+ =? + + ?null 大多数重要的差分方程可以求解 x通过式 0 l llj j xhu ? = = ∑ 和u联系起来 求解差分方程补充 若方程 1 01 0 kk k az az a ? +++=null 有不相等的根 12 ,,, k ? ??… 那么当 () 1 k l l j j j h γ ? = = ∑ 时， 1 l llj j xhu ? = = ∑ 为了理解h时怎么获得的，首先看一个简单的例子 假设 1lll xxu? ? =+且 () 010121212 1 0, , k l l jj j xxxuxxuuu x ??? γ? = = =+ =+=+ = ∑ 求解差分方程补充 三个重要的问题 若对于所有i，都有 1 i ? < ，那么 max lj j xC u≤ 其中C不依赖于l 若对于任何i，都有1 i ? > ， 那么存在一个有极限的 j u 使得 l x →∞ 若对于所有i，都有 1 i ? ≤ ，且当 1 i ? = 时， i ? 不等， 那么 max lj j xCl u≤ 多步法 多步法的稳定性 多步法差分方程 ( ) ( ) ( ) 1 00 11 ll lkl kk atEatE atEeλβ λβ λβ ?? ?? + ?? + + ?? =null 定义：对于任意 l e ，当且仅当 0, 0, 0 max max ll TT ll tt T tECe t ?? ?? ∈∈ ?? ?? ?? ?? ?? ?→ ≤ ? 时，多步法是稳定的。 定理：当且仅当方程 1 01 0 kk k az az a ? + ++=null 的根在数值上小于1或等于1但根不同时，多步法是稳定的。 差分方程 多步法 多步法的稳定性 稳定性定理“证明” 给定的多步法差分方程 ()( )() 1 00 11 ll lkl kk atEatE atEeλβ λβ λβ ?? ?? + ?? + + ?? =null 当方程 0 0 k kj j j az ? = = ∑ 的根为下列两者之一时 z数值上小于1 z数量上等于1，但是根不相同 那么由差分方程得 max ll l ECl e≤ 由此很容易证明稳定性 多步法 多步法的稳定性 稳定性定理“证明” 多步法 多步法的稳定性 最直接法 最直接2步法 “最”直接法非常不稳定 前欧拉法 多步法的稳定性 Dahlquist首次稳定性界限 对于一个稳定的，直接k步多步法，可以满足的精度约束的 最大值小于或等于k（注意由2k个系数）。对于非直接法，精 度约束在k为偶数时为k＋2，或k为奇数时精度约束为k＋1。 前欧拉法 收敛分析 收敛性、稳定性和连续性的条件 1）局部条件：一步误差较小（连续性） 对于 0 tt? <? ，精度约束达到 0 p（0 p 必须大于0） () 0 1 1 0, max p l T l t eCt + ?? ∈ ?? ? ?? ?≤? 2）全局条件：一步误差缓慢增加（稳定性） 方程的根在单元圆的内部或在圆上 2 0, 0, max max ll TT ll tt T EC e t ?? ?? ∈∈ ?? ?? ?? ?? ?? ?≤ ? 收敛结论： () 0 0, max p l T l t E CT t ?? ∈ ?? ? ?? ≤? 小结： 多步法的小时间步问题 局部切断误差和精度 差分方程的稳定性 稳定性＋连续性导出收敛 下一章要讨论的 两种时间刻度的绝对稳定性 振动 格朗－库拉方案