前言

我们知道差分里的CN格式是无条件稳定的。但是最近在学习有限元结合CN格式算长时间抛物问题的时候（在时间方向用差分空间方向用有限元）发现稳定性却不能保证，其数值解和真解误差会随着时间越来越大。并且实际能算的时间只有几秒，这几秒钟几乎是没有实际意义的。所以就想到一个问题，纯CN差分格式在长时间计算的时候，是否也是个理论看还行，实践臭弟弟的“花瓶”呢？
在以前做差分课本上的算例大多也是只有几秒，还真没注意到长时间稳定性的问题(¬_¬ )。所以本文找了最简单的一维抛物问题，把时间层增加至1000S，来看看CN格式实际是否真的如理论那样，能保证长时间的稳定性。

本文给出：
1、一个可以搬来用的CN格式一维、二维抛物问题python代码模板
2、长时间（1000s）稳定性结果展示。

注：稳定性是指前面所带来的小误差经过迭代、运算后对后面所得到结果误差的影响是可控的，相对的，如果不可控，则前面的误差会随着计算越来越大，这也叫数值溢出或者blow up。在不同的情况下有几种不同的稳定性定义，他们差别白话版理解：
1、长时间稳定：对给定的网格比，不管时间层n怎么取（哪怕时间层非常大），第n层的误差都不会出现数值溢出。
2、无条件稳定：对给定的时间层n，不管网格比怎么取，第n层的误差都都不会出现数值溢出。
3、绝对稳定：差分格式矩阵的增长因子可以被常数控制，并且这个常数与网格比无关。

可以看出，长时间稳定和无条件稳定之间不存在包含关系。

方程与问题

1、定义问题

定义:

f=(1+4tπ2)sin(2πx)f = (1+4t\pi ^2)sin(2\pi x)f=(1+4tπ2)sin(2πx)
初值: u(x,0)=0u(x,0)=0u(x,0)=0
真解：u=sin(2πx)tu =sin(2\pi x)tu=sin(2πx)t

2、CN格式离散
CN格式在抛物方程里具有蛮高的地位，本身的形式并不复杂，只用了两层，难得的一点是无条件稳定。所以在有限元结合差分问题分析时，一般都是选用CN格式。

原方程：

∂u∂t=∂2u∂x2+f\frac{\partial u}{\partial t} = \frac{\partial^2 u}{\partial x^2} +f∂t∂u=∂x2∂2u+f

现在对方程做CN差分格式离散，这个格式最关键的部分是定在时间层的t = n+1/2处。在时间方向，用一次中心差分格式，也就是
∂u∂t∣t=n+1/2=(ujn+1−ujn)/τ\frac{\partial u}{\partial t} |_{t = n+1/2}= (u_{j}^{n+1}-u^n_j)/\tau∂t∂u∣t=n+1/2=(ujn+1−ujn)/τ

在空间方向，我们实际在做剖分的时候，是没有n+1/2这个时间点的，所以就可以把实际存在但是剖分中不存在的这层的u用t = n和t = n+1层来表示，用大写的UUU来记作这一整层的x，即

U∣t=n+1/2=(U∣t=n+1+U∣t=n)/2U|_{t =n+1/2} = (U|_{t =n+1}+U|_{t =n})/2U∣t=n+1/2=(U∣t=n+1+U∣t=n)/2

然后再用二阶中心差分格式（见系列一）代替uxxu_{xx}uxx就好了。

经过上面所述离散后的方程为：

(ujk+1−ujk)/τ=12h2[uj+1k+1−2ujk+1+uj−1k+1+uj+1k−2ujk+uj−1k]+fj(u_j^{k+1}-u_j^k)/\tau =\frac 1 {2h^2}[u_{j+1}^{k+1}-2u_{j}^{k+1}+u_{j-1}^{k+1}+u_{j+1}^{k}-2u_{j}^{k}+u_{j-1}^{k}] +f_j(ujk+1−ujk)/τ=2h21[uj+1k+1−2ujk+1+uj−1k+1+uj+1k−2ujk+uj−1k]+fj

