在數學中,有限差分法(finite-difference methods,簡稱FDM),是一種微分方程數值方法,是通過有限差分來近似導數,從而尋求微分方程的近似解。
首先假設要近似函數的各級導數都有良好的性質,依照泰勒定理,可以形成以下的泰勒展開式:
![{\displaystyle f(x_{0}+h)=f(x_{0})+{\frac {f'(x_{0})}{1!}}h+{\frac {f^{(2)}(x_{0})}{2!}}h^{2}+\cdots +{\frac {f^{(n)}(x_{0})}{n!}}h^{n}+R_{n}(x),}](https://wikimedia.org/api/rest_v1/media/math/render/svg/914de73645542676cfad7bd9e16c93970a2ece1b)
其中n!表示是n的階乘,Rn(x)為餘數,表示泰勒多項式和原函數之間的差。可以推導函數f一階導數的近似值:
![{\displaystyle f(x_{0}+h)=f(x_{0})+f'(x_{0})h+R_{1}(x),}](https://wikimedia.org/api/rest_v1/media/math/render/svg/63cf9cba2c58c40088ff82f0898b015bd16b4480)
設定x0=a,可得:
![{\displaystyle f(a+h)=f(a)+f'(a)h+R_{1}(x),}](https://wikimedia.org/api/rest_v1/media/math/render/svg/a37e1e478731df6cb47927066529529fa09fb60f)
除以h可得:
![{\displaystyle {f(a+h) \over h}={f(a) \over h}+f'(a)+{R_{1}(x) \over h}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/2c0e1d411304ff08d7a960e1d532d664734ea7d9)
求解f'(a):
![{\displaystyle f'(a)={f(a+h)-f(a) \over h}-{R_{1}(x) \over h}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/e2d6682b32498b37de1ac5ec1c5c315f7567f6d3)
假設
相當小,因此可以將"f"的一階導數近似為:
![{\displaystyle f'(a)\approx {f(a+h)-f(a) \over h}.}](https://wikimedia.org/api/rest_v1/media/math/render/svg/73850165addec94b0dc55f6204a2b32b6371f522)
近似解的誤差定義為近似解及解析解之間的差值。有限差分法的兩個誤差來源分別是捨入誤差及截尾誤差(或稱為離散化誤差),前者是因為電腦計算小數時四捨五入造成的誤差,後者則是用有限階級數表示導數引起的誤差。
有限差分法是以在格點上函數的值為準
在運用有限差分法求解一問題(或是說找到問題的近似解)時,第一步需要將問題的定義域離散化。一般會將問題的定義域用均勻的網格分割(可參考右圖)。因此有限差分法會製造一組導數的離散數值近似值。
一般會關注近似解的局部截尾誤差,會用大O符號表示,局部截尾誤差是指應用有限差分法一次後產生的誤差,因此為
,此時
是實際值,而
為近似值。泰勒多項式的餘數項有助於分析局部截尾誤差。利用
泰勒多項式的餘數項,也就是
, 其中
,
可以找到局部截尾誤差的主控項,例如用前項差分法計算一階導數,已知
,
![{\displaystyle f(x_{0}+ih)=f(x_{0})+f'(x_{0})ih+{\frac {f''(\xi )}{2!}}(ih)^{2},}](https://wikimedia.org/api/rest_v1/media/math/render/svg/e3bd1b9a0b46bc10331724922cbafff4ef65114f)
利用一些代數的處理,可得
![{\displaystyle {\frac {f(x_{0}+ih)-f(x_{0})}{ih}}=f'(x_{0})+{\frac {f''(\xi )}{2!}}ih,}](https://wikimedia.org/api/rest_v1/media/math/render/svg/8a1c7b3020fc81e7e5103f29d5d113d3ea56c480)
注意到左邊的量是有限差分法的近似,右邊的量是待求解的量再加上一個餘數,因此餘數就是局部截尾誤差。上述範例可以用下式表示:
![{\displaystyle {\frac {f(x_{0}+ih)-f(x_{0})}{ih}}=f'(x_{0})+O(h).}](https://wikimedia.org/api/rest_v1/media/math/render/svg/3642c8eb151d313e6e37925c6d4ae9c29a655405)
在此例中,局部截尾誤差和時間格點的大小成正比。
例如考慮以下的常微分方程
![{\displaystyle u'(x)=3u(x)+2.\,}](https://wikimedia.org/api/rest_v1/media/math/render/svg/43fa12a893257b8930925ecafb4370f6cc8722ef)
利用數值方法中歐拉法求解,利用以下的有限差分式
![{\displaystyle {\frac {u(x+h)-u(x)}{h}}\approx u'(x)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/8f664c9609a75482c79d641f0a49ef2e62f60c6e)
來近似導數,並配合一些代數處理(等號兩側同乘以h,再加上u(x)),可得
![{\displaystyle u(x+h)=u(x)+h(3u(x)+2).\,}](https://wikimedia.org/api/rest_v1/media/math/render/svg/d47ab29816a0b6030b7a257e411339104a9f5375)
最後的方程式即為有限差分方程,求解此方程則可得到原方程的近似解。
考慮正規化的一維熱傳導方程式,為齊次的狄利克雷邊界條件
![{\displaystyle U_{t}=U_{xx}\,}](https://wikimedia.org/api/rest_v1/media/math/render/svg/e6750ab883a2f9e955ce6f518e997d2abaf2e0b2)
(邊界條件)
(初始條件)
對此問題求數值解的一種方式是用差分去近似所有的導數,可以將空間分割為
,將時間也分割為
。假設在時間及空間都是均勻的網格切割,空間中兩個連續位置的間隔為h,兩個連續時間之間的間隔為k。點
![{\displaystyle u(x_{j},t_{n})=u_{j}^{n}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/5a34fed61644bcbdbfedfbe1150911a9f5fe9dff)
表示
的數值近似解。
熱傳導方程最常用顯式方法的模版
利用在時間
的前向差分,以及在位置
的二階中央差分(FTCS 格式),可以得到以下的迭代方程:
![{\displaystyle {\frac {u_{j}^{n+1}-u_{j}^{n}}{k}}={\frac {u_{j+1}^{n}-2u_{j}^{n}+u_{j-1}^{n}}{h^{2}}}.\,}](https://wikimedia.org/api/rest_v1/media/math/render/svg/16537d53a7599049e4dcc8c13f749bdf3fc2b1b2)
這是用求解一維導熱傳導方程的顯式方法。
可以用以下的式子求解
![{\displaystyle u_{j}^{n+1}=(1-2r)u_{j}^{n}+ru_{j-1}^{n}+ru_{j+1}^{n}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/f3687ff4fab682f182cb4468868a39011ca2efd2)
其中
因此配合此迭代關係式,已知在時間n的數值,可以求得在時間n+1的數值。
及
的數值可以用邊界條件代入,在此例中為0。
此顯式方法在
時,為數值穩定且收斂[1]。其數值誤差和時間間隔成正比,和位置間隔的平方成正比:
![{\displaystyle \Delta u=O(k)+O(h^{2})\,}](https://wikimedia.org/api/rest_v1/media/math/render/svg/dc94dc45aea2ebd288b93cdd0b2a4a92eb4e77c0)
隱式方法的模版
若使用時間
的後向差分,及位置
的二階中央差分(BTCS 格式),可以得到以下的迭代方程:
![{\displaystyle {\frac {u_{j}^{n+1}-u_{j}^{n}}{k}}={\frac {u_{j+1}^{n+1}-2u_{j}^{n+1}+u_{j-1}^{n+1}}{h^{2}}}.\,}](https://wikimedia.org/api/rest_v1/media/math/render/svg/3016a002021b67daa7908a78f82fa2e8aeb90df6)
這是用求解一維導熱傳導方程的隱式方法。
在求解線性聯立方程後可以得到
:
![{\displaystyle (1+2r)u_{j}^{n+1}-ru_{j-1}^{n+1}-ru_{j+1}^{n+1}=u_{j}^{n}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/fa25e6866ec5a224a5f6462afab0ee10cbb03c60)
此方法不論
的大小,都數值穩定且收斂,但在計算量會較顯式方法要大,因為每前進一個時間間隔,就需要求解一個聯立的數值方程組。其數值誤差和時間間隔成正比,和位置間隔的平方成正比:
![{\displaystyle \Delta u=O(k)+O(h^{2})\,}](https://wikimedia.org/api/rest_v1/media/math/render/svg/dc94dc45aea2ebd288b93cdd0b2a4a92eb4e77c0)
若使用時間
的中間差分,及位置
的二階中央差分(CTCS 格式),可以得到以下的迭代方程:
![{\displaystyle {\frac {u_{j}^{n+1}-u_{j}^{n}}{k}}={\frac {1}{2}}\left({\frac {u_{j+1}^{n+1}-2u_{j}^{n+1}+u_{j-1}^{n+1}}{h^{2}}}+{\frac {u_{j+1}^{n}-2u_{j}^{n}+u_{j-1}^{n}}{h^{2}}}\right).\,}](https://wikimedia.org/api/rest_v1/media/math/render/svg/3bc7faf6e584d2c8aed8e3236a6dee1dc5ac8d62)
此公式為克蘭克-尼科爾森方法(Crank–Nicolson method)。
克蘭克-尼科爾森方法的模版
在求解線性聯立方程後可以得到
:
![{\displaystyle (2+2r)u_{j}^{n+1}-ru_{j-1}^{n+1}-ru_{j+1}^{n+1}=(2-2r)u_{j}^{n}+ru_{j-1}^{n}+ru_{j+1}^{n}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/46749522cb306348ad3673a3e4d49be5689cd3df)
此方法不論
的大小,都數值穩定且收斂,但在計算量會較顯式方法要大,因為每前進一個時間間隔,就需要求解一個聯立的數值方程組。其數值誤差和時間間隔的平方成正比,和位置間隔的平方成正比:
![{\displaystyle \Delta u=O(k^{2})+O(h^{2}).\,}](https://wikimedia.org/api/rest_v1/media/math/render/svg/63940bfe514171bf62ce12fe2e746fe8debfab1f)
若時間刻度較小時,克蘭克-尼科爾森方法是最精確的,而顯式方法是最不精確的,而且可能會不穩定,但是是最容易計算的,其數值計算量也最少。若時間刻度較大時,隱式方法的效果最好。
- ^ Crank, J. The Mathematics of Diffusion. 2nd Edition, Oxford, 1975, p. 143.