解双曲型非线性方程的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.*epsireturn phiyy
def Deltau(u2, u1):#\delta u(n,i+1/2),u2表示右边的return u2 - u1def Gi_0(u2,u1,u0):result = Minmod((u1-u0),(u2-u1))return resultdef Minmod(aa,bb):if aa*bb<=0.:result = 0.else:result = np.sign(aa)*min(abs(aa),abs(bb))return resultdef Hartenyee(u,x,t):corrant = (t[2] - t[1])/( x[2]-x[1])#corrant 数epsi=0.3for 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)#赋值x1print(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-E1beta_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-E1beta_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 resultdef 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 0x = 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')