常微分方程的數值解法的

收斂性、穩(wěn)定性

第7章--2

以上我們討論了求解問題（7-1）,（7-2）的單步法和多步法。應關注三個問題：、數值方法的局部截斷誤差和階二、在離散點tn處的數值解un是否收斂到精確解u(tn)

三、數值方法的穩(wěn)定性具體說，對于上述兩類方法求近似解（數值解）還誤差估計、收斂性和穩(wěn)定性。對于第一個問題前面我們已經討論過，而關于數值方法收斂性問題我們在這里不詳細討論，只給出一些基本結論性的結果，即：對單步法，當方法的階p≥1時，有整體誤差故有，因此方法是收斂的。對于多步法，若方法是k步p階法，那么（7-24）是

一個k階差分方程，引入多步法（7-24）的第一特征多項式和第二特征多項式：定義7.1

若（7-24）的第一特征多項式ρ(λ)的所有

根在單位圓內或圓上（︱λ︱≤1），且位于單位圓周上的根都是單根，稱多步法（7-24）滿足根條件。第二特征多項式第一特征多項式定理7.2

若線性多步法(7-24)的階p≥1，且滿足根條件，則方法是收斂的。對于常用的數值方法都是滿足收斂性條件的。下面我們著重討論第三個問題，即數值方法的穩(wěn)是有誤差的，且這些誤差將在計算中傳遞下去。定性問題。誤差積累無限增長，則會歪曲真解，這樣的算法是不如果能用的。用多步法計算時，各種因素如初值精確解為考慮二步三階顯式法：例如

初值問題取步長h=0.1，初值u0=1，附加值：

精確解數值解01.00000001.00000000.11.02010001.02010000.21.08160001.08120000.31.18810001.18923850.41.34560001.33886600.51.56250001.5929935………1.04.0000000-68.6398041.04.8841000+367.26392………2.025.0000000-6.96×108數值結果表在開始幾步數值解與精確解符合，但在再往后算，數值解的誤差則急劇增長，完全歪曲了真解.通常人們都是通過模型方程來討論方法的數值穩(wěn)定性。(7-32)而一般形式的一階微分方程總能化成（7-32）的形式。。因為實際計算時，h是固定的。

當某一步un有舍入誤差時，

若以后的計算中不會逐步擴大，稱這種穩(wěn)定性為絕對穩(wěn)定性。此后，若不做特殊說明，都是指絕對穩(wěn)定性。模型方程為：本書中數值方法的穩(wěn)定性也是如此。前提是求解好條件問題，其中Re(μ)＜0。另外，我們也不考慮h→0時方法的漸近穩(wěn)定性

例如，對最簡單的Euler法(7-33)用其求解模型方程（7-32）得到

當un有舍入誤差時，其近似解為，從而有取，得到誤差傳播方程記，只要都不會惡性發(fā)展，此時方法絕對穩(wěn)定。，則顯式Euler方法的解和誤差從可得

即時，(-1,0)為圓心，1為半徑的單位圓。又由于實數μ＜0，（7-33）絕對穩(wěn)定，若μ為復數，在的復平面上，則表示為以絕對穩(wěn)定區(qū)域絕對穩(wěn)定區(qū)間

定義7.2

一個數值方法用于求解模型問題（7-32），若在

平面中的某一區(qū)域D中方法都是絕對穩(wěn)定的，而在區(qū)域D外，方法是不穩(wěn)定的，則稱D是方法的絕對穩(wěn)定區(qū)域；絕對穩(wěn)定區(qū)間。它與實軸的交稱為例如，顯式Euler方法的絕對穩(wěn)定區(qū)域、區(qū)間。如圖現在考察多步法(7-24)，將它用于解模型方程(7-32)得到k階線性差分方程(7-34)

若取，則記（7-34）的特征方程為(7-35)

其中由k階線性差分方程的性質我們可以得到如下結論，區(qū)域：

例如，對于k=1時，考慮隱式方法中最簡單的后退Euler法方程（7-35）的根都在單位圓內（︱λ︱＜1)，則線性多步法（7-4）關于

絕對穩(wěn)定,其絕對穩(wěn)定域是復平面

上的其特征方程為：若特征得當時，故就是隱式Euler法的絕對穩(wěn)定區(qū)域。當μ＜0為實數時，絕對穩(wěn)定區(qū)間為

(-∞,0)。平面上以（1.0）為圓心的單位圓外區(qū)域。它是當Reμ＜0時，它位于平面上y軸左側區(qū)域。又如，梯形法其特征方程為：其根當Reμ＜0時，故梯形公式平面的左半平面。絕對穩(wěn)定區(qū)間為(-∞，0)。的絕對穩(wěn)定域是這樣檢驗絕對穩(wěn)定性歸結為檢驗特征方程(7-35)的根是否在單位圓內（︱λ︱＜1)。對此有很多判別法，如Schur準則、軌跡法。k=1～4的隱式Adams類方法的絕對穩(wěn)定區(qū)間（μ＜0為實數）。步階絕對穩(wěn)定區(qū)間12（-∞，0）23（-6.0，0）34（-3.0，0）45（-1.8，0）實系數二次方程λ2-bλ-c=0的根在單位圓內的充要條件為:

