偏微分方程数值解

抛物型方程的有限差分方法

最为常见的抛物型偏微分方程就是热传导方程
$\left\{\begin{array}{ll} u_{t}=Du_{x x}, & a \leqslant x \leqslant b \quad t \geqslant 0 \\ u(x, 0)=f(x), & a \leqslant x \leqslant b \\ u(a, t)=l(t), & t \geqslant 0 \\ u(b, t)=r(t), & t \geqslant 0 \end{array}\right.$
它本质上其实是一个初边值问题，直观来看是因为，它既需要初值条件，又需要边值条件。
对于偏导数的差分格式，一般分为隐式差分与显式差分两种，下面分别进行介绍。

显式差分格式

对于u的一阶导数，可以进行如下的近似：

$\frac{u_{n+1}-u_{n}}{h} \simeq u^{\prime}$
其中， $h$ 为设定的空间步长。如果设定h为单位长度，则

$u^{\prime}=\nabla u_{n+1}=u_{n+1}-u_{n}$
其中， $\nabla$ 是后向差分的符号，而不是梯度运算。
由此进而可以推导出 $u''$ 的表达式为：

$u^{\prime \prime}=\nabla u_{n+1}-\nabla u_{n}=u_{n+1}-2 u_{n}+u_{n-1}$
然后，再将 $h$ 代回，便可得到最终的差分格式：

$u_{x x}(x, t) \approx \frac{1}{h^{2}}(u(x+h, t)-2 u(x, t)+u(x-h, t))$
该式子中均是对空间的差分，所以 $t$ 并没有变化。
在对于时间的差分中，我们便有两种不同的选择，这也对应了显式和隐式差分格式。
对于显式差分，其时间的差分方式如下：
$u_{t}(x, t) \approx \frac{1}{k}(u(x, t+k)-u(x, t))$
其中， $k$ 为设定的时间步长。可以引入变量 $\sigma = \frac{ck}{h^2}$
我们定义 $w_{i,j}$ 就是在 $(x,t)$ 处的数值解，那么自然的，我们可以将它按照时间的顺序排列，得到下面的一个形式:
$w_{i, j+1}=w_{i j}+\sigma\left(w_{i+1, j}-2 w_{i j}+w_{i-1, j}\right)$

冯诺依曼分析法

它可以应用在各种常见的偏微分方程中。所谓的冯诺依曼分析法，其本质就是用向量，将所有的空间离散点拼起来。那么这样的话，这个矩阵所对应的表达式就仅仅与时间有关了。这样自然分析就会方便很多。
上述的差分形式，可以写成矩阵形式，如下：
$\left[\begin{array}{c} w_{1, j+1} \\ \vdots \\ w_{m, j+1} \end{array}\right]=\left[\begin{array}{ccccc} 1-2 \sigma & \sigma & 0 & \cdots & 0 \\ \sigma & 1-2 \sigma & \sigma & \ddots & \vdots \\ 0 & \sigma & 1-2 \sigma & \ddots & 0 \\ \vdots & \ddots & \ddots & \ddots & \sigma \\ 0 & \cdots & 0 & \sigma & 1-2 \sigma \end{array}\right]\left[\begin{array}{c} w_{1 j} \\ \vdots \\ w_{m j} \end{array}\right]+\sigma\left[\begin{array}{c} w_{0, j} \\ 0 \\ \vdots \\ 0 \\ w_{m+1, j} \end{array}\right]$
写的紧凑一点，就是公式 $W_{j+1} = AW_j + s_j$ ,那么为了让这个公式收敛，我们显然需要让 $A$ 的最大特征值小于1，最后一般需要 $\sigma<2$
实际的求解中，可以直接进行递推，然后对于首尾直接代入边界条件即可。

#一维热传导
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.animation as animation

a = 1
#第三类边界条件
c = 1
#dx太小或dt太大可能出错
dx = 0.01
dt = 0.00005

x = np.arange(0, 1, dx)
t = np.arange(0, 1, dt)

u = np.zeros((len(x), len(t)))
#初始条件
u[:, 0] = 0
#边界条件
m1 = lambda t : 1 + 0.0 * np.sin(t)
m2 = lambda t : 2 - 0.0 * np.sin(10 * t)
#构造对角矩阵1，-2，1
#其中diag为对角矩阵，1表示上移一格
A = -2 * np.eye((len(x))) + np.diag([1] * (len(x) - 1), 1) + np.diag([1] * (len(x) - 1), -1)

for n in range(len(t) - 1):
    u[:, n + 1] = u[:, n] + ((a**2 /(dx**2)) * np.dot(A, u[:, n]) + f.T) * dt
    #第一类边界条件
    #用边界条件确定第一个和最后一个
    u[0, n + 1] = m1(n + 1)
    u[-1, n + 1] = m2(n + 1)
figure = plt.figure()
ax = Axes3D(figure)
T, X = np.meshgrid(t, x)
ax.plot_surface(X, T, u, cmap='rainbow')
plt.xlabel('x')
plt.ylabel('t')
# plt.plot(x, u[:, -1])
# plt.grid()
plt.show()

若得到的解出现震荡，则可以适当增加dx或减小dt，最后得到的求解结果如下图：

隐式差分格式

与显式差分格式相比所不同的是对时间的差分格式，其差分格式如下
$w_{i j}-w_{i, j-1}=\sigma\left(w_{i+1, j}-2 w_{i j}+w_{i-1, j}\right)$
那么也容易通过化简得到我们的矩阵方程:
$\left[\begin{array}{ccccc} 1+2 \sigma & -\sigma & 0 & \cdots & 0 \\ -\sigma & 1+2 \sigma & -\sigma & \ddots & \vdots \\ 0 & -\sigma & 1+2 \sigma & \ddots & 0 \\ \vdots & \vdots & \vdots & \vdots & -\sigma \\ 0 & \cdots & 0 & -\sigma & 1+2 \sigma \end{array}\right]\left[\begin{array}{c} w_{1 j} \\ \vdots \\ w_{m j} \end{array}\right]=\left[\begin{array}{c} w_{1, j-1} \\ \vdots \\ w_{m j}-1 \end{array}\right]+\sigma\left[\begin{array}{c} w_{0 j} \\ 0 \\ \vdots \\ 0 \\ w_{m+1, j} \end{array}\right]$
写成紧凑形式便是， $Aw_j = w_{j-1}+s_j$ ，我们可以通过计算得到A的特征值为 $1+2 \sigma-2 \sigma \cos \frac{\pi j}{m+1}>1$ ，所以与显式差分相比该格式是无条件稳定的。但是并不能直接按时间进行递推求解，需要不断求解矩阵方程，而 $A$ 是以三对角矩阵，较为稀疏，直接对其进行求逆会造成效率低下，于是需要利用追赶法进行求解，可以大幅提升求解效率。

追赶法

对于一矩阵方程：
$AX = d$
其中 $A$ 为一非奇异三对角矩阵，可以对 $A$ 进行LU分解，得到一上三角矩阵和一下三角矩阵。之后可引入中间变量 $y = UX$ ,则有
$\begin{array}{l} A x=L U x=L y=d \\ L y=d \end{array}$
之后已知 $L,d$ 便可求得 $y$ ,已知 $U,y$ 后便可求得 $X$ 。具体代码实现如下：

import numpy as np
from scipy import linalg

#可解决大规模的三对角的矩阵方程
#如在偏微分方程的求解中
'''
b c 0 0
a b c 0
0 a b c
0 0 a b

其中b的绝对值
大于等于a与c的绝对值之和
'''
A = -2 * np.eye((5)) + np.diag([1] * (5 - 1), 1) + np.diag([1] * (5 - 1), -1)
b = np.array([1,1,1,1,1])
# print(b[1])
def Thomas(La, Mb, Uc, b):
    """
        La -- [lower item for tri-diagonal matrix] -- [a, a, a]
        Mb -- [mian item for tri-diagonal matrix] -- [b, b, b, b]
        Uc -- [upper item for tri-diagonal matrix] -- [c, c, c]
        b -- [AX = b, where A is the tri-diagonal matrix] [1,1,1,1,1]
    """
    n = len(Mb)
    Uc[0] = Uc[0] / Mb[0]
    for i in range(1, n-1):
        Uc[i] = Uc[i] / (Mb[i] - La[i - 1] * Uc[i - 1])
    b[0] = b[0] / Mb[0]
    for i in range(1, n):
        b[i] = (b[i] - La[i-1]*b[i-1]) / (Mb[i] - La[i-1] * Uc[i-1])
    ls = list(range(n-1))[::-1]
    for i in ls:
        b[i] = b[i] - Uc[i]*b[i+1]
    return b

之后便可利用追赶法和初始条件，遍历时间节点，逐次求解每个时间点的温度。

Crank-Nicolson方法

因为之前的两个方法它们的误差是 $O(n^2+k^2)$ ,这其实会有问题就是主要的误差其实都在时间步长上。我们需要把时间步长取得很小，可能才能够达到我们要的精度。因此 $Crank-Nicolson$ 方法可以很好地弥补这个缺陷。
这里我们的 $u_t$ 使用向后差分，和上面那一个一样。而我们的 $u_{xx}$ 使用的是混合差分。具体来说，就是用:
$u_{xx}(i,j)=\frac{1}{h^{2}}\left(\frac{1}{2}\left(w_{i+1, j}-2 w_{i j}+w_{i-1, j}\right)+\frac{1}{2}\left(w_{i+1, j-1}-2 w_{i, j-1}+w_{i-1, j-1}\right)\right)$
最后按时间进行整合后得到：
$-\sigma w_{i-1, j}+(2+2 \sigma) w_{i j}-\sigma w_{i+1, j}=\sigma w_{i-1, j-1}+(2-2 \sigma) w_{i, j-1}+\sigma w_{i+1, j-1}$
大体上的形式就是 $A \mathbf{w}_{j}=B \mathbf{w}_{j-1}+\sigma\left(\mathbf{s}_{j-1}+\mathbf{s}_{j}\right)$
其中：
$A=\left(\begin{array}{ccccc} 2+2 \sigma & -\sigma & 0 & \cdots & 0 \\ -\sigma & 2+2 \sigma & -\sigma & \ddots & \vdots \\ 0 & -\sigma & 2+2 \sigma & \ddots & 0 \\ \vdots & \ddots & \ddots & \ddots & -\sigma \\ 0 & \cdots & 0 & -\sigma & 2+2 \sigma \end{array}\right)$
$B=\left(\begin{array}{ccccc} 2-2 \sigma & \sigma & 0 & \cdots & 0 \\ \sigma & 2-2 \sigma & \sigma & \ddots & \vdots \\ 0 & \sigma & 2-2 \sigma & \ddots & 0 \\ \vdots & \ddots & \ddots & \ddots & \sigma \\ 0 & \cdots & 0 & \sigma & 2-2 \sigma \end{array}\right)$
可以计算得到 $A^{-1}B$ 的特征值为：
$\mu_{j}=\frac{2-2 \sigma+\sigma \lambda_{j}}{2+2 \sigma-\sigma \lambda_{j}}, \lambda_{j}=-2 \cos \pi j /(m+1)$
其特征值也是小于1的，故也是无条件稳定。同时其截断误差为 $O(k^2+h^2)$ ,这比上面两种方法的误差要小很多了