常微分方程初值問題數(shù)值解法

常微分方程初值問題數(shù)值解法9.1引言9.2簡單的數(shù)值方法與基本概念9.3龍格－庫塔方法9.4單步法的收斂性與穩(wěn)定性9.5線性多步法9.6方程組和高階方程9.1

引言

本章討論一階常微分方程的初值問題：只要函數(shù)適當(dāng)光滑—如滿足利普希茨條件：理論上就能保證初值問題的解存在并且唯一。所謂數(shù)值解法，就是尋求解在一系列離散點上的近似值，相鄰兩個點間的距離稱為步長。一般情況下我們?nèi)槌?shù)，這是節(jié)點為：9.1

引言

初值問題的求解有一個基本特點，它們都是采取“步進式”求解的，即，一步一步地求函數(shù)的值。求解的主要方法是：先對方程進行離散化，建立求數(shù)值解的遞推公式，一類在計算時只用到前面一步的，稱為單步法。另一類在計算時除了用到理論上就能保證初值問題的解存在并且唯一。

所謂數(shù)值解法，就是尋求解在一系列離散點上的近似值，相鄰兩個點間的距離稱為步長。一般情況下我們?nèi)槌?shù)，這是節(jié)點為：

初值問題的求解有一個基本特點，它們都是采取“步進式”求解的，即，一步一步地求函數(shù)的值。求解的主要方法是：先對方程進行離散化，建立求數(shù)值解的遞推公式，一類在計算時只用到前面一步的，稱為單步法。另一類在計算時除了用到前利用前面一步的，還要用到前面的等，這種方稱為多步法。其次，要研究公式的局部截斷誤差和階。數(shù)值解和精確解的誤差估計和收斂性，還有遞推迭代公式的數(shù)值穩(wěn)定性問題。顯式歐拉法在平面上，微分方程初值問題的解稱作方程的積分曲線。積分曲線上每一點的斜率等于該點的函數(shù)值。如果按函數(shù)在平面上建立一個方向場，那么，積分曲線上每一點的切線方向均與方向場在該點的方向一致。方稱為多步法。其次，要研究公式的局部截斷誤差和階。數(shù)值解和精確解的誤差估計和收斂性，還有遞推迭代公式的數(shù)值穩(wěn)定性問題。9.2

簡單的數(shù)值方法與基本概念顯式歐拉法在平面上，微分方程初值問題的解稱作方程的積分曲線。積分曲線上每一點的斜率等于該點的函數(shù)值。如果按函數(shù)在平面上建立一個方向場，那么，積分曲線上每一點的切線方向均與方向場在該點的方向一致?；谏厦娴膸缀谓忉專覀儚某跏键c出發(fā)，先依方向場在該點的方向推進到上一點，然后再從依方向場的方向推進到上一點，循此前進做出一條折線。9.2

簡單的數(shù)值方法與基本概念

一般地，設(shè)已做出該折線上的頂點，過依方向場的方向再推進到，顯然兩個頂點的坐標(biāo)有如下關(guān)系即

基于上面的幾何解釋，我們從初始點出發(fā)，先依方向場在該點的方向推進到上一點，然后再從依方向場的方向推進到上一點，循此前進做出一條折線。

一般地，設(shè)已做出該折線上的頂點，過依方向場的方向再推進到，顯然兩個頂點的坐標(biāo)有如下關(guān)系即這就是著名的歐拉公式。若初值已知，則依該公式可逐步算出：

例1求解初值問題：(其解為)解：根據(jù)歐拉方法，得到：這就是著名的歐拉公式。若初值已知，則依該公式可逐步算出：

