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

1、第九章 常微分方程的數(shù)值解法在自然科學(xué)的許多領(lǐng)域中，都會遇到常微分方程的求解問題。然而，我們知道，只有少數(shù)十分簡單的微分方程能夠用初等方法求得它們的解，多數(shù)情形只能利用近似方法求解。在常微分方程課中已經(jīng)講過的級數(shù)解法，逐步逼近法等就是近似解法。這些方法可以給出解的近似表達(dá)式，通常稱為近似解析方法。還有一類近似方法稱為數(shù)值方法，它可以給出解在一些離散點上的近似值。利用計算機解微分方程主要使用數(shù)值方法。我們考慮一階常微分方程初值問題在區(qū)間a, b上的解，其中f (x, y)為x, y的已知函數(shù)，y0為給定的初始值，將上述問題的精確解記為y(x)。數(shù)值方法的基本思想是：在解的存在區(qū)間上取n + 1個

2、節(jié)點這里差，i = 0,1, , n稱為由xi到xi+1的步長。這些hi可以不相等，但一般取成相等的，這時。在這些節(jié)點上采用離散化方法，（通常用數(shù)值積分、微分。泰勒展開等）將上述初值問題化成關(guān)于離散變量的相應(yīng)問題。把這個相應(yīng)問題的解yn作為y(xn)的近似值。這樣求得的yn就是上述初值問題在節(jié)點xn上的數(shù)值解。一般說來，不同的離散化導(dǎo)致不同的方法。1 歐拉法與改進(jìn)歐拉法1.歐拉法1對常微分方程初始問題用數(shù)值方法求解時，我們總是認(rèn)為(9.1)、(9.2)的解存在且唯一。歐拉法是解初值問題的最簡單的數(shù)值方法。從（9.2）式由于y (x0) = y0已給定，因而可以算出設(shè)x1 = h充分小，則近似地

3、有：(9.3)記 從而我們可以取作為y (x1)的近似值。利用y1及f (x1, y1)又可以算出y(x2)的近似值：一般地，在任意點xn+1 = (n + 1)h處y(x)的近似值由下式給出（9.4）這就是歐拉法的計算公式，h稱為步長。不難看出，近似解的誤差首先是由差商近似代替微商（見(9.3)）引起的，這種近似代替所產(chǎn)生的誤差稱為截斷誤差。還有一種誤差稱為舍入誤差，這種誤差是由于利用（9.4）進(jìn)行計算時數(shù)值舍入引起的。例9.1 用歐拉法求初值問題當(dāng)h = 0.02時在區(qū)間0, 0.10上的數(shù)值解。解 把代入歐拉法計算公式。就得具體計算結(jié)果如下表：nxnyny(xn)en = y(xn) -

4、 yn001.00001.0000010.020.98200.98250.000520.040.96500.96600.000530.060.94890.95030.001440.080.93360.93540.001850.100.91920.9230.0021在上表中y(xn)列，乃是初值問題（9.5）、（9.6）的真解在xn上的值。為近似值yn的誤差。從表中可以看出，隨著n的增大，誤差也在增大，所以說，歐拉法計算簡便，對一些問題有較大的使用價值，但是，它的誤差較大，所得的數(shù)值解精確度不高。2改進(jìn)歐拉法為了構(gòu)造比較精確的數(shù)值方法，我們從另一角度重新分析一下初值問題。一般說來，一階方程的初值

5、問題與積分方程（9.7）是等價的，當(dāng)x = x1時，（9.8）要得到y(tǒng)(x1)的值，就必須計算出（9.8）式右端的積分。但積分式中含有未知函數(shù)，無法直接計算,只好借助于數(shù)值積分。假如用矩形法進(jìn)行數(shù)值積分則因此有可見，用矩形法計算右端的積分與用歐拉法計出的結(jié)果完全相同。因此也可以說歐拉法的精度之所以很低是由于采用矩形法計算右端積分的結(jié)果。可以想象，用梯形公式計算（9.8）式右端的積分，可期望得到較高的精度。這時將這個結(jié)果代入（9.3）并將其中的y(x1)用y1近似代替，則得這里得到了一個含有y1的方程式，如果能從中解出y1，用它作為y(x1)的近似值，可以認(rèn)為比用歐拉法得出的結(jié)果要好些。仿照求y

