2D Poisson's Equation - Finite Difference

• 二维泊松方程边值问题:
  静电场中二维泊松方程如下,
  \begin{equation}
  \left( \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2} \right)u = -\frac{\rho(x, y)}{\varepsilon_0}, \quad (x, y)\in\Omega = \{(x, y) | \frac{x^2}{a^2} + \frac{y^2}{b^2} \leq 1 \}
  \label{eq_1}
  \end{equation}
  边界条件如下,
  \begin{equation}
  u(x, y) = \mu(\theta), \quad (x, y) \in \partial\Omega
  \label{eq_2}
  \end{equation}
  其中, $\theta$为边界方位角.
• 连续型系统的离散化:
  由于有限差分法适合处理矩形网格, 本文拟采用边界填充技术处理非矩形边界, 则二维泊松方程式$\eqref{eq_1}$及边界条件式$\eqref{eq_2}$分别转换如下,
  \begin{gather}
  \left( \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2} \right)u = -\frac{\rho(x, y)}{\varepsilon_0}, \quad (x, y)\in\Omega = [-a, a] \times [-b, b]
  \label{eq_3}   \\
  u(x, y) = \mu(\theta), \quad (x, y) \in \Omega \cap \{ (x, y) | \frac{x^2}{a^2} + \frac{y^2}{b^2} \geq 1 \}
  \label{eq_4}
  \end{gather}
  由于大型线性方程组求解的本质困难, 本文拟采用含时演化技术将椭圆型方程转化为抛物型方程进行求解. 如果含时演化的抛物型方程最终趋于稳定, 则收敛于相应椭圆型方程的解. 式$\eqref{eq_3}$泊松方程转换后的抛物型方程如下,
  \begin{equation}
  \frac{\partial u}{\partial t} = \left( \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2} \right)u + \frac{\rho(x, y)}{\varepsilon_0}, \quad (x, y)\in\Omega = [-a, a] \times [-b, b], t \in [0, +\infty]
  \label{eq_5}
  \end{equation}
  任意给定初始条件,
  \begin{equation}
  u(x, y, 0) = u_0(x, y), \quad (x, y) \in \Omega
  \label{eq_6}
  \end{equation}
  结合式$\eqref{eq_4}$给定边界条件,
  \begin{equation}
  u(x, y, t) = \mu(\theta), \quad (x, y) \in \Omega \cap \{ (x, y) | \frac{x^2}{a^2} + \frac{y^2}{b^2} \geq 1 \}, t \in [0, +\infty]
  \label{eq_7}
  \end{equation}
  式$\eqref{eq_1}$、$\eqref{eq_2}$静电场中泊松方程之求解转换为式$\eqref{eq_5}$、$\eqref{eq_6}$、$\eqref{eq_7}$抛物型方程之求解.
  本文拟采用中心差分离散化式$\eqref{eq_5}$,
  \begin{equation}
  \frac{\partial u}{\partial t} = \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) + \frac{\rho_{i, j}}{\varepsilon_0}
  \label{eq_8}
  \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}
  \rho_{0, 0} & \cdots & \rho_{0, n_x} \\
  \vdots & \ddots & \vdots \\
  \rho_{n_y, 0} & \cdots & \rho_{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_8}$转换如下,
  \begin{equation}
  \frac{\partial U}{\partial t} = \left( \frac{UA_x}{h_x^2} + \frac{A_yU}{h_y^2} \right) + \frac{V}{\varepsilon_0}
  \label{eq_9}
  \end{equation}
  注意, 上式$\eqref{eq_9}$中$U$之边界值需以边界条件式$\eqref{eq_7}$填充. 本文拟采用Runge-Kutta方法求解上式$\eqref{eq_9}$, 则有,
  \begin{equation}
  U^{k+1} = U^k + \frac{1}{6}(K_1 + 2K_2 + 2K_3 + K_4)h_t
  \label{eq_10}
  \end{equation}
  注意, 上式$\eqref{eq_10}$中$U$之初始值需以初始条件式$\eqref{eq_6}$填充. 其中,
  \begin{gather*}
  F(U) = \left( \frac{UA_x}{h_x^2} + \frac{A_yU}{h_y^2} \right) + \frac{V}{\varepsilon_0} \\
  K_1 = F(U) \\
  K_2 = F(U + \frac{h_t}{2}K_1) \\
  K_3 = F(U + \frac{h_t}{2}K_2) \\
  K_4 = F(U + h_tK_3)
  \end{gather*}
  由此, 根据式$\eqref{eq_10}$即可完成式$\eqref{eq_5}$抛物型方程的数值求解. 若,
  \begin{equation}
  \|U^{k+1} - U^k\| < \epsilon
  \end{equation}
  其中, $\epsilon$取值足够小正数, 则$U^k$趋于稳定并收敛于式$\eqref{eq_1}$的解.
  