例1求解初值問題：(其解為)解：根據(jù)歐拉方法，得到：取步長，計算得到如下結(jié)果：0.11.10001.09540.61.50901.48320.21.19181.18320.71.58031.54920.31.27741.26490.81.64981.61250.41.35821.34160.91.71781.67330.51.43151.41421.01.78481.7321隱式歐拉法在前面的討論中，近似計算公式可以看成是由在區(qū)間上積分得到，而右邊的積分是利用左矩形公式近似，再以代替得到，現(xiàn)在右端的積分用右矩形公式，則得到：取步長，計算得到如下結(jié)果：0.11.10001.09540.61.50901.48320.21.19181.18320.71.58031.54920.31.27741.26490.81.64981.61250.41.35821.34160.91.71781.67330.51.43151.41421.01.78481.7321隱式歐拉法在前面的討論中，近似計算公式可以看成是由在區(qū)間上積分得到，而右邊的積分是利用左矩形公式近似，再以代替得到，現(xiàn)在右端的積分用右矩形公式，則得到：此式的右端包含，是一種隱式的單步法，稱為隱式歐拉方法。利用此法，每一步都要把該式作為的方程來求解。從數(shù)值積分的誤差來分析，很難期望隱式歐拉方法比顯式歐拉方法精確。為了得到更精確的方法，我們用梯形公式來近似前面的積分，得到梯形方法：它也是一種隱式單步法。改進的歐拉法從梯形法的推導(dǎo)，可望它比歐拉法更精確，但它計算量較大，在實際計算中，可取歐拉法的結(jié)果為迭代計算的初值，然后再用梯形公式計算一次，得到：歐拉方法。利用此法，每一步都要把該式作為的方程來求解。從數(shù)值積分的誤差來分析，很難期望隱式歐拉方法比顯式歐拉方法精確。為了得到更精確的方法，我們用梯形公式來近似前面的積分，得到梯形方法：它也是一種隱式單步法。改進的歐拉法從梯形法的推導(dǎo)，可望它比歐拉法更精確，但它計算量較大，在實際計算中，可取歐拉法的結(jié)果為迭代計算的初值，然后再用梯形公式計算一次，得到：或?qū)懗桑焊倪M的歐拉法是一種顯式單步法，有時為了計算方便，也用下面的公式：其中：單步法的截斷誤差與階初值問題的單步法可用一般形式表示為：其中多元函數(shù)與有關(guān)，當(dāng)含有時，方法是隱式的，不含時則為顯式的，所以顯式單步法可表示為：稱為增量函數(shù)，對于歐拉法，或?qū)懗桑焊倪M的歐拉法是一種顯式單步法，有時為了計算方便，也用下面的公式：其中：單步法的截斷誤差與階初值問題的單步法可用一般形式表示為：其中多元函數(shù)與有關(guān)，當(dāng)含有時，方法是隱式的，不含時則為顯式的，所以顯式單步法可表示為：稱為增量函數(shù)，對于歐拉法，為了分析其截斷誤差，我們將在處作泰勒展開，得到：假設(shè)是精確的，于是可得歐拉法計算公式的誤差為：

定義1設(shè)是初值問題的準確解，稱為顯式單步法的局部截斷誤差。之所以稱為局部的，是假設(shè)在前各步?jīng)]有誤差，當(dāng)時，計算一步，則有開，得到：假設(shè)是精確的，于是可得歐拉法計算公式的誤差為：

定義1設(shè)是初值問題的準確解，稱為顯式單步法的局部截斷誤差。之所以稱為局部的，是假設(shè)在前各步?jīng)]有誤差，當(dāng)時，計算一步，則有所以，局部截斷誤差可以理解為用計算的方法計算一步的誤差，這樣，歐拉法的局部截斷誤差為：這里稱為局部截斷誤差主項，顯然一般情況下的定義如下：

定義2設(shè)是初值問題的準確解，若存在最大整數(shù)使顯式單步法的局部誤差滿足則稱該方法具有階精度。若將上式展開成則稱為局部截斷誤差主項。所以，局部截斷誤差可以理解為用計算的方法計算一步的誤差，這樣，歐拉法的局部截斷誤差為：這里稱為局部截斷誤差主項，顯然一般情況下的定義如下：

定義2設(shè)是初值問題的準確解，若存在最大整數(shù)使顯式單步法的局部誤差滿足則稱該方法具有階精度。若將上式展開成則稱為局部截斷誤差主項。以上定義對隱式單步法也適用，如對后退歐拉法，其局部截斷誤差為：這里，是1階方法，局部截斷誤差主項為同樣，對梯形方法有：所以梯形方法是2階的，其局部誤差主項為：

以上定義對隱式單步法也適用，如對后退歐拉法，其局部截斷誤差為：這里，是1階方法，局部截斷誤差主項為

