基于Python编写求解抛物型pde方程的经典数值格式模拟

前言

热方程的在很多领域都有所应用，熟知的在金融领域求解期权定价公式之Black-Scholes方程，就可以用数值格式求解此类方程，因为很多复杂的期权定价公式很难有显式解，数值方法在这方面就有很多优越性。本文将基于Python编写常见的三种数值格式来求解传统的初-边值一维热方程问题。

岁月如云，匪我思存，写作不易，望路过的朋友们点赞收藏加关注哈，在此表示感谢！

一：一维热传导方程简介

一维热传导方程如下：
$\frac{\vartheta u}{\vartheta t}=a\frac{\vartheta^2u}{\vartheta x^2}+f(x),0< t\leq T.（1）$
其中 a 是正常数， f(x) 是给定的连续函数，对于初-边值问题的描述为，对于有一定阶偏导微商函数 $u (x, t)$ ，满足方程（1）的初始条件： $u(x,0)=\Phi(x),0<x<l$ 以及边界条件 $u(0,t)=u(l,t)=0,0\leq t\leq T$ 。

特别提示：由于本文不是纯粹的数学推理文章，所以会涉及到一些条件和假设的简化，如果从纯粹数学角度考虑，一个条件的成立往往需要很多假设和前提，以及精准的定义，这一点跟工业中算法成立还是有很大区别的。

对于光滑 $f (x)$ 和 $\Phi(x)$ ，我们从初-边值考虑差分逼近。取空间步长 $h = l / J$ 和时间步长 $\tau=T/N$ ，其中 $J 、 N$ 都是自然数。用两族平行直线 $x=x_j=jh(j=0,1,...,J)$ 和 $t=t_n=n\tau(n=0,1,...N)$ 将如下矩形分割成网格，网格节点设为 $x_j,t_n)$ ，
在这里插入图片描述

我们用 $u_j^n$ 表示定义在网点 $x_j,t_n)$ 上的函数，接着用差商代替上述（1）方程，接下来介绍几种简便的经典差分格式。

二：差分格式

向前差分格式（显式格式）

即 $\frac{u_j^{n+1}-u_j^n}{\tau}=a\frac{u_{j+1}^{n}-2u_j^n+u_{j-1}^n}{h^2}+f_j$ ，其中 $j = 1, 2, . . ., J - 1, n = 0, 1, . . ., N - 1$ ，以 $\alpha=a\tau/h^2$ 表示网格比，上式我们变形可得： $u_j^{n+1}=\alpha u_{j+1}^n+(1-2\alpha)u_j^n+ru_{j-1}^n+\tau f_j$ ，取 $n = 0$ ，利用初值 $u_j^0=\varphi_j$ 和边值 $u_0^n=u_J^n=0$ 可以通过变形式算出 $u_j^1$ ，同理下去可以算出 $u_j^2....$ ，此格式不需要求解线性方程组，所以此格式视为显式格式，但显式格式并不是无条件稳定的，只有当 $\alpha\leq1/2$ 时，格式才是稳定的。通过泰勒公式与截断误差定义，我们能证明出此格式基于时间步长 $\tau$ 是一阶收敛的（由于此文不是存粹数学文章，这里涉及到的一些结论不再过多深入，有兴趣的朋友可以参考相关文献）。

向后差分格式（隐式格式）

即 $\frac{u_j^{n+1}-u_j^n}{\tau}=a\frac{u_{j+1}^{n+1}-2u_j^{n+1}+u_{j-1}^{n+1}}{h^2}+f_j$ ，其中 $j = 1, 2, . . ., J - 1, n = 0, 1, . . ., N - 1$ ，同理我们改写上式为： $-\alpha u_{j+1}^{n+1}+(1+2\alpha)u_j^{n+1}-\alpha u_{j-1}^{n+1}=u_j^n+\tau f_j$ ,令 $n=0,1,2,… $，则可利用 $u_j^0$ 和边值确定 $u_j^1$ ，利用 $u_j^1$ 和边值确定 $u_j^2$ 等，现在第 $n + 1$ 层的值不能用第 $n$ 层值明显表示，而是由线性方程组求解。由于变形式左端系数矩阵是严格对角占优的，所以方程肯定有解的，且该格式是无条件稳定的。同样此格式也是基于时间步长 $\tau$ 一阶收敛的。

六点对称格式（Crank-Nicolson 格式）

