2D Wave Equation

• 二维波动方程初边值问题:
  egin{equation}
  egin{cases}
  frac{partial^2u}{partial t^2} = D(frac{partial^2u}{partial x^2} + frac{partial^2u}{partial y^2}), quad (x, y) in Omega = (0, 2)^2, t in (0, 6)   \
  u(x, y, 0) = 0.1mathrm{e}^{-((x - x_c)^2 + (y - y_c)^2) / 0.001}, u_t(x, y, 0) = 0, quad (x, y) in Omega   \
  u(x, y, t) = 0, quad (x, y) in partialOmega, t in [0, 6]
  label{eq_1}
  end{cases}
  end{equation}
• 连续型系统的离散化:
  采用中心差分对波动方程离散化,
  egin{equation}
  frac{u_{i, j}^{k+1} + u_{i, j}^{k-1} - 2u_{i, j}^k}{h_t^2} = Dleft(frac{u_{i+1, j}^k + u_{i-1, j}^k - 2u_{i, j}^k}{h_x^2} + frac{u_{i, j+1}^k + u_{i, j-1}^k - 2u_{i, j}^k}{h_y^2} ight)
  label{eq_2}
  end{equation}
  由于该差分格式在时间上是多步法, 启动时需要给定前两层时间步上的数值才能构造迭代数列. 初值条件仅提供第0层时间步数值, 需要手动计算第1层时间步数值, 以填充迭代数列之启动数组.
  根据初值条件, 采用中心差分离散之,
  egin{equation}
  frac{u_{i, j}^1 - u_{i, j}^{-1}}{2h_t} = 0
  label{eq_3}
  end{equation}
  根据式( ef{eq_2}), 令$k = 0$有,
  egin{equation}
  frac{u_{i, j}^{1} + u_{i, j}^{-1} - 2u_{i, j}^0}{h_t^2} = Dleft(frac{u_{i+1, j}^0 + u_{i-1, j}^0 - 2u_{i, j}^0}{h_x^2} + frac{u_{i, j+1}^0 + u_{i, j-1}^0 - 2u_{i, j}^0}{h_y^2} ight)
  label{eq_4}
  end{equation}
  结合式( ef{eq_3})、式( ef{eq_4})可得第1层时间步数值,
  egin{equation}
  u_{i, j}^1 = frac{Dh_t^2}{2}left( frac{u_{i+1, j}^0 + u_{i-1, j}^0 - 2u_{i, j}^0}{h_x^2} + frac{u_{i, j+1}^0 + u_{i, j-1}^0 - 2u_{i, j}^0}{h_y^2} ight) + u_{i, j}^0
  label{eq_5}
  end{equation}
