常微分方程初值問題數(shù)值解法課件

1、第9章 常微分方程初值問題數(shù)值解法 9.1 引 言 本章要著重考察的一階方程的初值問題 （1.1）（1.2） 只要函數(shù) 適當光滑譬如關(guān)于 滿足利普希茨(Lipschitz）條件 （1.3）理論上就可以保證初值問題（1.1），（1.2）的解 存在并且唯一. 1 所謂數(shù)值解法，就是尋求解 在一系列離散節(jié)點 上的近似值 . 相鄰兩個節(jié)點的間距 稱為步長. 如不特別說明，總是假定 為定數(shù)，這時節(jié)點為 . 初值問題（1.1），（1.2）的數(shù)值解法的基本特點是采取“步進式”，即求解過程順著節(jié)點排列的次序一步一步地向前推進. 描述這類算法，只要給出用已知信息 計算 的遞推公式. 首先，要對方程（1.1）離散

2、化，建立求數(shù)值解的遞2推公式. 一類是計算 時只用到前一點的值 ，稱為單步法. 另一類是用到 前面 點的值 ，稱為 步法. 其次，要研究公式的局部截斷誤差和階，數(shù)值解 與精確解 的誤差估計及收斂性，還有遞推公式的計算穩(wěn)定性等問題. 3 9.2 簡單的數(shù)值方法與基本概念 9.2.1 歐拉法與后退歐拉法 在 平面上，微分方程（1.1）的解 稱作它的積分曲線. 積分曲線上一點 的切線斜率等于函數(shù) 的值. 如果按函數(shù) 在 平面上建立一個方向場，那么，積分曲線上每一點的切線方向均與方向場在該點的方向相一致. 基于上述幾何解釋，從初始點 出發(fā)，先依方向場在該點的方向推進到 上一點 ，然后再從 依方向場的方

3、向推進到 上一點 ，循此前進做出4一條折線 （圖9-1). 一般地，設已做出該折線的頂點 ，過 依方向場的方向再推進到 ，顯然兩個頂點 的坐標有關(guān)系 圖9-15即 （2.1）這就是著名的歐拉（Euler）公式. 若初值 已知，則依公式（2.1）可逐步算出 例1 求解初值問題 （2.2）6 解 歐拉公式的具體形式為 取步長 ，計算結(jié)果見表9-1. 7 初值問題（2.2）有解 ，按這個解析式子算出的準確值 同近似值 一起列在表9-1中，兩者相比較可以看出歐拉方法的精度很差. 還可以通過幾何直觀來考察歐拉方法的精度. 假設 ，即頂點 落在積分曲線 上，那么，按歐拉方法做出的折線 便是 過點 的切線（

4、圖9-2）. 圖9-28 從圖形上看，這樣定出的頂點 顯著地偏離了原來的積分曲線，可見歐拉方法是相當粗糙的. 為了分析計算公式的精度，通常可用泰勒展開將 在 處展開，則有 在 的前提下， .于是可得歐拉法（2.1）的公式誤差 （2.3）稱為此方法的局部截斷誤差. 9如果對方程（1.1）從 到 積分，得 （2.4）右端積分用左矩形公式 近似，再以 代替 代替 也得到（2.1），局部截斷誤差也是（2.3）. 如果在（2.4）中右端積分用右矩形公式 近似，則得另一個公式 （2.5）稱為后退的歐拉法. 后退的歐拉公式與歐拉公式有著本質(zhì)的區(qū)別，后者是關(guān)于 的一個直接的計算公式，這類公式稱作是顯式的；10

5、然而公式（2.5）的右端含有未知的 ，它實際上是關(guān)于 的一個函數(shù)方程，這類公式稱作是隱式的. 隱式方程（2.5）通常用迭代法求解，而迭代過程的實質(zhì)是逐步顯示化. 設用歐拉公式 給出迭代初值 ，用它代入（2.5）式的右端，使之轉(zhuǎn)化為顯式，直接計算得 然后再用 代入（2.5）式，又有 11如此反復進行，得 （2.6）由于 對 滿足利普希茨條件（1.3）. 由（2.6）減（2.5）得 由此可知，只要 迭代法（2.6）就收斂到解 . 關(guān)于后退歐拉方法的公式誤差，從積分公式看到它與歐拉法是相似的. 12 9.2.2 梯形方法 為得到比歐拉法精度高的計算公式，在等式（2.4）右端積分中若用梯形求積公式近似

