Pythonで攻略する流体力学の数値計算 第9回：常微分方程式の解法(2) ― 境界値問題と差分方程式

Python
この記事は約18分で読めます。

こんにちはJS2IIUです。Pythonで攻略する流体力学シミュレーション連載の第9回です。

前回は、時間とともに変化する現象を追いかける「初期値問題」を扱いました。オイラー法やルンゲ・クッタ法を用いて、「現在の状態」から「一歩先の未来」を予測する手法は、流体の動的な挙動をシミュレーションする際の核となる技術でした。

しかし、物理現象にはもう一つの側面があります。それは、時間が十分に経過して変化が落ち着いた「定常状態」や、空間全体のバランスが同時に決定されるようなケースです。このような問題を扱うのが、今回解説する境界値問題（Boundary Value Problem, BVP）です。

境界値問題では、時間発展のように一歩ずつ進めるのではなく、領域全体の値を「行列計算」によって一度に決定します。この考え方は、流体計算における圧力の決定や、熱伝導の定常解を求める際に不可欠な知識となります。今回もよろしくお願いします。

1. 境界値問題とは何か

初期値問題が「ある一点（初期時刻）から出発して未来へと進む」タイプの問題であるのに対し、境界値問題（Boundary Value Problem, BVP）は「領域の両端（境界）での条件があらかじめ決められており、その間の分布や状態を求める」タイプの問題です。

初期値問題との違い

初期値問題では、たとえば「時刻\(t_0\)での状態\(y(t_0)\)」が与えられ、そこから時間を進めて\(y(t)\)を計算していきます。これは「スタート地点からゴールに向かって一歩ずつ進む」イメージです。

一方、境界値問題では「空間領域の両端（\(x=0\)と\(x=L\)など）での値」が与えられ、その間の分布\(u(x)\)を一度に決定します。これは「スタートとゴールの両方が決まっていて、その間をどのように埋めるかを考える」イメージです。

物理的な例：1次元の棒の定常熱伝導

たとえば、1本の細長い棒の温度分布を考えます。

  • 左端（\(x=0\)）の温度が100度に固定されている。
  • 右端（\(x=L\)）の温度が0度に固定されている。
  • 棒の内部には熱源が分布しており、場所によって温度勾配が生じている。

このとき、棒の内部（\(0 < x < L\)）の温度分布\(u(x)\)がどのようになるかを求めるのが境界値問題です。

数学的な定式化

境界値問題は、一般に次のような微分方程式と境界条件の組み合わせで表されます。

$$
\frac{d^2 u}{dx^2} = f(x), \quad u(0) = \alpha, \quad u(L) = \beta
$$

ここで、

  • \(u(x)\)：未知の関数（例：温度分布）
  • \(f(x)\)：空間的な外力や熱源分布などの既知関数
  • \(u(0)=\alpha,\ u(L)=\beta\)：両端の境界条件（ディリクレ境界条件）

このような問題では、領域全体のバランスを一度に考える必要があり、初期値問題のように「一歩ずつ進める」ことはできません。

境界値問題の物理的・工学的な重要性

境界値問題は、熱伝導、弾性体の変形、静電場、流体の圧力分布など、空間全体の状態が境界条件によって決まる現象で広く現れます。

たとえば、

  • 棒や板の温度分布（定常熱伝導）
  • ビームや柱のたわみ（構造力学）
  • 電位分布（静電場）
  • 流体の圧力分布（非圧縮性流体のポアソン方程式）

などが典型例です。

流体力学における境界値問題

流体力学では、たとえば非圧縮性条件（流体の密度が変化しない条件）を満たすために、圧力のポアソン方程式という境界値問題を解く必要があります。これは、流体の速度場が「質量保存」を満たすように圧力分布を調整するための重要なステップです。

このように、境界値問題は「空間全体のバランスを一度に決定する」ための基本的な枠組みであり、数値流体力学や工学シミュレーションの根幹をなしています。

2. 差分方程式への変換

境界値問題をコンピュータで数値的に解くためには、連続的な微分方程式を離散化し、計算機で扱える差分方程式へと変換する必要があります。このとき活躍するのが、第6回で学んだ差分法（finite difference method, FDM）です。

差分法の基本的な考え方

連続空間をそのまま扱うことはできないため、領域全体を\(N\)個の格子点（グリッドポイント）に分割し、各格子点での値$u_i$を求める形にします。格子点間の間隔を$\Delta x$とします。

たとえば、\(x_i = i \Delta x\)（\(i=0,1,\ldots,N-1\)）のように等間隔で格子を並べます。

二階微分の差分近似

