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

最新推荐文章于 2023-03-10 09:46:08 发布

浪子膏丷

最新推荐文章于 2023-03-10 09:46:08 发布

阅读量8k

点赞数 18

分类专栏：微分方程数值解文章标签： python 有限差分法

本文链接：https://blog.csdn.net/qq_39482771/article/details/121587891

版权

微分方程数值解专栏收录该内容

5 篇文章 36 订阅

订阅专栏

前言

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

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

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

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

方程与问题

1、定义问题
在这里插入图片描述
定义:

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

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

原方程：

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

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

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

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

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

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

$(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$

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

$\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}$

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

写出矩阵的形式：
在这里插入图片描述

代码

1、定义方程并求解

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

import numpy as np

def grid(m,n):
    #将空间离散步数为m，时间离散步数为n，还需要定义L，R，t
    h = (R-L)/m
    tau = t/n
    X = 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)*t

def f(x,t):
    return np.sin(2*np.pi*x)+4*np.sin(2*np.pi*x)*t*np.pi**2

def right_solution(t):
    u = np.zeros(t.shape)
    return u

def left_solution(t):
    u = np.zeros(t.shape)
    return u


m = 100
n = 100
L = 0
R = 1
t = 1
h,tau,X,T=grid(m,n)
r = tau/h/h

U = 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*F
    b = A0@U[i,1:-1].T + RHS
    U[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不等式导致理论过不去，应该会用到三维问题的差分求解。

二维问题算例

$\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.$
求解区域： $\Omega = [0,1]\times[0,1]$

真解： $\exp (-t)*sin(\pi x)sin(\pi y)$

import numpy as np
from numpy import sin,cos,pi,exp

def grid(m,n):
    """
    #此处默认X轴和Y轴剖分步长一样
    """
    h = (R-L)/m  
    tau = t/n
    X1 = 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 u

m = 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_NEW
    U_NEW = np.linalg.solve(A1,b)  
    U[:,:,i+1] = U_NEW.reshape(11,-1)   #再用reshape改变回来
    
    #强制改变边界条件
    U[:,-1,i+1] = 0
    U[:,0,i+1] = 0
    U[-1,:,i+1] = 0
    U[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并行格式的稳定性

浪子膏丷

关注

18
点赞
踩
113

收藏

觉得还不错? 一键收藏
0
评论
【有限差分法】（三）一维和二维抛物方程CN格式以及长时间稳定性分析（附算例与Python代码）

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

扫一扫

专栏目录