6、，并用 代替 代替 ，則得 （2.7）稱為梯形方法. 梯形方法是隱式單步法，可用迭代法求解. 同后退的歐拉方法一樣，仍用歐拉方法提供迭代初值，則梯形法的迭代公式為 13（2.8） 為了分析迭代過程的收斂性，將（2.7）式與（2.8）式相減，得 于是有 式中 為 關(guān)于 的利普希茨常數(shù). 14如果選取 充分小，使得 則當 時有 ，這說明迭代過程（2.8）是收斂的. 15 9.2.3 單步法的局部截斷誤差與階 初值問題(1.1)，(1.2)的單步法可用一般形式表示為 （2.9）其中多元函數(shù) 與 有關(guān)，當 含有 時，方法是隱式的，若不含 則為顯式方法，所以顯式單步法可表示為 （2.10） 稱為增量函數(shù)

7、，例如對歐拉法（2.1）有 它的局部截斷誤差已由（2.3）給出.16 對一般顯式單步法則可如下定義. 定義1 設 是初值問題(1.1)，(1.2)的準確解，稱 （2.11）為顯式單步法（2.10）的局部截斷誤差. 之所以稱為局部的，是假設在 前各步?jīng)]有誤差. 當 時，計算一步，則有 17 根據(jù)定義，顯然歐拉法的局部截斷誤差 即為（2.3）的結(jié)果. 這里 稱為局部截斷誤差主項. 顯然 .的公式誤差. 所以，局部截斷誤差可理解為用方法（2.10）計算一步的誤差，也即公式（2.10）中用準確解 代替數(shù)值解產(chǎn)生18 定義2 設 是初值問題(1.1)，(1.2)的準確解，若存在最大整數(shù) 使顯式單步法(2

8、.10)的局部截斷誤差滿足 （2.12）則稱方法（2.10）具有 階精度. 若將（2.12）展開式寫成 則 稱為局部截斷誤差主項. 以上定義對隱式單步法（2.9）也是適用的. 19 對后退歐拉法（2.5）其局部截斷誤差為 這里 ，是1階方法，局部截斷誤差主項為 . 20 對梯形法（2.7）有 所以梯形方法（2.7）是二階的，其局部誤差主項為 . 21 9.2.4 改進的歐拉公式 梯形方法雖然提高了精度，但其算法復雜，在應用迭代公式（2.9）進行實際計算時，每迭代一次，都要重新計算函數(shù) 的值，而迭代又要反復進行若干次，計算量很大，而且往往難以預測. 為了控制計算量，通常只迭代一兩次就轉(zhuǎn)入下一步的

9、計算，這就簡化了算法. 具體地，先用歐拉公式求得一個初步的近似值 , 稱之為預測值，預測值 的精度可能很差，再用梯形公式（2.7）將它校正一次，即按（2.8）式迭代一次得，這個結(jié)果稱校正值，而這樣建立的預測-校正系統(tǒng)通常稱為改進的歐拉公式： 22預測校正（2.13）或表為下列平均化形式 23 例2 用改進的歐拉方法求解初值問題（2.2）. 解 改進的歐拉公式為 仍取 ，計算結(jié)果見表9-2. 同例1中歐拉法的計算結(jié)果比較，改進歐拉法明顯改善了精度. 2425 9.3 龍格-庫塔方法 9.3.1 顯式龍格-庫塔法的一般形式 上節(jié)給出了顯式單步法的表達式其局部截斷誤差為對歐拉法 ，即方法為 階，若用

10、改進歐拉法，它可表為 （3.1）26此時增量函數(shù) （3.2）它比歐拉法的 ，增加了計算一個右函數(shù) 的值，可望 . 若要使得到的公式階數(shù) 更大， 就必須包含更多的 值. 從方程（1.1）等價的積分形式（2.4），即 （3.3）若要使公式階數(shù)提高，就必須使右端積分的數(shù)值求積公式27精度提高，必然要增加求積節(jié)點，為此可將（3.3）的右端用求積公式表示為點數(shù) 越多，精度越高，上式右端相當于增量函數(shù) ，為得到便于計算的顯式方法，可類似于改進歐拉法（3.1），（3.2），將公式表示為 （3.4）其中 （3.5）28這里 均為常數(shù). (3.4)和(3.5)稱為 級顯式龍格-庫塔（Runge-Kutta）法，

