有限差分法

有限差分法（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}.$$

## 准确度及误差

近似解的误差定义为近似解及解析解之间的差值。有限差分法的两个误差来源分别是舍入误差及截尾误差（或称为离散化误差），前者是因为电脑计算小数时四舍五入造成的误差，后者则是用有限阶级数表示导数引起的误差。

<image>

在运用有限差分法求解一问题（或是说找到问题的近似解）时，第一步需要将问题的定义域离散化。一般会将问题的定义域用均匀的网格分割（可参考右图）。因此有限差分法会制造一组导数的离散数值近似值。

一般会关注近似解的局部截尾误差，会用大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)$的数值近似解。

### 显式方法

<image>

利用在时间$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) \,$$

### 隐式方法

<image>

若使用时间$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）。

<image>

在求解线性联立方程后可以得到$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