6、1的方法，可以逐個地求出y2, y3，。一般地當(dāng)求出yn以后，要求yn+1，則可歸結(jié)為解方程：這個方法稱為梯形法則。用梯形法則求解，需要解含有yn+1的方程式，這常常很不容易。為此，在實際計算時，可將歐拉法與梯形法則相結(jié)合，計算公式為（9.9）這就是先用歐拉法由(xn, yn)得出y(xn+1)的初始近似值，然后用（9.9）中第二式進(jìn)行迭代，反復(fù)改進(jìn)這個近似值，直到（e為所允許的誤差）為止，并把取作y(xn+1)的近似值yn+1。這個方法就叫改進(jìn)歐拉方法。顯然，應(yīng)用改進(jìn)歐拉法，如果序列收斂，它的極限便滿足方程即序列的極限可取作yn+1?？梢宰C明，如果有界，則只要h取得適當(dāng)小，上述序列必定收斂。

7、這樣當(dāng)h取得充分小，就可保證上述迭代過程收斂到一個解。當(dāng)步長h取得適當(dāng)時，歐拉方法算出的值已是較好的近似，因此改進(jìn)歐拉法收斂很快，通常只需二、三次迭代即可。如果迭代很多步仍不收斂，這表明表長h選的過大，應(yīng)縮小步長后再計算。通常把（9.9）叫做預(yù)報校正公式，其中第一式叫預(yù)報公式，第二式叫校正公式。這個公式還可寫為(9.9)3公式的截斷誤差現(xiàn)在來考察兩個公式的截斷誤差：有多大？這里假定前一步得的結(jié)果yn=y(xn)是準(zhǔn)確的。寫出y(xn+1)的泰勒展開式為由歐拉法得兩式相減得即歐拉法的截斷誤差為0(h2)，當(dāng)h 0時它與h2是同階無窮小量。對于改進(jìn)的歐拉方法，我們以迭代一次的預(yù)報校正格式(9.9)

8、為例來說明。因為用它們代入（9.9）式第二式即得這里末把含有h的三次冪以上的項寫出，因此有即迭代一次的預(yù)報校正格式（9.9）的截斷誤差為0(h3)?？梢姼倪M(jìn)的歐拉方法比歐拉法的階提高了。例9.2在區(qū)間0, 1.5上，取h = 0.1，求解。解：（1）用歐拉法計算公式如下：（2）用迭代一次的改進(jìn)歐拉法計算公式如下：本題的精確解為，可用來檢驗近似解的精確程度。計算結(jié)果如下表：xn歐拉法yn迭代一次改進(jìn)歐拉法yn準(zhǔn)確解01110.11.11.0959091.0954450.21.1918181.1840961.1832160.31.2774381.2602011.2649110.41.3582131

9、.3433601.3416410.51.4351331.4161021.4142140.61.5089661.4829561.4832400.71.5803381.5525151.5491930.81.6497831.6164761.6124520.91.7177791.6781681.6733201.01.7847701.7378691.7320511.11.851181.7958221.7888541.21.9174641.8522421.8439091.31.9840461.9073231.8973671.42.0514041.9612531.9493591.52.1200522.014

10、2072.0000002 龍格庫塔法由上節(jié)知道，截斷誤差的階是衡量一個方法精度高低的主要依據(jù)。能否用提高截斷誤差階來提高方法的精確度呢？回答是肯定的。本節(jié)介紹的泰勒級數(shù)法和龍格庫塔法就是基于這種思想構(gòu)造出來的。1泰勒級數(shù)法如果初值問題（9.10）的精確解y(x)在x0, x上存在k + 1階導(dǎo)數(shù)且連續(xù)，那么由泰勒公式（9.11）其中截斷誤差為（9.12）略去截斷誤差，并用近似值代替真值則得（9.13）用公式(9.13)解初值問題的數(shù)值方法稱為泰勒級數(shù)法，當(dāng) k = 3時，（9.13）變?yōu)椋?.14）這時的截斷誤差是（9.15）從截斷誤差的表示式中看出，如果微分方程（9.10）的真解y = y(

11、x)為次數(shù)不超過3的多項式時，公式（9.14）精確地成立，因此（9.14）是3階方法。例9.3 導(dǎo)出用三階泰勒級數(shù)法解方程的計算公式解:因故而其中表示f(x, y)對x的k階偏導(dǎo)數(shù)在x = xn點上的值。泰勒級數(shù)法只要初值問題的真解充分光滑，就可以獲得精確度較高的數(shù)值解。但是須計算y(x)的各階導(dǎo)數(shù)，這當(dāng)f (x, y)的表達(dá)式復(fù)雜時是很繁瑣的。因此泰勒級數(shù)法一般只用于求“表頭”（即開頭幾點的數(shù)值解，如y1, y2, y3, , y4等）。另外，用上述級數(shù)法計算表頭時，還可以得到選擇步長h的信息。假定我們要求計算誤差不超過e ，那么，當(dāng)h滿足條件（A）（B）時，應(yīng)該認(rèn)為是最好的。因為，當(dāng)條件（

