抛物线方程的有限差分法-Dirichlet条件下

偏微分方程求解问题中，是否有解析解很大程度上受到边界条件、初始条件以及方程形式的影响，在复杂边界条件下，寻求解析解是不可能完成的任务，因此，偏微分方程的数值解法具有很强的实用价值，对于物理场问题的描述，离不开数值解法。PDE的数值求解方法中，有限差分法是最基本的、也是求解思路最为清晰的一种数值方法，这一方法构造简洁，是在对场问题描述、理解的基础上得来的；在一维问题中，FDM要比FEM有着更高的精度和计算效率。对于PDE的数值求解，以FDM作为学习的开始，是一个不错的选择。

有限差分形式

PDE是是一种特殊的方程，在方程中往往包含着场函数的导数，而最终的求解是需要得到场函数本身，那么从结果出发来看，对于一个一维的非稳态问题来说，场函数的是一个二元函数，对于它的数值描述，实际可以这样理解，即在这个域上，离散后会得到许多个离散的点，每一个点的坐标实际上是已知的，这样，对于这个二元函数来说，它在每一个离散点处函数值，就是我们需要的数值解的描述，一般我们更愿意依据离散点的分布形式，将它的函数值放在一个同形状的矩阵中。数值方法是一种尽可能近似的、有限度的描述场函数分布的方法，理论上的真解是连续函数，而数值解描述的是区域上某点处的场函数值，因此，数值方法的实现往往第一步就是需要对问题的离散，有限差分法也不例外。

这里针对的是二阶线性偏微分方程，并且是各向同性的，即是一个均匀方程，那么不管怎样，PDE一定会有导数项，如果能将导数表示为场函数，那么方程就完全转化为了场函数本身，如果能够进行求解，那么所得结果正好就是我们所需要的数值解。对于导数差分形式的理解，从泰勒展开的角度来看，个人感觉这是一个比较清晰和简洁的方式。一个函数的泰勒展开形式为：

$$ f(x+\Delta x) = f(x)+\frac{1}{1!}f'(x)(\Delta x)^1+\frac{1}{2!}f''(x)(\Delta x)^2+o(x^3) $$

泰勒展开式提供了一种思路，即函数在$x + \Delta x$处的函数值可以利用函数在$x$处的函数值以及这个函数的导数去估计，当略去高阶项时，可以得到函数一阶导的表示方式，

$$ f'(x)=\frac{f(x+\Delta x) - f(x)}{\Delta x} $$

这就是一阶导数的差分形式，此时为两点前向差分公式，它将函数的一阶导数表示为两点处的函数值。

当采用点$x$和点$x - \Delta x$这两个点来构造时，得到的一阶导数的两点后向差分公式，后向差分形式是一种很重要的形式，虽然好像由此得到的方程的差分形式没有前向的简洁，但是它所得的差分迭代关系式在数值计算中通常是无条件稳定的，这一点复杂度的代价换来的求解优势还是很值得的。

$$ f'(x)=\frac{f(x)-f(x-\Delta x)}{\Delta x} $$

上述是一阶导数的差分公式，而PDE中二阶导数是最为常见的，通常需要三个点的函数值来描述二阶导数的值，当略去二阶以上的项时，$x+\Delta x$ 、$x-\Delta x$分别对应的泰勒展开为

$$ f(x + \Delta x) = f(x) + f'(x)\Delta x + \frac{1}{2}f''(x)(\Delta x)^2\\ f(x - \Delta x) = f(x) - f'(x)\Delta x + \frac{1}{2}f''(x)(\Delta x)^2 $$

上述两式相加，消去一阶导项，可以得到二阶导的三点中心差分公式：

$$ f''(x) = \frac{f(x-\Delta x) -2f(x)+f(x+\Delta x)}{(\Delta x)^2} $$

以上是最基本的差分公式，多数问题依据它们便可以得到解决，当然，可以通过泰勒展开来构造出更为复杂的形式。

抛物线方程的有限差分形式

