1. 引言

有限差分法（Finite Difference Method，FDM）是一种求解微分方程数值解的近似方法，其主要原理是对微分方程中的微分项进行直接差分近似，从而将微分方程转化为代数方程组求解。有限差分法的原理简单，粗暴有效，最早由远古数学大神欧拉（L. Euler 1707-1783）提出，他在1768年给出了一维问题的差分格式。1908年，龙格（C. Runge 1856-1927）将差分法扩展到了二维问题【对，就是龙格-库塔法中的那个龙格】。但是在那个年代，将微分方程的求解转化为大量代数方程组的求解无疑是将一个难题转化为另一个难题，因此并未得到大量的应用。随着计算机技术的发展，快速准确地求解庞大的代数方程组成为可能，因此逐渐得到大量的应用。发展至今，有限差分法已成为一个重要的数值求解方法，在工程领域有着广泛的应用背景。本文将从有限差分法的原理、基本差分公式、误差估计等方面进行概述，给出其基本的应用方法，对于一些深入的问题不做讨论。

2. 有限差分方法概述

首先，有限差分法是一种求解微分方程的数值方法，其面对的对象是微分方程，包括常微分方程和偏微分方程。此外，有限差分法需要对微分进行近似，这里的近似采取的是离散近似，使用某一点周围点的函数值近似表示该点的微分。下面将对该方法进行概述。

2.1. 有限差分法的基本原理

这里我们使用一个简单的例子来简述有限差分法的基本原理，考虑如下常微分方程

微分方程与代数方程最大的不同就是其包含微分项，这也是求解微分方程最难处理的地方。有限差分法的基本原理即使用近似方法处理微分方程中的微分项。为了得到微分的近似，我们最容易想到的即导数定义

上式后面的近似表示使用割线斜率近似替代切线斜率，Δx 即为步长，如图 1(a)所示。式(2)表明函数在某一点的微分可以由相邻点的函数值近似确定。显然，这里微分近似的精度与步长的选取有关，步长越短则越精确。

因此，这里首先需要对求解域进行离散，然后分别得到各离散点上的微分近似。对于示例的一维问题，将求解区间等分为 N 个区间，步长为 h，分别将包含首尾的各结点记为 x0,x1,x2,...,xN，对应的函数值为 u0,u1,u2,...,uN，对应各点的一阶微分记为 u0′,u1′,u2′,...,uN′，如图 1(b) 所示。这样，就把原问题的求解转化为了各结点函数值 ui 的求解。式(2)的离散形式表示为

将式（3）代入方程（1）可以得到

这里的结点坐标 xi=a+ih,(i=0,1,2,3,4)，步长 h=(b-a)/N 均为已知。记 ci=c(xi),fi=f(xi)， 将式(4)合并同类项可以得到如下递推关系

上式共有 N 个方程（i=0,1,2,...,N-1），包括 N 个未知数（u1,u2,...,uN），刚好可以求解得到各结点上的待求函数的值，从而得到原问题在求解域上的近似解。由于该问题中初值已给定，按照各结点依次迭代就可以得到该问题的解。此外，为了给出更一般化的求解方法，可以将(5)写成矩阵的形式，即

即

其中 u={u1,…,uN-1,uN}T 总共包含 N 个待求结点值。Ci=hci-1，Fi=hfi 。上述方程组的系数矩阵 A 显然可逆，则上述问题的解总是可以表示为 u=A-1F。为了验证上述结果的正确性，我们取 c(x)=1, f(x)=sin(x)+cos(x), 求解区间为 x∈[0,2π]，且满足边界条件 u(0)=0，则原问题(1)可以写为如下形式

则该问题对应方程组(6)中的 Ci=h-1，Fi=h[sin(xi)+cos(xi)]，d=0。将上述表达式代入式(7)并利用 u=A-1F 即可以得到该问题的近似解。此外，上述方程的解析解可以容易给出为 u(x)=sin(x)。分别取 N=5,10,20,100 时有限差分法结果与解析解对比如图(2)所示，可以看到网格划分越精细则求解结果越精确，与解析解的误差越小。当 N=100 时，两者的误差已经极小，表明有限差分法是有效的。

求解该问题时的Python计算代码如下：

from numpy import *