12、A）不滿足時，達(dá)不到指定精確度，而當(dāng)條件（B）不滿足則表明h過小。能否構(gòu)造一種格式，既保留泰勒級數(shù)法精確度較高的優(yōu)點，又避免過多的計算f (x, y)的各階偏導(dǎo)數(shù)呢？下面介紹的龍格庫塔方法就能辦到這一點。2龍格庫塔法從理論上講，只要函數(shù)y = y(x)在區(qū)間a, b上充分光滑，那么它的各階導(dǎo)數(shù)值y(k)(xn)與函數(shù)y(x)在區(qū)間a, b上某些點的值就相互有聯(lián)系，就是說，函數(shù)值可用各階導(dǎo)數(shù)值近似地表示出來，反之，各階導(dǎo)數(shù)值也可用函數(shù)在一些點上值的線性組合近似地表示出來。事實上，歐拉法和改進(jìn)歐拉法也可以看成是導(dǎo)數(shù)值用函數(shù)值的線性組合表示的特例，例如，改進(jìn)歐拉法可以寫成（9.16）此公式也可稱為二

13、階龍格庫塔公式。為了導(dǎo)出龍格庫塔法的一般公式，我們?nèi)∪缦碌木€性組合形式，（9.17）其中（9.18）即而w1, w2, wv；b1 =0, b2, b3, bv；a21, a31, , avv-1除b1=0外均為待定系數(shù)。適當(dāng)選取這些系數(shù)，使得局部截斷誤差的階盡可能高即可。顯然，當(dāng)g = 1時，（9.12）式就是歐拉公式。下面我們先導(dǎo)出g =2時的公式。將k1, k2在同一點（xn, yn）泰勒展開，則有（9.19）將(9.19)代入(9.17)并與y(xn+h)在xn點的泰勒展開式：逐項比較，令h、h2項的系數(shù)相等，便得到把b2作為自由參數(shù)來確定w1和w2，如取b2 = 1，則w1 = w2

14、 = ，a21 = 1，這時（9.17）正好就是改進(jìn)的歐拉方法，截斷誤差的階為0(h3)。對于g = 3的情形，我們也可以完全仿上述方法推導(dǎo)出三階龍格庫塔公式。這時參數(shù)滿足下列條件（9.20）（9.20）比較簡單的一組解為：b2 =，b3 = 1，a21 = ，a31 = -1，a32 = 2，w1 = ，w2 = ，w3 = 將它代入（9.17）得(9.21)這就是三階龍格庫塔公式。當(dāng)然，參數(shù)的不同選取公式（9.21）就有不一樣的形式，但它們的截斷誤差階都是0(h4)。通常人們所說的龍格庫塔法是指四階而言的。我們可以仿照二階的情形推導(dǎo)出此公式，不過太繁雜，此處從略，常用的四階公式是（9.22

15、）公式(9.22)的截斷誤差階為0(h5)。龍格庫塔法有精確度高、收斂、穩(wěn)定（在一定的條件下）計算過程中可以改變步長等優(yōu)點，但仍需計算f (x, y)在一些點的值，如四階龍格庫塔法每計算一步需要算四次f (x, y)的值，這就給實際計算帶來一定的復(fù)雜性。因此，它與泰勒級數(shù)法一樣，一般用于計算“表頭”。例9.4 用龍格庫塔法解初值問題y = x2 y (0x1) y(0) = 1（9.23）解 ： 取 h = 0.1，由（9.22）得把初始條件x0 = 0，y0 = 1，代入,得k1 = -1，k2 = -0.9475，k3 = -0.9501，k4 = 0.8950，將這些k值代（9.22），

16、得重復(fù)上述步驟可算出y2，y3，y10等。3 線性多步法前面介紹的方法，統(tǒng)稱為單步法，就是在計算yn+1時，只用到前面一步y(tǒng)n的值。本節(jié)我們將介紹利用前邊已經(jīng)算出來的若干個值yn-k，yn-k+1，yn-1，yn，來求得yn+1的高精度公式線性多步方法。我們已經(jīng)知道初值問題y = f (x, y) y (x0) = y0（9.24）與積分方程（9.25）等價。本節(jié)介紹的方法，其基本思想是用一個插值多項式P(x)來代替（9.25）中的被積函數(shù)f (x, y)，然后用（9.26）代替(9.25)。這種方法實際上要分兩步來做：（1）求出開頭幾個點上的近似值，即計算“表頭”；（2）利用(9.26)逐步