境界値問題でよく現れる二階微分$\frac{d^2u}{dx^2}$は、中心差分という手法で次のように近似できます：

$$
\frac{d^2u}{dx^2}\Big|_{x_i} \approx \frac{u_{i+1} – 2u_i + u_{i-1}}{\Delta x^2}
$$

この式は、「真ん中の点\(i\)の二階微分は、隣り合う点\(i-1\)と\(i+1\)の値を使って近似できる」という意味です。これは物理的にも「隣接点とのバランスが中央点の曲がり具合を決める」ことを表しています。

差分方程式への翻訳

もとの微分方程式

$$
\frac{d^2u}{dx^2} = f(x)
$$

を、各格子点\(x_i\)で差分近似すると、

$$
\frac{u_{i+1} – 2u_i + u_{i-1}}{\Delta x^2} = f_i
$$

となります。ここで\(f_i = f(x_i)\)です。

この式を整理すると、

$$
u_{i-1} – 2u_i + u_{i+1} = \Delta x^2 f_i
$$

という代数方程式になります。

境界条件の扱いと未知数の決定

この差分方程式は、内部のすべての格子点（\(i=1,2,\ldots,N-2\)）について成り立ちます。一方、両端の点（\(i=0\)と\(i=N-1\)）は境界条件（\(u_0=\alpha,\ u_{N-1}=\beta\)など）として既知です。

したがって、実際に未知数となるのは内部点\(u_1,\ldots,u_{N-2}\)だけです。これらについて、\(N-2\)本の連立一次方程式が得られます。

連立方程式化の意義

このようにして、もともと連続的だった微分方程式の問題が、

  • 「隣接点との関係式（差分方程式）」
  • 「境界条件」

という形で有限個の未知数に対する連立一次方程式に変換されます。これが差分法の本質です。

この連立方程式を解くことで、領域全体の分布\(u(x)\)を近似的に求めることができます。差分法は、シンプルかつ汎用性が高く、流体力学・熱伝導・構造力学など幅広い分野で使われています。

領域を \(N\) 個の格子点に分割し、格子間隔を \(\Delta x\) とします。地点 \(x_i\) における二階微分を二次の中心差分で近似すると、以下のような代数方程式が得られます。

$$ \frac{u_{i+1} – 2u_i + u_{i-1}}{\Delta x^2} = f_i $$

これを整理すると、隣接する格子点との関係式になります。

$$ u_{i-1} – 2u_i + u_{i+1} = \Delta x^2 f_i $$

この式は、内部のすべての格子点（\(i=1, 2, \dots, N-2\)）について成り立ちます。境界点（\(i=0\) と \(i=N-1\)）の値は境界条件として既知であるため、未知数 \(u_1, \dots, u_{N-2}\) に関する連立一次方程式が得られることになります。

3. 行列計算による一括解決

差分方程式によって内部格子点ごとに得られた\(N-2\)本の一次方程式は、行列形式でまとめて表現できます。これにより、領域全体の未知数を一度に効率よく計算できるようになります。

行列形式へのまとめ

各内部点\(i\)について、

$$
u_{i-1} – 2u_i + u_{i+1} = \Delta x^2 f_i
$$

という式が成り立ちます。これをすべての内部点について並べると、

$$
A\mathbf{u} = \mathbf{b}
$$

という連立一次方程式の形になります。

ここで、

  • \(A\)：係数行列（\(N-2\)次元の正方行列）
  • \(\mathbf{u}\)：内部点の未知数ベクトル（\(u_1, u_2, \ldots, u_{N-2}\)）
  • \(\mathbf{b}\)：右辺ベクトル（\(\Delta x^2 f_i\)と境界値の寄与を含む）

三対角行列の構造

例えば、内部格子点が3点（全体で5点）ある場合、行列は次のようになります：

$$
\begin{bmatrix}
-2 & 1 & 0 \\
1 & -2 & 1 \\
0 & 1 & -2
\end{bmatrix}
\begin{bmatrix}
u_1 \\
u_2 \\
u_3
\end{bmatrix}
=
\begin{bmatrix}
\Delta x^2 f_1 – u_0 \\
\Delta x^2 f_2 \\
\Delta x^2 f_3 – u_4
\end{bmatrix}
$$

この\(A\)は、三対角行列（tridiagonal matrix）と呼ばれます。すなわち、主対角成分（\(-2\)）と、そのすぐ上・下の成分（\(1\)）だけが非ゼロで、それ以外はすべてゼロです。

右辺ベクトル\(\mathbf{b}\)の意味

