【数値解析】１次元Poisson方程式：境界条件を考慮したGreen関数を用いた厳密解
このノートは以下のノートに関係するものです．
復習：Poisson方程式とGreen関数による解
ポアソン（Poisson）方程式
$$
\nabla^{2} \phi = - f
$$
について，Green関数
$$
\nabla^{2} G = - \delta(\bm{r})
$$
は1-3次元でそれぞれ以下のようになる：
ただし，$${x, \rho, r}$$は各次元でのユークリッドノルムである：
これを用いて，解は式（３）のように書ける：
復習： Poisson方程式の解における計算量
この Green 関数を使って，数値的に解を得る時を考える．
空間を $${N}$$ 点のサンプリングをしたとしよう．
式（３）をそのまま計算しようとすると，
各位置での$${\phi(\bm{r})}$$を，$${\bm{r}^{\prime}}$$に関する積分によって得るので，$${\mathscr{O}(N^2)}$$ の計算量になる．
一方で，解の積分を$${\phi = G * f}$$のように畳み込み積分の記号で書けば，Fourier変換 $${X = \mathcal{F}(x)}$$によって，以下の式（４）のように簡単に計算できるようになる：
畳み込み積分をFourier変換の積によって得る．
これらを使って，数値解を得てみよう．
境界条件を考慮する
先の議論では，境界条件が考慮されていないので，
境界条件を考慮しなくてはならない．
今回は以下の条件で解析してみよう．
1D Poisson 方程式 を扱う：
$$
\begin{equation*}
\begin{split}
\frac{d^2}{dx^2} \phi(x) &= - f(x), \phi(0) = \phi(L) = 0,
\end{split}
\end{equation*}
$$
境界条件条件付きGreen関数
1D Poisson 方程式についてGreen関数は，
$$
\begin{equation*}
\begin{split}
\frac{d^2}{dx^2} G(x, x') &= - \delta(x-x'), \\
\phi(x) &= \int_{0}^{L} dx' G(x,x') f(x')
\end{split}
\end{equation*}
$$
を満たす．
ポイント：Green関数に境界条件を適用
これより，境界条件から，以下の条件を得る．
さらに，$${x \neq x'}$$ において，
$$
\frac{d^2}{dx^2} G(x, x') = 0
$$
であるから，$${x = x'}$$ について，
下側$${x<x'}$$と上側$${x>x'}$$で分けて考えると
$$
\begin{equation*}
\begin{split}
G_{<}(x,x') &= a_{<}(x') x + b_{<}(x') \\
G_{>}(x,x') &= a_{>}(x') x + b_{>}(x')
\end{split}
\end{equation*}
$$
をそれぞれ得る．
式（５）の条件から，
$$
G_{<}(0,x') = 0 \iff b_{<}(x') = 0
$$
$$
G_{>}(L,x') = 0 \iff b_{>}(x') = - a_{>}(x') L
$$
が得られる．
ポイント：Green関数の支配式を積分，１階微分の式にする
Green 関数は，
$$
\int_{0}^{L} dx \frac{d^2}{dx^2} G(x, x')
= \frac{d}{dx} G_{>}(L, x') - \frac{d}{dx} G_{<}(0, x') ,
- \int_{0}^{L} dx \delta(x-x') = -1
$$
であるから，
$$
a_{>}(x') - a_{<}(x') = -1
$$
ポイント：Green関数の連続性を与える
さらに，Green関数が $${x = x'}$$で連続
$$
G_{<}(x', x') = G_{>}(x', x') \iff a_{<}(x') x' = a_{>}(x') (x' - L)
$$
したがって，
$$
a_{<}(x') = \frac{x' - L}{L}, a_{>}(x') = \frac{x'}{L}
$$
以下のGreen関数を得る：
これで数値計算が可能になるだろう．
畳み込み積分を直接計算
台形公式やシンプソン公式などを用いて直接積分してみるのが一つの方法だ．
この計算量はすでに議論したように，$${\mathscr{O}(N^2)}$$になる．
式（７）を用いて直接数値積分する．
Fourier変換を用いて計算できる？
境界条件をつけると，ほとんどの場合について，Green関数について，
$$
G(x, x') \neq G(x - x')
$$
となってしまい，畳み込み積分が使える状況にならない．
したがって，FFTの利用ができないことがわかる．
計算の効率化はできないようだ．
（筆者は他の方法をまだ知らない．）
1次元ポアソン方程式の厳密解（数値計算結果）
まずは，式（７）の Green 関数の設定は以下のようにできる：
import numpy as np
def green_1d_dirichlet(x, xi, L):
"""
1次元 Poisson 方程式 d^2/dx^2 phi = -f(x),
境界条件 phi(0)=0, phi(L)=0 のグリーン関数 G(x, xi) を返す関数
Parameters
----------
x : float または ndarray
観測点の座標
xi : float または ndarray
ソース点(デルタ関数の位置)の座標
L : float
問題の区間長 (0 < x < L)
Returns
-------
G : float または ndarray
グリーン関数 G(x, xi)
"""
# x, xi がともに配列の場合も扱えるように numpy の配列として扱う
x_arr = np.array(x, dtype=float, copy=False, ndmin=1)
xi_arr = np.array(xi, dtype=float, copy=False, ndmin=1)
# x_arr と xi_arr の形を揃える (ブロードキャストを利用)
# たとえば x がスカラー, xi がベクトルの場合などに対応
x_b, xi_b = np.broadcast_arrays(x_arr, xi_arr)
# G(x, xi) を要素ごとに計算する
# x <= xi のとき x*(L-xi)/L, それ以外は xi*(L-x)/L
# 条件: x_b <= xi_b
G_val = np.where(
x_b <= xi_b,
x_b*(L - xi_b)/L,
xi_b*(L - x_b)/L
)
# x, xi がスカラーだった場合はスカラーを返す
# さもなくば ndarray
if np.isscalar(x) and np.isscalar(xi):
return float(G_val)
else:
return G_val
# --- 以下，動作確認用の例 ---
if __name__ == "__main__":
# 例: x と xi をそれぞれ 0 から L=1 までの 50 分割点で作り，G(x, xi) を可視化
import matplotlib.pyplot as plt
L = 1.0
N = 50
x_list = np.linspace(0, L, N)
xi_list = np.linspace(0, L, N)
# 2次元格子上でグリーン関数を評価する (x を縦軸, xi を横軸とする)
X, XI = np.meshgrid(x_list, xi_list, indexing='ij')
G_2d = green_1d_dirichlet(X, XI, L) # shape = (N, N)
plt.figure(figsize=(6,5))
plt.pcolormesh(xi_list, x_list, G_2d, shading='auto')
plt.colorbar(label='G(x, xi)')
plt.xlabel('xi')
plt.ylabel('x')
plt.title('1D Poisson Green Function (Dirichlet BC)')
plt.show()
このコードを利用して直接，式（８）の積分を解こう．
台形公式による数値積分 np.trapz を使う．
def green_1d_dirichlet_array(N, dx, L):
"""
離散点 x_i = i*dx (i=0..N-1) における1次元Poisson方程式用の
グリーン関数離散値 G[i] = G(x_i), ただし x_i in [0, L], BC: Dirichlet
G(x) = x(L-x)/L ( 0 < x < L ), G(0)=G(L)=0
"""
G = np.zeros(N)
for i in range(N):
x_i = i*dx
G[i] = x_i*(L - x_i)/L # 端点含め外では 0
return G
# パラメータ設定
L = 10.0 # 計算領域の長さ
N = 256 # 格子点数
dx = L / (N - 1) # 格子間隔
# x軸の離散点
x = np.linspace(0, L, N)
# f(x) の設定（例: ガウス関数）
sigma = 0.5 # 幅
f_x = 10.0*np.exp(-((x - L/2.0) ** 2) / (2 * sigma**2))
f_x += -10.0*np.exp(-((x - L/2.0 + L / 5.0) ** 2) / (2 * sigma**2))
f_x += 8.0*np.exp(-((x - L/10.0) ** 2) / (2 * sigma**2))
f_x += 10.0*np.exp(-((x - L + L/10.0) ** 2) / (2 * sigma**2))
f_x /= np.trapz(f_x, x) # 規格化
phi_conv_direct = np.zeros(N)
for i in range(N):
# G(x_i, x_j)
G_ij = green_1d_dirichlet(x[i], x, L)
# Simpson/Trapezoidal などで数値積分
# ここでは簡単に台形公式 (np.trapz) を使用
phi_conv_direct[i] = np.trapz(G_ij * f_x, x)
# ----- 結果のプロット -----
plt.figure(figsize=(8, 4))
plt.plot(x, f_x, label="$f(x)$ (source term)")
plt.plot(x, phi_conv_direct, label="$\\phi(x)$ (Direct np.trapz)", linestyle='dotted')
plt.axhline(0, color='gray', linestyle='--', linewidth=0.8)
plt.xlabel("x")
plt.ylabel("Function values")
plt.legend()
plt.title("1D Poisson Equation Solution (Dirichlet BCs)")
plt.show()
式（７）と式（８）のGreen関数による解を利用．
まとめ
今回は，Dirichletの境界条件（$${\phi(0) = \phi(L) = 0}$$）が課された
１次元Poisson方程式を例にして，任意の $${f(x)}$$ の分布が与えられた時の厳密解を，直接数値計算によって得た．
次回は，これが，Jacobi 法などの反復法の結果と一致するのかをみてみよう．
そして，次はどこまで行けるか，まだわからないが，
multigrid法など反復法の収束を早めるための工夫まで説明できたら良いなと．