- 前言
- 一：一维热传导方程简介
- 二：差分格式
- 三：代码实现
- 四：数值结果
- 五：总结
热方程的在很多领域都有所应用，熟知的在金融领域求解期权定价公式之Black-Scholes方程，就可以用数值格式求解此类方程，因为很多复杂的期权定价公式很难有显式解，数值方法在这方面就有很多优越性。本文将基于Python编写常见的三种数值格式来求解传统的初-边值一维热方程问题。
岁月如云，匪我思存，写作不易，望路过的朋友们点赞收藏加关注哈，在此表示感谢！
一：一维热传导方程简介一维热传导方程如下： 特别提示：由于本文不是纯粹的数学推理文章，所以会涉及到一些条件和假设的简化，如果从纯粹数学角度考虑，一个条件的成立往往需要很多假设和前提，以及精准的定义，这一点跟工业中算法成立还是有很大区别的。 对于光滑
f
(
x
)
f(x)
f(x) 和
Φ
(
x
)
Phi(x)
Φ(x) ，我们从初-边值考虑差分逼近。取空间步长
h
=
l
/
J
h=l/J
h=l/J 和时间步长
τ
=
T
/
N
tau=T/N
τ=T/N ，其中
J
、
N
J、N
J、N 都是自然数。用两族平行直线
x
=
x
j
=
j
h
(
j
=
0
,
1
,
.
.
.
,
J
)
x=x_j=jh(j=0,1,...,J)
x=xj=jh(j=0,1,...,J) 和
t
=
t
n
=
n
τ
(
n
=
0
,
1
,
.
.
.
N
)
t=t_n=ntau(n=0,1,...N)
t=tn=nτ(n=0,1,...N) 将如下矩形分割成网格，网格节点设为
(
x
j
,
t
n
)
(x_j,t_n)
(xj,tn) ， 我们用
u
j
n
u_j^n
ujn 表示定义在网点
(
x
j
,
t
n
)
(x_j,t_n)
(xj,tn) 上的函数，接着用差商代替上述（1）方程，接下来介绍几种简便的经典差分格式。 即
u
j
n
+
1
−
u
j
n
τ
=
a
u
j
+
1
n
−
2
u
j
n
+
u
j
−
1
n
h
2
+
f
j
frac{u_j^{n+1}-u_j^n}{tau}=afrac{u_{j+1}^{n}-2u_j^n+u_{j-1}^n}{h^2}+f_j
τujn+1−ujn=ah2uj+1n−2ujn+uj−1n+fj ，其中
j
=
1
,
2
,
.
.
.
,
J
−
1
,
n
=
0
,
1
,
.
.
.
,
N
−
1
j=1,2,...,J-1,n=0,1,...,N-1
j=1,2,...,J−1,n=0,1,...,N−1 ，以
α
=
a
τ
/
h
2
alpha=atau/h^2
α=aτ/h2 表示网格比，上式我们变形可得：
u
j
n
+
1
=
α
u
j
+
1
n
+
(
1
−
2
α
)
u
j
n
+
r
u
j
−
1
n
+
τ
f
j
u_j^{n+1}=alpha u_{j+1}^n+(1-2alpha)u_j^n+ru_{j-1}^n+tau f_j
ujn+1=αuj+1n+(1−2α)ujn+ruj−1n+τfj ，取
n
=
0
n=0
n=0 ，利用初值
u
j
0
=
φ
j
u_j^0=varphi_j
uj0=φj 和边值
u
0
n
=
u
J
n
=
0
u_0^n=u_J^n=0
u0n=uJn=0 可以通过变形式算出
u
j
1
u_j^1
uj1 ，同理下去可以算出
u
j
2
.
.
.
.
u_j^2....
uj2....，此格式不需要求解线性方程组，所以此格式视为显式格式，但显式格式并不是无条件稳定的，只有当
α
≤
1
/
2
alphaleq1/2
α≤1/2 时，格式才是稳定的。通过泰勒公式与截断误差定义，我们能证明出此格式基于时间步长
τ
tau
τ 是一阶收敛的（由于此文不是存粹数学文章，这里涉及到的一些结论不再过多深入，有兴趣的朋友可以参考相关文献）。 即
u
j
n
+
1
−
u
j
n
τ
=
a
u
j
+
1
n
+
1
−
2
u
j
n
+
1
+
u
j
−
1
n
+
1
h
2
+
f
j
frac{u_j^{n+1}-u_j^n}{tau}=afrac{u_{j+1}^{n+1}-2u_j^{n+1}+u_{j-1}^{n+1}}{h^2}+f_j
τujn+1−ujn=ah2uj+1n+1−2ujn+1+uj−1n+1+fj ，其中
j
=
1
,
2
,
.
.
.
,
J
−
1
,
n
=
0
,
1
,
.
.
.
,
N
−
1
j=1,2,...,J-1,n=0,1,...,N-1
j=1,2,...,J−1,n=0,1,...,N−1 ，同理我们改写上式为：
−
α
u
j
+
1
n
+
1
+
(
1
+
2
α
)
u
j
n
+
1
−
α
u
j
−
1
n
+
1
=
u
j
n
+
τ
f
j
-alpha u_{j+1}^{n+1}+(1+2alpha)u_j^{n+1}-alpha u_{j-1}^{n+1}=u_j^n+tau f_j
−αuj+1n+1+(1+2α)ujn+1−αuj−1n+1=ujn+τfj ,令 $n=0,1,2,… $，则可利用
u
j
0
u_j^0
uj0 和边值确定
u
j
1
u_j^1
uj1 ，利用
u
j
1
u_j^1
uj1 和边值确定
u
j
2
u_j^2
uj2 等，现在第
n
+
1
n+1
n+1 层的值不能用第
n
n
n 层值明显表示，而是由线性方程组求解。由于变形式左端系数矩阵是严格对角占优的，所以方程肯定有解的，且该格式是无条件稳定的。同样此格式也是基于时间步长
τ
tau
τ 一阶收敛的。 即向前差分格式和向后差分格式做算术平均，即可得到CN格式如下： 接下来我们通过python的面向对象编程，逐一实现上面的三种数值格式。 以下代码经过本人调试，没任何bug，直接复制即可，有兴趣需要的朋友别忘记点赞关注加收藏哦！！谢谢！！ 编写差分格式类如下： 调用相关类中方法展示实验结果 从上面所有图中结果可知，三种格式最终在末端位置（时间末端和空间末端）的运算结果很接近（如果深入从收敛阶角度出发，其实CN格式的精度更高更接近真解，收敛速率最快）。 当我们网格比取值大于0.5时，显式格式符合理论上的结果，明显不稳定收敛了，但两种隐式格式还是依旧很稳定，从图7很明显能发掘到这一点。 本文主要是数值模拟实验类文章，也是学习数值求解偏微分方程的基础性文章。偏微分方程的理论是非常广且深，另外，它跟随机微分方程也有这千丝万缕的联系。如果有兴趣的小伙伴可以多多交流。 写作不易，望路过的朋友们点赞收藏加关注哈，在此表示感谢！
ϑ
u
ϑ
t
=
a
ϑ
2
u
ϑ
x
2
+
f
(
x
)
,
0
<
t
≤
T
.
（
1
）
frac{vartheta u}{vartheta t}=afrac{vartheta^2u}{vartheta x^2}+f(x),0< tleq T.（1）
ϑtϑu=aϑx2ϑ2u+f(x),0
u
j
n
+
1
−
u
j
n
τ
=
a
2
[
u
j
+
1
n
+
1
−
2
u
j
n
+
1
+
u
j
−
1
n
+
1
h
2
+
u
j
+
1
n
−
2
u
j
n
+
u
j
−
1
n
h
2
]
+
f
i
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
τujn+1−ujn=2a[h2uj+1n+1−2ujn+1+uj−1n+1+h2uj+1n−2ujn+uj−1n]+fi，
其中
j
=
1
,
2
,
.
.
.
,
J
−
1
,
n
=
0
,
1
,
.
.
.
,
N
−
1
j=1,2,...,J-1,n=0,1,...,N-1
j=1,2,...,J−1,n=0,1,...,N−1，同理上式我们可以改写为
−
α
2
u
j
+
1
n
+
1
+
(
1
+
α
)
u
j
n
+
1
−
α
2
u
j
−
1
n
+
1
=
r
2
u
j
+
1
n
+
(
1
−
α
)
u
j
n
+
α
2
u
j
−
1
n
+
τ
f
j
-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
−2αuj+1n+1+(1+α)ujn+1−2αuj−1n+1=2ruj+1n+(1−α)ujn+2αuj−1n+τfj ，利用
u
j
0
u_j^0
uj0 和边值便可以逐层求解到
u
j
n
u_j^n
ujn 。同样
C
N
CN
CN格式是隐式格式，无条件稳定的，由第 n 层计算第 n+1 层时，需要求解线性方程组，但此格式基于时间步长
τ
tau
τ 能达到二阶收敛。
##############导入相应模块
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)
四：数值结果
