問題設定

空間座標 $x,y$ を変数に持つ関数 $u(x,y)$ に関するラプラス方程式

\[\cfrac{\partial^2 u}{\partial x^2} + \cfrac{\partial^2 u}{\partial y^2} = 0 \tag{1}\]

を数値的に解きたい。

理論

2次元空間のラプラス方程式

連立方程式の導出

  • 座標空間に微小な間隔 $\Delta x, \Delta y$ で格子点 $x_1, \cdots, x_{N_x}$ および $y_1, \cdots, y_{N_y}$ を取る
    • $x_{m+1} = x_m + \Delta x$
    • $y_{l+1} = y_l + \Delta y$

ラプラス方程式はポアソン方程式の特殊ケースであるから、ポアソン方程式と同じ議論を経て、$\rho(x,y)=0,\,\Delta x = \Delta y$ とすれば、

\[-4U_m^l + U_{m+1}^l + U_{m-1}^l + U_m^{l+1} + U_m^{l-1} = 0 \qquad \left(U_m^l := u(x_m, y_l)\right) \tag{2}\]

微分方程式の解を得るには、全ての格子点に関して $(2)$ 式を立てて、連立方程式として解けば良い。

境界条件の設定

ラプラス方程式と同様、各格子点の連立方程式を立てるには隣接4点の格子点が必要となる。
下図の通り、考える領域内の格子点（図の青色）の1つ外側に仮想的な格子点（図のオレンジ色）を取り、これらの値を 境界条件 として与えることで、領域内の全ての格子点に隣接格子点を持たせる。

laplace-eq-grid

具体例

下図の例を考えてみる。

  • $1 \le x \le 3,\,1 \le y \le 4$
  • $N_x = 3,\,N_y = 4$
  • $\Delta x = \Delta y = 1$

laplace-eq-small

ここに境界条件

\[\begin{cases} u(x_0, y_l) &=& U_0^l &=& 0 \\ u(x_4, y_l) &=& U_4^l &=& 1 \\ u(x_m, y_0) &=& U_m^0 &=& 2 \\ u(x_m, y_5) &=& U_m^5 &=& 3 \end{cases}\]

を課す。
全ての格子点に関して $(2)$ 式を立てて、境界条件から得られた定数項を右辺に移項すると、

\[\begin{cases} -4U_1^1 + U_2^1 + U_1^2 &=& -2 \\ -4U_1^2 + U_2^2 + U_1^3 + U_1^1 &=& 0 \\ -4U_1^3 + U_2^3 + U_1^4 + U_1^2 &=& 0 \\ -4U_1^4 + U_2^4 + U_1^3 &=& -3 \\ -4U_2^1 + U_3^1 + U_1^1 + U_2^2 &=& -2 \\ -4U_2^2 + U_3^2 + U_1^2 + U_2^3 + U_2^1 &=& 0 \\ -4U_2^3 + U_3^3 + U_1^3 + U_2^4 + U_2^2 &=& 0 \\ -4U_2^4 + U_3^4 + U_1^4 + U_2^3 &=& -3 \\ -4U_3^1 + U_2^1 + U_3^2 &=& -3 \\ -4U_3^2 + U_2^2 + U_3^3 + U_3^1 &=& -1 \\ -4U_3^3 + U_2^3 + U_3^4 + U_3^2 &=& -1 \\ -4U_3^4 + U_2^4 + U_3^3 &=& -4 \end{cases}\]

これは $N_x \times N_y = 24$ 個の未知変数 $U_i^j$ に関する連立1次方程式となっている。
行列形式で書くと、

\[\begin{pmatrix} -4 & 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & -4 & 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & -4 & 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & -4 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 & -4 & 1 & 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 1 & -4 & 1 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 1 & -4 & 1 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 1 & -4 & 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & -4 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & -4 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & -4 & 1 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & -4 \end{pmatrix} \begin{pmatrix} U_1^1 \\ U_1^2 \\ U_1^3 \\ U_1^4 \\ U_2^1 \\ U_2^2 \\ U_2^3 \\ U_2^4 \\ U_3^1 \\ U_3^2 \\ U_3^3 \\ U_3^4 \end{pmatrix} = \begin{pmatrix} -2 \\ 0 \\ 0 \\ -3 \\ -2 \\ 0 \\ 0 \\ -3 \\ -3 \\ -1 \\ -1 \\ -4 \end{pmatrix}\]

ここからは手計算は辛いのでプログラムで計算：

A = np.matrix([
    [-4,1,0,0,1,0,0,0,0,0,0,0],
    [1,-4,1,0,0,1,0,0,0,0,0,0],
    [0,1,-4,1,0,0,1,0,0,0,0,0],
    [0,0,1,-4,0,0,0,1,0,0,0,0],
    [1,0,0,0,-4,1,0,0,1,0,0,0],
    [0,1,0,0,1,-4,1,0,0,1,0,0],
    [0,0,1,0,0,1,-4,1,0,0,1,0],
    [0,0,0,1,0,0,1,-4,0,0,0,1],
    [0,0,0,0,1,0,0,0,-4,1,0,0],
    [0,0,0,0,0,1,0,0,1,-4,1,0],
    [0,0,0,0,0,0,1,0,0,1,-4,1],
    [0,0,0,0,0,0,0,1,0,0,1,-4]
])
v = np.matrix([-2,0,0,-3,-2,0,0,-3,-3,-1,-1,-4]).T
U = (np.linalg.inv(A) * v).reshape(3,4)
"""
                               (x=0)
            0           0           0           0
       2  [[1.04883455, 0.78234411, 0.889871  , 1.44668401],   3
(y=0)  2   [1.41299409, 1.19067091, 1.33045585, 1.89686506],   3  (y=5)
       2   [1.41247092, 1.23688957, 1.34441645, 1.81032038]])  3
            1           1           1           1
                               (x=4)
"""

# 検算：隣接格子点4つの平均値になっているか？
print(U[0,0], (0+2+U[0,1]+U[1,0])/4)
# 1.0488345517877549 1.0488345517877549
print(U[1,2], (U[1,1]+U[1,3]+U[0,2]+U[2,2])/4)
# 1.3304558533999695 1.3304558533999695
print(U[2,3], (U[2,2]+U[1,3]+3+1)/4)
# 1.8103203777897097 1.81032037778971

実装

ラプラス方程式：

\[\cfrac{\partial^2 u}{\partial x^2} + \cfrac{\partial^2 u}{\partial y^2} = 0\]

境界条件：

\[\begin{cases} u(x_\mathrm{min}, y) &=& 0 \\ u(x_\mathrm{max}, y) &=& 0 \\ u(x, y_\mathrm{min}) &=& 0 \\ u(x, y_\mathrm{max}) &=& 100 \end{cases}\] 
solver = LaplaceEquationSolver()
solver.solve(x1_range=(2,4), x2_range=(0,1), dx=0.02)
solver.draw_result_heatmap()
solver.draw_result_3d()
ヒートマップ 3D
laplace-eq-heatmap laplace-eq-3d