數(shù)值分析第9章 常微分方程數(shù)值解

1、許多實(shí)際問題的數(shù)學(xué)模型是微分方程或微分方程的定解問題。如物體運(yùn)動(dòng)、電路振蕩、化學(xué)反映及生物群體的變化等。常微分方程可分為線性、非線性、高階方程與方程組等類；線性方程包含于非線性類中，高階方程可化為一階方程組。若方程組中的所有未知量視作一個(gè)向量，則方程組可寫成向量形式的單個(gè)方程。因此研究一階微分方程的初值問題 (9-1)的數(shù)值解法具有典型性。常微分方程的解能用初等函數(shù)、特殊函數(shù)或它們的級(jí)數(shù)與積分表達(dá)的很少。用解析方法只能求出線性常系數(shù)等特殊類型的方程的解。對(duì)非線性方程來說，解析方法一般是無能為力的，即使某些解具有解析表達(dá)式，這個(gè)表達(dá)式也可能非常復(fù)雜而不便計(jì)算。因此研究微分方程的數(shù)值解法是非常必要

2、的。只有保證問題(9-1)的解存在唯一的前提下，研究其數(shù)值解法或者說尋求其數(shù)值解才有意義。由常微分方程的理論知，如果(9-1)中的滿足條件（1）在區(qū)域上連續(xù)；（2）在上關(guān)于滿足Lipschitz條件，即存在常數(shù)，使得則初值問題(9-1)在區(qū)間上存在惟一的連續(xù)解。在下面的討論中，我們總假定方程滿足以上兩個(gè)條件。所謂數(shù)值解法，就是求問題(9-1)的解在若干點(diǎn)處的近似值的方法。稱為問題(9-1)的數(shù)值解，稱為由到的步長(zhǎng)。今后如無特別說明，我們總假定步長(zhǎng)為常量。建立數(shù)值解法，首先要將微分方程離散化，一般采用以下幾種方法：(1) 用差商近似導(dǎo)數(shù)在問題(9-1)中，若用向前差商代替，則得用其近似值代替，所

3、得結(jié)果作為的近似值，記為，則有這樣，問題(9-1)的近似解可通過求解下述問題(9-2)得到，按式(9-2)由初值經(jīng)過步迭代，可逐次算出。此方程稱為差分方程。需要說明的是，用不同的差商近似導(dǎo)數(shù)，將得到不同的計(jì)算公式。(2) 用數(shù)值積分法將問題(9-1)中的微分方程在區(qū)間上兩邊積分，可得(9-3)用,分別代替,，若對(duì)右端積分采用取左端點(diǎn)的矩形公式，即同樣可得出顯式公式(9-2)。類似地，對(duì)右端積分采用其它數(shù)值積分方法，又可得到不同的計(jì)算公式。(3) 用Taylor多項(xiàng)式近似。把在點(diǎn)處Taylor展開，取一次多項(xiàng)式近似，則得設(shè)，略去余項(xiàng)，并以代替，便得以上三種方法都是將微分方程離散化的常用方法，每一

