一维有限差分波动方程建模

Question

我试图用有限差分近似法模拟由高斯点源产生的一维波。下面是我的代码。在

import matplotlib.pyplot as plt
import numpy as np

########Pre-Defining Values########
# spacial extent
lox = -1000
upx = 1000

# space sampling interval (km)
dx = 2.0
dx2inv = 1/(dx*dx)

# temporal extent
lot = 0
upt = 60

# time sampling interval (s)
dt = 0.5
dt2 = dt*dt

x = np.arange(lox,upx,dx)
t = np.arange(lot,upt,dt)

# pressure source location
psx = 0

# velocity (km/s)
v = 2.0
v2 = v*v

# density change location
pcl = 500

# density
p1 = 1
p1inv = 1/p1
p2 = 0.2
p2inv = 1/p2
pinv = np.zeros_like(x)
p = np.zeros_like(x)
for i in range(0,(int)((upx+pcl)/dx),1):
  pinv[i] = p1inv
  p[i] = p1
for i in range((int)((upx+pcl)/dx),len(pinv),1):
  pinv[i] = p2inv
  p[i] = p2


# waveform
f = np.zeros((len(t),len(x)))

# source
amp = 20 
mu = 0
sig = 10/dx
s = np.zeros_like(f)
s[0] = 1/(sig*np.sqrt(2*np.pi)) * np.exp(-(x-mu)*(x-mu)/2/sig/sig)
maxinv = 1/np.amax(s[0])
for i in range(1,len(s[0])):
  s[0][i] *= amp*maxinv


########Calculating Waveform########
h = np.zeros_like(f)
n1 = len(f)
n2 = len(f[0])

def fdx(i1):
  for i2 in range(1,n2-1):
    gi  = f[i1][i2  ]
    gi -= f[i1][i2-1]
    gi *= pinv[i2]
    h[i1][i2-1] -= gi
    h[i1][i2  ]  = gi

#f[0] = s[0]
fdx(0)
for i2 in range(0,n2):
  f[1][i2] = 2*f[0][i2] + (s[0][i2] - h[0][i2] * dx2inv) * p[i2] * v2 * dt2
for i1 in range(1,n1-1):
  fdx(i1)
  for i2 in range(0,n2):
    f[i1+1][i2] = 2*f[i1][i2] - f[i1-1][i2] + (s[i1][i2] - h[i1][i2] * dx2inv) * p[i2] * v2 * dt2

########Plotting########
plt.plot(x,f[50])

maxf = 1.5*amp
minf = -1.5*amp
plt.axis([lox,upx,minf,maxf])
plt.xlabel('x')
plt.ylabel('f(x,t)')

# vertical colored bars representing density
plt.axvspan(lox, pcl, facecolor='g', alpha=0.1)
plt.axvspan(pcl, upx, facecolor='g', alpha=0.2)

# text with density values
plt.text(pcl-0.2*upx,0.8*maxf,r'$\rho = $%s'%(p1),fontsize=15)
plt.text(pcl+0.05*upx,0.8*maxf,r'$\rho = $%s'%(p2),fontsize=15)

plt.show()

不幸的是，这段代码没有产生正确的结果（两个高斯脉冲从x=0左右移动）。它产生一个随时间增长的高斯脉冲。有人知道我犯了什么错误吗？在

非常感谢。在

Answer 1

自从你发布这篇文章以来已经有一段时间了，但是如果有什么帮助的话，这里有一个生成高斯脉冲的代码。我不擅长编程，所以我很抱歉，如果这个代码是混淆。我用了电磁波的一维时域有限差分（FDTD）波传播方程（无单位）：

import numpy as np                      
import matplotlib.pyplot as plt
import matplotlib.animation as animation

#defining dimensions
xdim=720
time_tot = 500
xsource = xdim/2

#stability factor
S=1

#Speed of light
c=1
epsilon0=1
mu0=1

delta =1
deltat = S*delta/c

Ez = np.zeros(xdim)
Hy = np.zeros(xdim)

epsilon = epsilon0*np.ones(xdim)
mu = mu0*np.ones(xdim)

fig , axis = plt.subplots(1,1)
axis.set_xlim(len(Ez))
axis.set_ylim(-3,3)
axis.set_title("E Field")
line, = axis.plot([],[])

def init():
    line.set_data([],[])
    return line,

def animate(n, *args, **kwargs):
    Hy[0:xdim-1] = Hy[0:xdim-1]+(delta/(delta*mu[0:xdim-1]))*(Ez[1:xdim]-Ez[0:xdim-1])
    Ez[1:xdim]= Ez[1:xdim]+(delta/(delta*epsilon[1:xdim]))*(Hy[1:xdim]-Hy[0:xdim-1])
    Ez[xsource] = Ez[xsource] + 30.0*(1/np.sqrt(2*np.pi))*np.exp(-(n-80.0)**2/(100))
    ylims = axis.get_ylim()
    if (abs(np.amax(Ez))>ylims[1]):
        axis.set_ylim(-(np.amax(Ez)+2),np.amax(Ez)+2)
    line.set_data(np.arange(len(Ez)),Ez)
    return line,

ani = animation.FuncAnimation(fig, animate, init_func=init, frames=(time_tot), interval=10, blit=False, repeat =False)
fig.show()

我希望有帮助。：）