解双曲型非线性方程的Harden-Yee算法
先贴代码，教程后面有空再写

import matplotlib
import math
matplotlib.use('TkAgg')
import numpy as np  
import matplotlib.pyplot as plt
def Phiy(yy,epsi):#phi(y)
    if abs(yy) >= epsi:
        phiyy = abs(yy)
    else:
        phiyy = (yy*yy + epsi*epsi)/2.*epsi
    return phiyy
def Deltau(u2, u1):#\delta u(n,i+1/2),u2表示右边的
    return u2 - u1

def Gi_0(u2,u1,u0):
    result = Minmod((u1-u0),(u2-u1))
    return result

def Minmod(aa,bb):
    if aa*bb<=0.:
        result = 0.
    else:
        result = np.sign(aa)*min(abs(aa),abs(bb))
    return result

def Hartenyee(u,x,t):
    corrant = (t[2] - t[1])/( x[2]-x[1])#corrant 数
    epsi=0.3
    for j in range(0, t.size-1):# j 表示t方向，t从零开始
        u[j+1,1]=Lax_Wendroff(u[j,0],u[j,1],u[j,2],corrant)#赋值x1
        print(j)
        for i in range(2, x.size-2): # i 表示x 方向

            if Deltau(u[j,i+1],u[j,i]) == 0.:
                alphai_1plus2 = u[j,i]
                beta_1plus2 = 0.
            else:
                alphai_1plus2 = (u[j,i+1] + u[j,i])/2.#E2-E1
                beta_1plus2 = (Gi_0(u[j,i+2],u[j,i+1],u[j,i]) - Gi_0(u[j,i+1],u[j,i],u[j,i-1]))/(u[j,i+1]-u[j,i])


            if Deltau(u[j,i],u[j,i-1]) == 0.:
                alphai_1minus2 = u[j,i]
                beta_1minus2 = 0.
            else:
                alphai_1minus2 = (u[j,i] + u[j,i-1])/2.#E2-E1
                beta_1minus2 = (Gi_0(u[j,i+1],u[j,i],u[j,i-1]) - Gi_0(u[j,i],u[j,i-1],u[j,i-2]))/(u[j,i]-u[j,i-1])
            
            phini_plus12 = Gi_0(u[j,i+2],u[j,i+1],u[j,i])+Gi_0(u[j,i+1],u[j,i],u[j,i-1])-Phiy(alphai_1plus2+beta_1plus2,epsi)*(u[j,i+1]-u[j,i])
            phini_minus12 = Gi_0(u[j,i+1],u[j,i],u[j,i-1])+Gi_0(u[j,i],u[j,i-1],u[j,i-2])-Phiy(alphai_1minus2+beta_1minus2,epsi)*(u[j,i]-u[j,i-1])
            
            u[j+1,i] = u[j,i] - corrant/2.*(u[j,i+1]*u[j,i+1]/2.-u[j,i-1]*u[j,i-1]/2.)-corrant/2*(phini_plus12-phini_minus12)
                 
    return u
def Lax_Wendroff(u0,u1,u2,corrant):   #lax——Wendroff格式
    dE1 = u2*u2/4.-u0*u0/4.
    dE2 = u2*u2/4.-u1*u1/4.
    dE3 = u1*u1/4.-u0*u0/4.
    result = u1 - corrant/2.*dE1 + corrant*corrant/2.*((u1+u2)/2.*dE2-(u1+u0)/2.*dE3) 
    return result

def Plot(x, t, result,title):
    plt.figure()
    plt.plot(x, result[0,:])
    y = [t[int(t.size/5)], t[int(t.size/4)], t[int(t.size/3)], t[int(t.size/2)],t[int(t.size/1.1)]]

    labels = ['t='f'{num:.2f}' for num in y]#将变量转换为字符串

    plt.plot(x, result[int(t.size/5),:], label=labels[0])
    plt.plot(x, result[int(t.size/4),:],label=labels[1])
    plt.plot(x, result[int(t.size/3),:],label=labels[2])
    plt.plot(x, result[int(t.size/2),:],label=labels[3])
    plt.plot(x, result[int(t.size/1.1),:],label=labels[4])
    plt.legend()
    plt.xlabel('x')
    plt.ylabel('t')
    plt.title(title)
    plt.show()
    plt.close()

    plt.figure()
    plt.contourf(x, t, result, 50, cmap = 'jet')
    plt.colorbar()
    plt.savefig('CN.png')
    plt.xlabel('x')
    plt.ylabel('t')
    plt.title(title)
    plt.show()
    plt.close()
    return 0

x = np.linspace(0,4,400)
t = np.linspace(0,4.0,400)
u = np.zeros((t.size,x.size),dtype=float)#注意，这里u不是必须二维，我用二维主要为了测试和画图方便，当t比较大的时候会占用大量内存
u[0,0:int(x.size/2)] = 1.0
u[:,0] = 1.0
u[:,-1] = 0.0
corrant = (t[2] - t[1])/( x[2]-x[1])#corrant 数
'''#Lax_Wendroff 方法做对比
for j in range (1, t.size-1):
    print(j)
    for i in range (1,x.size-1):
        u[j,i] = Lax_Wendroff(u[j-1,i-1],u[j-1,i],u[j-1,i+1],corrant)
'''



u2=Hartenyee(u,x,t)
Plot(x,t,u2,'Harten-Yee solver')