同樣，對梯形方法有：所以梯形方法是2階的，其局部誤差主項為：改進歐拉法的例梯形法精度得到提高，但算法復(fù)雜，計算量大。為了控制計算量，通常只迭代一兩次就轉(zhuǎn)入下一步的計算，這樣就可以簡化計算。具體做法：先用歐拉公式求一個初值，稱為預(yù)測值，再用梯形公式校正一次，這樣建立的公式通常稱為改進的歐拉公式：預(yù)測：校正：或表示為：所以梯形方法是2階的，其局部誤差主項為：改進歐拉法的例梯形法精度得到提高，但算法復(fù)雜，計算量大。為了控制計算量，通常只迭代一兩次就轉(zhuǎn)入下一步的計算，這樣就可以簡化計算。具體做法：先用歐拉公式求一個初值，稱為預(yù)測值，再用梯形公式校正一次，這樣建立的公式通常稱為改進的歐拉公式：預(yù)測：校正：或表示為：例：用改進的歐拉方法求解初值問題：解：仍取，計算結(jié)果如表0.11.09591.09540.61.48601.48320.21.18411.18320.71.55251.54920.31.26621.26490.81.61531.61250.41.34341.34160.91.67821.67330.51.41641.41421.01.73791.7321顯式龍格－庫塔法的一般形式歐拉法的局部截斷誤差為，是1階的方法對于改進的歐拉法，它可表示為：此時增量函數(shù)例：用改進的歐拉方法求解初值問題：解：仍取，計算結(jié)果如表0.11.09591.09540.61.48601.48320.21.18411.18320.71.55251.54920.31.26621.26490.81.61531.61250.41.34341.34160.91.67821.67330.51.41641.41421.01.73791.73219.3

龍格－庫塔方法顯式龍格－庫塔法的一般形式歐拉法的局部截斷誤差為，是1階的方法對于改進的歐拉法，它可表示為：此時增量函數(shù)它比歐拉法的多計算了一個函數(shù)值，可望，若要使得到的公式階數(shù)更大，就必須包含更多的值，實際上從方程等價的積分形式，即

(*)9.3

龍格－庫塔方法要使公式的階數(shù)提高，就必須使右端積分的數(shù)值求積公式精度提高，它必然要增加求積節(jié)點，為此可將(*)的右端用求積公式表示為：一般來說，點數(shù)越多，精度越高，上式右端相當(dāng)于增量函數(shù)，為得到便于計算的顯式方法，可它比歐拉法的多計算了一個函數(shù)值，可望，若要使得到的公式階數(shù)更大，就必須包含更多的值，實際上從方程等價的積分形式，即

(*)要使公式的階數(shù)提高，就必須使右端積分的數(shù)值求積公式精度提高，它必然要增加求積節(jié)點，為此可將(*)的右端用求積公式表示為：一般來說，點數(shù)越多，精度越高，上式右端相當(dāng)于增量函數(shù)，為得到便于計算的顯式方法，可類似于改進歐拉法，將公式表示為：其中這里均為常數(shù)，方法稱為級顯式R-K方法

當(dāng)時，就是歐拉法，此時當(dāng)時，改進歐拉法就是其中的一種，下面證明，要使公式有更高的階，就要增加點數(shù)。下面我們僅就推導(dǎo)R-K方法，并給出時的常用公式，其推導(dǎo)方法與時類似，只是計算較為復(fù)雜。增量函數(shù)，為得到便于計算的顯式方法，可類似于改進歐拉法，將公式表示為：其中這里均為常數(shù)，方法稱為級顯式R-K方法

當(dāng)時，就是歐拉法，此時當(dāng)時，改進歐拉法就是其中的一種，下面證明，要使公式有更高的階，就要增加點數(shù)。下面我們僅就推導(dǎo)R-K方法，并給出時的常用公式，其推導(dǎo)方法與時類似，只是計算較為復(fù)雜。二階顯式R-K方法對于的R-K方法，可以得到如下的計算公式：這里均為待定常數(shù)，我們希望適當(dāng)選取這些系數(shù)，使公式階數(shù)盡量高。根據(jù)局部截斷誤差的定義，上式的局部截斷誤差為：

