【数値計算】 ヤコビ法(1)：ポアソン方程式の計算方法
ポアソン（Poisson）方程式の計算方法について，
シンプルな反復法による解法，ヤコビ（Jacobi）法による計算について考えよう．
Poisson's equation
$$
\Delta u(\bm{r}) = -f(\bm{r})
$$
この方程式は，電荷周りのクーロンポテンシャルや，質点周りの重力ポテンシャルの計算などで扱われるものだ．
もし，$${f(\bm{r}) = 0}$$なら，$${\Delta u(\bm{r}) = 0}$$ であり，これは，ラプラス方程式（Laplace's equation）として知られる．
この時，境界条件（Boundary condition）によって，解が決定される．
境界条件は，端点の値を直接与える，ディレクレ境界条件（Dirichlet boundary condition）や，端点の微分値を与える，ノイマン境界条件（Neumann boundary condition）が主に用いられる．
これを数値計算で解こう．
解くためには，ラプラシアン（ラプラス作用素，Laplace operator）
$${ \Delta = \nabla^{2} }$$ の離散化(Discretization)が必要だ．
Discretization
コンピュータで計算させるために，連続的なアナログ値を離散的なデジタル値にする，細かく測定点もしくは，予測点を置くようなものだ（標本化もしくは，samplingとも言われる）．まずはどのように点を置くかで，離散化の表現は変わってくる．
余談だが，人間の聴覚や視覚の解像度がそれぞれわかっているので，これの倍くらいで，等間隔に標本化しておけば，人間にはアナログな現実と，デジタル空間の差異はなくなる（視覚と聴覚に関してはね）．
もっと言えば，量子物理学によれば，現実の空間や時間はプランクスケールで離散化されている．
詳しくは，標本化定理（サンプリング定理，Nyquist-Shannon sampling theorem）あたりを参考にしてほしい．
今回行う，空間の離散化には，デカルト直交座標系で等間隔に標本化する．
簡単に，1次元格子もしくは，2次元の正方格子をここでは取ろう．
1次元格子
格子間隔$${h}$$を持つ，1次元格子において，離散化を行おう．
1次元の座標軸を$${x}$$軸と取れば，
原点$${x = 0}$$から$${i}$$番目の格子点の座標は$${x = ih}$$と表せる．
ここで，$${i}$$番目の格子点上に電荷もしくは質点があるとして考えてみよう．この条件は，$${f(ih) \neq 0}$$ のように書けるだろう．
簡単のために，関数の引数についてのみ，$${h = 1}$$のように$${h}$$を省略して表している．
数値計算では，オレンジのラインで描かれたような分布$${u(x)}$$を求めたい．これを実現するには，以下のように端点を除いた内部$${ i \in [1,n]}$$の格子点について離散化する：
$$
\frac{u(i - 1) - 2u(i) + u(i+1)}{h^{2}} + \mathscr{O}(h^4) = -f(i), i \in [1,n]
$$
$${ \mathscr{O}(h^4) }$$ は誤差のオーダーを表す．
（ランダウの記号，一番寄与の大きい部分を書くもの）
この関係は以下のように導出する：
１階微分は，前進差分（Forward Differentiation）
もしくは，後退差分（Backward Differentiation），
さらに，中央差分（Central Differentiation） のように離散化できる．
前進差分：
$${\frac{d}{dx} u(x) = \frac{u(i + 1) - u(i)}{h} + \mathscr{O}(h^2)}$$
後退差分：
$${\frac{d}{dx} u(x) = \frac{u(i) - u(i-1)}{h} + \mathscr{O}(h^2)}$$
中央差分：
$${\frac{d}{dx} u(x) = \frac{u(i+1) - u(i-1)}{2h} + \mathscr{O}(h^4)}$$
$${ \mathscr{O}(h^2), \mathscr{O}(h^4)}$$ はそれぞれ$${h^2, h^4}$$に比例する誤差が生じていることを意味する．
これを組み合わせれば，ラプラシアン（２階微分）の演算は，
$${\mathscr{O}(h^4)}$$の誤差を持つ形で，離散化されるだろう．
マクローリンもしくは，テイラー展開を知っているのであれば，証明は簡単だろう．
離散化した結果をもう一度書いておく：
$$
u(i - 1) - 2u(i) + u(i+1) + \mathscr{O}(h^4) = - f(i) h^{2}
$$
これをさらに変形して，
$${u(i)}$$がその周囲の$${u}$$とその点の$${f(i)}$$で決まるように表現しておこう：
$$
u(i) = \frac{u(i+1) + u(i - 1) + f(i) h^2}{2} + \mathscr{O}(h^4)
$$
$${f(i) = 0}$$なら，周囲$${u}$$の平均値によって$${u(i)}$$が得られるように見える．
しかし，$${u(i)}$$が決まれば，また，その周りの$${u(i+1)}$$なども変化させなくてはならない．先に触れておくと，ヤコビ法ではこの繰り返し，反復計算を収束するまで実行して，解を得る．
ちなみに，応用になるが，離散化したPoisson方程式を行列表示すると以下のようになる：
$$
A\bm{u} = \bm{f}
$$
行列$${A}$$は３重対角行列（tri-diagonal matrix）である．
(tri)diag(1, -2, 1)と書かれることもある．
2次元格子
1次元の場合と同じように，直交する二つの向きに関して等間隔$${h}$$で離散化し，2次元の正方格子を得る．
二次元の場合には，Poisson 方程式は以下のように離散化される：
$$
u(i+1, j) + u(i-1, j) + u(i, j+1) + u(i, j-1) - 4 u(i,j) + \mathscr{O}(h^4) = - f(i,j) h^2
$$
この場合も同様に，誤差は$${\mathscr{O}(h^4)}$$である．
これも同様に，$${u(i, j)}$$についての式に変形しておこう：
$$
u(i, j) = \frac{u(i+1, j) + u(i-1, j) + u(i, j+1) + u(i, j-1) + f(i, j) h^2}{4} + \mathscr{O}(h^4)
$$
やはり，2次元の場合も同様に，
隣接する格子点の$${u}$$の和と，その格子点の$${f(i,j)}$$の値から得られる．
$${f(i,j) = 0}$$なら，隣接する格子点の$${u}$$の平均値から得られる．
この計算は画像処理などで，ラプラシアンフィルタなどとして用いられたりもしている．この場合，格子点の近傍をどのように取るか（隣接する格子点はどれを取るか，どの格子点同士がつながっているとみるか）によって，フィルタとなる行列は異なる．これは，グラフ理論の一部として一般化されていて，ラプラシアン行列も合わせてみてほしい．
Jacobi method (a.k.a. Jacobi iteration method)
すでに少し触れたが，Jacobi 法では反復計算を使って，各格子点の$${u(i,j)}$$が収束する（ほとんど変化しなくなる）まで計算を行う．
反復する計算式を1次元と2次元のそれぞれについて，改めて見ておこう：
1次元
$$
u(i) = \frac{u(i+1) + u(i - 1) + f(i) h^2}{2} + \mathscr{O}(h^4)
$$
2次元
$$
u(i, j) = \frac{u(i+1, j) + u(i-1, j) + u(i, j+1) + u(i, j-1) + f(i, j) h^2}{4} + \mathscr{O}(h^4)
$$
3次元以降も同じような計算で行われると推論できるだろう．
反復計算を行うにあたって，計算を終了させる条件，収束条件が必要だ．
元の各格子点の$${u}$$について，
上記のラプラシアン行列の計算から得られた
新しい値を$${u_{\text{new}}(i)}$$とおいて，
次のような各格子点の誤差を表す式を使おう（例として１次元をとった）：
$$
e(i) = u_{\text{new}}(i) - u(i), ||e||_{\infty} = \text{max}|e(i)|
$$
各格子点のうち，最大の誤差を持つものを$${||e||_{\infty}}$$と呼んで定義している．これが，こちらで定めた，許容範囲（tolerance）以下になった際に計算を終了させる．
ヤコビ法のアルゴリズム
1次元の場合について，そのアルゴリズムを書いておこう：
// 収束するか，最大反復回数（iter_max）に至るまで計算
for iter in range(iter_max):
// 内部の各格子点について，$${u_{\text{new}}(i)}$$を計算
// 計算量 $${\mathscr{O}(n^{d})}$$
for i in range(1,n):
$${u_{new}(i) \leftarrow \frac{u(i+1) + u(i - 1) + f(i) h^2}{2}}$$
// 内部の各格子点について，誤差 $${e(i)}$$ を計算
// 計算量 $${\mathscr{O}(n^{d})}$$
for i in range(1,n):
$${e(i) = u_{\text{new}}(i) - u(i)}$$
// 誤差の最大値を計算
$${ ||e||_{\infty} = \text{max}|e(i)| }$$
// 誤差が許容範囲未満なら反復終了
if. $${ ||e||_{\infty} }$$ < tolerance:
break
全体の計算量は，Jacobi法がどの程度で収束するかにかかっている．
tolerance をどの程度に設定するかにもよるが，
ラプラシアンの計算は，１回の計算で，自身の持つ$${u(i)}$$の影響を隣接する格子点に伝播させるだけなので，端から端の一番遠い格子点同士の影響が届くまで，反復計算が必要だと考えられる．
つまり，反復計算は少なくとも，$${n^{d}}$$回は必要になるだろう．
計算量 $${\mathscr{O}(n^{d})}$$ のラプラシアンの計算を$${n^{d}}$$回必要になるので，
全体の計算量は，$${\mathscr{O}(n^{2d})}$$ 程度になるだろう．
行基本変形を用いた，ガウスの消去法（掃き出し法，Gauss elimination）は$${\mathscr{O}(n^{3d})}$$の計算回数が必要になるので，それよりは高速であることがわかる．
反復計算が収束する理由
なぜ，Jacobi 法などの反復計算が収束するのかについて，触れておこう．
例として，1次元の場合について考えながら，簡単に示しておく．
$$
A\bm{u} = \bm{f}
$$
ここで，行列$${A}$$は1次元の場合には，
以下の三重対角行列になる：
$$
A =
\begin{pmatrix}
-2 & 1 & & & &\\
1 & -2 & 1 & & & \\
& 1 & -2 & 1 & & \\
& & \ddots & & &\\
& & & 1 & -2 & 1 \\
& & & & 1 & -2
\end{pmatrix}
$$
$${ A = D + N }$$ のように対角行列$${D}$$と，非対角行列$${N}$$ に分解して，以下のように変形しよう：
$$
\bm{u} = - D^{-1} (N \bm{u} + \bm{f})
$$
ここで基本的な反復法では，これを以下のようにおいて，繰り返し計算を行う：
$$
\bm{u}^{(k)} = - D^{-1} (N \bm{u}^{(k-1)} + \bm{f})
$$
$${k \in \mathbb{N}}$$ は反復数のインデックスである．
先の等式は，1次元のJacobi法での反復計算では，以下の式のことである：
$$
u_{\text{new}}(i) = \frac{u(i+1) + u(i - 1) + f(i) h^2}{2} + \mathscr{O}(h^4)
$$
これを行列 $${ M = -D^{-1} N }$$を使って表現しておこう：
$$
\bm{u}_{\text{new}} = M \bm{u} - D^{-1} \bm{f} h^2 + \mathscr{O}(h^4)
$$
行列$${M}$$はある格子点について，周りの点の値を取ってきて，平均値を求める演算をしている．
1次元の場合は以下のような三重対角行列になる：
$$
M =
\frac{1}{2}
\begin{pmatrix}
0 & 1 & & & &\\
1 & 0 & 1 & & & \\
& 1 & 0 & 1 & & \\
& & \ddots & & &\\
& & & 1 & 0 & 1 \\
& & & & 1 & 0
\end{pmatrix}
$$
ここで，誤差ベクトル$${\bm{e}}$$について考えれば，
$$
\begin{equation*}
\begin{split}
\bm{e} &= \bm{u}_{\text{new}} - \bm{u} \\
&= M \bm{u} - D^{-1} \bm{f} h^2 + \mathscr{O}(h^4) - \bm{u} \\
&= (M - I) \bm{u} - D^{-1} \bm{f} h^2 + \mathscr{O}(h^4)
\end{split}
\end{equation*}
$$
ここで，$${I}$$は単位行列である．
$${k+1}$$回反復後の誤差$${\bm{e}^{(k+1)}}$$を考えよう．$${\bm{e}^{(k)}}$$や，それ以前の誤差との関係はどうなるだろうか？
$$
\begin{equation*}
\begin{split}
\bm{e}^{(k+1)} &= \bm{u}^{(k+1)} - \bm{u}^{(k)} \\
&=M \bm{u}^{(k)} -D^{-1} \bm{f} h^2 - (M \bm{u}^{(k-1)} -D^{-1} \bm{f} h^2) + \mathscr{O}(h^4) \\
&= M ( \bm{u}^{(k)} - \bm{u}^{(k-1)} ) + \mathscr{O}(h^4) \\
&= M \bm{e}^{(k)} + \mathscr{O}(h^4) \\
& \cdots \\
&= M^{k} \bm{e}^{(1)} + \mathscr{O}(h^4)
\end{split}
\end{equation*}
$$
十分大きな$${k}$$において，
$${||\bm{e}^{(k+1)}||_{2} < \|\bm{e}^{(0)}||_{2}}$$
になれば収束することがわかる．
ここで，これを不等式の関係を用いながら，変形していく：
$$
\begin{equation*}
\begin{split}
||\bm{e}^{(k+1)}||_{2}
&= || M \bm{e}^{(k)} + \mathscr{O}(h^4)||_{2} \\
&\leq || M \bm{e}^{(k)} ||_{2} \\
&= (M \bm{e}^{(k)})^T M \bm{e}^{(k)} \\
& \cdots \\
&\leq || M^{k} \bm{e}^{(1)} + \mathscr{O}(h^4)||_{2} \\
& \leq || M^{k} \bm{e}^{(1)} ||_{2} \\
&= (M^{k} \bm{e}^{(1)})^T M^{k} \bm{e}^{(1)} \\
\end{split}
\end{equation*}
$$
煩雑ですが，もうちょっと頑張りましょう．
もっと一般化できますが，今回は，ポアソン方程式を例に考えています．
そこで，行列$${M = -D^{-1} N}$$は対角成分がゼロとなっている，三重対角行列でであり，対称行列です．
実対称行列$${M}$$は対角化できる．
ここで，固有値$${\lambda_{p}}$$に対応する固有ベクトルを$${\bm{x}_{p}}$$として，
$$
M \bm{x}_{p} = \lambda_{p} \bm{x}_{p}
$$
さらに，対称行列の固有ベクトルは互いに直交する性質があります．
$$
(\bm{x}_{p})^{T} \bm{x}_{q} = \delta_{p,q}
$$
この固有ベクトル $${ \bm{x}_{p} }$$を用いて，誤差$${\bm{e}^{(k)}}$$ を展開しましょう：
$$
\bm{e}^{(k)} = \sum_{p = 1}^{n} a^{(k)}_{p} \bm{x}_{p}
$$
これより，
$$
M \bm{e}^{(k)} = \sum_{p = 1}^{n} a^{(k)}_{p} \lambda_{p} \bm{x}_{p}
$$
これを先の $${ \bm{e}^{(k + 1)} }$$ の関係式に代入すれば，
$$
\begin{equation*}
\begin{split}
||\bm{e}^{(k+1)}||_{2}
&\leq || M \bm{e}^{(k)} ||_{2} \\
&= (M \bm{e}^{(k)})^T M \bm{e}^{(k)} \\
&= \left(\sum_{p = 1}^{n} a^{(k)}_{p} \lambda_{p} \bm{x}_{p}\right)^{T}
\left(\sum_{q = 1}^{n} a^{(k)}_{q} \lambda_{q} \bm{x}_{q}\right)
\\
&= \sum_{p,q = 1}^{n} a^{(k)}_{p} a^{(k)}_{q} \lambda_{p} \lambda_{q} \delta_{p, q}
\\
&= \sum_{p = 1}^{n} (a^{(k)}_{p} \lambda_{p} )^2
\end{split}
\end{equation*}
$$
さらに，反復を考えれば，
$${ M^{k} \bm{x}_{p} = (\lambda_{p})^{k} \bm{x}_{p} }$$
になるから，容易に以下の関係が得られる：
$$
\begin{equation*}
\begin{split}
||\bm{e}^{(k+1)}||_{2}
&\leq \sum_{p = 1}^{n} (a^{(k)}_{p} \lambda_{p} )^2 \\
& \cdots \\
&\leq \sum_{p = 1}^{n} \left( a^{(1)}_{p} (\lambda_{p})^{k} \right)^2
\end{split}
\end{equation*}
$$
ここで，行列
$$
M =
\frac{1}{2}
\begin{pmatrix}
0 & 1 & & & &\\
1 & 0 & 1 & & & \\
& 1 & 0 & 1 & & \\
& & \ddots & & &\\
& & & 1 & 0 & 1 \\
& & & & 1 & 0
\end{pmatrix}
$$
の具体的な固有値を書いておきます：
$$
\lambda_{p} = \cos{ \left(\frac{p}{n+1} \pi \right) }, p \in [1, n]
$$
この導出については，例えば，以下のサイトを参考にしてください：
重要な条件は，$${|\lambda_{p}| < 1}$$ であることです．
この条件から以下のような計算を進めます．
まず，誤差ノルム$${||\bm{e}^{(k+1)}||_{2}}$$については，
$$
||\bm{e}^{(k+1)}||_{2}
\leq \sum_{p = 1}^{n} \left( a^{(1)}_{p} (\lambda_{p})^{k} \right)^2
$$
が成り立っています．
さらに，もう１回反復した後の誤差ノルム
$${||\bm{e}^{(k+1)}||_{2}}$$を考えれば，
$$
\begin{equation*}
\begin{split}
||\bm{e}^{(k+2)}||_{2}
&\leq
\sum_{p = 1}^{n} \left( a^{(1)}_{p} (\lambda_{p})^{k+1} \right)^2 \\
&=
\sum_{p = 1}^{n} \left( a^{(1)}_{p} (\lambda_{p})^{k} \right)^2 (\lambda_{p})^2 \\
&< ||\bm{e}^{(k+1)}||_{2} \leq\sum_{p = 1}^{n} \left( a^{(1)}_{p} (\lambda_{p})^{k} \right)^2
\end{split}
\end{equation*}
$$
３行目の不等号 $${<}$$ で先の重要な関係 $${|\lambda_{p}| < 1}$$ を用いた．
余談だが，一般の証明では，行列$${M}$$のスペクトル半径（固有値の中で絶対値が最大となるもの）が１未満であれば収束すると証明されている．
これより，誤差ノルムは反復する度に単調に減少していくことが示され，反復計算は収束することがわかる．
以上が，Jacobi法による数値計算方法である．
何か参考になったなら嬉しく思います．
Appendix: 格子点描画のスクリプト
格子を描くのに用いた python スクリプトをおいておこうと思う．
よく使うものだが，毎回，LLMで作成するのはちょっと面倒（プロンプトを書くのが）なので，ここにおいておく．（これだけのためにGIthubレポジトリ置くのも面倒というw）
import numpy as np
import matplotlib.pyplot as plt
def draw_1d_grid(n):
"""
1D の格子点（2^n 個）を等間隔で描画する。
"""
num_points = 2 ** n
x = np.linspace(0, 1, num_points)
y = np.zeros_like(x)
plt.figure(figsize=(8, 2))
plt.plot(x, y, 'ko-', markersize=4)
plt.axis('off')
plt.show()
def draw_2d_grid(n):
"""
2D の正方格子（各方向に 2^n 個の点、全体で 2^n × 2^n 個の点）を描画する。
"""
num_points = 2 ** n
x = np.linspace(0, 1, num_points)
y = np.linspace(0, 1, num_points)
X, Y = np.meshgrid(x, y)
plt.figure(figsize=(6, 6))
for i in range(num_points):
plt.plot([x[i]] * num_points, y, 'k-', linewidth=0.5) # 縦線
plt.plot(x, [y[i]] * num_points, 'k-', linewidth=0.5) # 横線
plt.scatter(X, Y, color='k', s=10)
plt.axis('off')
plt.show()
if __name__ == "__main__":
n = 3
draw_1d_grid(n)
draw_2d_grid(n)