4、類方法又可導(dǎo)出不同形式的計(jì)算公式。其中Taylor展開法，不僅可以得到求數(shù)值解的公式，而且容易估計(jì)截?cái)嗾`差。上面我們給出了求解初值問題(9-1)的一種最簡(jiǎn)單的數(shù)值公式(9-2)。雖然它的精度比較低，實(shí)踐中很少采用，但它的導(dǎo)出過程能較清楚地說明構(gòu)造數(shù)值解公式的基本思想，且?guī)缀我饬x明確，因此它在理論上仍占有一定的地位。1簡(jiǎn)單的數(shù)值方法和基本概念1.1 Euler法與向后Euler法一、Euler法Euler方法就是用差分方程初值問題(9-4)的解來近似微分方程初值問題(9-1)的解，即由公式(9-4)依次算出的近似值。從幾何上看，微分方程在平面上確定了一個(gè)向量場(chǎng)：點(diǎn)處的方向斜率為。問題(9-1)的

5、解代表一條過點(diǎn)的曲線，稱為積分曲線，且此曲線上每點(diǎn)的切向都與向量場(chǎng)在這點(diǎn)的方向一致。從點(diǎn)出發(fā)，以為斜率作一直線段，與直線交于點(diǎn)，顯然有，再?gòu)某霭l(fā)，以為斜率作直線段推進(jìn)到上一點(diǎn)，其余類推，這樣得到解曲線的一條近似曲線，它就是折線。因此Euler方法又稱為Euler折線法。二、向后Euler法在微分方程離散化時(shí)，用向后差商代替導(dǎo)數(shù)，即，則得到如下差分方程(9-5)用這組公式求問題(9-1)的數(shù)值解稱為向后Euler法。向后Euler法與Euler法形式上相似，但實(shí)際計(jì)算時(shí)卻復(fù)雜得多。Euler法計(jì)算的公式中不含有，這樣的公式稱為顯式公式；向后Euler法計(jì)算的公式中含有，稱為隱式公式。顯式公式與隱

6、式公式各有特點(diǎn)。顯式公式的優(yōu)點(diǎn)是使用方便，計(jì)算簡(jiǎn)單，效率高。其缺點(diǎn)是計(jì)算精度低，穩(wěn)定性差；隱式公式正好與它相反，它具有計(jì)算精度高，穩(wěn)定性好等優(yōu)點(diǎn)，但求解過程很復(fù)雜，一般采用迭代法。為了結(jié)合各自的優(yōu)點(diǎn)，通常將顯式公式與隱式公式配合使用，由顯式公式提供迭代初值，再經(jīng)隱式公式迭代校正。上面隱式公式中，在求解時(shí)，為已知，是方程的根。一般說來，這是一個(gè)非線性方程，因此我們通過構(gòu)造簡(jiǎn)單迭代法來求解。迭代格式為由于滿足Lipschitz條件，所以由此可知，只要，迭代法就收斂到解。1.2 梯形公式利用數(shù)值積分方法將微分方程離散化時(shí)，若用梯形公式計(jì)算式(9-3)中右端積分，即并用代替，則得計(jì)算公式(9-6)這就

7、是求解初值問題(9-1)的梯形公式。梯形公式也是隱式格式，一般需用迭代法求解，迭代公式為由于函數(shù)關(guān)于滿足Lipschitz條件，所以其中為L(zhǎng)ipschitz常數(shù)。因此，當(dāng)時(shí)，迭代法就收斂到解。1.3 局部截?cái)嗾`差與方法的精度為了刻畫近似解的準(zhǔn)確程度，引入局部截?cái)嗾`差與方法精度的概念。定義 假設(shè)在某一步的近似解是準(zhǔn)確的，即（這個(gè)假設(shè)稱為局部化假設(shè)）。在此前提下，用某公式推算所得，我們稱為該公式（即該方法）的局部截?cái)嗾`差。定義9.2 如果某種方法的局部截?cái)嗾`差是則稱該方法是階方法，或具有階精度。顯然越大，方法的精度越高。1）Euler法的截?cái)嗾`差假設(shè)問題的解充分光滑，且前步計(jì)算結(jié)果是精確的，即，于

8、是Euler法的截?cái)嗾`差是(9-7)這里稱為局部截?cái)嗾`差主項(xiàng)。顯然。2）向后Euler法的截?cái)嗾`差。計(jì)算公式是將對(duì)用微分中值定理，有（在與之間）將在處Taylor展開于是將方程的解作Taylor展開因此故(9-8)3）梯形法的計(jì)算公式是將對(duì)用微分中值定理，有（在與之間）將在處Taylor展開于是將方程的解作Taylor展開因此故(9-9)所以梯形法是二階方法。1.4 改進(jìn)的Euler法我們看到，雖然梯形方法提高了精度，但其算法復(fù)雜，在應(yīng)用迭代公式進(jìn)行實(shí)際計(jì)算時(shí)，每迭代一次，都要重新計(jì)算函數(shù)的值，而迭代又要反復(fù)進(jìn)行若干次，計(jì)算量很大。為了控制計(jì)算量，通常只迭代一兩次就轉(zhuǎn)入下一步的計(jì)算，這就簡(jiǎn)化了

9、算法。具體地說，我們先用歐拉公式求得一個(gè)初步的近似值，稱之為預(yù)測(cè)值，預(yù)測(cè)值的精度可能很差，再用梯形公式將它校正一次得，稱為校正值。這樣的預(yù)測(cè)校正系統(tǒng)通常稱為改進(jìn)的歐拉公式。即預(yù)測(cè)：校正：為了編制程序上機(jī)，將上式改寫成(9-10)算法(1) 輸入，整數(shù)，初值(2) 置 輸出(3) 計(jì)算 輸出(4) 若，置，轉(zhuǎn)3；否則停機(jī)。改進(jìn)歐拉法的截?cái)嗾`差。計(jì)算公式是在處作Taylor展開式，注意到，有于是(9-11)所以改進(jìn)Euler法是二階方法。2 龍格-庫(kù)塔（Runge-Kutta）法2.1 Runge-Kutta法的基本思想及一般形式設(shè)初值問題(9-1)的解，按微分中值定理，必存在，使(9-12)圖9

