Pythonで学ぶ数値計算のアルゴリズムと実装第13回：空間の熱と振動を解く：偏微分方程式の差分近似と反復解法

こんにちは、JS2IIUです。
全13回にわたる「Pythonで学ぶ数値計算のアルゴリズムと実装」の連載も、ついに最終回を迎えました。第1回で多項式計算の基礎から始まり、積分、線形代数、非線形方程式、そして前回の常微分方程式（ODE）へと進んできた私たちの旅は、最も応用範囲が広く、かつ計算負荷の高い「偏微分方程式（PDE）」の領域へと到達します。

常微分方程式が「時間とともに変化する一つの点」の挙動を追うものであったのに対し、偏微分方程式は「空間的な広がり（時間＋空間、あるいは多次元空間）」における現象を記述します。例えば、部屋の中の熱の伝わり方、波の伝播、流体の流れ、さらにはAIにおける拡散モデル（Diffusion Models）の数学的背景に至るまで、偏微分方程式は現代科学と工学のあらゆる場面で中心的な役割を果たしています。

今回は、偏微分方程式の中でも最も基本的かつ重要な「2次元ラプラス方程式」を題材に、それをコンピュータで解くための「有限差分近似」と、大規模な方程式を効率的に解くための「線形反復法（ガウス・サイデル法、SOR法）」について深く解説します。今回もよろしくお願いします。

1. 2次元ラプラス方程式と有限差分近似
1. コンピュータによる解法：格子分割と差分近似
2. 線形反復法の必要性
3. ガウス・サイデル法（Gauss-Seidel Method）
1. アルゴリズムの流れ
2. 物理的イメージとメリット
4. 加速の技術：SOR法（Successive Over-Relaxation）
1. SOR法のメリットと実用性
5. PyTorchによる実践：2次元ポテンシャルの計算
1. コードの解説
6. 線形反復法の収束定理
7. 連載の終わりに：数値計算の地平線
連載記事、各回へのリンク

1. 2次元ラプラス方程式と有限差分近似

今回扱う「2次元ラプラス方程式」は、物理や工学の多くの分野で現れる基本的な偏微分方程式です。式は次のようになります。

$$\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} = 0$$

この方程式は、例えば熱が外部から加えられず、時間的に安定した状態（定常状態）の温度分布や、静電場のポテンシャル分布などを記述します。意味としては「ある点の周囲の曲がり具合（2階微分）の合計がゼロ」、つまりその点は周囲の値の平均になっている、という物理的な調和状態を表します。

コンピュータによる解法：格子分割と差分近似

連続的な空間をコンピュータで扱うためには、空間を小さな格子（グリッド）に分割します。各格子点の間隔を $\Delta x, \Delta y$ とすると、連続的な微分を離散的な差分で近似できます。例えば、$x$方向の2階微分は次のように表せます。

$$\frac{\partial^2 u}{\partial x^2} \approx \frac{u_{i+1,j} – 2u_{i,j} + u_{i-1,j}}{(\Delta x)^2}$$

同様に$y$方向も差分で近似できます。これらをラプラス方程式に代入し、格子間隔が等しい（$\Delta x = \Delta y$）と仮定すると、各格子点の値は上下左右4点の平均値になるという非常にシンプルな関係式が得られます。

$$u_{i,j} = \frac{1}{4}(u_{i+1,j} + u_{i-1,j} + u_{i,j+1} + u_{i,j-1})$$

この式は「格子点 $(i, j)$ の値は、隣接する4点の値の平均になる」という物理的な意味を持ちます。これは、熱や電位が空間的に滑らかに分布し、急激な変化がない安定状態を反映しています。

ラプラス方程式は空間の調和状態を記述する。
コンピュータでは格子分割と差分近似で離散化する。
各格子点の値は、隣接点の平均値になるというシンプルな関係式に帰着する。

この考え方が、以降の反復法や数値解法の基礎となります。

2. 線形反復法の必要性

ラプラス方程式を格子状に離散化すると、各格子点の値が未知数となり、$N \times N$ の格子なら $N^2$ 個の未知数を持つ連立方程式が生じます。例えば $100 \times 100$ の格子では、10,000 個もの未知数を一度に解く必要があります。