当PDE中的导数项利用差分形式进行表示后，此时方程中是不含有导数项的，完全是由场函数值来组成的。而不论是几维问题，偏导数永远都是针对某一方向而言的，因此，前述的差分形式都是通用的，为了便于计算，可以这样定义：$x、x + \Delta、、x - \Delta$可以理解为中间点，后一个点、前一个点，一个PDE问题，其空间域往往是得到确定的，一维问题中x的长度是已知的，当采用均匀的划分后，$\Delta x$实际上已经得到了确定，如果说$x$方向上，区间左端点为$x_0$，那么下一个点记为$x_1$，其坐标必定为$x_0 + 1 \cdot \Delta x$，以此类推得到$x_i =x_0 + i \cdot \Delta x$，为了表达式的简洁，我们可以直接用点号来表示点，即$u(x_i)$记作$u(i)$。

以抛物线方程为例，通过上述的表示方式来构造出PDE对应的有限差分形式。PDE形式如下：

$$ D\frac{\partial ^ 2 u(x, t)}{\partial x ^ 2} = \frac{\partial u(x,t)}{\partial t} $$

这个方程是一维抛物线型方程的最简形式，它通常的物理意义可以表示一维的非稳态温度场、一维的固结问题等。对于其中的二阶导采用三点中心差分公式来替换，一阶导采用前向差分公式，且对于$x$方向来说点编号记作$i$，时间$t$方向上点编号记作$j$，得到

$$ D \cdot \frac{u(i-1, j) - 2u(i, j) + u(i+1, j)}{(\Delta x)^2} = \frac{u(i,j+1) - u(i,j)}{\Delta t} $$

该方程是一个非稳态问题控制方程，可以理解为场函数在前一时刻的分布情况下得到下一时刻场的分布，当位于初始时刻时，由于初始条件的存在，便可以得到下一个时刻的状态分布，而每一个时刻中，边界条件都是需要强制满足的，我们只需要将$j$时刻与$j+1$时刻的项分别放在等号左右两侧即可反映出这一种计算思路。令网格比为$\sigma = D \frac{\Delta t}{(\Delta x)^2}$

$$ u(i, j+1) = \sigma u(i-1,j)+(1-2\sigma)u(i,j)+\sigma u(i+1,j) $$

这是前向差分得到的迭代式，这里先不写出矩阵形式，因为矩阵形式需要考虑到边界条件的影响，确切地来说，边界条件的形式会决定自由点的个数，导致矩阵的大小不同。

对于时间的导数，采用后向差分法时，得到后向差分形式为：

$$ D \cdot \frac{u(i-1, j) - 2u(i, j) + u(i+1, j)}{(\Delta x)^2}=\frac{u(i,j) - u(i,j-1)}{\Delta t}\\ -\sigma u(i-1,j) +(2\sigma +1)u(i,j)-\sigma u(i+1,j) = u(i,j-1) $$