10、-3其中叫做在上的平均斜率。只要對(duì)平均斜率提供一種算法，公式(9-12)就給出一種數(shù)值解公式。例如，用代替，就得到前進(jìn)Euler公式；用替代，就得到后退Euler公式；如果用的算術(shù)平均值替代，則可得到2階精度的梯形公式，如圖9-3?？梢栽O(shè)想，如果在上能多預(yù)測(cè)幾個(gè)點(diǎn)的斜率值，用它們的加權(quán)平均值替代，就有望得到具有較高精度的數(shù)值解公式，這就是Runge-Kutta法的基本思想。Runge-Kutta公式的一般形式是(9-13)其中是在點(diǎn)的斜率預(yù)測(cè)值。均為常數(shù)。選取這些常數(shù)的原則是使公式(9-13)具有盡可能高的精度。公式(9-13)叫做級(jí)顯式Runge-Kutta公式，簡(jiǎn)稱RK方法。2.2 二階R

11、K法的推導(dǎo)的RK方法的近似公式為(9-14)這里均為待定常數(shù)，我們希望適當(dāng)選取這些系數(shù)，使公式階數(shù)盡量高。先展開，按照二元函數(shù)Taylor級(jí)數(shù)：，得(9-15)為了敘述方便，把及其偏導(dǎo)數(shù)中的省略不寫，將式(9-15)代入(9-14)的第一式得再展開，注意到公式得于是，局部截?cái)嗾`差為(9-16)要使式(9-16)的局部截?cái)嗾`差為，則應(yīng)要求的系數(shù)為零，于是有(9-17)方程有4個(gè)未知數(shù)，3個(gè)方程，所以有無窮多組解，它的每組解代入式(9-14)得到的近似公式，局部截?cái)嗾`差均為，故這些方法統(tǒng)稱為二階方法。例如，取，得，近似公式為(9-18)這就是改進(jìn)Euler公式。如果取，有，近似公式為(9-19)這

12、也是常用的二階公式，稱為中點(diǎn)公式。注：對(duì)于一般函數(shù)，由于，所以不論參數(shù)如何選取，只能有。這說明(9-14)至多是二階的方法。類似地，對(duì)和的情形，通過更復(fù)雜的計(jì)算，可以導(dǎo)出三階和四階RK公式，其中常用的三階和四階RK公式為(9-20)和(9-21)式(9-21)稱為經(jīng)典的四階RK公式，通常說四階RK方法就是指用式(9-21)求解。算法(1) 輸入(2) 置(3) 計(jì)算輸出(4) 若，則置轉(zhuǎn)3；否則停機(jī)。例4 用Euler法（），改進(jìn)Euler法（）和經(jīng)典R-K方法（）計(jì)算初值問題.其結(jié)果如下表9-1表9-1Euler方法()改進(jìn)Euler法()經(jīng)典R-K法()00000注：1）通過比較RK四階經(jīng)

13、典公式和前面用Euler法、改進(jìn)Euler法（它是一種二階RK方法）的計(jì)算結(jié)果，我們發(fā)現(xiàn)，四階RK方法的精度高得多。就計(jì)算量來說，雖然Euler法、改進(jìn)Euler法每步只需計(jì)算一個(gè)或二個(gè)函數(shù)值，而四階RK方法每步需計(jì)算四個(gè)函數(shù)值，但由于放大了步長(zhǎng)，計(jì)算量幾乎相同。可見，選擇有效的算法是非常重要的。2）RK方法的導(dǎo)出基于Taylor展開，故它要求問題的解具有較高的光滑度。當(dāng)解充分光滑時(shí)，四階RK方法確實(shí)優(yōu)于改進(jìn)Euler法。如果解的光滑性較差，則用四階RK方法求得的數(shù)值解，其精度可能反而比改進(jìn)Euler法的差。因此，在實(shí)際計(jì)算時(shí)，需根據(jù)問題的具體情況選擇合適的算法。3）當(dāng)時(shí)，RK公式的最高階數(shù)恰

14、是，當(dāng)時(shí)，RK公式的最高階數(shù)不是，如時(shí)，RK公式的最高階數(shù)仍為4，時(shí)，RK公式的最高階數(shù)仍為5。由此看來，四階RK公式是最經(jīng)濟(jì)的。2.3 變步長(zhǎng)的RK方法若問題(9-1)的解函數(shù)變化是不均勻的，用等步長(zhǎng)求數(shù)值解，則可能產(chǎn)生有些點(diǎn)處精度過高，有些點(diǎn)處精度過低的情況，為保證一定的精度，必須取較小的步長(zhǎng)。這樣做既增加了計(jì)算量，也導(dǎo)致誤差的嚴(yán)重積累，因而實(shí)際計(jì)算時(shí)，往往采用事后估計(jì)誤差、自動(dòng)調(diào)整步長(zhǎng)的RK方法。設(shè)用階方法計(jì)算，記從出發(fā)，以步長(zhǎng)計(jì)算一步所得的近似值，以步長(zhǎng)計(jì)算兩步得的近似值。由于階公式的局部截?cái)嗾`差為，且在小區(qū)間上變化不大，故有(9-22)(9-23)將式(9-24)乘以減(9-23)得

