第9章常微分方程初值問題數值解法19.1

引言本章要著重考察的一階方程的初值問題（1.1）（1.2）只要函數適當光滑——譬如關于滿足利普希茨(Lipschitz）條件（1.3）理論上就可以保證初值問題（1.1），（1.2）的解存在并且惟一.2所謂數值解法，就是尋求解在一系列離散節點上的近似值.相鄰兩個節點的間距稱為步長.如不特別說明，總是假定為定數，這時節點為.初值問題（1.1），（1.2）的數值解法的基本特點是采取“步進式”.即求解過程順著節點排列的次序一步一步地向前推進.（1.1）（1.2）3描述這類算法，只要給出用已知信息計算的遞推公式.一類是計算時只用到前一點的值，稱為單步法.另一類是用到前面點的值，稱為步法.其次，要研究公式的局部截斷誤差和階，數值解與精確解的誤差估計及收斂性，還有遞推公式的計算穩定性等問題.首先對方程離散化，建立求數值解的遞推公式.49.2

簡單的數值方法與基本概念

9.2.1

歐拉法與后退歐拉法積分曲線上一點的切線斜率等于函數的值.如果按函數在平面上建立一個方向場，那么，積分曲線上每一點的切線方向均與方向場在該點的方向一致.在平面上，微分方程的解稱作它的積分曲線.5基于上述幾何解釋，從初始點出發，先依方向場在該點的方向推進到上一點，然后再從依方向場的方向推進到上一點，循此前進做出一條折線（圖9-1).圖9-16一般地，設已作出該折線的頂點，過依方向場的方向再推進到，顯然兩個頂點的坐標有關系即（2.1）這就是著名的歐拉（Euler）公式.若初值已知，則依公式（2.1）可逐步算出7

例1（2.2）求解初值問題

解取步長，歐拉公式的具體形式為計算結果見表9-1.初值問題（2.2）的解為，按這個解析式子算出的準確值同近似值一起列在表9-1中，兩者相比較可以看出歐拉方法的精度很差.8還可以通過幾何直觀來考察歐拉方法的精度.假設，即頂點落在積分曲線上，那么，按歐拉方法作出的折線便是過點的切線（圖9-2）.9圖9-2從圖形上看，這樣定出的頂點明顯地偏離了原來的積分曲線，可見歐拉方法是相當粗糙的.為了分析計算公式的精度，通常可用泰勒展開將在處展開，則有10在的前提下，稱為此方法的局部截斷誤差.于是可得歐拉法（2.1）的公式誤差（2.3）（2.4）（2.1）如果對方程從到積分，得11右端積分用左矩形公式近似.再以代替如果在（2.4）中右端積分用右矩形公式（2.5）稱為后退的歐拉法.歐拉公式是關于的一個直接的計算公式，這類公式代替也得到（2.1），局部截斷誤差也是（2.3）.近似，則得另一個公式稱作是顯式的；（2.1）（2.3）（2.4）12公式（2.5）的右端含有未知的，它是關于的一個函數方程，隱式方程通常用迭代法求解，而迭代過程的實質是逐步顯示化.設用歐拉公式給出迭代初值，用它代入（2.5）式的右端，使之轉化為顯式，直接計算得這類公式稱作是隱式的.（2.5）（2.5）13然后再用代入（2.5）式，又有如此反復進行，得（2.6）由于對滿足利普希茨條件（1.3）.由（2.6）減（2.5）得由此可知，只要迭代法（2.6）就收斂到解.（2.5）（2.5）（1.3）14從積分公式可以看到后退歐拉方法的公式誤差與歐拉法是相似的.15

9.2.2

梯形方法若用梯形求積公式近似等式（2.4）右端的積分并分別用代替則可得到比歐拉法精度高的計算公式（2.7）稱為梯形方法.梯形方法是隱式單步法，可用迭代法求解.（2.4）16為了分析迭代過程的收斂性，將（2.7）與（2.8）式相減，得（2.8）同后退的歐拉方法一樣，仍用歐拉方法提供迭代初值，則梯形法的迭代公式為（2.7）17如果選取充分小，使得則當時有這說明迭代過程（2.8）是收斂的.于是有式中為關于的利普希茨常數.18