即向前差分格式和向后差分格式做算术平均，即可得到CN格式如下：
$\frac{u_j^{n+1}-u_j^n}{\tau}=\frac{a}{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]+f_i$ ，
其中 $j = 1, 2, . . ., J - 1, n = 0, 1, . . ., N - 1$ ，同理上式我们可以改写为 $-\frac{\alpha}{2}u_{j+1}^{n+1}+(1+\alpha)u_j^{n+1}-\frac{\alpha}{2}u_{j-1}^{n+1}=\frac{r}{2}u_{j+1}^n+(1-\alpha)u_j^n+\frac{\alpha}{2}u_{j-1}^n+\tau f_j$ ，利用 $u_j^0$ 和边值便可以逐层求解到 $u_j^n$ 。同样 $C N$ 格式是隐式格式，无条件稳定的，由第 n 层计算第 n+1 层时，需要求解线性方程组，但此格式基于时间步长 $\tau$ 能达到二阶收敛。

三：代码实现

接下来我们通过python的面向对象编程，逐一实现上面的三种数值格式。

以下代码经过本人调试，没任何bug，直接复制即可，有兴趣需要的朋友别忘记点赞关注加收藏哦！！谢谢！！

##############导入相应模块
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D

plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

编写差分格式类如下：

class diff_schemes:

    def __init__(self,dt,dx,r,x,t):##定义内置初始化函数
        self.dt = dt
        self.dx = dx
        self.r = r
        self.x = x
        self.t = t

    def make_figure(self,matx,title_msg):###作图
        fig = plt.figure()
        ax = fig.gca(projection='3d')
        x, y = np.meshgrid(self.x, self.t)
        z = matx
        ax.plot_surface(x, y, z, cmap=cm.coolwarm, linewidth=0, antialiased=False)
        ax.set_xlabel('X Label')
        ax.set_ylabel('Y Label')
        ax.set_zlabel('Z Label')
        plt.title(title_msg)
        plt.show()

    def init_condition(self):##定义初值函数，这里选择混个三角函数
        return np.sin(self.x) + np.cos(self.x)  ###也可以选择其他函数如指数函数exp(x)等等

    def forward_diff_scheme(self):##向前差分格式（显式）
        matx = np.zeros([len(self.t),len(self.x)])###默认边值都为0
        matx[0,:] = ic
        matx[:,0] = 0
        matx[:,-1] = 0
        for i in range(1,len(self.t)):
            for j in range(1,len(self.x)-1):
                matx[i,j] = self.r * matx[i - 1,j - 1] + (1 - 2 * self.r) * matx[i - 1,j] + self.r * matx[i - 1,j + 1]
        print(matx)
        return matx

    def backward_diff_scheme(self):
        matx = np.zeros([len(self.t), len(self.x)])
        matx[0, :] = ic
        matx[:, 0] = 0
        matx[:, -1] = 0
        matxa = np.zeros([len(self.x) - 2,len(self.x) - 2 ])
        for j in range(len(self.x) - 2):
            matxa[j, j] = 1 + 2 * self.r
            if j >= 1:
                matxa[j - 1, j] = - self.r
                matxa[j, j - 1] = - self.r
        iv = ic[1:-1]
        matxa_t = np.linalg.inv(matxa)##求解线性方程组
        for i in range(1, len(self.t)):
            res = np.dot(matxa_t,iv)
            matx[i,1:-1] = res
            iv = res
        print(matx)
        return matx

    def CN_diff_scheme(self):
        r1 = 1 + self.r
        r2 = 1 - self.r
        matx = np.zeros([len(self.t), len(self.x)])
        matx[0, :] = ic
        matx[:, 0] = 0
        matx[:, -1] = 0
        matxp = np.zeros([len(self.x) - 2, len(self.x) - 2])
        matxq = np.zeros([len(self.x) - 2, len(self.x) - 2])
        for j in range(len(self.x) - 2):
            matxp[j, j] = r1
            matxq[j, j] = r2
            if j >= 1:
                matxp[j - 1, j] = - self.r / 2
                matxp[j, j - 1] = - self.r / 2
                matxq[j - 1, j] = self.r / 2
                matxq[j, j - 1] = self.r / 2
        iv = np.dot(matxq,ic[1:-1])
        matxa_t = np.linalg.inv(matxp)
        for i in range(1, len(self.t)):
            res = np.dot(matxa_t, iv)
            matx[i, 1:-1] = res
            iv = np.dot(matxq,res)
        print(matx)
        return matx

调用相关类中方法展示实验结果

dt = 0.01##时间步长
dx = 0.01##空间步长
a = 0.5##网格比的取值，且显式格式时a的取值不能大于0.5
X_array = np.arange(-2.5,2.5 + dx,dx)
T_array = np.arange(0,1 + dt,dt)
res = diff_schemes(dt,dx,a,X_array,T_array)
ic = res.init_condition()###初始条件
rfd = res.forward_diff_scheme()
rbd = res.backward_diff_scheme()
rcnd = res.CN_diff_scheme()

res.make_figure(rfd,'网格比a=%s时，显式格式求解'%a)
res.make_figure(rbd,'网格比a=%s时，隐式格式求解'%a)
res.make_figure(rcnd,'网格比a=%s时，CN格式求解'%a)