15、從而有(9-24)這樣，可以事后誤差估計(jì)(9-24)中的大小來選擇合適的步長(zhǎng)。變步長(zhǎng)的RK方法與數(shù)值積分做法相同，可以從的值來選擇合適的步長(zhǎng)，從而得到合乎精度要求的。具體作法是：設(shè)要求精度是，如果，就反復(fù)加倍步長(zhǎng)進(jìn)行計(jì)算，直到為止。這時(shí)上一次的步長(zhǎng)就是合適的步長(zhǎng)，上一次計(jì)算所得的結(jié)果，就是合乎精度要求的；如果，則反復(fù)折半步長(zhǎng)，直到為止。這時(shí)，最后一次步長(zhǎng)就是合適的步長(zhǎng)，而這最后一次的計(jì)算結(jié)果就是滿足精度要求的。這種計(jì)算過程中自動(dòng)選擇步長(zhǎng)的方法，叫變步長(zhǎng)方法。以上所介紹的各種數(shù)值解法都是單步法，這時(shí)因?yàn)樗鼈冊(cè)谟?jì)算時(shí)，都只用到前一步的值，單步法的一般形式是 (9-25)其中為增量函數(shù)，例如Eule

16、r方法的增量函數(shù)為，改進(jìn)Euler法的增量函數(shù)為3 相容性、收斂性與穩(wěn)定性3.1 相容性與收斂性上面所介紹的方法都是用離散化的手段，將微分方程初值問題化為差分方程初值問題求解的。這些轉(zhuǎn)化是否合理？即當(dāng)時(shí)，差分方程是否能無限逼近微分方程，差分方程的解是否能無限逼近微分方程初值問題的準(zhǔn)確解，這就是相容性與收斂性問題。用單步法(9-25)求解初值問題(9-1)，即用差分方程初值問題 (9-26)的解作為問題(9-1)的近似解，如果近似是合理的，則應(yīng)有 (9-27)（體會(huì)：要使得差分方程能無限逼近微分方程，我們將微分方程的精確解代人差分方程中，當(dāng)時(shí)，應(yīng)該滿足）其中為問題(9-1)的精確解。因?yàn)楣视?9

