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

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

CN方法离散成差分方程组，注意对于边界条件旁边的未知数边界项需要移动到等式右边

这里外边界直接采用精确解作为边界条件

import matplotlib
import math
matplotlib.use('TkAgg')
import numpy as np  
import matplotlib.pyplot as plt


def 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 xx


def CN(x,t,a0,u):#Cranck-Nicolson格式
    n = (x.size-2)#  x方向未知数数量，x方向有两个边界作为已知条件，因此未知数是格点数-2
    m = t.size-1   #t方向未知数量，t方向未知数-1

    courant = 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<2
    print('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 u

def 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)