常微分方程數(shù)值解

1、第七章 常微分方程數(shù)值解§1 引言一 一階初值問題解的存在唯一性一階常微分方程初值問題 （*）其中是平面某一區(qū)域D上的連續(xù)函數(shù)，如果，滿足 存在，并滿足方程那么是初值問題（*），在上的解。對于（*）是否有唯一解？對還要附加一些條件。定義 如果存在正常數(shù)，使得對任意有則稱滿足Lipschitz條件，L稱為Lipschitz常數(shù)。 如果 ，那么有導數(shù)有界Þ滿足Lipschitz條件如果存在常數(shù)，使得對一切及有 則稱對滿足Lipschitz條件。同理，只要 f(x,y)對 y的偏導數(shù)有界，則f(x,y) 滿足對y的Lipschitz條件。定理（存在唯一性），設是在上的連續(xù)函數(shù)，而

2、且對滿足Lipschitz條件，則對任意，初值問題（*）在上存在唯一的連續(xù)可微解。二 本章研究的問題例1 ，滿足微分方程和初始條件的解是 y(x)=，無法給出具體表達式。為了對初值問題進行求解，一些簡單問題有解析解，大量非線性問題沒有解析表達式，就是線性問題也不一定有解析解，因此，近似求解和數(shù)值求解常微分方程是非常必要的。 1 初值問題數(shù)值解基本概念已知 (1.1)滿足（1.1）的解是過點（的一條曲線y(x)首先對連續(xù)區(qū)間離散化為常數(shù)。離散點是等矩的也可以是不等距的，下面僅討論等距情況。為方便起見：把在處的精確值記為，其近似值用表示 (1.2) 方程（1.2）叫差分方程 這種求微分方程近似解的

3、過程稱為步進式的方法，計算若用到前面不止一個信息量，叫多步法。 ，即這一步的公式誤差。綜上所述，求微分方程數(shù)值解需要處理以下幾個問題：1 把微分方程（連續(xù)的）離散化為差分方程（離散的）2 用差分方程和初始條件計算出微分方程數(shù)值解3 有關理論（1） 誤差，局部截斷誤差，主局部截斷誤差，階（2） 收斂性（3） 穩(wěn)定性（絕對穩(wěn)定性）本章講授的內(nèi)容：1 單步法：顯示Euler方法，隱式Euler方法，梯形方法，改進Euler方法，Rung-kutta方法2 單步法的收斂性，穩(wěn)定性，相容性3 線性多步法三 預備知識1 一元Taylor多項式 y 2 二元Taylor多項式k)f(x,y)+ , =3 數(shù)

4、值積分4 Lagrang插值§2 簡單數(shù)值方法已知 (2.1) （I）顯式Euler方法（也稱為Euler方法）Euler方法是求常微分方程初值問題的最簡單辦法。上節(jié)的公式(1.2)就是顯示Euler方法 (2.2)問題： 要從(2.1) 離散化得到(2.2)還有什么辦法？1數(shù)值積分方法 把（2.1）寫成如下形式： (*)用左矩形公式近似左邊積分,用得到結(jié)果 (2.2)2 Taylor展開方法略去高階項 即 設是初值問題的解，那么有從而有用近似值由此得出公式（2.2）3數(shù)值微分方法稱為差商，即用差商近似微商，也得出：由此得出公式（2.2） 當已知時，可由公式（2.2）簡單地求出，方法

5、為顯式的，由上的近似值可求出上的近似值，稱為單步公式，（2.2）稱為顯式Euler方法也稱為Euler方法。（II）隱式Euler方法（也叫后退Euler方法） 數(shù)值積分方法對(*) 公式右邊積分用右矩形積分公式得到 (2.3) 也可用Taylor展開方法 略去高階項，并設是初值問題的解，則有即寫為同樣用 的近似值代入有 (2.3)的右端含有，一般，不能直接由（2.3）得出，這種方法是隱式的。（2 .3）稱為隱式Euler方法。該方法也可用數(shù)值微分方法得出。（III）梯形方法在上對上式進行積分有等式右邊積分采用梯形公式近似有 用來代替就得到 （2.4）（2.4）稱為梯形公式(IV) 用公式做計