這里我們給出一種簡單的、常用的判別法：例證明求解一階常微分方程初值問題：的差分格式收斂并求其局部截斷誤差主項、絕對穩(wěn)定區(qū)間。解：由差分格式可知，則其特征值滿足根條件。令得λ1=0,λ2=1。故此為隱式二步三階法，其局部截斷誤差主項為：注意，從而由定理7.2

可知，此方法收斂。而自然成立。得即有可得其絕對穩(wěn)定區(qū)間：又其特征方程為而使得

︱λ︱＜1的充要條件為：現在再由進一步而自然成立。顯式Runge-Kutta法

第7章--37.1.4

四階顯式Runge-Kutta法

通過觀察我們發(fā)現顯式Euler法和隱Euler法各用到了u(t)在[t,t+h]上的一個一階導數值，它們都是一階方法。改進的Euler法用到了u(t)在[t,t+h]上的兩個一階導數值，它

梯形法和們都是二階方法。我們要研究的Runge-Kutta方法是一種高階單步法，它使用u(t)在[t,t+h]上的斜率f在一些點的值非線性表示使得其局部截斷誤差的階和Taylor展開法相等。

Euler是最簡單的單步法。單步法不需要附加初值，所需的存儲量小，改變步長靈活，但線性單步法的階最高為2，Taylor展開法，用在同一點(tn,un)的高階導數表示，這不便于計算。先引進若干記號，首先[t,t+h]取上的m個點:令

Runge-Kutta矩陣B為嚴格下三角矩陣：滿足顯式Runge-Kutta

公式假設三組系數已給定，則求解（7-1），（7-2）的一般（7-12）其中(7-13)(7-14)顯式Runge-Kutta法的計算過程如下：現在推導一些常用的計算方案，特別地，給出m=3顯式首先將u(t+h)在t處展開到h的三次冪,即：(7-15)其中(7-16)

Runge-Kutta法的推導。其次，由二元函數f(t,u(t))在(t,u)點處的Taylor展開式可得：

于是，將k1,k2,k3代入（7-13）中，即(7-17)由(7-16)已得其中合并f(t,u,h)展開式中的各階hl(l=0,1,2)的系數，得比較和的同次冪系數，可得（一）m=1此時c2=c3=0,f(t,u,h)=c1f,

比較h的零次冪，知方法（7-21）為一級一階Runge-Kutta法，實際上為Euler法。

（二）m=2,此時

c3=0，則與比較1,h的系數，則

它有無窮多組解，從而有無窮多個二級二階方法。（1）稱為中點法。此時(7-18)三個常見的方法是：（2）稱為改進的Euler法。此時(7-19)

（3）此時（三）m=3比較（7-16）和（7-17），令

f,h,h2的系數相等，并注意的任意性，得四個方程不能完全確定六個系數，因此這是含兩個參數的三級三階方法類。常見方案有：

Heun三階方法。

此時取(7-20)（2）Kutta三階方法，(7-21)此時（四）m=4將（7-16）和（7-17）展開到h3,比較的系數，則含13個待定系數的11個方程，由此得到含兩個參數的四級四階Runge-Kutta方法類，其中最常用的有以下兩個方法：經典四階Runge-Kutta方法：（7-22）Butcher表分別為：

以上討論的是m級Runge-Kutta法在m=1,2,3,4時，可分別得到最高階級一、二、三、四階，但是，通常m級Runge-Kutta

方法最高階不一定是m階。若設p(m)是m級Runge-Kutta方法可達到的最高階，可證：改進的Euler法計算公式為：經典Runge-Kutta法計算公式為：例1

分別用Euler法，改進的Euler法(7-27)和經典Runge-Kutta法(7-30)求解初值問題:解：Euler法計算公式為：

三個方法計算結果比較表n

tn精確解u(tn)Euler法改進Euler法經典Runge-Kutta法數值解un誤差數值解un誤差數值解un誤差00000000010.50.4333330.5000000.0666670.4000000.0333330.4332180.00011521.00.6666670.8000000.1333330.6350000.0316670.6663120.00035531.50.8076920.9000000.0923080.7875960.0200960.8074230.00026942.00.9333530.9856150.0512820.9210250.0123080.9331560.000171作比較

，計算結果見下表：取步長h=0.5,tn=0.5n,n=0,1,2,3。并與精確解：下面考察Runge-Kutta法的絕對穩(wěn)定性。根據定義，對m級p階Runge-Kutta法（7-12）取f=mu，則（其中