17、-27)得如果增量函數(shù)關(guān)于連續(xù)，則有 (9-28)如果單步法的增量函數(shù)滿足條件(9-28)，則稱單步法(9-25)與初值問題(9-1)相容。通常稱(9-28)為單步法的相容條件。滿足相容條件(9-28)是可以用單步法求解初值問題(9-1)的必要條件。容易驗(yàn)證Euler法和改進(jìn)Euler法均滿足相容性條件。事實(shí)上，對(duì)Euler法，增量函數(shù)為自然滿足條件(9-28)，改進(jìn)Euler法的增量函數(shù)為因?yàn)檫B續(xù)，從而有所以Euler法和改進(jìn)Euler法均與初值問題(9-1)相容。一般地，如果單步法有階精度（），則其截?cái)嗾`差為上式兩端同除以，得令，如果連續(xù)，則有所以的單步法均與問題(9-1)相容。由此即得各

18、階RK方法與初值問題(9-1)相容。定義9.4 一種數(shù)值方法稱為是收斂的，如果對(duì)于任意初值及任意固定的，都有其中為初值問題(9-1)的精確解。如果我們?nèi)∠植炕俣?，使用某單步法公式，從出發(fā)，一步一步地推算到處的近似值。若不計(jì)各步的舍人誤差，而每步都有局部截?cái)嗾`差，這些局部截?cái)嗟姆e累就是整體截?cái)嗾`差。定義9.5我們稱為某數(shù)值方法的整體誤差。其中為初值問題(9-1)的精確解，為不計(jì)舍人誤差時(shí)用某數(shù)值方法從開始，逐步得到的在處的近似值（不考慮舍人誤差的情況下，截?cái)嗾`差的積累）。定理設(shè)單步法(9-25)具有階精度，其增量函數(shù)關(guān)于滿足Lipschitz條件，問題(9-1)的初值是精確的，即，則單步法的

19、整體截?cái)嗾`差為證明 由已知，關(guān)于滿足Lipschitz條件，故存在，使得對(duì)任意的及，都有記，因?yàn)閱尾椒ň哂须A精度，故存在，使得從而有反復(fù)遞推得因?yàn)椋?，又，于是所以推? 設(shè)單步法具有()階精度，增量函數(shù)在區(qū)域：上連續(xù)，且關(guān)于滿足Lipschitz條件，則單步法是收斂的。當(dāng)在區(qū)域上連續(xù)，且關(guān)于滿足Lipschitz條件時(shí)，改進(jìn)Euler方法，各階RK方法的增量函數(shù)在區(qū)域上連續(xù)，且關(guān)于滿足Lipschitz條件，因而它們都是收斂的。關(guān)于單步法收斂的一般結(jié)果是：定理設(shè)增量函數(shù)在區(qū)域上連續(xù)，且關(guān)于滿足Lipschitz條件，則單步法收斂的充分必要條件是相容性條件(9-28)。3.2 穩(wěn)定性穩(wěn)定性與收

20、斂性是兩個(gè)不同的概念，收斂性是在假定每一步計(jì)算都準(zhǔn)確的前提下，討論當(dāng)步長(zhǎng)時(shí)，方法的整體截?cái)嗾`差是否趨于零的問題。而穩(wěn)定性則是討論舍人誤差的積累能否對(duì)計(jì)算結(jié)果有嚴(yán)重影響的問題。 若一種數(shù)值方法在節(jié)點(diǎn)值上有一個(gè)大小為的擾動(dòng)，于以后各節(jié)點(diǎn)上產(chǎn)生的偏差均不超過，則稱該方法是穩(wěn)定的。我們以Euler方法為例進(jìn)行討論。假設(shè)由于舍人誤差，實(shí)際得到的不是而是，其中是誤差。由此再計(jì)算一步，得到把它與不考慮舍人誤差的Euler計(jì)算公式相減，并記，就有其中。如果滿足條件， (9-29)則從到的計(jì)算，誤差是不增的，可以認(rèn)為計(jì)算是穩(wěn)定的。如果條件(9-29)不滿足，則每步誤差將增大。當(dāng)時(shí)，顯然條件(9-29)不可能滿足