9.2.3

單步法的局部截斷誤差與階初值問題(1.1)，(1.2)的單步法可用一般形式表示為（2.9）其中多元函數與有關，當含有時，方法是隱式的，若不含則為顯式方法，（2.10）稱為增量函數.所以顯式單步法可表示為例如對歐拉法（2.1）有它的局部截斷誤差已由（2.3）給出.（1.1）（1.2）（2.1）（2.3）19對一般顯式單步法則可如下定義.

定義1設是初值問題(1.1)，(1.2)的準確解，（2.11）為顯式單步法（2.10）的局部截斷誤差.

之所以稱為局部的，是假設在前各步沒有誤差.當時，計算一步，則有稱（1.1）（1.2）（2.10）20在前一步精確的情況下用公式（2.10）計算產生的公式誤差.根據定義，歐拉法的局部截斷誤差即為（2.3）的結果.這里稱為局部截斷誤差主項.局部截斷誤差可理解為用方法（2.10）計算一步的誤差，即顯然（2.10）（2.3）21

定義2設是初值問題(1.1)，(1.2)的準確解，若存在最大整數使顯式單步法(2.10)的局部截斷誤差滿足（2.12）則稱方法（2.10）具有階精度.若將（2.12）展開式寫成則稱為局部截斷誤差主項.以上定義對隱式單步法（2.9）也是適用的.（1.1）（1.2）（2.10）（2.10）（2.9）22對后退歐拉法（2.5）其局部截斷誤差為這里，是1階方法，局部截斷誤差主項為.（2.5）23對梯形法（2.7）有所以梯形方法是2階的，其局部誤差主項為（2.7）24

9.2.4

改進的歐拉公式梯形方法雖然提高了精度，但其算法復雜.在應用迭代公式（2.9）進行實際計算時，每迭代一次，都要重新計算函數的值.為了控制計算量，通常只迭代一兩次就轉入下一步的計算，這就簡化了算法.具體地，先用歐拉公式求得一個初步的近似值,而迭代又要反復進行若干次，計算量很大，而且往往難以預測.稱之為預測值，（2.9）25這樣建立的預測-校正系統通常稱為改進的歐拉公式：預測值的精度可能很差，再用梯形公式（2.7）將它校正一次，即按（2.8）式迭代一次得，這個結果稱校正值.預測校正（2.13）也可以表為下列平均化形式（2.7）（2.8）26

例2

解用改進的歐拉方法求解初值問題（2.2）.（2.2）這里改進的歐拉公式為27仍取，計算結果見表9-2.同例1中歐拉法的計算結果比較，改進歐拉法明顯改善了精度.289.3

龍格-庫塔方法29

9.3.1

顯式龍格-庫塔法的一般形式上節給出了顯式單步法的表達式其局部截斷誤差為對歐拉法，即方法為階.（3.1）若用改進歐拉法，它可表示為30此時增量函數（3.2）與歐拉法的相比，增加了計算一個右函數的值，可望.若要使得到的公式階數更大，就必須包含更多的值.（3.3）從方程等價的積分形式（2.4），即31若要使公式階數提高，就必須使右端積分的數值求積公式精度提高，必然要增加求積節點.為此可將（3.3）的右端用求積公式表示為點數越多，精度越高，上式右端相當于增量函數，為得到便于計算的顯式方法，可類似于改進歐拉法，將公式表示為（3.4）其中（3.3）32（3.5）這里均為常數.

(3.4)和(3.5)稱為級顯式龍格-庫塔（Runge-Kutta）法，簡稱R-K方法.當時，就是歐拉法，此時方法的階為.當時，改進歐拉法(3.1)，(3.2)也是其中的一種.（3.4）33下面將證明階.要使公式(3.4),(3.5)具有更高的階,就要增加點數.下面就推導R-K方法.（3.4）（3.5）34

9.3.2

二階顯式R-K方法