右辺\(\mathbf{b}\)は、各点での\(\Delta x^2 f_i\)に加え、両端の境界値\(u_0\)や\(u_{N-1}\)の影響が含まれます。たとえば最初の式には\(-u_0\)、最後の式には\(-u_{N-1}\)が現れます。これにより、境界条件が内部点の解に正しく反映されます。

行列計算の物理的・計算的な意義

このように行列形式にまとめることで、

  • 領域全体のバランスを一度に計算できる
  • 数値的に安定で効率的なアルゴリズム（ガウス消去法、トーマス法など）が使える
  • 大規模な問題にも拡張しやすい
    といった利点があります。

特に三対角行列は、計算量が\(O(N)\)で済む専用アルゴリズム（トーマス法）が存在し、数千～数百万次元の問題でも高速に解くことができます。

また、現代のPythonや数値計算ライブラリ（NumPy, SciPyなど）では、こうした行列計算が非常に効率よく実装されており、複雑な物理現象のシミュレーションにも応用できます。

大規模・高次元への拡張性

1次元の場合は三対角行列ですが、2次元・3次元の境界値問題でも同様に「疎な（ほとんどがゼロの）行列」が現れます。これらも疎行列用のアルゴリズムやライブラリを使うことで、現実的な計算時間で解くことが可能です。

このように、行列計算による一括解決は、数値シミュレーションの根幹を支える強力な道具となっています。

4. Pythonによる実装：定常熱伝導問題

それでは、PythonとNumPyを用いて、1次元の境界値問題を解いてみましょう。
今回は、\(\frac{d^2 u}{dx^2} = \sin(x)\) という熱源分布を持つ領域において、両端の温度が固定された状態の定常解を求めます。

Python
import numpy as np
import matplotlib.pyplot as plt

def solve_bvp_1d(n_points, x_max, bc_left, bc_right, source_func):
    """
    1次元境界値問題 (d^2u/dx^2 = f(x)) を行列計算で解く。

    Args:
        n_points (int): 格子点数 (境界を含む)
        x_max (float): 領域の長さ
        bc_left (float): 左端の境界値 u(0)
        bc_right (float): 右端の境界値 u(L)
        source_func (callable): 右辺の関数 f(x)

    Returns:
        x (ndarray): 座標
        u (ndarray): 計算された解
    """
    x = np.linspace(0, x_max, n_points)
    dx = x[1] - x[0]

    # 未知数の数 (境界を除いた内部点)
    n_internal = n_points - 2

    # 係数行列 A の作成 (三対角行列)
    # A * u = b の形式にする
    # 各行は [..., 1, -2, 1, ...] という構造を持つ
    A = np.zeros((n_internal, n_internal))
    for i in range(n_internal):
        A[i, i] = -2.0
        if i > 0:
            A[i, i-1] = 1.0
        if i < n_internal - 1:
            A[i, i+1] = 1.0

    # 右辺ベクトル b の作成
    # dx^2 * f(x) の計算
    b = (dx**2) * source_func(x[1:-1])

    # 境界条件の反映
    # i=0 の式には u_left, i=n_internal-1 の式には u_right が含まれるため、
    # これを既知数として右辺に移動させる
    b[0] -= bc_left
    b[-1] -= bc_right

    # 連立一次方程式を解く (NumPyの線形代数ソルバーを利用)
    u_internal = np.linalg.solve(A, b)

    # 全体の解を作成 (境界値を追加)
    u = np.concatenate(([bc_left], u_internal, [bc_right]))

    return x, u

# --- パラメータ設定 ---
N = 50           # 格子点数
L = np.pi * 2    # 領域の長さ
u_left = 1.0     # 左端の温度
u_right = 0.0    # 右端の温度

# 熱源分布関数 (例: sin(x))
def f(x):
    return -np.sin(x)

# 境界値問題の解決
x_res, u_res = solve_bvp_1d(N, L, u_left, u_right, f)

# 理論解 (u'' = -sin(x) の一般解 u = sin(x) + ax + b に境界条件を適用)
# u(x) = sin(x) + ( (u_right - u_left)/L ) * x + u_left
u_exact = np.sin(x_res) + ((u_right - u_left) / L) * x_res + u_left

# --- 可視化 ---
plt.figure(figsize=(10, 6))
plt.plot(x_res, u_res, 'bo', label='Numerical Solution (FDM)')
plt.plot(x_res, u_exact, 'r-', label='Exact Solution', alpha=0.7)
plt.title("Steady State Solution of 1D BVP (Stationary Heat Equation)")
plt.xlabel("x")
plt.ylabel("u(x)")
plt.legend()
plt.grid(True)
plt.show()