h = 2*pi/N
x = linspace(0, 2*pi, N+1)    
A = zeros((N, N))
F = zeros((N, 1))
u = zeros((N+1, 1))
for i in range(N):
    A[i, i] = 1
    if i > 0: A[i, i-1] = h-1 
    F[i] = h*(sin(x[i]) + cos(x[i]))
u[1:] = dot(linalg.inv(A), F)

2.2. 微分的近似表示

上一节以一个一维一阶微分方程为例简单描述了有限差分的基本原理，下面我们对其更广泛的使用进行扩展。注意到式（3）中的一阶微分是使用相邻结点的函数值来表示，一般地，假设函数 u(x) 在求解域内连续可导，则由泰勒级数我们可以得到

基于泰勒级数并做适当的截断，我们就可以得到各阶微分的近似表达式，或者记做“差分公式”。下面主要对一阶和二阶微分的差分公式进行讨论，更高阶的微分可以同理推导得到。

2.2.1. 一阶微分

我们将(8)式另写作

这里 ξ∈(0,h) 。对上式变形即可以得到一阶微分的向前差分公式：

将(10)式中的 h 用 -h 替代，则可以得到一阶微分的向后差分公式：

联立式(11)与(12)，则可以得到一阶微分的中心差分公式：

式(11)-(13)是几种常用的一阶微分的差分公式，式(11)-(13)中右边的最后一项即为截断误差。由于这里的误差与步长相关，因此可以看到向前和向后差分公式具有一阶近似精度 [O(h)]，而中心差分则具有二阶近似精度 [O(h2)]，三种差分公式的区别如图(3)所示。

2.2.2. 二阶微分

类似地，将式(8)多写几项表示为

类似的，在式(14)中令 h→-h，并联立原式可以得到二阶微分的中心差分公式

这里列出的差分公式仅仅只是一些常用的形式，不同的差分公式具有不同的求解误差，会产生不一样的求解效果。对于一般问题，由于中心差分公式的精度较高，因此使用较多。

2.2.3. 差分公式的推导

前面给出了常用一阶和二阶差分公式，而实际问题中经常需要给定差分项的差分公式，这时可以采用待定系数法来给出。例如，对于一阶微分，我们想要在差分公式中同时包含 u(x)，u(x+h)，u(x+2h) 三项，即需要同时使用 ui，ui+1 和 ui+2 来表示 ui′，则令

由泰勒公式

将式(17)代入(16)并对比系数可以得到

即

这样就得到了指定差分项的差分公式。通过计算可知，上式给出的差分公式具有二阶近似精度，比向前差分与向后差分公式的精度都要高。同样地，其它差分公式也可以通过类似的方法推导得到。通过上述方法我们可以根据实际求解问题的不同给出不同的差分公式以达到最高的求解效率。

2.3. 偏微分方程和时域问题的有限差分法

通常工程中遇到的许多问题都是偏微分方程，即包含两个及以上自变量变化的问题。对于偏微分方程，有限差分法也可以进行类似处理。对于二维、三维或者时域问题，其处理思路与方法类似，这里为了简便，我们以二维问题为例来进行讨论，三维问题和时域问题可以同理进行推广。对于二维偏微分方程，其包含两个自变量的变化，求解域为面域。我们还是以一个实例来演示其基本方法，比如，考虑如下二维泊松方程

首先，将求解域进行离散化。由于求解域为矩形，分别将 x 和 y 方向等分为 M 与 N 份，从而构建网格，两个方向的离散坐标分别记为 xi(i=0,1,2,..,M) 和 yj(j=0,1,2,..,N)， 各结点的函数值简记为 u(xi,yj)=ui,j，如下图所示。

进一步，对方程(20)中的二阶微分使用中心差分公式有

其中 hx=a/M 与 hy=b/N 分别为 x 方向与 y 方向的步长。记 p=1/hx2，q=1/hy2，r=-2(p+q)，f(xi,yj)=fi,j，将式(21)代入原方程得到

递推式(22)表明每一次迭代中有4个结点值是相互独立的，即ui+1,j,ui-1,j,ui,j,ui,j+1,ui,j-1，这5个点组成一个菱形，根据其中任意4个结点值可以算出最后一个结点的值。由此，不断迭代可以求得所有结点的值。根据所划分的网格，除去所求矩形区域的四个顶点（由边界条件给出），图4网格中的所有结点将会被菱形覆盖，根据边界条件可以依次向内求解得到全域的解。类似地，可以将上述差分方程写作矩阵形式

