2D Wave Equation (2) - Finite Difference

• 2维波动方程初边值问题:
  2维波动方程如下,
  \begin{equation}
  \frac{\partial^2u}{\partial t^2} = D\left(\frac{\partial^2u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}\right), \quad (x, y) \in \Omega = [-L_x. L_x] \times[-L_y, L_y], t\in[0, T]
  \label{eq_1}
  \end{equation}
  初始条件如下,
  \begin{equation}
  \left\{
  \begin{split}
  &u(x, y, 0) = u_0(x, y),\quad &(x, y) \in \Omega \\
  &u_t(x, y, 0) = 0,\quad &(x, y) \in \Omega
  \end{split}
  \right.
  \label{eq_2}
  \end{equation}
  边界条件如下,
  \begin{equation}
  u(x, y, t) = 0, \quad (x, y) \in \partial \Omega, t \in [0, T]
  \label{eq_3}
  \end{equation}
• 连续型系统的离散化:
  本文拟降阶时间微分项, 令
  \begin{equation}
  \frac{\partial u}{\partial t} = v
  \label{eq_4}
  \end{equation}
  式$\eqref{eq_1}$转化如下,
  \begin{equation}
  \left\{
  \begin{split}
  \frac{\partial u}{\partial t} &= v \\
  \frac{\partial v}{\partial t} &= D\left( \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} \right) \\
  \end{split}
  \right .
  \label{eq_5}
  \end{equation}
  式$\eqref{eq_2}$转化如下,
  \begin{equation}
  \left\{
  \begin{split}
  u(x, y, 0) &= u_0(x, y) \\
  v(x, y, 0) &= 0 \\
  \end{split}
  \right.
  \label{eq_6}
  \end{equation}
  本文拟采用中心差分处理空间微分项, 式$\eqref{eq_5}$转换如下,
  \begin{equation}
  \left\{
  \begin{split}
  \frac{\partial u_{i, j}}{\partial t} &= v_{i, j} \\
  \frac{\partial v_{i, j}}{\partial t} &= D\left( \frac{u_{i+1,j} + u_{i-1, j} - 2u_{i,j}}{h_x^2} + \frac{u_{i,j+1} + u_{i,j-1} - 2u_{i,j}}{h_y^2} \right)
  \end{split}
  \right .
  \label{eq_7}
  \end{equation}
  令,
  \begin{equation*}
  U = \begin{bmatrix}
  u_{0,0} & \cdots & u_{0,n_x} \\
  \vdots & \ddots & \vdots \\
  u_{n_y,0} & \cdots & u_{n_y,n_x} \\
  \end{bmatrix}\quad
  V = \begin{bmatrix}
  v_{0,0} & \cdots & v_{0,n_x} \\
  \vdots & \ddots & \vdots \\
  v_{n_y,0} & \cdots & v_{n_y,n_x} \\
  \end{bmatrix}\quad
  A_x = \begin{bmatrix}
  -2 & 1 & & & \\
  1 & -2 & 1 & & \\
  & \ddots & \ddots & \ddots & \\
  & & 1 & -2 & 1 \\
  & & & 1 & -2 \\
  \end{bmatrix}\quad
  A_y = \begin{bmatrix}
  -2 & 1 & & & \\
  1 & -2 & 1 & & \\
  & \ddots & \ddots & \ddots & \\
  & & 1 & -2 & 1 \\
  & & & 1 & -2 \\
  \end{bmatrix}
  \end{equation*}
  式$\eqref{eq_7}$转换如下,
  \begin{equation}
  \left\{
  \begin{split}
  \frac{\partial U}{\partial t} &= V \\
  \frac{\partial V}{\partial t} &= D\left( \frac{UA_x}{h_x^2} + \frac{A_yU}{h_y^2} \right)
  \end{split}
  \right.
  \label{eq_8}
  \end{equation}
  注意, 上式$\eqref{eq_8}$中$U$之边界值需以边界条件式$\eqref{eq_3}$填充. 本文拟采用Runge-Kutta方法求解上式$\eqref{eq_8}$, 则有,
  \begin{equation}
  \begin{bmatrix}U \\ V\end{bmatrix}^{k+1}
  = \begin{bmatrix}U \\ V\end{bmatrix}^k + \frac{1}{6}\left(\vec{K}_1 + 2\vec{K}_2 + 2\vec{K}_3 + \vec{K}_4\right)h_t
  \label{eq_9}
  \end{equation}
  注意, 上式$\eqref{eq_9}$中$U$、$V$之初始值需以初始条件式$\eqref{eq_6}$填充. 其中,
  \begin{gather*}
  \vec{F}(U, V) = \begin{bmatrix} V \\ D\left( \frac{UA_x}{h_x^2} + \frac{A_yU}{h_y^2} \right) \end{bmatrix} =
  \begin{bmatrix} F_1 \\ F_2 \end{bmatrix}\\
  \vec{K}_1 = \vec{F}(U, V) \\
  \vec{K}_2 = \vec{F}(U + \frac{h_t}{2}F_1^{(1)}, V+\frac{h_t}{2}F_2^{(1)}) \\
  \vec{K}_3 = \vec{F}(U + \frac{h_t}{2}F_1^{(2)}, V+\frac{h_t}{2}F_2^{(2)}) \\
  \vec{K}_4 = \vec{F}(U + h_tF_1^{(3)}, V + h_tF_2^{(3)}) \\
  \end{gather*}
  其中, 上标$^{(1)}$、$^{(2)}$、$^{(3)}$分别代表相应分量来自于$\vec{K}_1$、$\vec{K}_2$、$\vec{K}_3$. 由此, 根据式$\eqref{eq_9}$即可完成式$\eqref{eq_1}$2维波动方程的数值求解.
• 代码实现:
  Part Ⅰ:
  现以如下初始条件为例进行算法实施,
  \begin{equation*}
  u_0(x,y) = 0.1\mathrm{e}^{-((x-x_c)^2 + (y-y_c)^2)/0.0005}
  \end{equation*}
  Part Ⅱ:
  利用finite difference求解2D Wave equation, 实现如下,
  

    1 # 2D Wave equation之求解 - 有限差分法
    2 
    3 import numpy
    4 from matplotlib import pyplot as plt
    5 from matplotlib import animation
    6 
    7 
    8 class WaveEq(object):
    9     
   10     def __init__(self, nx=300, ny=300, nt=1000, Lx=1, Ly=1, Lt=6):
   11         self.__nx = nx                        # x轴子区间划分数
   12         self.__ny = ny                        # y轴子区间划分数
   13         self.__nt = nt                        # t轴子区间划分数
   14         self.__Lx = Lx                        # x轴半长
   15         self.__Ly = Ly                        # y轴半长
   16         self.__Lt = Lt                        # t轴全长
   17         self.__D = 0.1
   18         
   19         self.__hx = 2 * Lx / nx
   20         self.__hy = 2 * Ly / ny
   21         self.__ht = Lt / nt
   22         
   23         self.__X, self.__Y = self.__build_gridPoints()
   24         self.__T = numpy.linspace(0, Lt, nt + 1)
   25         self.__Ax, self.__Ay = self.__build_2ndOrdMat()
   26         
   27         
   28     def get_solu(self):
   29         '''
   30         数值求解
   31         '''
   32         UList = list()
   33         
   34         U0 = self.__get_initial_U0()
   35         V0 = self.__get_initial_V0()
   36         self.__fill_boundary_U(U0)
   37         UList.append(U0)
   38         for t in self.__T[:-1]:
   39             # print(t, numpy.max(U0))
   40             Uk, Vk = self.__calc_Uk_Vk(t, U0, V0)
   41             UList.append(Uk)
   42             U0, V0 = Uk, Vk
   43         
   44         return self.__X, self.__Y, self.__T, UList
   45             
   46     
   47     def __calc_Uk_Vk(self, t, U0, V0):
   48         K1 = self.__calc_F(U0, V0)
   49         tmpU, tmpV = U0 + self.__ht / 2 * K1[0], V0 + self.__ht / 2 * K1[1]
   50         self.__fill_boundary_U(tmpU)
   51         
   52         K2 = self.__calc_F(tmpU, tmpV)
   53         tmpU, tmpV = U0 + self.__ht / 2 * K2[0], V0 + self.__ht / 2 * K2[1]
   54         self.__fill_boundary_U(tmpU)
   55         
   56         K3 = self.__calc_F(tmpU, tmpV)
   57         tmpU, tmpV = U0 + self.__ht * K3[0], V0 + self.__ht * K3[1]
   58         self.__fill_boundary_U(tmpU)
   59         
   60         K4 = self.__calc_F(tmpU, tmpV)
   61         
   62         Uk = U0 + (K1[0] + 2 * K2[0] + 2 * K3[0] + K4[0]) / 6 * self.__ht
   63         Vk = V0 + (K1[1] + 2 * K2[1] + 2 * K3[1] + K4[1]) / 6 * self.__ht
   64         self.__fill_boundary_U(Uk)
   65         return Uk, Vk
   66         
   67         
   68     def __calc_F(self, U0, V0):
   69         F1 = V0
   70         term0 = numpy.matmul(U0, self.__Ax) / self.__hx ** 2
   71         term1 = numpy.matmul(self.__Ay, U0) / self.__hy ** 2
   72         F2 = self.__D * (term0 + term1)
   73         return F1, F2
   74     
   75         
   76     def __fill_boundary_U(self, U):
   77         '''
   78         填充边界条件
   79         '''
   80         U[0, :] = 0
   81         U[-1, :] = 0
   82         U[:, 0] = 0
   83         U[:, -1] = 0
   84         
   85         
   86     def __get_initial_V0(self):
   87         '''
   88         获取V之初始条件
   89         '''
   90         V0 = numpy.zeros(self.__X.shape)
   91         return V0
   92         
   93         
   94     def __get_initial_U0(self):
   95         '''
   96         获取U之初始条件
   97         '''
   98         xc = 0
   99         yc = 0
  100         U0 = 0.1 * numpy.exp(-((self.__X - xc) ** 2 + (self.__Y - yc) ** 2) / 0.0005)
  101         return U0
  102         
  103         
  104     def __build_2ndOrdMat(self):
  105         '''
  106         构造2阶微分算子的矩阵形式
  107         '''
  108         Ax = self.__build_AMat(self.__nx + 1)
  109         Ay = self.__build_AMat(self.__ny + 1)
  110         return Ax, Ay
  111         
  112         
  113     def __build_AMat(self, n):
  114         term0 = numpy.identity(n) * -2
  115         term1 = numpy.eye(n, n, 1)
  116         term2 = numpy.eye(n, n, -1)
  117         AMat = term0 + term1 + term2
  118         return AMat
  119         
  120         
  121     def __build_gridPoints(self):
  122         '''
  123         构造网格节点
  124         '''
  125         X = numpy.linspace(-self.__Lx, self.__Lx, self.__nx + 1)
  126         Y = numpy.linspace(-self.__Ly, self.__Ly, self.__ny + 1)
  127         X, Y = numpy.meshgrid(X, Y)
  128         return X, Y
  129         
  130         
  131 class WavePlot(object):
  132     
  133     fig = None
  134     ax1 = None
  135     line = None
  136     X, Y, T, Z = None, None, None, None
  137     
  138     @classmethod
  139     def plot_ani(cls, waveObj):
  140         cls.X, cls.Y, cls.T, cls.Z = waveObj.get_solu()
  141         
  142         cls.fig = plt.figure(figsize=(10, 10), dpi=40)
  143         cls.ax1 = plt.subplot()
  144         cls.line = cls.ax1.pcolormesh(cls.X, cls.Y, cls.Z[0][:-1, :-1], cmap="jet")
  145         
  146         ani = animation.FuncAnimation(cls.fig, cls.update, frames=range(1, len(cls.T), 25), blit=True, interval=5, repeat=True)
  147         
  148         ani.save("plot_ani.gif", writer="PillowWriter", fps=200)
  149         plt.close()
  150         
  151     
  152     @classmethod
  153     def update(cls, frame):
  154         cls.line.set_array(cls.Z[frame][:-1, :-1])
  155         cls.line.set_norm(norm=None)
  156         return cls.line,
  157         
  158         
  159 if __name__ == "__main__":
  160     nx = 300
  161     ny = 300
  162     nt = 4000
  163     
  164     Lx = 0.5
  165     Ly = 0.5
  166     Lt = 10
  167     waveObj = WaveEq(nx, ny, nt, Lx, Ly, Lt)
  168     # waveObj.get_solu()
  169     
  170     WavePlot.plot_ani(waveObj)

  View Code
• 结果展示:
  注意, 以上为幅值分布情况, 为加强比对, 每一帧之色彩分布均是相对的.
• 使用建议:
  ①. 利用变数变换降阶处理高阶PDE, 转换为低阶PDE系统便利离散化.
• 参考文档:
  吴一东. 数值方法7:偏微分方程5:双曲型微分方程. bilibili, 2020.

标签：begin,.__,partial,Equation,equation,2D,end,Finite,self
来源： https://www.cnblogs.com/xxhbdk/p/14614327.html