21、，我們認(rèn)為問題本身具有先天的不穩(wěn)定性。當(dāng)時(shí)，為了滿足收斂性要求(9-29)，有時(shí)要很小，這樣步數(shù)就相當(dāng)多，這時(shí)誤差的累積可能是十分嚴(yán)重的，出現(xiàn)數(shù)值不穩(wěn)定現(xiàn)象?，F(xiàn)在我們考慮的初值問題 (9-30)其解是。今設(shè)在的初值有誤差，實(shí)際求解的問題是 (9-31)它的準(zhǔn)確解是。因此如果用很精確的方法，求出(9-31)的解，則它和(9-30)的解在處誤差為，而在處的誤差則是。它隨著的增大而增大，與所選的數(shù)值方法無關(guān)，是問題本身固有的特性。再看一個(gè)的例子，初值問題為其準(zhǔn)確解是。用歐拉法求解有 (9-32)若取，則歐拉公式的具體形式為同樣討論有初始誤差，則在的誤差是。數(shù)值不穩(wěn)定。式(9-32)中取，則歐拉公式的

22、具體形式為。對(duì)初始誤差，在的誤差是。數(shù)值穩(wěn)定。以上討論表明穩(wěn)定性不但與方法有關(guān)，也與步長(zhǎng)的大小有關(guān)，當(dāng)然也與方程中的有關(guān)。為了只考察數(shù)值方法本身，一般只檢驗(yàn)數(shù)值方法用于求解模型方程的穩(wěn)定性，模型方程為 (9-33)其中為復(fù)數(shù)。這個(gè)方程的分析比較簡(jiǎn)單，一般的方程可以通過局部線性化為這種形式，例如在的鄰域略去高階項(xiàng)，再作變量替換就得到的形式。事實(shí)上，方程可簡(jiǎn)寫為，作變換即可得。對(duì)于個(gè)方程的方程組，可以線性化為，其中為的Jacobi矩陣。若有個(gè)不同的特征值，則可對(duì)角化為，再作變換，得到個(gè)非耦合的方程，其中可以復(fù)數(shù)，所以一般討論(9-33)中的為復(fù)數(shù)。對(duì)于方程(9-33)，若，類似以上分析，認(rèn)為方程是