這里為得到的階，要將上式各項在處做泰勒展開，由于是二元函數(shù)，故要用到二元泰勒展開，各項展開為：二階顯式R-K方法對于的R-K方法，可以得到如下的計算公式：這里均為待定常數(shù)，我們希望適當(dāng)選取這些系數(shù)，使公式階數(shù)盡量高。根據(jù)局部截斷誤差的定義，上式的局部截斷誤差為：這里為得到的階，要將上式各項在處做泰勒展開，由于是二元函數(shù)，故要用到二元泰勒展開，各項展開為：其中將以上結(jié)果代入：得到：其中將以上結(jié)果代入：得到：要使公式成為2階的，必須有：即：上式的解不唯一，可以令，則得這樣得到的公式稱為二階R-K方法，如取，則，這就是改進的歐拉法。若取，則，得到：要使公式成為2階的，必須有：即：上式的解不唯一，可以令，則得這樣得到的公式稱為二階R-K方法，如取，則，這就是改進的歐拉法。若取，則，得到：稱為中點公式，相當(dāng)于數(shù)值積分的中矩形公式，也即

對的R-K公式能否使局部誤差提高到？為此需要把多展開一項，通過分析可以得到不可能的結(jié)論。因此時只能得到2階的公式。三階與四階顯式R-K方法要得到三階的顯式R-K方法，必須，此時計算若取，則，得到：稱為中點公式，相當(dāng)于數(shù)值積分的中矩形公式，也即

對的R-K公式能否使局部誤差提高到？為此需要把多展開一項，通過分析可以得到不可能的結(jié)論。因此時只能得到2階的公式。三階與四階顯式R-K方法要得到三階的顯式R-K方法，必須，此時計算公式為：其中及均為待定參數(shù)，誤差為只要將按二元函數(shù)的泰勒公式展開，使可得待定參數(shù)需滿足的方程為：

要得到三階的顯式R-K方法，必須，此時計算公式為：其中及均為待定參數(shù)，誤差為只要將按二元函數(shù)的泰勒公式展開，使可得待定參數(shù)需滿足的方程為：這里為8個未知數(shù)，6個方程，解也是不唯一的，可以得到很多公式，下面是其中稱為庫塔的公式這里為8個未知數(shù)，6個方程，解也是不唯一的，可以得到很多公式，下面是其中稱為庫塔的公式繼續(xù)上述過程，經(jīng)過較復(fù)雜的演算，可以導(dǎo)出各種四階龍格－庫塔公式，下面是常用的經(jīng)典公式：例3：用四階龍格－庫塔方法求解定解問題：取步長，從到。解：寫出具體的四階龍格－庫塔的具體格式為：繼續(xù)上述過程，經(jīng)過較復(fù)雜的演算，可以導(dǎo)出各種四階龍格－庫塔公式，下面是常用的經(jīng)典公式：例3：用四階龍格－庫塔方法求解定解問題：取步長，從到。解：寫出具體的四階龍格－庫塔的具體格式為：

列出計算結(jié)果如下表：0.21.18321.18320.41.34171.34160.61.48331.48320.81.61251.61251.01.73211.7321列出計算結(jié)果如下表：可見計算精度較以前大為提高。變步長的龍格－庫塔方法單從每一步看，步長越小，截斷誤差就越小，但隨著步長的縮小，在求解區(qū)間內(nèi)的步數(shù)就增加了，這在增加計算量的同時也增加了計算的舍入誤差，因此也有一個選擇步長的問題。0.21.18321.18320.41.34171.34160.61.48331.48320.81.61251.61251.01.73211.7321選擇步長需要考慮的兩個問題：怎樣衡量和檢驗計算結(jié)果的精度？如何依據(jù)所獲得的精度處理步長？我們考慮經(jīng)典的四階龍格－庫塔公式，從節(jié)點出發(fā)，先以為步長求出一個近似值，記為，由于公式的局部截斷誤差為，故變步長的龍格－庫塔方法單從每一步看，步長越小，截斷誤差就越小，但隨著步長的縮小，在求解區(qū)間內(nèi)的步數(shù)就增加了，這在增加計算量的同時也增加了計算的舍入誤差，因此也有一個選擇步長的問題。選擇步長需要考慮的兩個問題：怎樣衡量和檢驗計算結(jié)果的精度？如何依據(jù)所獲得的精度處理步長？我們考慮經(jīng)典的四階龍格－庫塔公式，從節(jié)點出發(fā)，先以為步長求出一個近似值，記為，由于公式的局部截斷誤差為，故然后將步長折半，即取為步長，從跨兩步到再來求得一個近似值，每一步的截斷誤差是因此有比較前后兩式，我們得到，步長折半后，誤差大約減少到，即有由此易得下列事后估計式：這樣，我們可以通過檢查步長，折半前后兩次計算結(jié)果的偏差來判定所選的步長是否合適，公式的局部截斷誤差為，故然后將步長折半，即取為步長，從跨兩步到再來求得一個近似值，每一步的截斷誤差是因此有

