常微分方程初值問題數值解法

(1)

利用數值方法解問題 (1)時，通常假定解存在且惟一，解函數y（x）及右端函數?(x,y)具有所需的光滑程度。數值解法的基本思想是:先取自變量一系列離散點,把微分問題(1)離散化，求出離散問題的數值解,並以此作為微分問題解y(x)的近似。例如取步長h＞0,以h剖分區間【α,b】，令xi=α+ih,把微分方程離散化成一個差分方程。以y(x)表微分方程初值問題的解,以yi表差分問題的解，就是近似解的誤差，稱為全局誤差。因此，設計各種離散化模型,求出近似解,估計誤差以及研究數值方法的穩定性和收斂性等構成了數值解法的基本內容。

離散化方法 　常用的有三種：

① 基於數值微分的方法 將方程(1)左端的導數用某個一階數值微分公式代替,例如在xn點以(yn+1-yn)/h代替yń即得到歐拉向前公式

(2)

若在xn+1點以(yn+1-yn)/h代替則得到歐拉向後公式

(3)

取(2)、(3)的平均，可導出二階精度的梯形公式

(4)

② 基於泰勒展開的方法　設計一個算法,假定公式中含有某些待定常數，在函數光滑的假定下，將其按泰勒展開並與微分方程解y(xn+h)的展式中h的同冪次項相比較，按照給定的精度階得到待定常數應滿足的一些方程，通過這些方程確定待定常數，即可得到所要的差分公式。由此法可導出龍格－庫塔公式。

設計算公式有下列形式

(5)

αi、βij為待定常數。取定N值，可按上述泰勒展開的方法確定它們。最常用的顯式4階龍格-庫塔公式為

(6)

式中　③ 基於函數數值積分的方法 將微分方程的解y(x)代入方程(1),在子區間【xn-ih,xn+jh】上積分得到公式

(7)

積分號下是x的函數，若用某些結點上?的值的數值積分公式近似這個積分，便得到各種差分公式，特別地，若取i=0,j=1,並用?在結點xn,xn-1,xn-2,…上的插值代替(7)式中的被積函數，便得亞當斯外推公式。4階亞當斯外推公式為

(8)

若取?在結點xn+1,xn,xn-1,…上的插值代替(7)式的被積函數，則得亞當斯內插公式。4階亞當斯內插公式為

(9)

解法可按計算yn+1時用多少個結點上的值分為單步法和多步法,又可以按yn+1出現的形式分為顯式法和隱式法。

單步法是指已知結點xn上yn的值便可計算yn+1的值的解法,如(2)、(3)、(4)。單步法是可以自己起步的,即可從方程的初值y0一步步算出y1,y2,…的值。

多步法是指已知yn,yn-1,…,yn-k+1(k≥2)的值才能計算yn+1的值的解法,又稱k步法。例如，(8)是四步法，(9)是三步法。多步法不能自己起步,即給了初值y0以後,還要用其他解法(如單步法)，算出y1,y2,…,yk-1後,才能使用多步法，繼續往下計算。多步法公式若對yi和?i都是線性的,則稱作線性多步法，k步線性多步法的一般形式為　顯式法的公式中,未知的yn+1明顯地被表示，即公式中除yn+1一項外，其他的項中不再含有yn+1,如公式(2)。

隱式法的公式不顯含yn+1,求未知的yn+1時一般需要解方程，如公式(3)或(4)。通常用各種迭代方法解隱式差分方程,也可採用較簡單的預估-校正方法，如使用梯形公式(4)時，可先用顯式公式(2)求得yn+1的預估值,代入式(4)的函數?n+1中, 再求得yn+1的值。此法又稱改進的歐拉折線法。

數值解法滿足相容的、收斂的、數值穩定的條件時，才有實用價值。為此要研究以下的一些問題。

相容性 　將微分方程離散化所帶來的誤差叫截斷誤差。當h→0時，截斷誤差趨於零，則稱離散化後的方程與微分方程具有相容性，表示離散化後的方程是微分方程的近似。若截斷誤差的主要項為Chp+1,則稱截斷誤差的階是p+1,而稱該解法是p階的。p越大表示離散化後的方程與微分方程近似程度越高。

收斂性 　是指當h→0時,全局誤差εi→0，即離散問題的解yn收斂於微分問題的解y(x)，這是離散解可用的理論基礎。p階的解法,即是當h→0時,εi以hp的速度收斂。

誤差估計 　對全局誤差εi的估計，是應用數值解法時最關心的問題。先驗估計通常只能給出誤差的階，即誤差的主要項中步長h的冪次。一般採用事後估計,即在計算的過程中估計誤差，例如用理查森外推法估計誤差。外推法也是提高解的精確度的有效方法。

數值穩定性 　是指計算過程中，某一步上產生的誤差一步一步地傳遞下去,是衰減、不增或有界,使得傳遞下來的誤差不致於影響數值解的精度，至少是不會湮沒數值解。數值穩定性是常微分方程數值積分時必須考慮的問題。

1956年G.達赫爾斯特證明:存在2k階k步線性多步法，但數值穩定的k步線性多步法，當k為偶數時，其階不能超過k+2，當k為奇數時,其階不能超過k+1。稱為限制性定理。

判別一個數值方法的穩定性時，微分方程

(10)

有較廣的代表性，這裏λ=α+iβ,αh複平面上使方法數值穩定的區域。歐拉公式（2）的穩定條件為|1+λh|h取得使所有λh（λ是A的本徵值）都落在穩定區域內，用這個步長積分時數值穩定，理論上穩定區域可以用來選擇可用的步長h。由圖1、圖2可知，對於αh取得足夠小，λh就可落在穩定區域內。對於圖3，h沒有限制。穩定區域包含 λh平面的整個左半平面的方法叫做A穩定的。如公式(3)、(4)是A穩定的。

剛性方程組、常微分方程組的初值問題為：式中皆為n維向量。假定微分方程右端函數?的雅可比矩陣的特徵值為且 Reλj(j=1,2,…,n)，當比值時,用通常顯式公式計算,│hλj│(j=1,2,…,n)不能超過某個量，也就是h必須很小,從而大大增加計算時間。這類問題稱為剛性問題。因為A穩定的方法從數值穩定來説對步長h沒有限制,適用於剛性方程組，如歐拉向後公式(3),梯形公式(4)，隱式龍格-庫塔公式等都是有效的方法。

參考書目

P. Henrici,Discrete variable Methods in Ordinary Differential Equations, John Wiley & Sons, New York, 1962.

C. W.吉爾著,費景高、劉德貴、高永春譯：《常微分方程初值問題的數值解法》，科學出版社，北京，1978。(C. W. Gear,Numerical Initial Value Problems in Ordinary Differential Equations,Prentice-Hall, Englewood Cliffs, New Jersey,1971.)