23、不穩(wěn)定的。所以我們考慮的情形，這時(shí)不同的數(shù)值方法可能是數(shù)值穩(wěn)定或不穩(wěn)定的。當(dāng)一個(gè)單步法用于試驗(yàn)方程，從計(jì)算一步得 (9-34)其中依賴于所選的方法。因?yàn)橥ㄟ^點(diǎn)試驗(yàn)方程的解曲線（它滿足）為，而一個(gè)階單步法的局部截?cái)嗾`差在時(shí)有，所以有 (9-35)這樣可以看出是的一個(gè)近似值。由(9-34)可以看到，若計(jì)算中有誤差，則計(jì)算時(shí)將產(chǎn)生誤差，所以有下面定義。定義9.6 如果(9-34)式中，則稱單步法(9-25)是絕對(duì)穩(wěn)定的。在復(fù)平面上復(fù)變量滿足的區(qū)域，稱為方法(9-25)的絕對(duì)穩(wěn)定區(qū)域，它與實(shí)軸的交稱為絕對(duì)穩(wěn)定區(qū)間。在上述定義中，規(guī)定嚴(yán)格不等式成立，是為了和線性多步法的絕對(duì)穩(wěn)定性定義一致。事實(shí)上，時(shí)也可

24、以認(rèn)為誤差不增長(zhǎng)。(1) Euler方法的穩(wěn)定性Euler方法用于模型方程(9-33)，得，所以有。所以絕對(duì)穩(wěn)定條件是，它的絕對(duì)穩(wěn)定區(qū)域是復(fù)平面上以為中心的單位圓。而為實(shí)數(shù)時(shí)，絕對(duì)穩(wěn)定區(qū)間是。(2) 梯形公式的穩(wěn)定性對(duì)模型方程，梯形公式的具體表達(dá)式，即，所以梯形公式的絕對(duì)穩(wěn)定條件為?；?jiǎn)得，因此梯形公式的絕對(duì)穩(wěn)定區(qū)域?yàn)槠矫娴淖笃矫妗Ｌ貏e地，當(dāng)為負(fù)實(shí)數(shù)時(shí)，對(duì)任意的，梯形公式都是穩(wěn)定的。圖9-4 圖9-5(3) RK方法的穩(wěn)定性與前面的討論相仿，將RK方法用于模型方程(9-33)，可得二、三、四階RK方法的絕對(duì)穩(wěn)定區(qū)域分別為當(dāng)為實(shí)數(shù)時(shí)，三、四階顯式RK方法的絕對(duì)穩(wěn)定區(qū)域分別為，。表9-2 常用幾個(gè)

25、公式的絕對(duì)穩(wěn)定區(qū)間方法階數(shù)區(qū)間Euler公式（顯式）1后退Euler公式（隱式）1梯形公式（隱式）2二級(jí)R-K公式（顯式）2二級(jí)R-K公式（隱式）2三級(jí)R-K公式（顯式）3經(jīng)典R-K公式（顯式）4例5 分別取和0.2，用經(jīng)典R-K方法計(jì)算解 本例，分別為和，前者在絕對(duì)穩(wěn)定區(qū)間內(nèi)，后者則不在，這里列出計(jì)算誤差表9-3115251256253125從表9-3可以看出，當(dāng)時(shí)，雖然計(jì)算結(jié)果的相對(duì)誤差很大，但計(jì)算過程是穩(wěn)定的。而當(dāng)時(shí)，計(jì)算過程不穩(wěn)定。這是因?yàn)闀r(shí)，屬于穩(wěn)定區(qū)間，當(dāng)時(shí)，不屬于穩(wěn)定區(qū)間。對(duì)于一般方程，在考慮穩(wěn)定性對(duì)步長(zhǎng)的限制時(shí)，常用代替（因?yàn)榍懊孀儞Q中），只要步長(zhǎng)的選取使得在所用方法的絕對(duì)穩(wěn)定

26、區(qū)域內(nèi)即可。例如若用Euler方法求解初值問題因?yàn)榍?，要使，只? 線性多步法前面所講的各種數(shù)值解法，都是單步法。因?yàn)樗鼈冊(cè)谟?jì)算時(shí)，都僅僅用到前面的一個(gè)信息；如果在計(jì)算值時(shí)，能夠比較充分地利用前面的已知信息，如，那么就可望使所得到的更加精確。這就是多步法的基本思想。多步法中最常用的是線性多步法，其一般公式是 (9-36)其中均為常數(shù)，。當(dāng)時(shí)，上式稱為顯式，否則稱為隱式。線性多步法與Runge-Kutta法都是高精度方法，不同的是線性多步法是利用前面已求得的在各節(jié)點(diǎn)處的近似值及導(dǎo)數(shù)近似值的線性組合來逼近的，而Runge-Kutta法則是用在內(nèi)若干點(diǎn)處的一階導(dǎo)數(shù)預(yù)測(cè)值的線性組合，來逼近在區(qū)間內(nèi)的平

27、均斜率，從而求得的近似值。構(gòu)造線性多步法的公式有多種途徑，最常用的是Taylor展開方法和數(shù)值積分的方法.定義設(shè)是初值問題(9-1)的精確解，線性多步法(9-36)在上的局部截?cái)嗾`差為若，則稱方法(9-36)是階的，則稱方法(9-36)與問題(9-1)的第一個(gè)方程相容的。4.1 線性多步法公式的推導(dǎo)利用Taylor展開導(dǎo)出線性多步公式(9-36)的基本方法是：將線性多步公式(9-36)在處進(jìn)行Taylor展開，然后與在處的Taylor展開式相比較，要求它們前面的項(xiàng)重合，由此確定參數(shù)。設(shè)初值問題(9-1)的解充分光滑，記，把它在處展成Taylor級(jí)數(shù)則有 (9-37)及 (9-38)用替代式(9

28、-37)、(9-38)中的，則得及把這兩個(gè)式子代入式(9-36)，得(9-39)為了使式(9-36)具有階精度，只須使式(9-39)的前項(xiàng)與在處的Taylor展開式(9-40)的前項(xiàng)系數(shù)對(duì)應(yīng)相等即可。對(duì)比關(guān)于的同次項(xiàng)系數(shù)，得到確定的方程組： (9-41)可見，只要式(9-36)的系數(shù)滿足式(9-41)，方法(9-36)就具有階精度。式(9-40)減去(9-39)，可得局部截?cái)嗾`差為 (9-42)4.2常用的線性多步法公式式(9-41)共有個(gè)方程，個(gè)待定常數(shù)：，只要，即，就可以由式(9-41)，定出具有階精度的公式(9-36)的系數(shù)。下面介紹幾個(gè)常用的線性多步法公式。一、阿達(dá)姆斯（Adams）公

29、式取，有(9-43)此時(shí)方程有9個(gè)未知數(shù)，5個(gè)方程，有四個(gè)自由變量。特別取，可解得相應(yīng)的線性多步法公式為(9-44)因?yàn)?，?9-44)稱為Adams顯式公式，它是四階公式，局部截?cái)嗾`差為 (9-45)如果令，由式(9-43)可得解相應(yīng)的線性多步法公式為(9-46)因?yàn)?，?9-46)稱為Adams隱式公式，它是四階公式，局部截?cái)嗾`差為 (9-47)二、米爾恩（Milne）公式對(duì)方程組(9-42)，令，解出相應(yīng)的線性多步法公式 (9-48)稱為Milne公式，其局部截?cái)嗾`差為 (9-49)Milne公式是四階四步顯式公式。三、海明（Hamming）公式對(duì)方程組(9-41)，取，并令，可得Ham

