1D Poisson's Equation

• 算法特征:
  ①. 选定基函数; ②. 作用测试函数; ③. 求解弱形式
• 算法推导:
  Part Ⅰ: 投影之引入
  本文拟采用$L_2$范数衡量函数距离. 在给定函数空间$X$中, 获取函数$f$在区间$[a, b]$的最优逼近, 即
  egin{equation}
  min_{g in X}quad mathrm{dist}(f, g) = left ( int_a^b(f - g)^2 mathrm{d}x ight )^{1 / 2} quadRightarrowquad min_{g in X}quad int_a^b(f - g)^2 mathrm{d}x
  label{eq_1}
  end{equation}
  现令$g = g^ast + g^prime$, 则
  egin{equation}
  egin{split}
  int_a^b (f - g)^2 mathrm{d}x &= int_a^b (f - g^ast - g^prime)^2 mathrm{d}x \
  &= int_a^b (f - g^ast)^2 - 2g^prime (f - g^ast) + g^{prime 2}mathrm{d}x
  end{split}
  label{eq_2}
  end{equation}
  如果函数$g^ast$为函数$f$在空间$X$的投影, 则
  egin{equation}
  int_a^b g^prime (f - g^ast) mathrm{d}x = 0
  label{eq_3}
  end{equation}
  结合式$eqref{eq_2}$, 有
  egin{equation}
  int_a^b (f - g)^2 mathrm{d}x = int_a^b (f - g^ast)^2 mathrm{d}x + int_a^b g^{prime 2}mathrm{d}x
  label{eq_4}
  end{equation}
  由于$g$是任意的, $g^ast$是$f$的投影, 则$g^prime$是任意的, 结合式$eqref{eq_1}$与式$eqref{eq_4}$,
  egin{equation}
  g^ast = mathop{argmin}_{g in X} mathrm{dist}(f, g)
  label{eq_5}
  end{equation}
  即函数$f$在给定函数空间$X$中的最优逼近即为它的投影. 由此也将式$eqref{eq_1}$之优化问题等效为式$eqref{eq_3}$之方程求解.
  Part Ⅱ: 投影之计算
  在给定函数空间$X$(假定为$n$维线性函数空间)中选择一组完备的基函数${ phi_1, phi_2, cdots, phi_n }$, 则
  egin{equation}
  g^ast = sum_{i=1}^n a_iphi_i
  label{eq_6}
  end{equation}
  代入式$eqref{eq_3}$, 有
  egin{equation}
  int_a^b g^prime (f - sum_{i=1}^n a_iphi_i)mathrm{d}x = 0
  label{eq_7}
  end{equation}
  由于$g^prime$在投影空间中是任意的, 可取任意基函数作为测试函数, 上式转换如下,
  egin{equation}
  int_a^b phi_j (f - sum_{i=1}^n a_iphi_i)mathrm{d}x = int_a^b phi_j f mathrm{d}x - sum_{i=1}^n a_iint_a^bphi_jphi_imathrm{d}x = 0, quad mathrm{where} j = 1, cdots, n
  label{eq_8}
  end{equation}
  令,
  egin{equation*}
  A = egin{bmatrix}
  int_a^bphi_1phi_1mathrm{d}x & cdots & int_a^bphi_1phi_nmathrm{d}x   \
  vdots & ddots & vdots   \
  int_a^bphi_nphi_1mathrm{d}x & cdots & int_a^bphi_nphi_nmathrm{d}x
  end{bmatrix} quad
  b = egin{bmatrix}
  int_a^bphi_1 fmathrm{d}x   \
  vdots   \
  int_a^bphi_n fmathrm{d}x
  end{bmatrix} quad
  alpha = egin{bmatrix}
  a_1   \
  vdots   \
  a_n
  end{bmatrix}
  end{equation*}
  式$eqref{eq_8}$转换如下,
  egin{equation}
  Aalpha = b
  label{eq_9}
  end{equation}
  求解上式中$alpha$, 即完成投影之计算.
  Part Ⅲ: 1维泊松方程求解
  给定如下1维泊松方程及其边界条件,
  egin{equation}
  frac{partial^2 u}{partial x^2} + f = 0, quad u(a) = u(b) = 0
  label{eq_10}
  end{equation}
  作用任意测试函数$g(x)$, 并积分之,
  egin{equation}
  int_a^b g(frac{partial^2 u}{partial x^2} + f)mathrm{d}x = 0
  label{eq_11}
  end{equation}
  对上式$eqref{eq_11}$分部积分, 引入弱形式,
  egin{equation}
  gfrac{partial u}{partial x}igg|_a^b - int_a^bfrac{partial g}{partial x}cdot frac{partial u}{partial x}mathrm{d}x + int_a^b gfmathrm{d}x = 0
  label{eq_12}
  end{equation}
  本文拟选在分段线性空间中获取$u$的数值解. 分段线性空间常用基函数如下,
  egin{equation}
  phi_i(x) = egin{cases}
  1, quad mathrm{where} x = x_i   \
  0, quad mathrm{where} x eq x_i
  end{cases}
  label{eq_13}
  end{equation}
  其中, $x_i$表示网格节点. 本文将区间$[a, b]$剖为$n$份, 则$i = 0, 1, cdots, n$, 且$x_0 = a$, $x_n = b$.
  根据基函数的特征, 令
  egin{equation}
  u = sum_{i=0}^n u_iphi_i = u(a)phi_0 + u(b)phi_n + sum_{i=1}^{n-1} u_iphi_i = sum_{i=1}^{n-1} u_iphi_i
  label{eq_14}
  end{equation}
  取测试函数$g(x)$为基函数$phi_1, phi_2, cdots, phi_{n-1}$, 结合上式$eqref{eq_14}$, 弱形式$eqref{eq_12}$转换如下,
  egin{equation}
  -sum_{i=1}^{n-1}u_iint_a^bfrac{partial phi_j}{partial x}cdotfrac{partial phi_i}{partial x}mathrm{d}x + int_a^bphi_j fmathrm{d}x = 0, quadmathrm{where} j = 1, cdots, n-1
  label{eq_15}
  end{equation}
  令,
  egin{equation*}
  A = egin{bmatrix}
  int_a^b frac{partialphi_1}{partial x}frac{partialphi_1}{partial x}mathrm{d}x & cdots & int_a^bfrac{partial phi_1}{partial x}frac{partial phi_{n-1}}{partial x}mathrm{d}x   \
  vdots & ddots & vdots   \
  int_a^b frac{partialphi_{n-1}}{partial x}frac{partialphi_1}{partial x}mathrm{d}x & cdots & int_a^bfrac{partial phi_{n-1}}{partial x}frac{partial phi_{n-1}}{partial x}mathrm{d}x   \
  end{bmatrix}quad
  b = egin{bmatrix}
  int_a^b phi_1 fmathrm{d}x   \
  vdots   \
  int_a^b phi_{n-1} fmathrm{d}x   \
  end{bmatrix}quad
  alpha = egin{bmatrix}
  u_1   \
  vdots   \
  u_{n-1}   \
  end{bmatrix}
  end{equation*}
  式$eqref{eq_15}$转换如下,
  egin{equation}
  Aalpha = b
  label{eq_16}
  end{equation}
  求解上式中$alpha$, 代入式$eqref{eq_14}$即可获得泊松方程解$u$.
