有限差分法

有限差分法（finite-difference methods，FDM），是在數學中的一種微分方程數值方法，是通過有限差分來近似導數，從而尋求微分方程的近似解。

## 由泰勒展開式的推導

首先假設要近似函數的各級導數都有良好的性質，依照泰勒定理，可以形成以下的泰勒展開式：

$$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的階乘，Rn(x)為餘數，表示泰勒多項式和原函數之間的差。可以推導函數f一階導數的近似值：

$$f(x_0 + h) = f(x_0) + f'(x_0)h + R_1(x),$$

設定x0=a，可得：

$$f(a+h) = f(a) + f'(a)h + R_1(x),$$

除以h可得：

$${f(a+h)\over h} = {f(a)\over h} + f'(a)+{R_1(x)\over h}$$

求解f'(a)：

$$f'(a) = {f(a+h)-f(a)\over h} - {R_1(x)\over h}$$

假設$R_1(x)$相當小，因此可以將"f"的一階導數近似為：

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

## 準確度及誤差

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

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

一般會關注近似解的局部截尾誤差，會用大O符號表示，局部截尾誤差是指應用有限差分法一次後產生的誤差，因此為$f'(x_i) - f'_i$，此時$f'(x_i)$是實際值，而$f'_i$為近似值。泰勒多項式的餘數項有助於分析局部截尾誤差。利用$f(x_0 + h)$泰勒多項式的餘數項，也就是

$<pre> R_n(x_0 + h) = \frac{f^{(n+1)}(\xi)}{(n+1)!} (h)^{n+1} </pre>$, 其中$x_0 < \xi < x_0 + h$,

可以找到局部截尾誤差的主控項，例如用前項差分法計算一階導數，已知$f(x_i)=f(x_0+i h)$,

$$f(x_0 + i h) = f(x_0) + f'(x_0)i h + \frac{f''(\xi)}{2!} (i h)^{2},$$

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

$$\frac{f(x_0 + i h) - f(x_0)}{i h} = f'(x_0) + \frac{f''(\xi)}{2!} i h,$$

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

$$\frac{f(x_0 + i h) - f(x_0)}{i h} = f'(x_0) + O(h).$$

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

## 範例：常微分方程

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

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

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

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

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

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

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

## 範例：熱傳導方程

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

$$U_t=U_{xx} \,$$

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

$U(x,0) =U_0(x) \,$（初始條件）

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

$$u(x_j,t_n) = u_{j}^n$$

表示$u(x_j, t_n)$的數值近似解。

### 顯式方法

利用在時間$t_n$的前向差分，以及在位置$x_j$的二階中央差分（FTCS格式），可以得到以下的迭代方程：

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

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

可以用以下的式子求解$u_j^{n+1}$

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

其中$r=k/h^2.$

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

此顯式方法在$r\le 1/2$時，為數值穩定且收斂。其數值誤差和時間間隔成正比，和位置間隔的平方成正比：

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

### 隱式方法

若使用時間$t_{n+1}$的後向差分，及位置$x_j$的二階中央差分（BTCS 格式），可以得到以下的迭代方程：

$$\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}. \,$$

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

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

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

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

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

### 克蘭克－尼科爾森方法

若使用時間$t_{n+1/2}$的中間差分，及位置$x_j$的二階中央差分（CTCS 格式），可以得到以下的迭代方程：

$$\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）。

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

$$(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$$

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

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

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

## 參考資料

- K.W. Morton and D.F. Mayers, Numerical Solution of Partial Differential Equations, An Introduction. Cambridge University Press, 2005.
- Oliver Rübenkönig, The Finite Difference Method (FDM) - An introduction, (2006) Albert Ludwigs University of Freiburg
- Autar Kaw and E. Eric Kalu, Numerical Methods with Applications, (2008)

## 外部連結

- Finite-Difference Method (see and listen to lecture 9)
- List of Internet Resources for the Finite Difference Method for PDEs
- Finite Difference Method of Solving ODEs (Boundary Value Problems) Notes, PPT, Maple, Mathcad, Matlab, Mathematica
- Lecture Notes Shih-Hung Chen, National Central University
- , Finite Difference Methods for Ordinary and Partial Differential Equations , SIAM, 2007.
- Finite Difference Method
- Finite Difference Method for Boundary Value Problems
- Finite Difference Methodology in Materials Science
- Numerical Methods for time-dependent Partial Differential Equations