在前面的笔记里孤光一点萤：数值常微分方程-欧拉法与龙格-库塔法​zhuanlan.zhihu.com

整理了常微分方程的一些数值解法。类似的方法可以拓展到解偏微分方程的问题，这里整理有限差分法的相关笔记。

计算机无法处理连续介质的无限多点，因此只能通过网格划分的方式，将带求解区域内连续分布的函数进行离散化，这样只用求解有限多的格点即可。同样地，在计算机里无法取极限求导数，只能通过离散化之后，用有限差商代替导数。由此，可以将带求解的偏微分方程转化为一组离散格点上的差分方程。解这组方程的问题，最终又归化为一个矩阵的求解问题，而矩阵求解这种“体力活”是计算机所擅长的。

由于我们使用划分的离散网格来近似连续的待求解空间，因此网格划分越密，求解越精确。同时，网格划分越密集，意味着最后得到的矩阵方程维度越高，计算量会非常庞大。这时候常用的高斯消元法(一个高斯消元的例子)往往会吃不消，需要使用松弛迭代法等。

首先，还是从单变量的方程开始。考虑函数

，其中

。将区间

分为N份，每一份就是步长

，则离散化后的格点坐标为：

。对于任意格点

：

一阶向前差分

一阶向后差分

一阶中心差分

由此可得：

一阶向前差商

一阶向后差商

一阶中心差商

二阶中心差商

现在拓展到二维情况：

例如对于方程

，待求解区域为正方形域

将待求解区域划分为矩形网格，每一个小矩形的长和宽(步长)分别为

和

，使用中心差商，可得：

对于这里的正方形区域，自然而然可以选取

，代入方程，得：

以此类推，可以得到对于点

，有：

这便是点

满足的差分方程，由于该方程包含点

及其上下左右一共5个点，也称为“五点差分格式”。

五点差分格式的截断误差为

阶。如果把左上、左下、右上、右下四个点也用上，为“九点差分格式”，误差可以降低到

阶。

为简便起见，这里以拉普拉斯方程

为例，则五点差分格式简化为：

在边界

上，

在划分网格时候要注意，如果划分网格为

个正方形格子，那么待求解的点有

个，因为区域边界上的点(即边界条件)是给定的，不用求。因此，给网格编号时候，从0开始编号再合适不过了，如下：

考虑第一个点

，根据五点差分格式，它的值由上下左右四个点确定，很显然，上、左分别为边界点

和

，右、下分为为点

和

，写为差分格式为：

第一列其他点同理，以此类推，如果将网格第

列上每个点的解记为：

，则最终可以写为：

其中

为

阶单位矩阵，

这里出现了一个迭代方程，想解

需要

，解

需要

......

然后把这一系列矩阵方程写成一个更大的矩阵方程：

把整个网格的每一列纵向连接为一个列向量：

则有：

，其中：

，

终于，我们把要求解的问题化简为了矩阵方程

的形式，可以方便地使用计算机求解，这里A是系数矩阵，列向量

便是问题的解。

对于上面的拉普拉斯方程，用Python求解：

为了方便，求解时直接使用了线性代数库

，这样计算的是矩阵方程的精确解，耗时有点长。可以用松弛迭代法来求数值解，有空再更。

import numpy as np

import matplotlib.pyplot as plt

from mpl_toolkits.mplot3d import Axes3D

N=79 # 矩阵维数

Imat = -np.eye(N)/4

def Gmatrix(N): # 生成G矩阵

mtr = np.eye(N)

for i in range(0, N - 1):

mtr[i][i + 1] = -1 / 4

for i in range(1, N):

mtr[i][i - 1] = -1 / 4

return mtr

def zero_mat(N,m): # 生成m个横向连接N阶零矩阵,m>=1

zero_matrix = np.zeros((N,N))

for i in range(m-1):

zero_matrix = np.block([zero_matrix,np.zeros((N,N))])

return zero_matrix

def zero_arr(N, m): # 生成m个纵向连接N阶列向量,m>=1

zero_array = np.zeros((N,1))

for i in range(m-1):

zero_array = np.block([[zero_array],[np.zeros((N,1))]])

return zero_array

def Kmatrix(N,Gmat,Imat):

Kmar_list = []

for i in range(N):

if i == 0:

Kmat_line = np.block([np.block([Gmat, Imat]),zero_mat(N,N-2)])

Kmar_list.append(Kmat_line)

elif i == 1:

Kmat_line = np.block([Imat, Kmar_list[-1][:,:-N]])

Kmar_list.append(Kmat_line)

else:

Kmat_line = np.block([np.zeros((N,N)), Kmar_list[-1][:,:-N]])

Kmar_list.append(Kmat_line)

Kmat = Kmar_list[0]

for j in range(N-1):

Kmat = np.block([[Kmat],[Kmar_list[j+1]]])

return Kmat

def Barray(N): # 方程右边

barr = np.ones((N,1))/4

barr = np.block([[barr],[zero_arr(N,N-2)]])

barr = np.block([[barr],[np.ones((N,1))/4]])

return barr

Gmat = Gmatrix(N)

Kmat = Kmatrix(N,Gmat,Imat)

Barr = Barray(N)

Phi = np.linalg.solve(Kmat,Barr)

def trans(N,Phi):

# 把线性方程组解得的列矢量分割，从左到右排列转化为方矩阵，方便可视化

phi = Phi[0:N]

for i in range(N-1):

j = i+1

phi = np.block([phi,Phi[(j*N):((j+1)*N)]])

return phi

Phi = trans(N,Phi)

x = np.linspace(0,1,N)

y = np.linspace(0,1,N)

X,Y = np.meshgrid(x,y)

fig = plt.figure()

ax = Axes3D(fig)

ax.plot_surface(X,Y,Phi,rstride=1,cstride=1,cmap=plt.get_cmap('rainbow'))

ax.contourf(X,Y,Phi,zdir='z',offset=1,cmap='rainbow')

plt.title('Difference methods for Laplace equation',fontsize='12')

plt.show()

作图如下：

通过控制作图参数里的offset，可以调整等高线的位置：

差分法求解的误差主要来自于差分公式对偏微分的近似替代，为了提高求解精度，可以把网格划分得更密，同时计算量也会显著增大。

Reference：

1、马文淦，《计算物理学》

2、李荣华，刘播，《微分方程数值解法》