• 代码实现:
  Part Ⅰ:
  现以如下电荷密度、初始条件及边界条件为例进行算法实施,
  \begin{gather*}
  \rho(x, y) = 100\mathrm{e}^{-100[(x-c)^2+y^2]}-100\mathrm{e}^{-100[(x+c)^2+y^2]},\quad \text{where }c = \sqrt{a^2 - b^2} \\
  u_0(x, y) = 0 \\
  \mu(\theta) = \cos(5\theta)
  \end{gather*}
  Part Ⅱ:
  利用finite difference求解2D Poisson‘s equation, 实现如下,
  
  

    1 # Poisson‘s equation之求解 - 有限差分法
    2 
    3 import numpy
    4 from matplotlib import pyplot as plt
    5 
    6 
    7 class PoissonEq(object):
    8     
    9     def __init__(self, nx=150, ny=100, a=3, b=2):
   10         self.__nx = nx                        # x轴子区间划分数
   11         self.__ny = ny                        # y轴子区间划分数
   12         self.__Lx = a                         # x轴椭圆半长
   13         self.__Ly = b                         # y轴椭圆半长
   14         self.__epsilon = 1                    # 真空介电常数取值
   15         
   16         assert a > b, "a must larger than b"
   17         
   18         self.__hx = 2 * a / nx
   19         self.__hy = 2 * b / ny
   20         self.__ht = min([self.__hx, self.__hy]) ** 3
   21         
   22         self.__X, self.__Y = self.__build_gridPoints()
   23         self.__Ax, self.__Ay = self.__build_2ndOrdMat()
   24         self.__V = self.__get_V()
   25         
   26         self.__OuterLoc = self.__X ** 2 / a ** 2 + self.__Y ** 2 / b ** 2 >= 1
   27         self.__InnerLoc = self.__X ** 2 / a ** 2 + self.__Y ** 2 / b ** 2 <= 1
   28         self.__Angle = numpy.angle(self.__X + 1j * self.__Y)
   29         
   30         
   31     def get_solu(self, maxIter=1000000, varepsilon=1.e-9):
   32         ‘‘‘
   33         数值求解
   34         maxIter: 最大迭代次数
   35         varepsilon: 收敛判据
   36         ‘‘‘
   37         U0 = self.__get_initial_U()
   38         self.__fill_boundary_U(U0)
   39         for i in range(maxIter):
   40             t = i * self.__ht
   41             Uk = self.__calc_Uk(t, U0)
   42             
   43             # print(i, numpy.linalg.norm(Uk - U0, numpy.inf))
   44             if self.__isConverged(U0, Uk, varepsilon):
   45                 break
   46             
   47             U0 = Uk
   48         else:
   49             raise Exception(">>> Not converged after {} iterations!<<<".format(maxIter))
   50         
   51         return numpy.ma.array(self.__X, mask=~self.__InnerLoc), numpy.ma.array(self.__Y, mask=~self.__InnerLoc), numpy.ma.array(Uk, mask=~self.__InnerLoc)
   52         
   53         
   54     def get_Efield(self, U):
   55         ‘‘‘
   56         获取电场
   57         ‘‘‘
   58         Ey, Ex = numpy.gradient(U, self.__hy, self.__hx)
   59         return numpy.linspace(-self.__Lx, self.__Lx, self.__nx + 1), numpy.linspace(-self.__Ly, self.__Ly, self.__ny + 1), -Ex, -Ey
   60                         
   61         
   62     def __isConverged(self, U0, Uk, varepsilon):
   63         normVal = numpy.linalg.norm(Uk - U0, numpy.inf)
   64         if normVal < varepsilon:
   65             return True
   66         return False
   67             
   68             
   69     def __calc_Uk(self, t, U0):
   70         K1 = self.__calc_F(U0)
   71         Uk = U0 + K1 * self.__ht
   72         self.__fill_boundary_U(Uk)
   73         return Uk
   74         
   75         
   76     def __calc_F(self, U):
   77         term0 = numpy.matmul(U, self.__Ax) / self.__hx ** 2
   78         term1 = numpy.matmul(self.__Ay, U) / self.__hy ** 2
   79         term2 = self.__V / self.__epsilon
   80         FVal = term0 + term1 + term2
   81         return FVal
   82         
   83         
   84     def __fill_boundary_U(self, U):
   85         ‘‘‘
   86         填充边界条件
   87         ‘‘‘
   88         U[self.__OuterLoc] = numpy.cos(self.__Angle[self.__OuterLoc])
   89         
   90         
   91     def __get_initial_U(self):
   92         ‘‘‘
   93         获取初始条件
   94         ‘‘‘
   95         U0 = numpy.zeros(self.__X.shape)
   96         return U0
   97         
   98         
   99     def __get_V(self):
  100         ‘‘‘
  101         获取电荷密度
  102         ‘‘‘
  103         c = numpy.math.sqrt(self.__Lx ** 2 - self.__Ly ** 2)
  104         term0 = -100 * ((self.__X - c) ** 2 + self.__Y ** 2)
  105         term1 = -100 * ((self.__X + c) ** 2 + self.__Y ** 2)
  106         V = 100 * numpy.exp(term0) - 100 * numpy.exp(term1)
  107         return V
  108         
  109         
  110     def __build_2ndOrdMat(self):
  111         ‘‘‘
  112         构造2阶微分算子的矩阵形式
  113         ‘‘‘
  114         Ax = self.__build_AMat(self.__nx + 1)
  115         Ay = self.__build_AMat(self.__ny + 1)
  116         return Ax, Ay
  117         
  118         
  119     def __build_AMat(self, n):
  120         term0 = numpy.identity(n) * -2
  121         term1 = numpy.eye(n, n, 1)
  122         term2 = numpy.eye(n, n, -1)
  123         AMat = term0 + term1 + term2
  124         return AMat
  125         
  126         
  127     def __build_gridPoints(self):
  128         ‘‘‘
  129         构造网格节点
  130         ‘‘‘
  131         X = numpy.linspace(-self.__Lx, self.__Lx, self.__nx + 1)
  132         Y = numpy.linspace(-self.__Ly, self.__Ly, self.__ny + 1)
  133         X, Y = numpy.meshgrid(X, Y)
  134         return X, Y
  135         
  136         
  137 class PoissonPlot(object):
  138     
  139     @staticmethod
  140     def plot_fig(poiObj):
  141         maxIter = 1000000
  142         varepsilon = 1.e-9
  143         
  144         X, Y, U = poiObj.get_solu(maxIter, varepsilon)
  145         X1, Y1, Ex, Ey = poiObj.get_Efield(U)
  146         
  147         fig = plt.figure(figsize=(6, 4))
  148         ax1 = plt.subplot()
  149         
  150         cs1 = ax1.contour(X, Y, U, levels=50, cmap="jet")
  151         ax1.streamplot(X1, Y1, Ex, Ey, density=1, linewidth=1, color="black")
  152         fig.colorbar(cs1)
  153         fig.tight_layout()
  154         fig.savefig("plot_fig.png", dpi=100)
  155         
  156 
  157 if __name__ == "__main__":
  158     nx = 150
  159     ny = 100
  160     a = 1.5
  161     b = 1
  162     poiObj = PoissonEq(nx, ny, a, b)
  163     
  164     PoissonPlot.plot_fig(poiObj)

  View Code
• 结果展示:
  注意, 图中等高线代表等势线, 流线代表电场走向.
• 使用建议:
  ①. 利用边界填充技术处理非矩形边界;
  ②. 利用含时演化技术求解椭圆型方程;
  ③. 使用Euler向前差分替代Runge-Kutta方法获取含时演化抛物型方程的稳态解.
• 参考文档:
  吴一东. 数值方法6: 偏微分方程4: 二维拉普拉斯方程，泊松方程. bilibili, 2020.

2D Poisson's Equation - Finite Difference

原文：https://www.cnblogs.com/xxhbdk/p/14534909.html