• 代码实现:
  Part Ⅰ:
  现以如下Poisson's equation为例进行算法实施,
  egin{equation*}
  frac{partial^2 u}{partial x^2} + sinpi x = 0, quad u(0) = u(1) = 0
  end{equation*}
  解析解如下,
  egin{equation*}
  u = frac{sinpi x}{pi^2}
  end{equation*}
  Part Ⅱ:
  利用finite element求解Poisson's equation, 实现如下:
  

  # Poisson's equation之实现 - 有限元法
  # 注意, 采用Gauss-Legendre Quadrature进行数值积分
  
  import numpy
  from matplotlib import pyplot as plt
  
  
  numpy.random.seed(0)
  
  
  # Gauss-Legendre Quadrature之实现
  def getGaussParams(num=6):
      if num == 1:
          X = numpy.array([0])
          A = numpy.array([2])
      elif num == 2:
          X = numpy.array([numpy.math.sqrt(1/3), -numpy.math.sqrt(1/3)])
          A = numpy.array([1, 1])
      elif num == 3:
          X = numpy.array([numpy.math.sqrt(3/5), -numpy.math.sqrt(3/5), 0])
          A = numpy.array([5/9, 5/9, 8/9])
      elif num == 6:
          X = numpy.array([0.238619186081526, -0.238619186081526, 0.661209386472165, -0.661209386472165, 0.932469514199394, -0.932469514199394])
          A = numpy.array([0.467913934574257, 0.467913934574257, 0.360761573028128, 0.360761573028128, 0.171324492415988, 0.171324492415988])
      elif num == 10:
          X = numpy.array([0.973906528517240, -0.973906528517240, 0.433395394129334, -0.433395394129334, 0.865063366688893, -0.865063366688893, 
              0.148874338981367, -0.148874338981367, 0.679409568299053, -0.679409568299053])
          A = numpy.array([0.066671344307672, 0.066671344307672, 0.269266719309847, 0.269266719309847, 0.149451349151147, 0.149451349151147, 
              0.295524224714896, 0.295524224714896, 0.219086362515885, 0.219086362515885])
      else:
          raise Exception(">>> Unsupported num = {} <<<".format(num))
      return X, A
  
  
  def getGaussQuadrature(func, a, b, X, A):
      '''
      func: 原始被积函数
      a: 积分区间下边界
      b: 积分区间上边界
      X: 求积节点数组
      A: 求积系数数组
      '''
      term1 = (b - a) / 2
      term2 = (a + b) / 2
      term3 = func(term1 * X + term2)
      term4 = numpy.sum(A * term3)
      val = term1 * term4
      return val
  
  
  class PoissonEq(object):
      
      def __init__(self, n, order=6):
          self.__n = n                                              # 子区间划分数
          self.__order = order                                      # Legendre多项式之阶数
          
          self.__X = self.__build_gridPoints(n)                     # 构造网格节点
          self.__XQuad, self.__AQuad = getGaussParams(order)        # 获取数值积分之求积节点、求积系数
          
          
      def get_solu(self):
          A = self.__build_A()
          b = self.__build_b()
          alpha = numpy.matmul(numpy.linalg.inv(A), b)
          
          U = numpy.zeros(self.__X.shape)
          U[1:-1] = alpha.flatten()
          
          return self.__X, U
          
          
      def __build_b(self):
          '''
          构造列向量b
          '''
          self.__xiMinusOne, self.__xi, self.__xiPlusOne = None, None, None
          
          b = numpy.zeros((self.__n-1, 1))
          for i in range(b.shape[0]):
              self.__xiMinusOne = self.__X[i]
              self.__xi = self.__X[i+1]
              self.__xiPlusOne = self.__X[i+2]
              quadValLeft = getGaussQuadrature(self.__funcLeft, self.__xiMinusOne, self.__xi, self.__XQuad, self.__AQuad)
              quadValRight = getGaussQuadrature(self.__funcRight, self.__xi, self.__xiPlusOne, self.__XQuad, self.__AQuad)
              b[i, 0] = quadValLeft + quadValRight
          return b
      
      
      def __funcLeft(self, x):
          val = (x - self.__xiMinusOne) / (self.__xi - self.__xiMinusOne) * numpy.sin(numpy.pi * x)
          return val
      
      
      def __funcRight(self, x):
          val = (self.__xiPlusOne - x) / (self.__xiPlusOne - self.__xi) * numpy.sin(numpy.pi * x)
          return val
          
          
      def __build_A(self):
          '''
          构造矩阵A
          '''
          dPhiLeft = 1 / (self.__X[1:-1] - self.__X[:-2])
          dPhiRight = -1 / (self.__X[2:] - self.__X[1:-1])
          interval = self.__X[1:] - self.__X[:-1]
          
          A = numpy.zeros((self.__n-1, self.__n-1))
          for i in range(A.shape[0]):
              for j in range(A.shape[1]):
                  if j == i-1:
                      A[i, j] = dPhiLeft[i] * dPhiRight[j] * interval[i]
                  elif j == i:
                      A[i, j] = dPhiLeft[i] * dPhiLeft[j] * interval[i] + dPhiRight[i] * dPhiRight[j] * interval[i+1]
                  elif j == i+1:
                      A[i, j] = dPhiRight[i] * dPhiLeft[j] * interval[i+1]
          return A
                      
          
      def __build_gridPoints(self, n):
          '''
          随机初始化网格节点
          '''
          xMin, xMax = 0, 1
          X = numpy.sort(numpy.random.uniform(xMin, xMax, n+1))
          X[0] = xMin
          X[-1] = xMax
          return X
          
          
  class PoissonPlot(object):
      
      @staticmethod
      def plot_fig(poissonObj):
          X, U = poissonObj.get_solu()
          
          xMin, xMax = numpy.min(X), numpy.max(X)
          X_ = numpy.linspace(xMin, xMax, 1000)
          U_ = numpy.sin(numpy.pi * X_) / numpy.pi ** 2
          
          fig = plt.figure(figsize=(6, 4))
          ax1 = plt.subplot()
          
          ax1.plot(X, U, marker='.', ls="--", label="numerical solution")
          ax1.plot(X_, U_, label="analytical solution")
          ax1.set(xlabel="$x$", ylabel="$u$")
          ax1.legend()
          fig.savefig("plot_fig.png", dpi=100)
          
          
          
  if __name__ == "__main__":
      n = 20
      order = 6
      poissonObj = PoissonEq(n, order)
      
      PoissonPlot.plot_fig(poissonObj)

• 结果展示:
  可以看到, 数值解与解析解是足够逼近的.
• 使用建议:
  ①. 微分方程转积分方程;
  ②. 分部积分求解弱形式.
• 参考文档:
  Numerical Methods for PDE 2016 by Professor Qiqi Wang & Jacob White from MIT