CN格式求解二维抛物线方程（python代码）

针对之前的CN格式求解二维抛物线方程的MATLAB代码，现在给出python代码，至于算例以及格式可见之前的文章，代码如下：

import numpy as np


def CN_pw_two(hx, hy, tau):
    Lx = 1
    Ly = 1
    T = 1
    x = np.linspace(0, Lx, int(1 / hx + 1))
    y = np.linspace(0, Ly, int(1 / hy + 1))
    t = np.linspace(0, T, int(1 / tau + 1))
    m = len(x)
    n = len(y)
    l = len(t)
    r1 = tau / (hx *hx)
    r2 = tau / (hy *hy)
    # 计算精确解
    ue = np.zeros((m, n, l), dtype=float)
    uc = np.zeros((m, n, l), dtype=float)
    for i in range(m):
        for j in range(n):
            for k in range(l):
                ue[i][j][k] = np.exp((x[i] + y[j])*0.5 - t[k])
    # 定义初始条件和边界条件
    for i in range(m):
        for j in range(n):
            uc[i][j][0] = np.exp((x[i] + y[j])*0.5)
    for j in range(n):
        for k in range(l):
            uc[0][j][k] = np.exp(y[j]*0.5 - t[k])
            uc[m - 1][j][k] = np.exp((1 + y[j]) / 2 - t[k])
    for i in range(m):
        for k in range(l):
            uc[i][0][k] = np.exp(x[i] / 2 - t[k])
            uc[i][n - 1][k] = np.exp((1 + x[i]) / 2 - t[k])
    # 构造系数矩阵
    A = np.zeros((m - 2, m - 2), dtype=float)
    B = np.zeros((n - 2, n - 2), dtype=float)
    for i in range(m - 2):
        A[i][i] = 1 + r1
    for i in range(m - 3):
        A[i][i + 1] = -r1 / 2
        A[i + 1][i] = -r1 / 2
    for i in range(n - 2):
        B[i][i] = 1 + r2
    for i in range(n - 3):
        B[i][i + 1] = -r2 / 2
        B[i + 1][i] = -r2 / 2
    d1 = np.zeros((m - 2, 1), dtype=float)
    v = np.zeros((m , n - 1), dtype=float)
    d2 = np.zeros((n - 2, 1), dtype=float)
    # 开始时间层循环，k
    for k in range(l - 1):
    # 求解v,固定j
        for j in range(n-2):
            for i in range(m-2):
                d1[i][0]=r2/2*(uc[i+1][j][k]+uc[i+1][j+2][k])+(1-r2)*uc[i+1][j+1][k]+tau/2*(-3/2)*(np.exp((x[i+1]+y[j+1])/2-(t[k]+t[k+1])/2))
            v[0][j+1]=(1-r2)/2*uc[0][j+1][k]+(1+r2)/2*uc[0][j+1][k+1]+r2/4*(uc[0][j][k]+uc[0][j+2][k]-uc[0][j][k+1]-uc[0][j+2][k+1])
            v[m-1][j+1]=(1-r2)/2*uc[m-1][j+1][k]+(1+r2)/2*uc[m-1][j+1][k+1]+r2/4*(uc[m-1][j][k]+uc[m-1][j+2][k]-uc[m-1][j][k+1]-uc[m-1][j+2][k+1])
            d1[0][0]=d1[0][0]+r1/2*v[0][j+1]
            d1[m-3][0]=d1[m-3][0]+r1/2*v[m-1][j+1]
            v1=np.linalg.inv(A)@d1
            for i in range(m-2):
                v[i+1][j+1]=v1[i]
        for i in range(m-2):
            for j in range(n-2):
                d2[j][0]=r1/2*(v[i][j+1]+v[i+2][j+1])+(1-r1)*v[i+1][j+1]+tau/2*(-3/2)*(np.exp((x[i+1]+y[j+1])/2-(t[k]+t[k+1])/2))
            d2[0][0]=d2[0][0]+r2/2*uc[i+1][0][k+1]
            d2[n-3][0]=d2[n-3][0]+r2/2*uc[i+1][n-1][k+1]
            v2=np.linalg.inv(B)@d2
            for j in range(n-2):
                uc[i+1][j+1][k+1]=v2[j]
    #计算误差
    err=np.abs(ue-uc)
    return err

#计算时间收敛阶，空间收敛阶亦是如此，读者可自行验证
dt=[1/5,1/10,1/20,1/40,1/80]
m=len(dt)
error=[]
for i in range(m):
    err = np.max(np.max(np.max(CN_pw_two(1 / 200, 1 / 400, dt[i]))))
    error.append(err)
for i in range(m-1):
    rate=np.log2(error[i]/error[i+1])
    print(rate)

总结：上述代码运行结果与MATLAB保持一致，但是所花费时间较长，其次，在编写过程中如遇到问题，若采用pycharm或者IDLE编辑器很难发现其中的错误，因此这里特别推荐使用spyder，它可以直接显示出各个中间变量的值，方便与MATLAB做比较，发现其中哪个环节出了问题。