對的R-K方法，計算公式如下（3.6）這里均為待定常數.希望適當選取這些系數，使公式階數盡量高.根據局部截斷誤差的定義，（3.6）的局部截斷誤差為（3.7）35這里.為得到的階，要將上式各項在處做泰勒展開，由于是二元函數，故要用到二元泰勒展開，其中各項展開式為（3.8）36將以上結果代入局部截斷誤差公式則有要使公式（3.6）具有階，必須使（3.6）37即（3.9）的解是不惟一的.令，則得這樣得到的公式稱為二階R-K方法，如取，則這就是改進歐拉法（3.1）.（3.9）38若取則.稱為中點公式，（3.10）也可表示為得計算公式（3.10）相當于數值積分的中矩形公式.的R-K公式(3.6)的局部誤差不可能提高到.（3.6）39把多展開一項，從（3.8）的

看到展開式中的項是不能通過選擇參數消掉的.實際上要使的項為零，需增加3個方程，要確定4個參數，這是不可能的.故的顯式R-K方法的階只能是，而不能得到三階公式.40

9.3.3

三階與四階顯式R-K方法要得到三階顯式R-K方法，必須.（3.11）其中及均為待定參數.此時（3.4）,（3.5）的公式表示為公式（3.11）的局部截斷誤差為（3.4）（3.5）41只要將按二元函數泰勒展開，使，可得待定參數滿足方程（3.12）42這是8個未知數6個方程的方程組，解也不是惟一的.所以這是一簇公式.滿足條件（3.12）的公式（3.11）統稱為三階R-K公式.一個常見的公式為此公式稱為庫塔三階方法.43繼續上述過程，經過較復雜的數學演算，可以導出各種四階龍格-庫塔公式，下列經典公式是其中常用的一個：可以證明其截斷誤差為.四階龍格-庫塔方法的每一步需要計算4次函數值，（3.13）44

例3

解設取步長，從到用四階龍格-庫塔方法求解初值問題這里，經典的四階龍格-庫塔公式為45表9-3列出了計算結果，同時列出了相應的精確解.比較例3和例2的計算結果，顯然龍格-庫塔方法的精度高.46但由于這里放大了步長，所以表9-3和表9-2所耗費的計算量幾乎相同.龍格-庫塔方法的推導基于泰勒展開方法，因而它要求所求的解具有較好的光滑性質.反之，如果解的光滑性差，那么，使用四階龍格-庫塔方法求得的數值解，其精度可能反而不如改進的歐拉方法.四階龍格-庫塔方法每一步要4次計算函數,其計算量比每一步只要2次計算函數的二階龍格-庫塔方法----改進的歐拉方法大一倍.47

9.3.4

變步長的龍格-庫塔方法單從每一步看，步長越小，截斷誤差就越小，但隨著步長的縮小，在一定求解范圍內所要完成的步數就增加了.步數的增加不但引起計算量的增大，而且可能導致舍入誤差的嚴重積累.因此同積分的數值計算一樣，微分方程的數值解法也有個選擇步長的問題.在選擇步長時，需要考慮兩個問題：

1.怎樣衡量和檢驗計算結果的精度？48

2.如何依據所獲得的精度處理步長？考察經典的四階龍格-庫塔公式:（3.13）從節點出發，先以為步長求出一個近似值，49（3.14）然后將步長折半，即取為步長從跨兩步到，再求得一個近似值，每跨一步的截斷誤差是（3.15）比較（3.14）式和（3.15）式我們看到，步長折半后，由于公式的局部截斷誤差為，故有因此有誤差大約減少到.即有50由此易得下列事后估計式這樣，可以通過檢查步長、折半前后兩次計算結果的偏差來判定所選的步長是否合適.具體地說，將區分以下兩種情況處理：51

1.對于給定的精度，如果，反復將步長折半進行計算，直至為止.這時取最終得到的作為結果；

2.如果，反復將步長加倍，直到為止，這種通過加倍或折半處理步長的方法稱為變步長方法.這時再將步長折半一次，就得到所要的結果.表面上看，為了選擇步長，每一步的計算量增加了，但總體考慮往往是合算的.529.4單步法的收斂性與穩定性

9.4.1

收斂性與相容性