现在把τ/h2\tau/h^2τ/h2记作rrr，经过简单的变换后最后的结果：

−r2uj+1k+1+(1+r)2ujk+1−r22uj−1k+1=r2uj+1k−(1−r)2ujk+r22uj−1k- \frac{r}{2}u_{j+1}^{k+1}+(1+r)2u_{j}^{k+1}-\frac{r}{2}2u_{j-1}^{k+1}= \frac{r}{2}u_{j+1}^{k}-(1-r)2u_{j}^{k}+\frac{r}{2}2u_{j-1}^{k}−2ruj+1k+1+(1+r)2ujk+1−2r2uj−1k+1=2ruj+1k−(1−r)2ujk+2r2uj−1k

这样可以直观的看出上一层的已知项如何求解下一层的未知项。

写出矩阵的形式：

代码

1、定义方程并求解

先算一个t =1时的来验证CN方法大致的准确性，在第三步再加长时间。如果一开始就把时间加的过长，就无法从图像上比较好的验证方法正确性。

import numpy as npdef grid(m,n):#将空间离散步数为m，时间离散步数为n，还需要定义L，R，th = (R-L)/mtau = t/nX = np.linspace(L,R,m+1)T = np.linspace(0,t,n+1)return h,tau,X,T
def init_solution(x):u0 = np.zeros(x.shape)return u0
def ture_solution(x,t):return np.sin(2*np.pi*x)*tdef f(x,t):return np.sin(2*np.pi*x)+4*np.sin(2*np.pi*x)*t*np.pi**2def right_solution(t):u = np.zeros(t.shape)return udef left_solution(t):u = np.zeros(t.shape)return um = 100
n = 100
L = 0
R = 1
t = 1
h,tau,X,T=grid(m,n)
r = tau/h/hU = np.zeros((len(T),len(X)))
#U[0,:] = init_solution(X)  #此处因为初值为0 ，所以不需要额外定义，如果初值为其他非0的时候就需要定义
U[:,0] = left_solution(T)
U[:,-1] = right_solution(T)#crank_nicholson
d1 = 1 + np.ones((len(X) -2 ,))*r
d2 = 1 - np.ones((len(X) -2 ,)) *r
c = 0.5* np.ones((len(X) -3 ,)) *r
A1 = np.diag (-c , -1) + np.diag ( -c ,1) + np.diag ( d1 )
A0 = np.diag (c , -1) + np.diag (c ,1) +  np.diag( d2 )
for i in range(len(T)-1):F=np.zeros((len(X)-2,))for j in range(len(X)-2):F[j] = f((j+1)*h,i*tau)RHS = tau*Fb = A0@U[i,1:-1].T + RHSU[i+1,1:-1] =np.linalg.solve(A1,b)
#print(U)

2、绘制真解和数值解图像

#把真解表示出来
X,Y = np.meshgrid(X,T)
ture_U = ture_solution(X,Y)
#print(ture_U )## 表面图
import matplotlib.pyplot as plt
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D
plt.rcParams['font.sans-serif'] = ['SimHei']#可以plt绘图过程中中文无法显示的问题fig =plt.figure(figsize=(15,7))ax = fig.add_subplot(1, 2, 1, projection='3d')
ax.plot_surface( X,Y, U, cmap=cm.coolwarm, linewidth=0, antialiased=False)
ax.set_zlim(-1.01, 1)
plt.title("数值解")ax = fig.add_subplot(1, 2, 2, projection='3d')
ax.plot_surface( X,Y,ture_U, cmap=cm.coolwarm, linewidth=0, antialiased=False)
ax.set_zlim(-1, 1)
plt.title("真解")
plt.show()

由图像可以看到，两者基本接近。也可以通过最大模误差验证两者的“接近度”（此处可以用下面代码验证）

3、现在完成本文的第二项任务，分析长时间稳定性

为了达到分析的结果，只要将前面的t=1用t=1000替换，并且将时间层的网格加密，其他地方不变。这里只给出最后分析的代码