这里，根据边界条件，四个顶点的值均为常数，此外，向量 u→i,0，u→i,M，u→0,j 和 u→M,j 均由边界条件给定。可以看出，相比于常微分方程，由于求解结点的增多，偏微分方程的求解复杂度大大增加。同样地，为了验证计算结果，我们取 f(x,y)=-2π2sin(πx)sin(πy)，c(y)=d(y)=g(x)=h(x)=0，(x,y)∈[0,1]×[0,1]，使用有限差分方法对其进行求解。同时，其对应的解析解为 u(x,y)=sin(πx)sin(πy)，在求解网格为 100×100 时两者的结果对比如图5所示。

可以看到有限差分法求解得到的结果与解析解基本一致。我们选取 x 与 y 方向相同的网格密度，即令 M=N=S，得到其最大相对误差随网格密度 S 的变化如图6所示。可以看到在 50×50 计算网格下的计算误差已经可以忽略不计，因此对该问题的有限差分法求解是可靠的。

对应该问题的Python求解代码如下：

from numpy import *

M, N = 100, 100
a, b = 1, 1  
hx, hy = a/M, b/N
p, q = 1/hx**2, 1/hy**2
r = -2*(p + q)

U = zeros((M-1, M-1))
for i in range(M-1):
    U[i, i] = r
    if i < M-2: U[i, i+1] = p
    if i > 0:   U[i, i-1] = p 
V = diag([q]*(M-1))
Zero_mat = zeros((M-1, M-1))

A_blc = empty((N-1, N-1), dtype=object) # 矩阵A的分块形式
for i in range(N-1):
    for j in range(N-1):
        if i == j:
            A_blc[i, j] = U
        elif abs(i-j) == 1:
            A_blc[i, j] = V
        else:
            A_blc[i, j] = Zero_mat
            
A = vstack([hstack(A_i) for A_i in A_blc]) # 组装得到矩阵A

x_i = linspace(0, a, M+1)
y_i = linspace(0, b, N+1)
F = vstack([-2*pi**2 * sin(pi*x_i[1:M].reshape((M-1, 1))) * sin(pi*j) for j in y_i[1:N]])

u = dot(linalg.inv(A), F).reshape(M-1, N-1)
u_f = vstack([zeros((1,M+1)),  # 最后组装边界条件得到全域的解
              hstack([zeros((N-1,1)), u, zeros((N-1,1))]), 
              zeros((1,M+1))])

3. 应用实例
根据上面给出的有限差分法的基本方法，这里给出一个具体的应用实例，即枝晶固化生长的相场模拟，使用的也是经典的 Kobayashi 模型[1]。该问题也是一个典型的相变问题，枝晶固化生长即由液态变为固态的过程。使用变量 p(r,t)∈[0,1] 来表示不同相，这里 r,t 分别为位置变量和时间变量。 并且，p=0 表示液相（未相变）， p=1 表示固相（已相变）。其中变量 p 满足 Ginzburg–Landau 方程

这里 τ 为小常数，δ 为变分符号；Φ 为自由能，表达式为：

其中 ε 为梯度能系数，与界面的厚度相关。为了表征界面的各向异性，假设 ε 为界面外法线方向 v 的函数，即 ε(v)=ε(-p)，将式 (26) 带入 (25) 可以得到

这里的推导用到了如下泛函的变分公式：若

则

这里我们以二维的情况为例进行讨论，引入偏移角 θ 表示 v 与 x 轴正向的夹角，即

此时的梯度能系数 ε 可以进一步表示为

其中 σ 为界面各项异性的强度（strength of anisotropy）；J 为界面的各向异性振型数（mode number of anisotropy）； θ0 为初始偏移角。由此，式 (27) 可以进一步变形为

此外，自由能中的参数 m 是界面运动的驱动力，与局部温度 T 相关，

其中 α 为正常数，Teq 为平衡温度。同时，局部温度分布满足方程

这样，式(27)-(34)就给出了该问题的控制方程，包含两个变量 p 和 T 需要求解。需要说明的是，上述给出的控制方程均已经归一化，相关的参数选择为 τ=3×10-4， ε=0.01， σ=0.02， J=4， θ0=0.2， α=0.9， γ=10， Teq=1， κ=1.8。这里我们假设空间求解域为方形区域 r=[x,y]，x,y 方向上的分别划分为 M,N 段；时间域 t 上划分为 K 段，在各方向上的步长分别对应记为 Δx,Δy,Δt。初始 t=0 时，在求解域中间取一个区域作为固化的初始点，记为 SN，如图 7 所示；则在固化点 SN 区域上 p=1，在其余区域 p=0,T=0，并随着时间增长枝晶从固化点开始生长。