數值解法的基本思想是通過某種離散化手段將微分方程轉化為差分方程，（4.1）它在處的解為，而初值問題（1.1）,（1.2）在處的精確解為記，稱為整體截斷誤差.如單步法（2.10），即（1.1）（1.2）53收斂性就是討論當固定且時的問題.

定義3若一種數值方法對于固定的,當時有，其中是(1.1),(1.2)的準確解，則稱該方法是收斂的.顯然數值方法收斂是指.對單步法（4.1）有下述收斂性定理.（1.1）（1.2）（4.1）54

定理1假設單步法（4.1）具有階精度，且增量函數關于滿足利普希茨條件（4.2）又設初值是準確的，即，（4.3）則其整體截斷誤差

注：依據這一定理，判斷單步法（4.1）的收斂性，歸結為驗證增量函數能否滿足利普希茨條件（4.2）且p≥1.對于歐拉方法，由于其增量函數就是，故當關于滿足利普希茨條件且p=1,所以它是收斂的.55再考察改進的歐拉方法，其增量函數由給出，這時有假設關于滿足利普希茨條件，記利普希茨常數為,設為定數），上式表明關于的利普希茨常數則由上式推得又因為p=2,因此改進的歐拉方法也是收斂的.56類似地，也可驗證其他龍格-庫塔方法的收斂性.定理1表明時單步法收斂，并且當是初值問題（1.1），（1.2）的解，（4.1）具有階精度時，有展開式所以的充要條件是，（4.1）（1.1）（1.2）57而，于是可給出如下定義.

定義4若單步法（4.1）的增量函數滿足則稱單步法（4.1）與初值問題（1.1），（1.2）相容.以上討論表明階方法（4.1）當時與(1.1),(1.2)相容，反之相容方法至少是1階的.于是由定理1可知若增量函數方法（4.1）收斂的充分必要條件是此方法是相容的.關于滿足利普希茨條件58

穩定性引例/*Stability*/例：考察初值問題在區間[0,0.5]上的解。分別用歐拉顯、隱式格式和改進的歐拉格式計算數值解。0.00.10.20.30.40.5精確解改進歐拉法

歐拉隱式歐拉顯式

節點xi

1.0000

2.00004.0000

8.00001.6000101

3.2000101

1.00002.5000101

6.25001021.56251023.90631039.76561041.00002.50006.25001.56261013.90631019.76561011.00004.97871022.47881031.23411046.14421063.0590107Whatiswrong??!

AnEngineercomplains:"Maththeoremsaresounstablethatasmallperturbationontheconditionswillcauseacrashontheconclusions!"599.4.2絕對穩定性與絕對穩定域

定義5以歐拉法為例考察計算穩定性.

例4其準確解是一個按指數曲線衰減得很快的函數，若一種數值方法在節點值上大小為的擾動于以后各節點值上產生的偏差均不超過，則稱該方法是穩定的.考察初值問題如圖9-3所示.60若取，則歐拉公式的具體形式為計算結果列于表9-4的第2列.可以看到，歐拉方法的解（圖×號標出）在準確值的上下波動，計算過程明顯地不穩定.圖9-3用歐拉法解方程得61再考察后退的歐拉方法，取時計算公式為計算結果列于表第3列（圖9-3中標以·號），這時計算過程穩定.但若取則計算過程穩定.這表明穩定性不但與方法有關，也與步長h的大小有關，當然也與方程中的f(x,y)有關.62為了只考察數值方法本身，通常只檢驗將數值方法用于解模型方程的穩定性，（4.8）其中為復數.模型方程為模型方程的穩定性歐拉方法的穩定性.模型方程的歐拉公式為（4.9）設在節點值上有一擾動值，它的傳播使節點值產生大小為的擾動值，63假設用按歐拉公式得出的計算過程不再有新的誤差，則擾動值滿足可見擾動值滿足原來的差分方程（4.9）.如果差分方程的解是不增長的，即有則它就是穩定的.顯然，為要保證差分方程（4.9）的解不增長，只要選取充分小，（4.10）在的復平面上，這是以為圓心，1為半徑的單位圓域.即使64到的解，滿足，則稱方法單步法（4.1）用于解模型方程（4.8），若得（4.1）是絕對穩定