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, 实现如下,
    2D Wave Equation (2) - Finite Difference
      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
  • 结果展示:
    2D Wave Equation (2) - Finite Difference注意, 以上为幅值分布情况, 为加强比对, 每一帧之色彩分布均是相对的.
  • 使用建议:
    ①. 利用变数变换降阶处理高阶PDE, 转换为低阶PDE系统便利离散化.
  • 参考文档:
    吴一东. 数值方法7:偏微分方程5:双曲型微分方程. bilibili, 2020.
上一篇:Luogu P6042 「ACOI2020」学园祭 题解


下一篇:CSS3-2D变换与3D变换、过渡与动画一文全搞定