30、ming公式 (9-50)其局部截?cái)嗾`差為 (9-51)Hamming公式是四階三步隱式公式。由于線性多步法公式(9-36)，只有在知道了前面的之后，才能使用。所以開頭的的值必須先用同階的單步法（Taylor級(jí)數(shù)法或Runge-Kutta法）求出，然后才能用線性多步法公式（例如，若使用四階線性多步法公式，必須用同階單步法，算出）。例6 分別用四階Adams顯式和隱式公式求初值問題的數(shù)值解，取。解 根據(jù)題意，由四階Adams顯式公式(9-43)有由四階Adams隱式公式(9-44)有由上式可解出利用精確解求出起步值（一般可用四階RK公式計(jì)算起步值）后，按上面的公式計(jì)算，結(jié)果見下表表9.4Adam

31、s顯式法Adams隱式法00從表9.4可以看出，四階Adams隱式法比顯式法的精度高，比較這兩種方法的局部截?cái)嗾`差公式(9-45)與(9-47)，可以說明這一現(xiàn)象。一般地，同階的隱式法比顯式法精確，而且數(shù)值穩(wěn)定性也好。但在隱式公式中，通常很難用迭代法求解，這樣又增加了計(jì)算量。因此實(shí)際計(jì)算時(shí)，很少單獨(dú)使用顯式公式或隱式公式，而是將它們聯(lián)合使用：先用顯式公式求出的預(yù)測(cè)值，記作，再用隱式公式對(duì)預(yù)測(cè)值進(jìn)行校正，求出的近似值。4.3 預(yù)測(cè)-校正系統(tǒng)一般地，采用同階的顯式公式與隱式公式配對(duì)使用，即把由顯式求出的（記）作為的預(yù)測(cè)值，然后再代入隱式公式進(jìn)行校正，求出更接近的值。這樣就構(gòu)成了預(yù)測(cè)-校正系統(tǒng)。采用的預(yù)測(cè)校正系統(tǒng)有兩種：Adams顯式-隱式公式 (9-52)Milne-Hamming預(yù)測(cè)校正系統(tǒng) (9-53)說明：(1) 以上兩種預(yù)測(cè)校正公式均為四階公式，其起步值通常用四階RK公式計(jì)算。(2) 有時(shí)為提高精度，校正公式可迭代進(jìn)行多次，但迭代次數(shù)一般超過3次。為減少一次迭代所產(chǎn)生的誤差，常常用局部截?cái)嗾`差進(jìn)一步修正預(yù)測(cè)值與校正值