有限差分法

在数学中，**有限差分法**（**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*的階乘，*R**n*(*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)泰勒多項式的餘數項，也就是

:
  R_n(x_0 + h) = \frac{f^{(n+1)}(\xi)}{(n+1)!} (h)^{n+1}
, 其中 x_0 ,

可以找到局部截尾誤差的主控項，例如用前項差分法計算一階導數，已知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}, </math>

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

:<math> \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, *[https://web.archive.org/web/20100111163529/http://www.imtek.de/simulation/mathematica/IMSweb/imsTOC/Lectures%20and%20Tips/Simulation%20I/FDM_introDocu.html The Finite Difference Method (FDM) - An introduction]*, (2006) Albert Ludwigs University of Freiburg
* Autar Kaw and E. Eric Kalu, *Numerical Methods with Applications*, (2008) http://www.autarkaw.com/books/numericalmethods/index.html 

## 外部連結
* 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
* , *[http://faculty.washington.edu/rjl/fdmbook/ 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