Cranck-Nicolson隐式方法解线性双曲型方程

Cranck-Nicolson隐式方法解线性双曲型方程
Cranck-Nicolson方法在抛物型方程里面比较常用，双曲型方程例子不多，该方法是二阶精度，无条件稳定，然而，数值震荡比较明显，特别是时间演化比较大以及courant数比较大的时候。
对于这类方程的隐式解法，系数矩阵是三对角矩阵，每次固定时间t，通过解方程的方法解出来下一个时间步上的函数值，比如X方向分了N个格点，去掉边界条件，每个时间步N-2个未知数，构成三对角矩阵。如果采用追赶法求解，courant数需要小于2
网上的例子不多，我这里写了一下，供大家参考
在这里插入图片描述
CN方法离散成差分方程组，注意对于边界条件旁边的未知数边界项需要移动到等式右边

在这里插入图片描述
这里外边界直接采用精确解作为边界条件

python">import matplotlib
import math
matplotlib.use('TkAgg')
import numpy as np  
import matplotlib.pyplot as pltdef catchup(n,a0,b0,c0,F0):             #输入三个对角线上的值a,b,c，F0是等号右边的向量
#注意，追赶法要求：|b1|>|c1|，|bn|>|an|, |bi|>|ai|+|ci|m = np.zeros(F0.size, dtype=float)#中间变量xx = np.zeros(F0.size, dtype=float)for i in range(1, n): #正序，追m[i] = a0[i]/b0[i-1]#          a0是左副对角线，c0是右副对角线b0[i] = np.copy(b0[i]) - m[i]*c0[i-1]#b是主对角线F0[i] = np.copy(F0[i]) - m[i]*F0[i-1]xx[-1] = F0[-1]/b0[-1]for j in range (n-2, -1, -1):#倒序，赶xx[j] = (F0[j]-c0[j]*xx[j+1])/b0[j]return xxdef CN(x,t,a0,u):#Cranck-Nicolson格式n = (x.size-2)#  x方向未知数数量，x方向有两个边界作为已知条件，因此未知数是格点数-2m = t.size-1   #t方向未知数量，t方向未知数-1courant = a0*(t[2]-t[1])/(x[2]-x[1])#courant数，虽然CN算法对courant数没有限制，但追赶法有限制if courant > 2:print('courant number is too large',courant )#由于追赶法要求三条对角线|bi|>|ai|+|ci|，对应从courant<2print('courant=',courant)b = np.zeros(n, dtype=float)#方程组的系数矩阵，b是主对角线，ac是两条副对角线，这是一个三对角矩阵a = np.zeros(n, dtype=float)c = np.zeros(n, dtype=float)f = np.zeros(n, dtype=float)#差分方程等号的右边,和对角线上元素数量一致for j in range (0,m):for i in range(0,n): #对三条对角线分别赋值，对方程组等号右边f也赋值，注意第一行和最后一行需要额外+-courant/4.*u[m+1,0]b[i] = 1.#     主对角线if i >0:   a[i] = -1.*courant/4.if i<(n-1):c[i] =courant/4.f[i] = u[j,i+1]-courant/4.*(u[j,i+2]-u[j,i])'''#这里是系数矩阵，实际上不需要，可以打印出来看一下A = np.zeros((n,n),dtype=float)for i in range(0,n):         A[i,i] = 1.#     主对角线if i >0:A[i-1,i] = -1.*courant/2.if i<(n-1):A[i,i+1] =courant/2.f[i] = u[j,i+1]-courant/4.*(u[j,i+2]-u[j,i])f[0] = np.copy(f[0])+courant/4.*u[j+1,0]f[-1] = np.copy(f[-1])-courant/4.*u[j+1,-1]aa,bb,cc,=get_tridiagonal2(A,f)dd=np.copy(f)'''#u[m+1,-1] = u[m,-1]#+((t[2]-t[1])/(x[2]-x[1]))*(u[m,-1]-u[m,-2])f[0] =  u[j,1]-courant/4.*(u[j,2]-u[j,0])+courant/4.*u[j+1,0]#对方程组等号右边f首元素单独赋值（见公式）f[-1] = np.copy(f[-1])-courant/4.*u[j+1,-1]#对方程组等号右边f末元素单独赋值（见公式）result = catchup(n, a, b, c, f)#采用追赶法，只需输入三条对角线和方程右边u[j+1, 1:(result.size+1)] = np.copy(result)return udef plot(x, t, result):plt.figure()plt.plot(x, result[-1,:])plt.show()plt.close()plt.figure()plt.contourf(x, t, result, 50, cmap = 'jet')plt.colorbar()plt.savefig('CN.png')plt.show()plt.close()return 0#初始化
a0 = 250.
x = np.linspace(0., 400., 200)
t = np.linspace(0., 0.5, 200)u = np.zeros((t.size,x.size), dtype=float)
u[0,:] = 100.*np.cos(math.pi*x/60.)
u[:,0] = 100.*np.cos(math.pi*(0.- a0*t)/60.)#100.*np.cos(-1.*math.pi*a0*t/60.)
u[:,-1] = 100.*np.cos(math.pi*(x[-1] - a0*t)/60.)result = CN(x, t, a0, u)plot(x, t, result)