このような大規模な連立方程式を、ガウス消去法やLU分解などの「直接法」で解こうとすると、計算量は $O((N^2)^3) = O(N^6)$ となり、格子を細かくするほど計算時間やメモリ消費が急激に増大します。現実的なサイズの問題では、直接法は計算資源の制約から実用的ではありません。

一方、ラプラス方程式の差分化によって得られる連立方程式の係数行列は、ほとんどの要素がゼロである「スパース（疎）行列」となります。例えば、各格子点は上下左右の4点としか関係しないため、1行あたりの非ゼロ要素はごくわずかです。この疎な構造を活かすことで、メモリ効率よく計算を進めることができます。

そこで登場するのが「線形反復法」です。反復法では、まず全ての未知数に適当な初期値（例えばゼロや境界条件に基づく値）を与えます。その後、各格子点の値を周囲の点の値を使って順次更新し、全体の解が徐々に正しい値へと収束するまでこの操作を繰り返します。

具体的な流れは以下の通りです：

初期値を設定する（例：内部はゼロ、境界は物理条件に従う）。
各格子点について、周囲の値を使って新しい値を計算する。
全体の変化量が十分小さくなるまで、2の操作を繰り返す。

この方法により、膨大な計算量を必要とする直接法を使わず、疎な構造を活かして効率的に大規模な問題を解くことが可能となります。線形反復法は、現代の数値計算やシミュレーションの現場で不可欠な手法です。

3. ガウス・サイデル法（Gauss-Seidel Method）

ガウス・サイデル法は、ラプラス方程式の差分近似で得られる連立方程式を効率的に解くための代表的な線形反復法です。基本的なアイデアは「各格子点の値を、周囲の最新の値の平均で順次更新していく」ことです。

この手法は、ヤコビ法（Jacobi法）とよく比較されます。ヤコビ法では、全ての格子点の値を一度に更新するため、各ステップで使う値は前回の値のみです。一方、ガウス・サイデル法では、格子をスキャンする際に、すでに更新した新しい値をすぐに次の計算に利用します。これにより、収束が大幅に速くなり、メモリ効率も向上します。

アルゴリズムの流れ

格子全体に初期値を設定する（内部はゼロ、境界は物理条件など）。
格子点 $(i, j)$ を順番に走査し、以下の式で新しい値を計算する：
$$u_{i,j}^{new} = \frac{1}{4}(u_{i+1,j}^{old/new} + u_{i-1,j}^{old/new} + u_{i,j+1}^{old/new} + u_{i,j-1}^{old/new})$$
※ すでに更新済みの点は新しい値、未更新の点は前回の値を使う。
全体の変化量が十分小さくなるまで、2の操作を繰り返す。

物理的イメージとメリット

ガウス・サイデル法は、格子点の値が「周囲の最新の値の平均」にどんどん近づいていくことで、空間全体が調和状態（ラプラス方程式の解）に収束します。ヤコビ法よりも収束が速く、計算資源の節約にもつながるため、実用的な数値計算で広く使われています。

また、格子の走査順序や並列化の工夫（レッド・ブラック法など）によって、さらに効率を高めることも可能です。

4. 加速の技術：SOR法（Successive Over-Relaxation）

ガウス・サイデル法は収束が速い反復法ですが、格子サイズが大きくなると収束までに非常に多くの反復回数が必要になることがあります。そこで、さらに収束を加速するために考案されたのが「逐次加速緩和法（SOR法：Successive Over-Relaxation）」です。

SOR法では、ガウス・サイデル法で得られる新しい値（平均値）に対して、現在の値との差分を強調して更新します。具体的には、次のような更新式を用います：

$$u_{i,j}^{new} = (1 – \omega)u_{i,j}^{old} + \omega \bar{u}_{GS}$$

ここで $u_{i,j}^{old}$ は前回の値、$\bar{u}_{GS}$ はガウス・サイデル法で計算した周囲の平均値、$\omega$ は「加速パラメータ（緩和係数）」です。