比較前后兩式，我們得到，步長折半后，誤差大約減少到，即有由此易得下列事后估計式：這樣，我們可以通過檢查步長，折半前后兩次計算結(jié)果的偏差來判定所選的步長是否合適，分下列兩種情況處理，稱為變步長方法：

1)對于給定的精度，如果，我們反復(fù)將步長折半進行計算，直到為止，這時取最終得到的作為結(jié)果；

2)如果，我們將反復(fù)將步長加倍，直到為止，這時再將步長折半一次，得到結(jié)果。9.4

單步法的收斂性與穩(wěn)定性收斂性與相容性數(shù)值解法的基本思想是通過離散將微分方程轉(zhuǎn)化為差分方程，如單步法

(1)它在處的解為，而初值問題的解為，記分下列兩種情況處理，稱為變步長方法：

1)對于給定的精度，如果，我們反復(fù)將步長折半進行計算，直到為止，這時取最終得到的作為結(jié)果；

2)如果，我們將反復(fù)將步長加倍，直到為止，這時再將步長折半一次，得到結(jié)果。9.4

單步法的收斂性與穩(wěn)定性收斂性與相容性數(shù)值解法的基本思想是通過離散將微分方程轉(zhuǎn)化為差分方程，如單步法

(1)它在處的解為，而初值問題的解為，記，稱為整體截斷誤差，收斂性就是討論當(dāng)固定，且時的問題。定義3若一種數(shù)值方法對于固定的，當(dāng)時有，則稱該方法是收斂的。顯然，數(shù)值方法收斂是指，對單步法(1)有下述收斂性定理：9.4

單步法的收斂性與穩(wěn)定性

定理1假設(shè)單步法(1)具有階精度，且增量函數(shù)關(guān)于滿足利普希茨條件

(2)又設(shè)初值準確，即，則其整體截斷誤差為

(3)證明：設(shè)以表示取用公式(1)求得的結(jié)果,

固定，且時的問題。定義3若一種數(shù)值方法對于固定的，當(dāng)時有，則稱該方法是收斂的。顯然，數(shù)值方法收斂是指，對單步法(1)有下述收斂性定理：

定理1假設(shè)單步法(1)具有階精度，且增量函數(shù)關(guān)于滿足利普希茨條件

(2)又設(shè)初值準確，即，則其整體截斷誤差為

(3)證明：設(shè)以表示取用公式(1)求得的結(jié)果,即(4)則為局部截斷誤差，由于所給的方法具有階精度，按定義2，存在定數(shù)，使又由(4)和(1)，有利用條件(2)，有：從而有即對整體截斷誤差成立下列遞推關(guān)系：反復(fù)遞推，可得：證明：設(shè)以表示取用公式(1)求得的結(jié)果,即(4)則為局部截斷誤差，由于所給的方法具有階精度，按定義2，存在定數(shù)，使又由(4)和(1)，有利用條件(2)，有：從而有即對整體截斷誤差成立下列遞推關(guān)系：反復(fù)遞推，可得：再注意到當(dāng)時，，最終得到下列估計式：由此可以斷定，如果初值是準確的，則(3)成立。

推論：對于一個階的顯式單步法，若微分方程的右端函數(shù)關(guān)于滿足利普希茨條件，且初值是精確的，則顯式歐拉法，改進歐拉法和龍格－庫塔方法是收斂的。定理1表明，時單步法收斂，并且當(dāng)是初值問題的解，且具有階精度時，有展開式：再注意到當(dāng)時，，最終得到下列估計式：由此可以斷定，如果初值是準確的，則(3)成立。

