有限差分法

在数学中，有限差分法（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),}$ 

其中n!表示是n的階乘，R_n(x)為餘數，表示泰勒多項式和原函數之間的差。可以推導函數f一階導數的近似值：

${\displaystyle f(x_{0}+h)=f(x_{0})+f'(x_{0})h+R_{1}(x),}$ 

設定x₀=a，可得：

${\displaystyle f(a+h)=f(a)+f'(a)h+R_{1}(x),}$ 

除以h可得：

${\displaystyle {f(a+h) \over h}={f(a) \over h}+f'(a)+{R_{1}(x) \over h}}$ 

求解f'(a)：

${\displaystyle f'(a)={f(a+h)-f(a) \over h}-{R_{1}(x) \over h}}$ 

假設${\displaystyle R_{1}(x)}$ 相當小，因此可以將"f"的一階導數近似為：

${\displaystyle f'(a)\approx {f(a+h)-f(a) \over h}.}$ 

準確度及誤差

近似解的誤差定義為近似解及解析解之間的差值。有限差分法的兩個誤差來源分別是捨入誤差及截尾誤差（英语：truncation error）（或稱為離散化誤差），前者是因為電腦計算小數時四捨五入造成的誤差，後者則是用有限阶级数表示导数引起的误差。

 

有限差分法是以在格點上函數的值為準

在運用有限差分法求解一問題（或是說找到問題的近似解）時，第一步需要將問題的定義域離散化。一般會將問題的定義域用均勻的網格分割（可參考右圖）。因此有限差分法會制造一組導數的離散數值近似值。

一般會關注近似解的局部截尾誤差（英语：local truncation error），會用大O符號表示，局部截尾誤差是指應用有限差分法一次後產生的誤差，因此為${\displaystyle f'(x_{i})-f'_{i}}$ ，此時${\displaystyle f'(x_{i})}$ 是實際值，而${\displaystyle f'_{i}}$ 為近似值。泰勒多項式的餘數項有助於分析局部截尾誤差。利用${\displaystyle f(x_{0}+h)}$ 泰勒多項式的餘數項，也就是

${\displaystyle R_{n}(x_{0}+h)={\frac {f^{(n+1)}(\xi )}{(n+1)!}}(h)^{n+1}}$ , 其中${\displaystyle x_{0}<\xi <x_{0}+h}$ ,

可以找到局部截尾誤差的主控項，例如用前項差分法計算一階導數，已知${\displaystyle f(x_{i})=f(x_{0}+ih)}$ ,

${\displaystyle f(x_{0}+ih)=f(x_{0})+f'(x_{0})ih+{\frac {f''(\xi )}{2!}}(ih)^{2},}$ 

利用一些代數的處理，可得

${\displaystyle {\frac {f(x_{0}+ih)-f(x_{0})}{ih}}=f'(x_{0})+{\frac {f''(\xi )}{2!}}ih,}$ 

注意到左邊的量是有限差分法的近似，右邊的量是待求解的量再加上一個餘數，因此餘數就是局部截尾誤差。上述範例可以用下式表示：

${\displaystyle {\frac {f(x_{0}+ih)-f(x_{0})}{ih}}=f'(x_{0})+O(h).}$ 

在此例中，局部截尾誤差和時間格點的大小成正比。

範例：常微分方程

例如考慮以下的常微分方程

${\displaystyle u'(x)=3u(x)+2.\,}$ 

利用數值方法中歐拉法求解，利用以下的有限差分式

${\displaystyle {\frac {u(x+h)-u(x)}{h}}\approx u'(x)}$ 

來近似導數，並配合一些代數處理（等號兩側同乘以h，再加上u(x)），可得

${\displaystyle u(x+h)=u(x)+h(3u(x)+2).\,}$ 

最後的方程式即為有限差分方程，求解此方程則可得到原方程的近似解。

範例：熱傳導方程

考慮正規化的一維熱傳導方程式，為齊次的狄利克雷邊界條件

${\displaystyle U_{t}=U_{xx}\,}$ 

${\displaystyle U(0,t)=U(1,t)=0\,}$ （邊界條件）

${\displaystyle U(x,0)=U_{0}(x)\,}$ （初始條件）

對此問題求數值解的一種方式是用差分去近似所有的導數，可以將空間分割為${\displaystyle x_{0},...,x_{J}}$ ，將時間也分割為${\displaystyle t_{0},....,t_{N}}$ 。假設在時間及空間都是均勻的網格切割，空間中兩個連續位置的間隔為h，兩個連續時間之間的間隔為k。點

${\displaystyle u(x_{j},t_{n})=u_{j}^{n}}$ 

表示${\displaystyle u(x_{j},t_{n})}$ 的數值近似解。

顯式方法

 

熱傳導方程最常用顯式方法的模版（英语：Stencil (numerical analysis)）

利用在時間${\displaystyle t_{n}}$ 的前向差分，以及在位置${\displaystyle x_{j}}$ 的二階中央差分（FTCS 格式（英语：FTCS scheme）），可以得到以下的迭代方程：

${\displaystyle {\frac {u_{j}^{n+1}-u_{j}^{n}}{k}}={\frac {u_{j+1}^{n}-2u_{j}^{n}+u_{j-1}^{n}}{h^{2}}}.\,}$ 

這是用求解一維導熱傳導方程的顯式方法。

可以用以下的式子求解${\displaystyle u_{j}^{n+1}}$ 

${\displaystyle u_{j}^{n+1}=(1-2r)u_{j}^{n}+ru_{j-1}^{n}+ru_{j+1}^{n}}$ 

其中${\displaystyle r=k/h^{2}.}$ 

因此配合此迭代關係式，已知在時間n的數值，可以求得在時間n+1的數值。${\displaystyle u_{0}^{n}}$ 及${\displaystyle u_{J}^{n}}$ 的數值可以用邊界條件代入，在此例中為0。

此顯式方法在${\displaystyle r\leq 1/2}$ 時，為數值稳定且收斂^[1]。其數值誤差和時間間隔成正比，和位置間隔的平方成正比：

${\displaystyle \Delta u=O(k)+O(h^{2})\,}$ 

隱式方法

 

隱式方法的模版

若使用時間${\displaystyle t_{n+1}}$ 的後向差分，及位置${\displaystyle x_{j}}$ 的二階中央差分（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}}}.\,}$ 

這是用求解一維導熱傳導方程的隱式方法。

在求解線性聯立方程後可以得到${\displaystyle u_{j}^{n+1}}$ ：

${\displaystyle (1+2r)u_{j}^{n+1}-ru_{j-1}^{n+1}-ru_{j+1}^{n+1}=u_{j}^{n}}$ 

此方法不論${\displaystyle r}$ 的大小，都數值稳定且收斂，但在計算量會較顯式方法要大，因為每前進一個時間間隔，就需要求解一個聯立的數值方程組。其數值誤差和時間間隔成正比，和位置間隔的平方成正比：

${\displaystyle \Delta u=O(k)+O(h^{2})\,}$ 

克兰克－尼科尔森方法

若使用時間${\displaystyle t_{n+1/2}}$ 的中間差分，及位置${\displaystyle x_{j}}$ 的二階中央差分（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).\,}$ 

此公式為克兰克－尼科尔森方法（Crank–Nicolson method）。

 

克兰克－尼科尔森方法的模版

在求解線性聯立方程後可以得到${\displaystyle u_{j}^{n+1}}$ ：

${\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}}$ 

此方法不論${\displaystyle r}$ 的大小，都數值稳定且收斂，但在計算量會較顯式方法要大，因為每前進一個時間間隔，就需要求解一個聯立的數值方程組。其數值誤差和時間間隔的平方成正比，和位置間隔的平方成正比：

${\displaystyle \Delta u=O(k^{2})+O(h^{2}).\,}$ 

若時間刻度較小時，克兰克－尼科尔森方法是最精確的，而顯式方法是最不精確的，而且可能會不穩定，但是是最容易計算的，其數值計算量也最少。若時間刻度較大時，隱式方法的效果最好。

相關條目

• 有限元分析
• 差分
• 時域有限差分
• 模版 (數值分析)（英语：Stencil (numerical analysis)）
• 有限差分係數
• 五點Stencil（英语：Five-point stencil）
• Lax等價定理（英语：Lax equivalence theorem）
• 期权定价的有限差分法（英语：finite difference methods for option pricing）

參考資料

1. Crank, J. The Mathematics of Diffusion. 2nd Edition, Oxford, 1975, p. 143.

• K.W. Morton and D.F. Mayers, Numerical Solution of Partial Differential Equations, An Introduction. Cambridge University Press, 2005.
• Oliver Rübenkönig, , (2006) Albert Ludwigs University of Freiburg
• Autar Kaw and E. Eric Kalu, Numerical Methods with Applications, (2008) [1] （页面存档备份，存于互联网档案馆）

外部連結

• Finite Difference Method of Solving ODEs (Boundary Value Problems) Notes, PPT, Maple, Mathcad, Matlab, Mathematica （页面存档备份，存于互联网档案馆）
• Shih-Hung Chen, National Central University
• Randall J. LeVeque（英语：Randall J. LeVeque）, Finite Difference Methods for Ordinary and Partial Differential Equations （页面存档备份，存于互联网档案馆）, SIAM, 2007.
• Finite Difference Method for Boundary Value Problems （页面存档备份，存于互联网档案馆）
• Numerical Methods for time-dependent Partial Differential Equations （页面存档备份，存于互联网档案馆）