求解微分方程的線性多步法

摘 要: 微分方程初值問題的求解一直是人們所關注的熱點問題，本文針對微分方程初值問題的線性多步法，利用它們的局部截斷誤差的主項系數的某種組合對它們進行修正，提出了預測—校正算法，并給出了幾種常用的預測—校正算法。

關鍵詞: 微分方程 初值問題 線性多步法

1.微分方程初值問題的基本解法

=f(t，u)u(t)=u(1)

其中f為t和u的已知函數，u為給定的初值。我們假設函數f(t，u)在區域:t≤t≤T，|u|<∞內連續，并且u滿足Lipschitz條件，即存在常數L，對所有t∈[t，T]和u，u，有|f(t，u)-f(t，u)|≤L|u-u|。在上述基本假設下，初值問題(1)在區間[t，T]上有唯一解，并且u(t)為連續可微的。

將區間[0，T]作N等分，小區間的長度h=T/N稱為步長，點列t=nh(n=0，1，…，N)稱為節點，t=0。由已知初值u(t)=u，可算出u(t)在t=t的導數值u′(t)=f(t，u(t))=f(t，u)。利用Taylor展式得u=u+hf(t，u)，u就是u(t)的近似值。利用u又可算出u，如此下去可算出u在所有結點上的值，一般遞推公式為u=u+hf(t，u)，n=0，1，…，N-1。

由文獻[1]知，求解該類方程分為單步法和多步法。只要利用h，t和u即可算出u的算法稱為單步方法，單步法包括Euler法、梯形法、Runge-kutta法，其一般形式為u=u+hφ(t，u;h)，函數φ稱為增量函數，它依賴于所給的微分方程。用單步法求出{φ}，m=1，2，…，N只需要一個初值u。

需要用到h，t，t，…，t和h，u，u，…，u才能求出u(k>1)的算法稱為多步方法，線性多步法包括數值積分法和待定系數法，其一般形式為αu=hβf，這里α，β(j=0，1，K，k)是常數，它需要有k個初值u，u，…，u才能得到整個序列{u}，m=1，2，…，N。

2.線性多步法的改進

常微分方程初值問題只給我們提供了一個初值，但線性k步方法需要k個初始值u，u，…，u，所以這其余的k-1個初始值就要通過其他方法來先期得到。最容易想到的是采用簡單的單步方法，如Euler方法，或改進的Euler法，但它們的收斂階較低。為保證整個計算的精度，我們應使初始值的計算精度至少與所用的線性多步方法同階。

因此，若使用的線性多步方法為q階[2]，[3]的，可以用q階Taylor展開方法或Runge-kutta方法來求u，u，…，u，然后進入正式的求解過程。用Taylor展開確定初始值的同時還能得到關于選擇h的信息。如要求計算誤差不超過ε，則應使

h|u(t)|≤εh|u(t)|>ε

同時成立。當微商次數增大時，計算量也將急速增大，因此利用低階微商構造高精度的方法對實際計算是很有意義的，對u(t)=u(t)+?蘩u′(t)dt利用u′(t)在點t，t處的函數值及其各階導數值構造Hernite型插值多項式，再代積分，舍去誤差項即可得到帶導數值的計算公式。

對隱式的計算格式一般采用迭代解法，將隱式線性多步方法改寫成

u-hβf(t，u)=-(αu-hβf)(2)

由于F是已知的，于是(2)式可用迭代法求解。

迭代法都需要有u，即u的迭代初值。選取一個好的u無疑是十分重要的，很自然的想法是用顯格式算出u——這稱為預測，再用隱格式進行迭代求u——這稱為校正，因此整個過程稱為預測—校正算法[4]，簡稱PC算法。顯然，預測的P式與校正的C式的階一般應取成相同的。

預測—校正算法求u的過程一般如下:

P∶u=αu-hβf，

E∶f=f(t，u)C∶u=-αu+h(∑βf+βfs=1，2，…，n，

u=u，f=f(t，u)。

由于使用了同階的預測公式求得初始近似u，因此一般說來校正次數不會太多，大約兩三次就夠了。若校正步數太多則往往提示步長過大，應適當減小步長后再進行計算。實際計算中經常利用算式與算式同階的特點，對(3)式再作一些修正。用線性組合的方法消去每個格式中的截斷誤差主項。僅花費很小的代價就使方法的精度再提高一階，從而不需要再通過迭代來改善計算的近似值。

設預測算式和校正算式的局部截斷誤差分別為

Δ[u(t)，h]=chu(t)+o(h)(3)

Δ[u(t)，h]=chu(t)+o(h)(4)

整理后分別得到:

u(t)-[u+(u-u)]=o(h)(5)

u(t)-[u+(u-u)]=o(h)(6)

這兩個式子告訴我們，在求(u和u)后，再利用它們的局部截斷誤差[5]，[6]的主項系數的某種組合對它們進行修正，可使得修正后的局部截斷誤差比原來高一階，即達到o(h)。這里和被稱為修正系數。

直接按照(5)式對u進行修正是不行的，因為此時的u還沒有算出。但是，在計算u的近似值時，u已經算出來了，因此u和u也都已經有了，這時可以用u-u去代替(5)中的u-u。

綜上所述，我們可以將帶修正的預測—校正算法的一般形式表示如下:

P∶u=-α+hβfM∶+(u-u)E∶=f(t，)C∶u=-α+hβ+hβfM∶u=u+(u-u)E∶f=f(t，u)(7)

在多步方法中，改善步長將帶來一個新問題，就是結點變了。在預期結點處已經算出的近似值，對于計算當前結點來說，很多都用不上了，而在以新步長為單位的結點處，往往并沒有計算過相應的近似函數值。在使用的方法中，一般采用插值的方法去補上所需的這些點上的近似值再進行計算。

下面給出幾種常用的帶修正的預測—校正算法。

(1)Milne預校算法。Hamming針對Milne算法絕對不穩定，提出了如下修正:

P∶ρ(λ)=λ-1，σ(λ)=(2λ-λ-λ);

C∶ρ(λ)=λ-1，σ(λ)=(λ+4λ+1)。

P和C均為四階方法，M和M中的修正系數分別為和-。

(2)Hamming預校算法。該預校算法的P式和C式分別為:P∶與Milne預校算法相同。C∶ρ(λ)=λ-λ+，σ(λ)=(λ+2λ+λ)。

P和C也均為四階方法，M和M中的修正系數分別為和-。

(3)Adams四階預校算法。此算法中的P算式和C算式分別為:P∶Adams四步四階外插公式。P∶Adams三步四階內插公式。

M和M中的修正系數分別為和-。理論研究表明，這個算法的穩定性比前兩個算法好，是常用的預校算法。

3.結語

針對微分方程初值問題的線性多步法，我們利用局部截斷誤差的主項系數的某種組合對其進行修正，使修正后的局部截斷誤差比原來高一階，即達到o(h)，并給出了帶修正的預測—校正算法的一般形式，見文中(7)式。

參考文獻:

[1]李立康，於崇華，朱政華編著.微分方程數值解法[M].復旦大學出版社，2005.

[2]胡健偉.微分方程數值方法[M].科學出版社，1999.

[3]李榮華.微分方程數值解法(第三版)[M].高等教育出版社，1996.

[4]Dennis G.Zill編著.微分方程與邊界值問題[M].機械工業出版社，2003.

[5]黃振侃.編著.數值計算——微分方程數值解[M].北京工業大學出版社，2006.

[6][美]茲爾著.陳啟宏等譯.微分方程與邊界值問題[M].華章數學譯叢.機械工業出版社，2005.