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

前言

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

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

一：一维热传导方程简介

一维热传导方程如下：
$\frac{\vartheta u}{\vartheta t}=a\frac{{\vartheta }^{2}u}{\vartheta x^2}+f(x),0<t\le 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\le t\le 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+r{u}_{j-1}^n+\tau f_j$ ，取 $n=0$ ，利用初值 $u_j^0={\phi }_j$ 和边值 $u_0^n=u_J^n=0$ 可以通过变形式算出 $u_j^1$ ，同理下去可以算出 $u_j^2....$，此格式不需要求解线性方程组，所以此格式视为显式格式，但显式格式并不是无条件稳定的，只有当 $\alpha \le 1/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$ 。同样$CN$格式是隐式格式，无条件稳定的，由第 n 层计算第 n+1 层时，需要求解线性方程组，但此格式基于时间步长 $\tau$ 能达到二阶收敛。

三：代码实现

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

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

##############导入相应模块
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3Dplt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

编写差分格式类如下：

class diff_schemes:def __init__(self,dt,dx,r,x,t):##定义内置初始化函数self.dt = dtself.dx = dxself.r = rself.x = xself.t = tdef make_figure(self,matx,title_msg):###作图fig = plt.figure()ax = fig.gca(projection='3d')x, y = np.meshgrid(self.x, self.t)z = matxax.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)])###默认边值都为0matx[0,:] = icmatx[:,0] = 0matx[:,-1] = 0for 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 matxdef backward_diff_scheme(self):matx = np.zeros([len(self.t), len(self.x)])matx[0, :] = icmatx[:, 0] = 0matx[:, -1] = 0matxa = np.zeros([len(self.x) - 2,len(self.x) - 2 ])for j in range(len(self.x) - 2):matxa[j, j] = 1 + 2 * self.rif j >= 1:matxa[j - 1, j] = - self.rmatxa[j, j - 1] = - self.riv = 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] = resiv = resprint(matx)return matxdef CN_diff_scheme(self):r1 = 1 + self.rr2 = 1 - self.rmatx = np.zeros([len(self.t), len(self.x)])matx[0, :] = icmatx[:, 0] = 0matx[:, -1] = 0matxp = 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] = r1matxq[j, j] = r2if j >= 1:matxp[j - 1, j] = - self.r / 2matxp[j, j - 1] = - self.r / 2matxq[j - 1, j] = self.r / 2matxq[j, j - 1] = self.r / 2iv = 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] = resiv = 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)

四：数值结果

从上面所有图中结果可知，三种格式最终在末端位置（时间末端和空间末端）的运算结果很接近（如果深入从收敛阶角度出发，其实CN格式的精度更高更接近真解，收敛速率最快）。

当我们网格比取值大于0.5时，显式格式符合理论上的结果，明显不稳定收敛了，但两种隐式格式还是依旧很稳定，从图7很明显能发掘到这一点。

五：总结

本文主要是数值模拟实验类文章，也是学习数值求解偏微分方程的基础性文章。偏微分方程的理论是非常广且深，另外，它跟随机微分方程也有这千丝万缕的联系。如果有兴趣的小伙伴可以多多交流。

写作不易，望路过的朋友们点赞收藏加关注哈，在此表示感谢！