このパラメータ $\omega$ の選び方によって、収束速度が大きく変わります。

$\omega = 1$ の場合は通常のガウス・サイデル法と同じ。
$1 < \omega < 2$ の範囲で選ぶと、解の更新量が増幅され、収束が加速されます（Over-Relaxation）。
$\omega$ が大きすぎると逆に発散することもあるため、最適値の探索が重要です。

SOR法のメリットと実用性

SOR法を用いることで、同じ精度の解を得るまでの反復回数を大幅に減らすことができます。特に大規模な格子問題では、計算時間が数分の一から数十分の一に短縮されることもあり、実用的な数値シミュレーションにおいて「魔法の加速法」と呼ばれるほど重宝されています。

物理的なイメージとしては、「解の変化を積極的に反映させることで、空間全体の調和状態への到達を早める」手法です。

SOR法は、ガウス・サイデル法の拡張として非常に広く使われており、有限要素法や流体計算など多くの分野で標準的な加速技術となっています。

5. PyTorchによる実践：2次元ポテンシャルの計算

それでは、正方形の領域の境界に異なる温度（電位）を設定し、内部がどのように平衡状態に達するかを PyTorch でシミュレーションしてみましょう。

Python

import torch
import matplotlib.pyplot as plt

def solve_laplace_sor(n_grid, omega, tol=1e-5, max_iter=5000):
    """
    SOR法を用いて2次元ラプラス方程式を解く
    """
    # グリッドの初期化 (内部は0, 境界条件を設定)
    u = torch.zeros((n_grid, n_grid))

    # 境界条件 (Dirichlet境界条件)
    u[0, :] = 100.0   # 上端を100度
    u[-1, :] = 0.0    # 下端を0度
    u[:, 0] = 50.0    # 左端を50度
    u[:, -1] = 50.0   # 右端を50度

    for it in range(max_iter):
        u_old = u.clone()

        # ガウス・サイデル的な更新 (ベクトル化は難しいのでループ、
        # またはチェッカーボード法が使われるがここでは基本に忠実に)
        # 境界を除く内部点を更新
        for i in range(1, n_grid - 1):
            for j in range(1, n_grid - 1):
                # 4点平均
                avg = 0.25 * (u[i+1, j] + u[i-1, j] + u[i, j+1] + u[i, j-1])
                # SORの更新式
                u[i, j] = (1 - omega) * u[i, j] + omega * avg

        # 収束判定 (最大変化量が許容値以下か)
        diff = torch.max(torch.abs(u - u_old))
        if diff < tol:
            print(f"SOR converged at iteration {it} (omega={omega})")
            return u, it

    return u, max_iter

# --- パラメータの設定と比較 ---
n = 40 # 40x40の格子
u_gs, it_gs = solve_laplace_sor(n, omega=1.0)  # ガウス・サイデル
u_sor, it_sor = solve_laplace_sor(n, omega=1.8) # SOR法 (加速あり)

print(f"Gauss-Seidel iterations: {it_gs}")
print(f"SOR iterations:          {it_sor}")

# 結果の可視化
plt.figure(figsize=(8, 6))
plt.imshow(u_sor.numpy(), cmap='inferno', origin='lower')
plt.colorbar(label='Potential / Temperature')
plt.title("Steady State distribution (Laplace Equation)")
plt.show()

import torch
import matplotlib.pyplot as plt

def solve_laplace_sor(n_grid, omega, tol=1e-5, max_iter=5000):
    """
    SOR法を用いて2次元ラプラス方程式を解く
    """
    # グリッドの初期化 (内部は0, 境界条件を設定)
    u = torch.zeros((n_grid, n_grid))

    # 境界条件 (Dirichlet境界条件)
    u[0, :] = 100.0   # 上端を100度
    u[-1, :] = 0.0    # 下端を0度
    u[:, 0] = 50.0    # 左端を50度
    u[:, -1] = 50.0   # 右端を50度

    for it in range(max_iter):
        u_old = u.clone()

        # ガウス・サイデル的な更新 (ベクトル化は難しいのでループ、
        # またはチェッカーボード法が使われるがここでは基本に忠実に)
        # 境界を除く内部点を更新
        for i in range(1, n_grid - 1):
            for j in range(1, n_grid - 1):
                # 4点平均
                avg = 0.25 * (u[i+1, j] + u[i-1, j] + u[i, j+1] + u[i, j-1])
                # SORの更新式
                u[i, j] = (1 - omega) * u[i, j] + omega * avg

        # 収束判定 (最大変化量が許容値以下か)
        diff = torch.max(torch.abs(u - u_old))
        if diff < tol:
            print(f"SOR converged at iteration {it} (omega={omega})")
            return u, it

    return u, max_iter

