有限差分法

在数学中，有限差分法（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, The Finite Difference Method (FDM) - An introduction, (2006) Albert Ludwigs University of Freiburg
• Autar Kaw and E. Eric Kalu, Numerical Methods with Applications, (2008) [1] （页面存档备份，存于互联网档案馆）

外部链接

• 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
• Randall J. LeVeque（英语：Randall J. LeVeque）, 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 （页面存档备份，存于互联网档案馆）