err =[]
for i in range(len(T)):err.append(max(abs(U[i,:]-ture_U[i,:])))
import matplotlib.pyplot as plt
t = np.arange(len(err))
plt.scatter(t,err)
# 横坐标t为离散的时间层，在此算例里计算时间为1000秒，离散层数为10000层
# 纵坐标err为对应每一个时间层的最大模误差

并没有出现blow up。即使算到了1000秒，最大模误差也都控制在0.3以内，这是一个不错的结果。

一点想法

差分方法虽然是整个PDE计算方法中最早发展完善、理论最简单的方法，但他现在依然能广泛使用不是没有原因的。对具体的问题应该具体分析，有限元长时间不能解决的问题，回到差分不失为一种好方法。
下一次将会更新三维的抛物问题求解，最近老板接了个锂电池热传导的项目，在用有限元分析时，因为gronwall不等式导致理论过不去，应该会用到三维问题的差分求解。

二维问题算例

{2π2∂u∂t=∂2u∂x2+∂2u∂y2u=0on∂Ω\left\{ \begin{aligned} &2\pi^2 \frac{\partial u}{\partial t} = \frac{\partial^2 u}{\partial x^2} +\frac{\partial^2 u}{\partial y^2} \\ &u = 0 \qquad on \partial \Omega \end{aligned} \right. ⎩⎪⎨⎪⎧2π2∂t∂u=∂x2∂2u+∂y2∂2uu=0on∂Ω
求解区域：Ω=[0,1]×[0,1]\Omega = [0,1]\times[0,1]Ω=[0,1]×[0,1]

真解：u=exp⁡(−t)∗sin(πx)sin(πy)u = \exp (-t)*sin(\pi x)sin(\pi y)u=exp(−t)∗sin(πx)sin(πy)

import numpy as np
from numpy import sin,cos,pi,expdef grid(m,n):"""#此处默认X轴和Y轴剖分步长一样"""h = (R-L)/m  tau = t/nX1 = np.linspace(L,R,m+1)  #X轴X2 = np.linspace(U,D,m+1)  #Y轴T = np.linspace(0,t,n+1)return h,tau,X1,X2,T  #X为空间方向的节点，T为时间方向的节点def ture_u(x,y,t):return exp(-t)*sin(pi*x)*sin(pi*y)def init_u(x,y):return sin(pi*x)*sin(pi*y)def right_solution(t):u = np.zeros(t.shape)return um = 10
n = 10
U = 0
D = 1
L = 0
R = 1
t = 1
h,tau,X1,X2,T=grid(m,n)
r = tau/h/h/(2*pi**2)U = np.zeros((len(X1),len(X2),len(T)))#crank_nicholson
node = (len(X1))*(len(X1))   #自由度个数
d1 = np.ones((node,))*(-2*r+1)
d2 = np.ones((node,))*(2*r+1)d3 = np.ones((node - 1,))*(r/2)
d4 = np.ones((node - len(X1),))*(r/2)A1 = np.diag(d2) - np.diag(d3,-1) - np.diag(d3,1) - np.diag(d4,- len(X1)) - np.diag(d4, len(X1))
A0 =  np.diag(d1) + np.diag(d3,-1) + np.diag(d3,1) + np.diag(d4,- len(X1)) + np.diag(d4, len(X1))#定义初值
X1,X2 = np.meshgrid(X1,X2)
U[:,:,0] = init_u(X1,X2)for i in range(len(T)-1):U_NEW = U[:,:,i].reshape(-1)   #展开为一个一维的向量b = A0@U_NEWU_NEW = np.linalg.solve(A1,b)  U[:,:,i+1] = U_NEW.reshape(11,-1)   #再用reshape改变回来#强制改变边界条件U[:,-1,i+1] = 0U[:,0,i+1] = 0U[-1,:,i+1] = 0U[0,:,i+1] = 0
#print(U[:,:,-1])
TURE_U = ture_u(X1,X2,1)   #t =1 时刻真解
#print(TURE_U)

计算出U的数值解后，可以把两者的图像画出来对比一下

## 表面图
import matplotlib.pyplot as plt
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D
plt.rcParams['font.sans-serif'] = ['SimHei']#可以plt绘图过程中中文无法显示的问题
plt.rcParams['axes.unicode_minus'] = False   #解决负号为方块的问题fig =plt.figure(figsize=(15,7))ax = fig.add_subplot(1, 2, 1, projection='3d')
ax.plot_surface( X1,X2, U[:,:,-1], cmap=cm.coolwarm, linewidth=0, antialiased=False)
ax.set_zlim(-1, 1)
plt.title("数值解")ax = fig.add_subplot(1, 2, 2, projection='3d')
ax.plot_surface( X1,X2,TURE_U, cmap=cm.coolwarm, linewidth=0, antialiased=False)
ax.set_zlim(-1, 1)
plt.title("真解")
plt.show()

参考文献：
【1】：Du Fort—Frankel格式及DFF—I并行格式的稳定性

【有限差分法】（三）一维和二维抛物方程CN格式以及长时间稳定性分析（附算例与Python代码）相关推荐

二维声波方程的有限差分法数值模拟
二维声波方程的有限差分法数值模拟文章目录二维声波方程的有限差分法数值模拟一.实现效果二.matlab代码分享三.python代码分享一.实现效果二.matlab代码分享 close al ...
二维有限元方程matlab,有限元法求解二维Poisson方程的MATLAB实现
有限元法求解二维 Poisson 方程的 MATLAB 实现陈莲a ,郭元辉b ,邹叶童a ( 西华师范大学 a．数学与信息学院; b．教育信息技术中心,四川南充 6437009) 摘要: ...
RBF的一维和二维逼近
RBF在函数逼近上有非常好的表现,这里给出了RBF在单自变量单因变量,以及二自变量单因变量的逼近方法,进行Python实现,多维问题的逼近同样很方便,只需增加相应的维度即可,只是需要更多的数据. 一. ...
二维Poisson方程五点差分格式及简单求解方法Python实现
二维Poisson方程简介给出二维 PoissonPoissonPoisson 方程 DirichletDirichletDirichlet 边值问题: {−Δu=f(x,y)(x,y)∈Ωu=φ ...
python二维数组去重_np.unique()对一维和二维数组去重
一维数组对一维数组或列表,unique()函数去除其中重复元素,并按元素大小返回一个新的无重复元组或列表. import numpy as np A = [1, 2, 2, 5,3, 4, 3] a ...
galerkin有限元法matlab实现,有限元法求解二维Poisson方程的MATLAB实现
有限元法求解二维Poisson方程的MATLAB实现陈莲a,郭元辉b,邹叶童a [摘要]文章讨论了圆形区域上的三角形单元剖分.有限元空间,通过变分形式离散得到有限元方程. 用MATLAB编程求得数值 ...
mysql第三章关系模型_一个MySQL关系模型只有三个关系(二维表)组成。_学小易找答案...
[判断题]DELETE语句功能是对表中所有记录或满足条件的记录进行批量删除. [填空题]The computer's entire ____ was on a single board. [单选题]下 ...
二维burgers方程_二维Burgers方程的RKDG有限元解法
二维 Burgers 方程的 RKDG 有限元解法 ∗ 马艳春 1, 张寅虎 2, 冯新龙 1 [摘要] 摘要 : 本文应用 RKDG 有限元方法求解具有周期边界条件的二维非粘性 Burgers ...
高斯低通滤波 matlab_一维和二维高斯函数及其一阶和二阶导数
二维高斯函数高斯函数在图像滤波.边缘检测等中发挥着重要的作用.高斯滤波是典型的低通滤波,对图像有平滑作用.同时,高斯函数的一阶.二阶导数也可以用于高通滤波,比如canny算子中用到的是高斯函数的一阶 ...
二维burgers方程_用格子Boltzmann方法研究二维Burgers方程
用格子 Boltzmann 方法研究二维 Burgers 方程张伟 ; 李文杰 [期刊名称] <天津城市建设学院学报> [年 ( 卷 ), 期] 2012(018)001 [ 摘要 ] ...

【有限差分法】（三）一维和二维抛物方程CN格式以及长时间稳定性分析（附算例与Python代码）

前言