对于这个问题，由于全域的初始值已给定，因此使用迭代法来求解各时间步内的解比较方便。为了方便，使用如下记号：pi,jk=p(xi,yj,tk)，其中 i,j,k 分别对应 x,y,t 方向上的结点。下面使用有限差分法对控制方程 (27) 与 (34) 进行离散，其中时域使用经典欧拉显示格式来离散，即

空间域则使用前面讲到的一般性离散方法，例如

则原控制方程的迭代方程可以推导为：

取 M=N=300,K=4000;Δx=Δy=0.03,Δt=1×10-4，求解得到枝晶的生长如图 8 所示，该结果与文献中的实验结果一致，表明有限差分法在该问题的求解中是有效的。

可以看到，枝晶边缘的变化曲线比较模糊，这是由于我们划分的网格相对比较粗糙。这里为了计算速度，我们并未取较细的网格。枝晶生长过程的动画如图 9 所示。需要注意的是，在枝晶生长的前10计算步中出现了一些 p<0 的情况，这里为了方便直接取这些元素为零。

该问题的Python求解代码如下：

from numpy import *

# Space and time domain
M,  N,  K   = 300,  300,  4000
Dx, Dy, Dt  = 0.03, 0.03, 1e-4

# Material Properties
tau = 3e-4
eps_bar = 0.01
sigma = 0.02
J = 4.
theta_0 = 0.2
alpha = 0.9
gamma = 10.
T_eq = 1.
kappa = 1.8

#%% Evolution
p = zeros((M, N))
T = zeros((M, N))
# Initial Solidification area
for i in range(M):
    for j in range(N):
        if (i - M/2)**2 + (j - N/2)**2 < 5.0:
            p[i, j] = 1.0

# Define Laplacian operator
def Lap(p):
    p_i_j  = delete(delete(p, [0, -1], axis=0), [0, -1], axis=1)
    p_im_j = delete(delete(p, [0, -1], axis=0), [-1,-2], axis=1)
    p_ip_j = delete(delete(p, [0, -1], axis=0), [0,  1], axis=1)
    p_i_jm = delete(delete(p, [0, -1], axis=1), [0,  1], axis=0)
    p_i_jp = delete(delete(p, [0, -1], axis=1), [-1,-2], axis=0)
    Lap_p  = (p_im_j + p_ip_j + p_i_jm + p_i_jp - 4*p_i_j)/Dx**2
    Lap_pj = vstack((Lap_p[0,:], Lap_p, Lap_p[-1,:]))
    return hstack((Lap_pj[:,0].reshape(N,1), Lap_pj, Lap_pj[:,-1].reshape(N,1)))

# Phase field evolution
def Phase_field(p, T):
    theta = arctan2(gradient(p, Dy, axis=1), gradient(p, Dx, axis=0))
    eps = eps_bar * (1. + sigma * cos(J * (theta - theta_0)))
    g = -eps * eps_bar * sigma * J * sin(J * (theta - theta_0)) * gradient(p, Dy, axis=1)
    h = -eps * eps_bar * sigma * J * sin(J * (theta - theta_0)) * gradient(p, Dx, axis=0)
    m = alpha/pi * arctan(gamma * (T_eq - T))
    term_1 = - p*(p - 1.)*(p - 0.5 + m)
    term_2 = - gradient(g, Dx, axis=0)
    term_3 = gradient(h, Dy, axis=1)
    term_4 = eps**2 * Lap(p)
    p_ev = Dt / tau * (term_1 + term_2 + term_3 + term_4)
    return p + p_ev

# Temperature evolution
def Temp(T, p_new, p_old):
    T_ev = Dt*Lap(T) + kappa*(p_new - p_old)
    return T + T_ev

# Evolution process
p_hist = []
T_hist = []
p_old = p; T_old = T
for t_step in range(K):
    p_new = Phase_field(p_old, T_old)
    T_new = Temp(T_old, p_new, p_old)
    p_old = p_new
    T_old = T_new
    if t_step % 50 == 0:
        p_hist.append(p_new)
        T_hist.append(T_new)
        print('step finished:', t_step,'/',str(K))