6、算顯示Euler公式由，但是，對于（2.3）、（2.4）是隱式公式，不能由直接計算出，而是要解方程，一般用迭代方法，以（2.3）為例，取，或用顯式公式求出作為 ， 即 。 當 時，取這樣方法稱為迭代法下面考慮迭代收斂性： 當收斂；稱為迭代收斂條件對于梯形公式（2.4），由于等式右邊含有，因而是隱式方法。用它們來求時必須解方程，一般用迭代求解。取，迭代公式為同樣，當時，??；仿隱式Euler方法推導，梯形公式迭代收斂條件為例1 已知 取 計算到解：用Euler方法把 代入有 梯形方法 把代入有 由于是線性方程，可把隱式方法顯式化，因此不用進行迭代Euler方法與梯形方法計算結(jié)果比較：從數(shù)值結(jié)果看出

7、，梯形公式比Euler公式好，但一般梯形公式需要進行迭代，因此做一步費時；對于“計算效率”的比較，應從精度，耗機時等方面進行比較，還應從實際對精度要求來考慮。（V）予估一校正方法為了消除迭代，出現(xiàn)了予估一校正的方法，先給出粗糙估計，然后再給出稍精確的求解，這是微分方程數(shù)值解常用方法。改進Euler方法 予估 校正或?qū)懗?（2.5）這公式稱為改進的Euler公式，其精度比Euler公式好，比梯形公式稍差些，但是，它是顯示方法。例2 用Euler方法和改進Euler方法解初值問題。 步長取 ；由0計算到3。方程準確解是 解：Euler方法 改進Euler方法 (可以看出改進Euler方法較為精確（

8、VI）顯式單步方法基本概念（局部截斷誤差，主局部截斷誤差，階）已經(jīng)引入了四種方法；Euler方法，隱式Euler方法，梯形方法，以及改進Euler方法，其中二個為隱式方法，Euler方法與改進Euler方法為顯式方法。下面重點討論顯式方法.其中，一般顯式方法可以統(tǒng)一寫成如下形式 （2.6）Euler方法 是典型單步方法，從開始進行計算， 在處微分方程初值問題 的精確解。 與之差 稱為方法在處的整體截斷誤差，當然這與整個計算中每步情況有關，但一般求得較為困難，因此先考慮一步的誤差。1 局部截斷誤差定義2.1 （2.7）稱為顯式單步法（2.6）在處的局部截斷誤差，其中是微分方程初值問題在處的精確解

9、。顯然對于一般單步方法 如果每步是精確的，即有 那么 即 為是精確值，用顯式單步法計算一步的誤差因此稱為局部截斷誤差。關于誤差僅討論截斷誤差，舍入誤差暫不討論。2 單步法的階 定義2.2 設y(x)是（2.1）的精確解，（2.6）為顯示單步法，若 y(x+h)-y(x)-h);h)= (2.8) 則稱整數(shù)p是單步法（2.6）的階，稱（2.6）是p 階方法，且局部截斷誤差 3 主局部截斷誤差 定義2.3 若（2.6）是p階方法，其局部截斷誤差 = (2.9)則稱是單步法（2.6）的主局部截斷誤差或稱局部截斷誤差主項。例3 寫出顯式Euler方法的局部截斷誤差，并求階，主局部截斷誤差解： 按定義：

10、 = =顯示Euler方法是一階方法，主局部截斷誤差是同理可得隱式Euler方法也是一階方法，主局部截斷誤差是梯形方法 梯形方法是二階方法，主局部截斷誤差是討論改進Euler方法的局部截斷誤差 要用到二元Taylor展開， 注意到：， 可以看出 此方法是二階的§3 Rung-Kutta方法下面我們要解決的問題是如何構(gòu)造高階單步法？（I）用Taylor展開構(gòu)造高階方法利用方程 把這些導數(shù)代入的有限項展式，就得到Taylor展開方法。例如 Euler公式Euler方法是1階方法，而上面是2階方法。 這是三階方法。 還可以類似推導下去，得到高階方法，但是要計算很多偏導數(shù)。特別是f復雜時，這

11、方法不實用。我們希望仍采用這樣思想，又不計算導數(shù)而構(gòu)造出高階方法，由此導出Runge-kutta方法。（II）Runge-Kutta方法Runge首先提出間接采用Taylor展開方法，即用在幾個節(jié)點上函數(shù)值的線性組合來代替f的導數(shù)，然后按Taylor展開，確定其系數(shù)，以期提高方法精度，這樣既避免f的導數(shù)計算，同時又保證了精度。一 RK方法的一般形式 （3.1）其中為待定權(quán)因子，R為使用的f值的個數(shù)，表示如下： （3.2） 具體寫出有： （3.3） 其中參數(shù)為提高精度創(chuàng)造了條件；用R個稱為R級R-K方法。二 構(gòu)造RK方法例1 取R=1 此時只有一個K1并有 （3.4）利用Taylor公式來確定C

