title: 1D Finite Difference Wave Equation Modeling
body: <p>I am attempting to model a 1D wave created by a Gaussian point source using the finite difference approximation method. Below is my code.</p>
  
  <pre><code>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()
  </code></pre>
  
  <p>Unfortunately this code does not produce the correct result (two Gaussian pulses traveling left and right away from x=0). It instead produces one Gaussian pulse that grows with time. Does anyone know what error I am making?</p>
  
  <p>Thank you very much.</p>