11、簡稱R-K方法. 當 時，就是歐拉法，此時方法的階為 . 當 時，改進歐拉法(3.1)，(3.2)也是其中的一種，下面將證明階 . 要使公式(3.4)，(3.5)具有更高的階 , 就要增加點數(shù) . 下面只就 推導R-K方法. 29 9.3.2 二階顯式R-K方法 對 的R-K方法，由（3.4），（3.5）可得到如下的計算公式 （3.6）這里 均為待定常數(shù)，希望適當選取這些系數(shù)，使公式階數(shù) 盡量高. 根據(jù)局部截斷誤差定義，（3.6）的局部截斷誤差為 （3.7）30這里 . 為得到 的階 ，要將上式各項在 處做泰勒展開，由于 是二元函數(shù)，故要用到二元泰勒展開，各項展開式為 其中 （3.8）31將以

12、上結(jié)果代入（3.7）則有 要使公式（3.6）具有 階，必須使 32（3.9）即 （3.9）的解是不唯一的. 令 ，則得 這樣得到的公式稱為二階R-K方法，如取 ，則這就是改進歐拉法（3.1）.33若取 ，則 . 得計算公式 （3.10）稱為中點公式，相當于數(shù)值積分的中矩形公式. （3.10）也可表示為 34 的R-K公式(3.6)的局部誤差不可能提高到 . 把 多展開一項，從（3.8）的 看到展開式中 的項是不能通過選擇參數(shù)消掉的. 實際上要使 的項為零，需增加3個方程，要確定4個參數(shù) ，這是不可能的. 故 的顯式R-K方法的階只能是 ，而不能得到三階公式. 35 9.3.3 三階與四階顯式R

13、-K方法 要得到三階顯式R-K方法，必須 . 此時（3.4），（3.5）的公式表示為 （3.11）其中 及 均為待定參數(shù)，公式（3.11）的局部截斷誤差為 只要將 按二元函數(shù)泰勒展開，使 ，可36得待定參數(shù)滿足方程 （3.12）這是8個未知數(shù)6個方程的方程組，解也不是唯一的. 可以得到很多公式. 37滿足條件（3.12）的公式（3.11）統(tǒng)稱為三階R-K公式. 一個常見的公式為 此公式稱為庫塔三階方法. 繼續(xù)上述過程，經(jīng)過較復雜的數(shù)學演算，可以導出各種四階龍格-庫塔公式，下列經(jīng)典公式是其中常用的一個： 38 四階龍格-庫塔方法的每一步需要計算四次函數(shù)值 ，可以證明其截斷誤差為 . 39 例3

14、設取步長 ，從 直到 用四階龍格-庫塔方法求解初值問題 解 這里，經(jīng)典的四階龍格-庫塔公式（3.13）具有形式 4041 比較例3和例2的計算結(jié)果，顯然以龍格-庫塔方法的精度為高. 雖然四階龍格-庫塔方法的計算量（每一步要4次計算函數(shù) ）比改進的歐拉方法（它是一種二階龍格-庫塔方法，每一步只要2次計算函數(shù) ）大一倍，但由于這里放大了步長 ，表9-3和表9-2 所耗費的計算量幾乎相同. 龍格-庫塔方法的推導基于泰勒展開方法，因而它要求所求的解具有較好的光滑性質(zhì). 反之，如果解的光滑性差，那么，使用四階龍格-庫塔方法求得的數(shù)值解，其精度可能反而不如改進的歐拉方法. 42 9.3.4 變步長的龍格-

15、庫塔方法 單從每一步看，步長越小，截斷誤差就越小，但隨著步長的縮小，在一定求解范圍內(nèi)所要完成的步數(shù)就增加了. 步數(shù)的增加不但引起計算量的增大，而且可能導致舍入誤差的嚴重積累. 因此同積分的數(shù)值計算一樣，微分方程的數(shù)值解法也有個選擇步長的問題. 在選擇步長時，需要考慮兩個問題： 1 怎樣衡量和檢驗計算結(jié)果的精度？ 2 如何依據(jù)所獲得的精度處理步長？ 考察經(jīng)典的四階龍格-庫塔公式（3.13）. 43 從節(jié)點 出發(fā)，先以 為步長求出一個近似值，記為 ，由于公式的局部截斷誤差為 ，故有 （3.14）然后將步長折半，即取 為步長從 跨兩步到 ，再求得一個近似值 ，每跨一步的截斷誤差是 ，因此有 （3.15）比較（3.14）式和（3.15）式我們看到，步長折半后