12、1設為微方程初值問題的光滑解。比較得 是顯示Euler公式，局部截斷誤差為 下面考慮則 （3.5）仍假定，對在處作Taylor展開把代入下式， 有= =+=得出 （3.6） 有四個未知數(shù)，但僅有三個方程。把四個參數(shù)滿足（3.6）的一組方法叫做二階RK方法。取為自由參數(shù)， ，取 ，得 （3.7） 中點公式取 ，得 改進Euler公式 （3.8） Heun二階公式。對于三級方法 （3.9）利用Taylor展開，可以構(gòu)造三階方法。八個參數(shù)滿足如下方程 (3.10)這是8個未知數(shù)，6個方程，有二個自由參數(shù)可選取。 對于R=4，可導出4階方法。此時有13個參數(shù)，11個方程，推導更為復雜，常用的有經(jīng)典Ru

13、nge-Kutta方法。三 四階經(jīng)典RK方法 (3.11)其中 這是經(jīng)常使用的方法，一般都有標準程序可調(diào)用例1 用Euler方法，h=0.025，改進Euler方法（二階Runge-Kutta方法），h=0.05以及4階Runge-Kutta方法，h=0.1解初值問題：比較計算結(jié)果利用可以得出精確解為Euler方法 改進Euler方法 4階Runge-Kutta方法 計算結(jié)果： Euler 改進Euler 4階R-K 精確解§4 單步法的收斂性和絕對穩(wěn)定性(I)收斂性已知 (4.1) (4.2)記=y( 差分方程精確解 - - 1 定義 定義4.1 設初值問題（4.1）的解存在唯一，

14、若單步法（4.2）產(chǎn)生的近似解對任一固定的x 則稱該單步法是收斂的。例1 用顯示Euler方法求解，討論收斂性解：設 方法收斂2 收斂定理 用定義只能針對具體的方程來討論，對于一般方程很難討論，下面定理給出了一般的判斷收斂的辦法。定理4.1 若初值問題的一個單步方法 的局部截斷誤差為精確成立，并且對y滿足Lipschitz條件，即存在正常數(shù)L，使 有 成立則單步法收斂并有整體截斷誤差 （4.3）證明: 根據(jù)收斂定義因此必須估計 事實上，即估計 局部截斷誤差由定理條件 這樣遞推下去有取，并且 并有 定理4.1 告訴我們，初值沒有誤差，只要單步法的階p1，增量函數(shù)滿足關于y的Li氏條件，則單步法收

15、斂，且有誤差估計式（4.3） 例2 討論Euler方法的收斂性解： 是一階方法 由于對滿足Lipschitz條件對也滿足Lipschitz條件。應用定理知Euler方法收斂。例3 討論Runge-Kutta方法，對R=2的改進Euler方法。的收斂性解： p=2 假定步長，取為 關于的Lipschitz常數(shù) 所以，改進Euler方法收斂對于一般Runge-Kutta方法： ， 記 對滿足Lipschitz條件 (為權(quán)應大于0) 于是，存在，當時有 對于 還可以取，使當時有同樣可得：由此，取有（收斂性證明中，因此可取h充分?。↖I）相容性1 相容的定義收斂性定理中要求局部截斷誤差 若按變量在處

16、作Taylor展開，那么有，而，也就是說是有界量這相當于含h的項必為零，即 滿足微分方程，由此有 定義4.2 單步法 滿足條件 （4.4）則稱單步法（4.2）與微分方程初值問題（4.1）是相容的。事實上，相容的方法必有 相容方法至少是一階的。，則單步法是相容的2 相容方法與方法收斂的關系對于相容單步法，若對滿足Lipschitz條件方法收斂即 即有 令 時，有 即計算格式趨于微分方程。相容本質(zhì)的意義在于“差分格式”收斂于“微分方程”。（III）穩(wěn)定性（絕對穩(wěn)定性） 研究收斂性時，沒有考慮舍人誤差的影響，實際上每步計算都有舍人誤差，穩(wěn)定性就是研究舍人誤差傳播問題，當求解過程中舍人誤差不增長，則稱

17、該數(shù)值方法是穩(wěn)定的。有單步法 設是帶舍入誤差的值，而是單步法精確計算而得的準確值。 即 在之間為使舍入誤差不增長，則應有 由于與有關，所以穩(wěn)定性與方程（微分）右端項有關，而右端項各式各樣，如何來測試某個方法的穩(wěn)定性？一般把該數(shù)值方法用于一個“模型方程”（試驗方程），（為復數(shù)， （4.5）來考察方法的穩(wěn)定性，模型方程的解析解是y(x)=c。選擇模型方程的原因： 討論方便，如果對這樣簡單方程不穩(wěn)定，那么復雜方程也不穩(wěn)定 一般方程可簡化為模型方程形式 模型方程的解是穩(wěn)定的例4 討論顯示Euler方法的穩(wěn)定性用顯示Euler方法解模型方程 有 （4.6）相應誤差方程 （4.7）可以看出，誤差方程（4.

18、7）與原來單步法（4.6）一致，這是因為模型方程是常系數(shù)線性方程而得到的，用于模型方程，考慮的增長與誤差增長是一樣的。因此，討論穩(wěn)定性直接用（4.6）就行了。 即 顯然， 時，誤差不增長，Euler方法絕對穩(wěn)定。對于隱式Euler方法對于梯形方法 對于改進Euler方法 對于一般單步法用于模型方程有 先考慮一下的性態(tài)令 有誤差 有誤差有誤差。 有誤差 有誤差如果很大，產(chǎn)生不穩(wěn)定，依賴于方法選取，Euler方法； 方法絕對穩(wěn)定下面給出絕對穩(wěn)定定義：定義4.3 單步方法 解模型問題 ，得到的解，若，則稱單步法是絕對穩(wěn)定的。在復平面中，滿足的區(qū)域，稱為單步法的絕穩(wěn)定區(qū)域，它與實軸的交集稱為絕對穩(wěn)定區(qū)

19、間。 試寫出隱式Euler方法，梯形方法，改進Euler方法的E(對于隱式Euler方法對于梯形方法 對于改進Euler方法 例5 求Euler方法的絕對穩(wěn)定區(qū)域和絕對穩(wěn)定區(qū)間解：，對于Euler方法有 為復數(shù)，令，那么，這是以（-1，0）為圓心1為半徑的單位圓的內(nèi)部，此為Euler法的絕對穩(wěn)定性區(qū)域。區(qū)域與實軸交集是區(qū)間，此為絕對穩(wěn)定性區(qū)間。要使方法絕對穩(wěn)定步長例6 考察隱式Euler方法的絕對穩(wěn)定性解：E(1,0)為圓心單位園的外部，絕對穩(wěn)定區(qū)間是（，要使方法絕對穩(wěn)定對步長 h沒有限制例7 考察改進Euler方法的絕對穩(wěn)定性解： 絕對穩(wěn)定性條件為 此等價于 改進Euler方法的絕對穩(wěn)定區(qū)間

20、為。要使方法絕對穩(wěn)定要求步長 可以證明梯形方法的絕對穩(wěn)定區(qū)間為，要使方法絕對穩(wěn)定對步長h沒有限制。對于4階經(jīng)典R-K方法 應用到方程 代入有： ，絕對穩(wěn)定當為實數(shù)，絕對穩(wěn)定區(qū)間是 -2.785例8 用改進Euler方法，梯形方法解初值問題 1 寫出計算公式2 討論絕對穩(wěn)定性3 用定義討論收斂性解： 1 改進Euler方法 梯形方法 2 若計算所以要使方法絕對穩(wěn)定，有<1h<0.4 同理可得梯形方法的絕對穩(wěn)定性，對于任意的h>0,因為E(分母的絕對值大于分子的絕對值，所以，方法絕對穩(wěn)定區(qū)間是 -，對步長h沒有限制,3 由改進Euler方法計算公式遞推得到設 x是任意固定點，且x

21、>0, h=近似解序列收斂到方程的精確解。定義4.4 絕對穩(wěn)定性區(qū)域包含整個左半平面，這種方法稱為A-穩(wěn)定的。向后Euler方法和梯形方法都是A穩(wěn)定的。§5 線性多步法（I）基本概念顯式單步法一般為 即由上近似值上近似值，隱式也是由解方程（迭代）求出，即求上的近似值，僅與前面一個點的近似值相聯(lián)系。提高精度的Runge-Kutta方法，一般也不很簡單。問題：還有沒有提高精度的辦法？回答是肯定的，就是采用前面多個信息，比如：上的近似來求，這樣的數(shù)值方法稱為多步方法。一 線性多步法的一般形式初值問題 (5.1) 先看一個簡單例子，對方程（5.1）兩端做積分有:右邊采用Simpson求

22、積公式，這樣有用來表示得計算，要用到以及相應的，并且公式中對是線性的。這樣方法稱為線性二步法。已知 求處y(x)的近似值線性多步法的一般形式為：其中：k 若 ，可以顯式計算，（5.2）為顯式方法。若，為隱式方法。為求解，必須進行迭代 其中： 可由相應顯式給定，迭代收斂條件 ， L為關于的Lipschitz常數(shù)線性單步法是線性步方法的特例，例如 有 不同可得各種顯、隱單步法。二 局部截斷誤差，階，主局部截斷誤差定義5.1 （局部截斷誤差） 設 是 的解， 稱為線性步法（5.2）在的局部截斷誤差。 定義5.2 把按h展開的首項稱為主局部截斷誤差,即 （5.4）為(5.2)的主局部截斷誤差。相應的多

23、步法稱為P階方法，這個定義包含了線性步法（當然包含線性單步法），特別包含了線性隱式單步方法。例1 Simpson公式 用微分方程充分光滑解代入，寫出局部截斷誤差表達式，并用Taylor展開有 局部截斷誤差主項為Simpsion方法是四階方法。（II）構(gòu)造線性多步法 一 基于數(shù)值積分的方法1 Adams方法 1.1 顯式Adams方法Adams方法是基于數(shù)值積分的方法，但積分區(qū)間為，對微分方程在區(qū)間積分有為求近似積分，采用插值多項式來近似這樣可以求得k-1次Lagrange插值多次式其中為對應點上的k-1次插值多項式基函數(shù) 從而有 （5.5）其中 這是顯式Adams方法，稱Adams-Bashf

24、orth方法。例2 k=2 用 做線性插值其中 容易得到：得多步法 （二步法） 其局部截斷誤差 1.2 隱式Adams方法在顯式Adams方法中，采用上的來插值求，這相當于外推，精度受到影響，改進辦法是把作為一個插值的節(jié)點，即共有k+1個節(jié)點，可得插值多項式取，那么 直接用局部截斷誤差定義求是一樣的。一般形式有：（隱式Adams方法or Adams-Moulton方法） ）下面對常用顯式Adams和隱式Adams列表如下 k為步數(shù)，P為方法的階，Cp+1是局部截斷誤差主項的系數(shù)（主局部截斷誤差的系數(shù)） 顯式Adams（Adams-Bashforth）（書P303表9.5）隱式Adams (Ad

25、ams-Moulton方法)（書P304表9.6）二 待定系數(shù)方法（基于Taylor展開的方法）設y是微分方程的充分光滑解，已知線性多步法（5.2），恰當選擇參數(shù)使其為P階方法。即把（5.4）中右端各項在處做Taylor展開，要使方法為P階的，必須的系數(shù)全為零，于是得到這些參數(shù)滿足的方程組，解此方程組求出,就得到了P階方法。 例4 已知確定參數(shù)使方法階盡量高。解： 把上式右端各項在處做Taylor展開，于是得到解上方程組得出 主局部截斷誤差 四階方法，稱為Milne方法例5 考慮4步方法即 要求方法是4階的. 解：寫出局部截斷誤差表達式 把上式右端各項在處做Taylor展開，令得到5個方程，有

26、未知數(shù)9個取 得由此得到 的Adams-Bashforth方法如果取得：稱為Milne公式 局部截斷誤差為 推導如下： 例6 考慮三步方法要求方法是4階的。解：用同樣的方法得到有7個未知數(shù)的5個方程，取 ，可以解得得到Hamming方法三 預估一校正方法對于隱式方法，每一節(jié)點上近似值用迭代方法得到，這必須大大增加計算量。若用一個恰當?shù)娘@式方法，求出作為隱式方法的預估值，然后用隱式方法對預估值作校正，并以這個校正值作為所求節(jié)點上的，那么將克服迭代法缺點。用Eular方法作預估，用梯形方法作校正的預估校正方法。 預估值 計算函數(shù)值 相當于修正Euler方法 通常使用預估校正方法有： Adams-Bashforth-Moulton方法4步4階顯式Adams方法作為預估，3步4階隱式Adams方法 作校正 預估 求值 校正這是使用相當廣泛的方法如果h充分小，由可知，當h充分小，從而有 其中 可以看出，誤差來估計而不用 The milne-Simpson Me