# --- パラメータの設定と比較 ---
n = 40 # 40x40の格子
u_gs, it_gs = solve_laplace_sor(n, omega=1.0)  # ガウス・サイデル
u_sor, it_sor = solve_laplace_sor(n, omega=1.8) # SOR法 (加速あり)

print(f"Gauss-Seidel iterations: {it_gs}")
print(f"SOR iterations:          {it_sor}")

# 結果の可視化
plt.figure(figsize=(8, 6))
plt.imshow(u_sor.numpy(), cmap='inferno', origin='lower')
plt.colorbar(label='Potential / Temperature')
plt.title("Steady State distribution (Laplace Equation)")
plt.show()

Plaintext

SOR converged at iteration 1803 (omega=1.0)
Gauss-Seidel iterations: 1803
SOR iterations:          5000

SOR converged at iteration 1803 (omega=1.0)
Gauss-Seidel iterations: 1803
SOR iterations:          5000

コードの解説

境界条件: u[0, :] = 100.0 などで領域の縁の値を固定しています。これは「ディリクレ境界条件」と呼ばれます。
二重ループ: 反復法では、新しい値をすぐに次の計算に使う性質上、単純な行列演算（ベクトル化）が難しいという側面があります。GPUを活用する場合は、格子を市松模様に塗り分けて並列計算する「レッド・ブラック法」などが使われます。
加速の効果: omega=1.0 と 1.8 の収束回数を比較すると、SOR法の圧倒的な速さが確認できるはずです。

6. 線形反復法の収束定理

線形反復法（ガウス・サイデル法やSOR法など）は、必ずしも全ての問題で収束するわけではありません。収束するためには、数学的にいくつかの条件が必要です。

収束の数学的条件

反復法の収束性は、反復操作を行う「反復行列」のスペクトル半径（最大固有値の絶対値）$\rho(G)$ が 1 未満であることが必要です。

$$
\rho(G) < 1
$$

この条件が満たされていれば、反復を繰り返すごとに誤差が指数的に減少し、最終的に真の解に収束します。

実用的な収束判定：対角優位

実際の数値計算では、係数行列 $A$ が「対角優位（diagonally dominant）」であるかどうかが重要な収束判定基準となります。

行列 $A$ の各行について、対角成分 $|a_{ii}|$ が他の成分の絶対値の和より大きい場合、

$$
|a_{ii}| > \sum_{j \neq i} |a_{ij}|
$$

この条件が全ての行で満たされていれば、ガウス・サイデル法は必ず収束することが数学的に保証されています。

ラプラス方程式の安定性

ラプラス方程式の差分化によって得られる係数行列は、内部点では必ず対角優位となる構造を持っています（境界条件も含めて）。そのため、ガウス・サイデル法やSOR法は非常に安定して収束しやすいという特徴があります。

反復法の収束には「スペクトル半径 $\rho(G) < 1$」が必要。
実用的には「対角優位」であれば収束が保証される。
ラプラス方程式の差分問題はこの条件を満たしやすく、安定して解ける。

7. 連載の終わりに：数値計算の地平線

全13回にわたるこの連載を通して、私たちはコンピュータがどのようにして「数」を扱い、「数式」を「現実的な解」へと変えていくのかを見てきました。