对比前向差分形式与后向差分形式可以发现，前向差分形式实际上是一个显式形式，因为前一时刻的状态都是已知的，所以下一时刻的$i$处的值由前一时刻的3个点完全可以得到确定；而对于后向形式而言，却恰好相反，已知的一个点与下一时刻的3个节点值有关，这是一种隐式。前向形式简洁明了，但是存在一个较大的隐患，即对网格比有要求，为了迭代的收敛，$\sigma$必须小于等于0.5才行，这样使得时间步长$\Delta t$与空间步长$\Delta x$的取值都提出了限制，这会极大地限制数值解的描述范围。而后向差分形式得到的迭代关系式是一个隐式，它是一种无条件稳定的，所以对于问题的解的描述自由度很大；造成这一种结果的根本原因是迭代过程中矩阵的谱半径必须要小于1才能够保证误差不会被放大。从求解精度来看，不论是前向还是后向，它们的精度都是一样的，为$O(\Delta t + (\Delta x)^2$。

除了上述的两种形式外还有很多中差分形式，但是最为著名和常用的应该是Crank-Nicolson法，它是显式与隐式的一种结合体，对于空间项而言，该方法采用了平均的思想，即取$j-1$时刻与$j$时刻的二阶导的平均值来表示二阶导项，而时间项依然采用后向差分形式。

$$ D(\frac{u(i-1,j-1) - 2 u(i, j - 1) + u(i+1, j-1)}{(\Delta x)^2} + \frac{u(i-1,j) - 2u(i, j)+u(i+1, j)}{(\Delta x)^2})\\ =\frac{u(i, j) - u(i, j-1)}{\Delta t} $$

将$j-1$与$j$项分开，得到Crank-Nicolson法的迭代关系式为：

$$ -\sigma u(i-1,j)+(2 \sigma +2) u(i, j) - \sigma u(i+1, j) = \sigma u(i-1,j-1)+ (2-2 \sigma)u(i, j-1)+ \sigma u(i+1, j-1) $$

C-N法也是一种隐式形式，所以需要通过方程组的求解才能确定下一个时刻的状态值，该形式下精度为$O((\Delta t)^2+ (\Delta x)^2)$，精度比后向、前向差分形式要高，并且也是无条件稳定的。

边界条件的引入

PDE的求解问题一般也称为边值问题BVP（稳态）或者IBVP初边值问题（非稳态），实际上这反映了PDE的求解的核心是边界条件、初始条件。PDE反映了一类具有共性的物理场问题，而边界条件、初始条件的引入使得问题转为某一特定的问题或者研究对象。

常见的边界条件有三种。第一类边界条件，也称为本质边界条件，亦称为Dirichlet条件，这一种边界条件反映了物理场状态量函数本身的取值情况，这一种边界条件的描述使得边界处的场函数值是已知的；第二类边界条件，即自然边界条件、Nuemann条件，它描述了场函数在边界处的导数值，但是边界处场函数的值来说是未知的。而有限差分在得到差分形式后，所要求解的就是各个差分点处的场函数值，而Nuemann条件使得在边界处的差分点场函数值与内部的点一样，都是未知的，所以都是自由点，这样从每一个时刻要求解的自由点的个数要比第一类边界条件下的要多1个。

后向差分在Dirichlet条件下的矩阵形式

对于一维空间的抛物线型PDE来说，当两个端点的边界条件均为Dirichlet条件时，如果空间步数为$n$那么，那么实际$x$方向的离散点的个数为$n+1$个，因为第一类边界条件的含义就是指定边界点的场函数取值，所以就整个过程而言，这些点的值都是已知的。由于后向差分、C-N形式都是隐式差分形式，无法显式地表示出下一时刻与上一时刻的关系，将他们的迭代关系写为矩阵形式，那么最终的求解就是反复地线性方程组求解。本质上，都是利用上一时刻的各点处的函数值去计算得到下一时刻的点的函数值。

后向差分形式时，迭代矩阵可以表示为：

$$ A \cdot u_{j} = u_{j - 1} + \sigma \cdot s_{j} $$

其中：

$$ A = \begin{bmatrix} 1 + 2 \sigma & -\sigma & 0 & \dots & 0\\ -\sigma & 1 + 2 \sigma & -\sigma & \ddots & \vdots\\ 0 & -\sigma & 1 + 2 \sigma & \ddots & 0\\ \vdots & \ddots & \ddots & \ddots & -\sigma\\ 0 & \dots & 0 & -\sigma & 1 + 2 \sigma \end{bmatrix} $$

$$ s_{j}= \begin{bmatrix} u_{0j}\\ 0\\ \vdots\\ 0\\ u_{nj} \end{bmatrix} $$

注意到，右侧的一维向量就是第一类边界条件的体现，采用$n$段划分时，第一个节点为0号节点，那么一共具有$n + 1$个节点，最后一个端点（右端点）为$n$号节点，那么自由点的点号范围时$1$~$n-1$，共$n - 1$个自由点。

Crank-Nicolson法的矩阵形式:

C-N法的迭代关系要比后向差分形式复杂，前一时刻的状态量本身也会成为一个矩阵，但是最终的计算都是归结于线性方程组的求解。

$$ A \cdot u_{j} = B \cdot u_{j-1} + \sigma \cdot s $$

其中：

$$ A = \begin{bmatrix} 2 \sigma + 2 & -\sigma & 0 & \dots & 0\\ -\sigma & 2 + 2\sigma & -\sigma & \ddots & \vdots\\ 0 & -\sigma & 2 + 2\sigma & \ddots & 0\\ \vdots & \ddots & \ddots & \ddots & -\sigma\\ 0 & \dots & 0 & -\sigma & 2 + 2\sigma \end{bmatrix} $$

$$ B = \begin{bmatrix} 2 - 2\sigma & \sigma & 0 & \dots & 0\\ \sigma & 2 - 2\sigma & \sigma & \ddots & \vdots\\ 0 & \sigma & 2 - 2\sigma &\ddots & 0\\ \vdots & \ddots & \ddots & \ddots & \sigma\\ 0 & \dots & 0 & \sigma & 2 - 2\sigma \end{bmatrix} $$

$s$向量用于表示第一类边界条件的影响。

$$ s = \begin{bmatrix} u_{0, j-1}\\ 0\\ \vdots\\ 0\\ u_{n, j-1} \end{bmatrix} + \begin{bmatrix} u_{0, j}\\ 0\\ \vdots\\ 0\\ u_{n, j} \end{bmatrix} $$

第一类边界条件也成为本质条件，它完完全全地针对场状态变量进行限定，因此，从求解问题的角度看，第一类边界条件使得边界处的点是已知的，它们会作为系数的形式参与到计算中。正式因为这一种特性，使得在不同类型的边界条件下，隐式差分形式的计算会些许有点不同，主要是影响到自由点的个数，Nuemann条件下的计算会在以后进行介绍。

实现

import numpy as np
# 边界条件为完全的Dirichlet条件
class Parabolic_CnMethod(object):
    def __init__(self, D, x_range, tmax, xn, tn):
        self.D = D
        self.x_range = x_range
        self.tmax = tmax
        self.xn = np.int32(xn)
        self.tn = np.int64(tn)
        self.dx = (x_range[-1] - x_range[0]) / self.xn
        self.dt = tmax / self.tn
        self.sigma = (self.D * self.dt) / self.dx**2
        self.x = np.arange(x_range[0], x_range[-1] + self.dx, self.dx)
        self.t = np.arange(0, tmax + self.dt, self.dt)
        self.left_matrix()
        self.right_matrix()
        
    def left_matrix(self):
        # 因为边界全部为第一类边界，所以两端边界处的点是已知的
        x_number = self.x.shape[0] - 2  # 自由点个数
        # 自由点的个数决定了矩阵的形状大小
        l_m0 = np.diag(np.ones(x_number) * (self.sigma*2 + 2), 0)
        l_ml = np.diag(np.ones(x_number - 1) * -self.sigma, -1)
        l_mu = np.diag(np.ones(x_number - 1) * -self.sigma, 1)
        self.l_m = l_m0 + l_ml + l_mu
        return self.l_m
    
    def right_matrix(self):
        r_m0 = np.diag(np.ones(self.x.shape[0] - 2) * (2 - 2*self.sigma), 0)
        r_ml = np.diag(np.ones(self.x.shape[0] - 2 - 1) * self.sigma, -1)
        r_mu = np.diag(np.ones(self.x.shape[0] - 2 - 1) * self.sigma, 1)
        self.r_m = r_m0 + r_ml + r_mu
        return self.r_m
    
    def solve(self, lbc, rbc, ic):
        # 初始化一个自由点的结果矩阵w
        w = np.zeros([self.t.shape[0],self.x.shape[0]-2])
        w[0, :] = ic
        s = np.zeros(self.x.shape[0] - 2)
        for iters in np.arange(self.tn):
            s[0] = lbc[iters] + rbc[iters + 1]
            s[-1] = lbc[iters] + rbc[iters + 1]
            b = np.dot(self.r_m, w[iters, :]) + self.sigma * s
            w[iters + 1, :] = np.linalg.solve(self.l_m, b)
        ww = np.c_[lbc, w, rbc]    
        
        return ww

针对的方程为：

$$ \frac{\partial ^2 u}{\partial x ^ 2}=\frac{\partial u}{\partial t} $$

其中$D=1$，边界条件为

$$ u(0, t) = 0\\ u(1, t)=0 $$

初始条件：$u(x, 0)=100, (0<x<1)$。

上述创建的这个类，是Crank-Nicolson法，针对的是两个边界均为第一类条件，且边界值为一个常数，如果边界条件为一个函数，那么需要稍微改造下即可。这里的计算总时间取0.3s。$x$方向划分的段数为100，总时间步取1000。

demo = Parabolic_CnMethod(1, [0, 1], 0.3, 100, 1000)
u = demo.solve(np.zeros(demo.t.shape[0]), np.zeros(demo.t.shape[0]), np.ones(demo.x.shape[0] - 2) * 100)

import matplotlib.pyplot as plt
cb = plt.pcolormesh(res.x, res.t, ww, cmap='Spectral_r', shading='gouraud')
plt.contour(res.x, res.t, ww, colors='gray', levels = 10)
plt.colorbar(cb)
plt.show()

2020！！新年快乐！

本站文章如非特别说明，均为原创。未经本人同意，请勿转载！转载请务必注明出处！