コードの解説と動作のポイント

  1. 三対角行列の構築:
    A = np.zeros((n_internal, n_internal)) で行列を初期化し、ループで対角成分に -2.0、隣接成分に 1.0 を代入しています。これが二階差分演算子 \(\frac{d^2}{dx^2}\) の離散化表現です。
  2. 境界条件の反映:
    行列計算を行う際、境界点は「変数」ではなく「定数」として扱われます。そのため、連立方程式の右辺ベクトル b から境界値の寄与を引き算する必要があります。
  3. np.linalg.solve の活用:
    NumPyのソルバーは内部で最適化されたLAPACKライブラリを呼び出しており、今回のような小さな行列から、数千次元の行列まで非常に効率的に解くことができます。
  4. 精度の確認:
    グラフを見ると、数値解（青い点）が理論解（赤い線）の上に完璧に乗っていることがわかります。境界値問題は、領域全体のバランスを一度に解くため、適切に離散化すれば非常に高い信頼性で解を得ることができます。

5. 流体シミュレーションにおける役割

この「境界値問題」と「行列計算」の組み合わせは、流体シミュレーションの中で極めて重要な役割を果たします。

圧力ポアソン方程式の物理的意味

流体力学の基礎方程式であるナビエ・ストークス方程式を数値的に解く際、非圧縮性（流体の密度が変化しない）を保つためには、速度場が「質量保存（連続の式）」を満たす必要があります。

この条件を満たすために、各タイムステップごとに「圧力場」を調整する必要があり、その調整は圧力のポアソン方程式という境界値問題を解くことに帰着します。

圧力ポアソン方程式の一般的な形は次の通りです：

$$
\nabla^2 p = f(\mathbf{x})
$$

ここで\(\nabla^2\)はラプラシアン（空間二階微分）、\(p\)は圧力、\(f(\mathbf{x})\)は速度場や外力から計算される既知項です。

数値計算での流れ

この方程式を差分法で離散化すると、前節で説明したような大規模な連立一次方程式（行列方程式）になります。2次元や3次元の場合、未知数の数は格子点数に比例して爆発的に増加します。

たとえば、\(100 \times 100\)の格子なら\(10,000\)個、\(100 \times 100 \times 100\)なら\(1,000,000\)個もの未知数を一度に解く必要があります。

なぜ行列計算が不可欠か

圧力場の決定は、速度場の「発散（divergence）」をゼロにするための調整であり、領域全体のバランスを一度に考慮しなければなりません。したがって、初期値問題のように一歩ずつ進めることはできず、行列計算による一括解決が不可欠です。

疎行列・三対角行列・アルゴリズムの重要性

このとき現れる行列は、1次元なら三対角行列、2次元・3次元なら「疎行列（ほとんどがゼロ）」となります。これらを効率よく解くためには、

  • 三対角行列用のトーマス法
  • 疎行列用の反復法（共役勾配法、SOR法、マルチグリッド法など）
    といった専用アルゴリズムやライブラリの知識が不可欠です。

実際の計算負荷と現代的な工夫

実際の流体シミュレーションでは、圧力ポアソン方程式の解決が計算時間の大部分を占めます。特に高解像度・大規模な3次元計算では、効率的な行列計算がシミュレーターの性能を大きく左右します。

そのため、

  • 疎行列のデータ構造
  • 並列計算やGPU活用
  • 反復法の収束性向上
    など、現代的な数値計算技術が積極的に導入されています。

まとめ

このように、境界値問題と行列計算の知識は、流体シミュレーションの「心臓部」とも言える圧力場の計算に直結しています。効率的なアルゴリズムや疎行列の扱いを理解することは、現代的な流体数値計算を行う上で不可欠なスキルです。

6. まとめと次回の予告

第9回では、時間発展ではなく空間全体のバランスを決定する手法を学びました。

  • 両端の条件から内部を決定する境界値問題。
  • 微分方程式を格子点ごとの関係式に変換した差分方程式。
  • 全体を一括で解くための行列計算と三対角行列。
  • Pythonでの具体的なソルバーの実装と理論解との比較。

これで、初期値問題（時間）と境界値問題（空間）という、数値計算の二大巨頭を攻略しました。

しかし、現実の物理現象は、これらが複雑に絡み合った「偏微分方程式」によって記述されます。熱が伝わる現象、波が伝播する現象、そして流体が流れる現象。これらは数学的にどのように分類され、どのような数値戦略を立てるべきなのでしょうか。

次回、第10回では、「偏微分方程式の分類と物理的特徴」について解説します。楕円型、放物線型、双曲型という3つの型を理解することで、シミュレーションの設計図を正しく描けるようになります。

最後まで読んでいただきありがとうございます。

コメント

タイトルとURLをコピーしました