近似の精神: テイラー展開や数値積分で学んだように、複雑なものは単純なものの積み重ねで近似できます。
誤差への敬意: ノルム、条件数、丸め誤差。計算の「正しさ」を疑う視点は、エンジニアにとっての良心です。
効率の追求: LU分解、帯行列、そして今回のSOR法。限られたリソースでいかに速く答えを出すか、その工夫こそがアルゴリズムの真髄です。
現象の再現: 微分方程式を解く力は、AIが現実世界の物理法則を理解し、共創するための基礎となります。

現代のAI技術は目覚ましい進化を遂げていますが、その計算の底流には、今回私たちが学んだような100年以上前から磨き上げられてきた数値計算の知恵が脈々と流れています。Pythonという強力な道具を使いこなし、ライブラリの裏側にある「原理」を知った今のあなたなら、ブラックボックスに頼るだけでなく、真に効率的で信頼性の高いシステムを構築できるはずです。

数値計算の世界は、ここで終わりではありません。有限要素法（FEM）、モンテカルロ法、最適化理論など、まだまだ深い森が広がっています。この連載が、皆さんがその広大な知の領域を探索していくための一助となれば幸いです。

最後までお読みいただき、本当にありがとうございました。

連載記事、各回へのリンク

Pythonで学ぶ数値計算のアルゴリズムと実装第1回：多項式計算とテイラー展開 | アマチュア無線局JS2IIU
https://js2iiu.com/2025/12/31/python_numerical_01/

Pythonで学ぶ数値計算のアルゴリズムと実装第2回：多項式補間とチェビシェフ補間の戦略 | アマチュア無線局JS2IIU
https://js2iiu.com/2026/01/01/python_numerical_02/

Pythonで学ぶ数値計算のアルゴリズムと実装第3回：積分の基礎：中点則、台形則からシンプソン則まで | アマチュア無線局JS2IIU
https://js2iiu.com/2026/01/01/python_numerical_03/

Pythonで学ぶ数値計算のアルゴリズムと実装第4回：高度な数値積分：二重指数型（DE）公式と広義積分 | アマチュア無線局JS2IIU
https://js2iiu.com/2026/01/01/python_numerical_04/

Pythonで学ぶ数値計算のアルゴリズムと実装第5回：計算の「正しさ」を評価する：ノルムと条件数 | アマチュア無線局JS2IIU
https://js2iiu.com/2026/01/01/python_numeric_05/

Pythonで学ぶ数値計算のアルゴリズムと実装第6回：連立一次方程式の直截解法：ガウス消去法とピボット選択 | アマチュア無線局JS2IIU
https://js2iiu.com/2026/01/01/python_numerical_06/

Pythonで学ぶ数値計算のアルゴリズムと実装第8回：スパース性を活かす：帯行列の高速計算アルゴリズム | アマチュア無線局JS2IIU
https://js2iiu.com/2026/01/02/python_numerical_08/

Pythonで学ぶ数値計算のアルゴリズムと実装第9回：最小二乗法とQR分解：ハウスホルダー変換による安定化 | アマチュア無線局JS2IIU
https://js2iiu.com/2026/01/02/python_numerial_09/

Pythonで学ぶ数値計算のアルゴリズムと実装第10回：非線形方程式の探索：2分法からニュートン法まで | アマチュア無線局JS2IIU
https://js2iiu.com/2026/01/03/python_numerical_10/

Pythonで学ぶ数値計算のアルゴリズムと実装第11回：行列の本質を捉える：固有値問題の数値解法 | アマチュア無線局JS2IIU
https://js2iiu.com/2026/01/03/python_numerical_11/

Pythonで学ぶ数値計算のアルゴリズムと実装第12回：現象の変化を追う：ルンゲ・クッタ法による常微分方程式の解法 | アマチュア無線局JS2IIU
https://js2iiu.com/2026/01/03/python_numerical_12/

Pythonで学ぶ数値計算のアルゴリズムと実装第13回：空間の熱と振動を解く：偏微分方程式の差分近似と反復解法 | アマチュア無線局JS2IIU
https://js2iiu.com/2026/01/03/python_numerical_13/

月	火	水	木	金	土	日
			1	2	3	4
5	6	7	8	9	10	11
12	13	14	15	16	17	18
19	20	21	22	23	24	25
26	27	28	29	30	31