• 离散型方程的数值求解(代码实现):
  

  # 2D Wave Equation - Finite Difference之实现
  
  import os
  import shutil
  import numpy
  from matplotlib import pyplot as plt
  from PIL import Image
  
  
  # 波动方程求解
  class WaveEq(object):
      
      def __init__(self, nx, ny, nt, D=0.1, xc=1, yc=1):
          self.__nx = nx                 # x轴网格数
          self.__ny = ny                 # y轴网格数
          self.__nt = nt                 # t轴网格数
          self.__D = D
          self.__xc = xc
          self.__yc = yc
          
          self.__init_grid()             # 网格初始化
          
      
      def __init_grid(self):
          xMin, xMax = 0, 2
          yMin, yMax = 0, 2
          tMin, tMax = 0, 6
          self.__hx = (xMax - xMin) / self.__nx
          self.__hy = (yMax - yMin) / self.__ny
          self.__ht = (tMax - tMin) / self.__nt
          self.__X = numpy.linspace(xMin, xMax, self.__nx + 1)
          self.__Y = numpy.linspace(yMin, yMax, self.__ny + 1)
          self.__T = numpy.linspace(tMin, tMax, self.__nt + 1)
          self.__X, self.__Y = numpy.meshgrid(self.__X, self.__Y)
  
  
      def get_solu(self):
          U0 = self.__calc_U0()
          yield self.__X, self.__Y, U0
          U1 = self.__calc_U1(U0)
          yield self.__X, self.__Y, U1
          Uk_1, Uk_2 = U1, U0
          for i in range(2, self.__nt+1):
              Uk = self.__calc_Uk(Uk_1, Uk_2)
              yield self.__X, self.__Y, Uk
              Uk_1, Uk_2 = Uk, Uk_1
          
          
      def __calc_Uk(self, Uk_1, Uk_2):
          '''
          计算第k层时间步数值
          '''
          Uk = numpy.zeros(Uk_1.shape)
          for i in range(1, self.__nx):
              for j in range(1, self.__ny):
                  term1 = (Uk_1[j, i+1] + Uk_1[j, i-1] - 2 * Uk_1[j, i]) / self.__hx ** 2
                  term2 = (Uk_1[j+1, i] + Uk_1[j-1, i] - 2 * Uk_1[j, i]) / self.__hy ** 2
                  Uk[j, i] = (term1 + term2) * self.__D * self.__ht ** 2 + 2 * Uk_1[j, i] - Uk_2[j, i]
          return Uk        
      
      
      def __calc_U1(self, U0):
          '''
          计算第1层时间步数值
          '''
          U1 = numpy.zeros(U0.shape)
          for i in range(1, self.__nx):
              for j in range(1, self.__ny):
                  term1 = (U0[j, i+1] + U0[j, i-1] - 2 * U0[j, i]) / self.__hx ** 2
                  term2 = (U0[j+1, i] + U0[j-1, i] - 2 * U0[j, i]) / self.__hy ** 2
                  U1[j, i] = (term1 + term2) * self.__D * self.__ht ** 2 / 2 + U0[j, i]
          return U1
      
          
      def __calc_U0(self):
          '''
          计算第0层时间步数值
          '''
          U0 = self.__calc_u0(self.__X, self.__Y)
          self.__fill_boundary(U0)                                 # 填充边界条件
          return U0
          
          
      def __fill_boundary(self, mat):
          mat[0, :] = 0
          mat[-1, :] = 0
          mat[:, 0] = 0
          mat[:, -1] = 0
                  
          
      def __calc_u0(self, x, y):
          '''
          计算初始profile
          '''
          u0 = 0.1 * numpy.exp(-((x - self.__xc) ** 2 + (y - self.__yc) ** 2) / 0.001)
          return u0
          
          
      def __calc_v0(self, x, y):
          '''
          计算初始velocity
          '''
          v0 = 0 * x
          return v0
          
          
  # 动态图绘制
  class WavePlot(object):
      
      def __init__(self, waveObj):
          self.__waveObj = waveObj
          
      
      def ani_plot(self, aniPath):
          if not os.path.exists(aniPath):
              os.mkdir(aniPath)
          
          imgs = list()
          for idx, solu in enumerate(self.__waveObj.get_solu()):
              print(idx)
              if (idx % 5 == 0):
                  X, Y, Uk = solu
                  fig = plt.figure(figsize=(8, 8))
                  ax1 = plt.subplot()
                  ax1.pcolor(X, Y, Uk[:-1, :-1], cmap="jet", vmin=0)
                  filename = "{}/{}.png".format(aniPath, idx)
                  fig.savefig("{}".format(filename), dpi=80)
                  plt.close()
                  img = Image.open(filename)
                  imgs.append(img)
          img.save("wave_plot.gif", save_all=True, append_images=imgs, duration=5)
          shutil.rmtree(aniPath)
  
  
  
  if __name__ == "__main__":
      waveObj = WaveEq(300, 300, 1000, 0.1)
      wpltObj = WavePlot(waveObj)
      wpltObj.ani_plot("./wave_plot")

• 结果展示:
  
• 参考文档:
  https://zhuanlan.zhihu.com/p/111640958