推論：對于一個階的顯式單步法，若微分方程的右端函數(shù)關(guān)于滿足利普希茨條件，且初值是精確的，則顯式歐拉法，改進歐拉法和龍格－庫塔方法是收斂的。定理1表明，時單步法收斂，并且當(dāng)是初值問題的解，且具有階精度時，有展開式：所以的充要條件是，而，于是可給出如下的定義：定義4單步法的增量函數(shù)稱為該單步法與初值問題相容。

以上討論表明階方法當(dāng)與初值問題相容，反之，相容的方法至少是一階的。于是，由定理1可知，線性單步方法收斂的充分必要條件是該方法是相容的。絕對穩(wěn)定性與絕對穩(wěn)定域

所以的充要條件是，而，于是可給出如下的定義：定義4單步法的增量函數(shù)稱為該單步法與初值問題相容。

以上討論表明階方法當(dāng)與初值問題相容，反之，相容的方法至少是一階的。于是，由定理1可知，線性單步方法收斂的充分必要條件是該方法是相容的。絕對穩(wěn)定性與絕對穩(wěn)定域定義5若一種數(shù)值方法在節(jié)點值上有大小為的擾動，而于以后各節(jié)點上產(chǎn)生的偏差均不超過，則稱該方法是穩(wěn)定的。例4考察初值問題其準確解是一個按指數(shù)曲線衰減的很快的函數(shù)，如圖：絕對穩(wěn)定性與絕對穩(wěn)定域

定義5若一種數(shù)值方法在節(jié)點值上有大小為的擾動，而于以后各節(jié)點上產(chǎn)生的偏差均不超過，則稱該方法是穩(wěn)定的。例4考察初值問題其準確解是一個按指數(shù)曲線衰減的很快的函數(shù)，如圖：

用歐拉方法解該方程，得到：若取，則歐拉法的具體形式為：若取，則歐拉法的具體形式為：顯然，前一形式計算不穩(wěn)定，后一形式計算穩(wěn)定。再考察后退歐拉方法，可以得到，當(dāng)時，計算是穩(wěn)定的。具體數(shù)據(jù)見書上356頁表9－4。例題表明：穩(wěn)定性不但與方法有關(guān)，也與步長有關(guān)。當(dāng)然也與有關(guān)，為了只考察數(shù)值方法本身，通常只檢驗將數(shù)值方法用于解模型方程的穩(wěn)定性，模型方程為：，其中為復(fù)數(shù)。

對于一般的方程，將方程中的在解域內(nèi)某點作泰勒展開并局部線性化：忽略高階項，令和對上式顯然，前一形式計算不穩(wěn)定，后一形式計算穩(wěn)定。再考察后退歐拉方法，可以得到，當(dāng)時，計算是穩(wěn)定的。具體數(shù)據(jù)見書上356頁表9－4。例題表明：穩(wěn)定性不但與方法有關(guān)，也與步長有關(guān)。當(dāng)然也與有關(guān)，為了只考察數(shù)值方法本身，通常只檢驗將數(shù)值方法用于解模型方程的穩(wěn)定性，模型方程為：，其中為復(fù)數(shù)。

對于一般的方程，將方程中的在解域內(nèi)某點作泰勒展開并局部線性化：忽略高階項，令和對上式作變量代換，得到：

這就是模型方程。首先研究歐拉方程的穩(wěn)定性。模型方程的歐拉公式為：假設(shè)在節(jié)點上有一擾動，它的傳播使節(jié)點值產(chǎn)生大小為的擾動值，假設(shè)用按歐拉公式得出的計算過程不再有新的誤差，則擾動值滿足可見擾動值滿足原來的差分方程，這樣如果差分方程的解是不增長的，即有：

這就是模型方程。首先研究歐拉方程的穩(wěn)定性。模型方程的歐拉公式為：假設(shè)在節(jié)點上有一擾動，它的傳播使節(jié)點值產(chǎn)生大小為的擾動值，假設(shè)用按歐拉公式得出的計算過程不再有新的誤差，則擾動值滿足可見擾動值滿足原來的差分方程，這樣如果差分方程的解是不增長的，即有：則它就是穩(wěn)定的。這一結(jié)論對于以下研究的其它方法也適用。顯然，為了保證，只要，這在的復(fù)平面上是一個以(-1,0)為圓心，1為半徑的單位圓我們稱其為歐拉法的絕對穩(wěn)定域，一般情形可有如下的定義。

定義6單步法用于模型，若得到的解，