17、求后面點xk上的值yk。選取不同的插值點作插值多項式，就會得出不同的數(shù)值解法，我們這里只討論其中的一種，叫阿當(dāng)姆斯（Adams）方法。（一）阿當(dāng)姆斯外推公式 為了簡單起見，我們以k = 2為例，導(dǎo)出阿當(dāng)姆斯外推公式。于是以xn-2，xn-1，xn為節(jié)點作牛頓向后插值多項式P2(x)。（9.27）其中而余項為（9.28）將（9.27）代入（9.26）式并經(jīng)變量替換x = xn + th，便得（9.29）在上述公式中被插值點xnxxn+1不包含在插值節(jié)點所決定的最大區(qū)間（xn-2, xn）內(nèi)，故（9.29）稱阿當(dāng)姆斯外推公式(或稱外插公式)，顯然，它是顯式的且每前進(jìn)一步只計算一次f (x, y)的

18、值即可。公式（9.26）的截斷誤差，可由（9.28）及積分第二中值定理得 上式表明阿當(dāng)姆斯外推法(當(dāng)k = 2時)的截斷誤差的階為0(h3)。同理，對k = 3的情形即可求得外推公式（9.31）而余項為（9.32）公式（9.29）和（9.31）要輸入一個差分表。為此，我們將差分表示成函數(shù)值的和的形式：（9.33）其中是組合數(shù)，即這樣（9.29）可改寫為（9.34）而（9.31）可改寫為（9.35）公式（9.34）和（9.35）是常用的外推公式。例9.5 求y = x + y, y (0) = 1當(dāng)x = 0.1到0.5，步長h = 0.1時的數(shù)值解。解 先用前面講過的方法計算出表頭y0 = 1

19、； y1 = 1.11034；y2 = 1.24281； y3 = 1.39972將上述值代入（9.35）式，計算得：y4 = 1.58364； y5 = 1.79742 （二）阿當(dāng)姆斯內(nèi)插法根據(jù)插值理論知道，節(jié)點的選擇對于精度有直接的影響。同樣次數(shù)的內(nèi)插公式比外插公式更精確。如果（9.26）中的被積函數(shù)是以xn-1，xn，xn+1為插值節(jié)點的內(nèi)插多項式，即（9.36）將（9.36）代入（9.26）便得（9.37）把它改寫成便于在計算機上實現(xiàn)的形式，就有（9.38）上式便是k = 1時的阿當(dāng)姆斯內(nèi)插公式，它的截斷誤差的階為（9.39）同理，對k = 2的情形可導(dǎo)出公式（9.40）或（9.41）

20、而（9.42）這就是常用的阿當(dāng)姆斯內(nèi)插公式。比較外推法和內(nèi)插法的截斷誤差可知，在同樣利用k + 1個已知值時，阿當(dāng)姆斯內(nèi)插法的截斷誤差階比外推法高一階，因而內(nèi)插法更精確。但內(nèi)插法（9.38）及（9.41）是隱式的，需要解方程，通常用迭代法求解。如用外推法（9.34）算出的值作為初始近似，然后相應(yīng)地用內(nèi)插法公式（9.41）進(jìn)行迭代，即（9.43）若將（9.41）記成的形式，容易算出因此，當(dāng)h充分小時，可有即迭代程序（9.43）收斂，并且h越小，收斂越快。當(dāng)h選得適當(dāng)時，一般迭代二、三次即可得到滿足精度要求的yn+1值。同理可導(dǎo)出一般的阿當(dāng)姆斯內(nèi)插公式。為了提高精度，經(jīng)常把阿當(dāng)姆斯外推法與內(nèi)插聯(lián)合

21、起來交替使用。例如第一個方程的精度為0(h4)，用第二個方程迭代一次精度仍達(dá)0(h5)。這樣兩個方程交替使用可達(dá)較好的效果。第一個方程是預(yù)報方程，第二個方程是較正方程。例96 對例9.5用內(nèi)插法求解解 先用前面講過的方法計算出表頭y0 = 1； y1 = 1.11034； y2 = 1.24281再用（9.43）的第一個式子計算出，最后用（9.43）的第二個式子進(jìn)行迭代，得y3 = 1.39972同樣，可算出y4及y5。4 解二階常微分方程邊值問題的差分法考慮常微分方程的邊值問題：（9.44）其中p(x)，q(x)和f (x)均為a, b上給定的函數(shù)，a，b為已知數(shù)。邊值問題的主要特點是在兩個端點上給出了定解條件。在以下討論中我們總假定p(x)、q(x)及f (x)均為a, b上充分光滑的函數(shù)，且q(x)0，這時，邊值問題（9.44）存在連續(xù)可微的解，且唯一。用差分法解邊值問題的主要步驟是：（1）將區(qū)間a, b離散化；（2）在這些節(jié)點上，將導(dǎo)數(shù)差商化，從而把微分方程化