始めに

熱伝導の数値シミュレーションを通して，陽解法と陰解法の特性を比較します．併せて，陰解法に登場する幾つかの反復解法を勉強してみます．かなり基礎的な内容ですが，何事も基礎からということでお付き合いください．

以下の内容は主にこちらを参考に，同記事に則した内容となります．

数値解法

自然現象はとても難しく，解析解が与えられる問題は非常に少ないです．しかし，コンピュータの力を借りて数値解を逐次求めていくことで，多くの問題に対して現象を追跡することができます．

本稿では，1次元の線形な熱方程式：

\frac{\partial u}{\partial t} = D \frac{\partial^2 u}{\partial x^2}

を対象とした数値計算を実施します．ここでは，$u$が温度を代表する物理量としています（$x \in \Omega, t \in [0, T], u: \Omega \times [0, T] \rightarrow \mathbb{R}$）．

陽解法 (Explicit method)

陽解法は以下のような漸化式に基づいた時間積分を実行します：

x_{n+1} = x_n + h f(x_n, t_0 + nh)

ここで，$h$は離散化パラメータと呼ばれるものです¹．上式より，$x_0$を与えれば，$x_n (n=1, 2, \dots, N)$が次々と計算できます．

FTCS (Forward-Time Centered-Space) method

FTCS法についてです．
上記の熱方程式は，時間に関する1階の微分と空間に関する2階の微分を含みます．まず，1階微分について考えます．任意の関数$v (\xi)$について，$(\xi + \Delta \xi)$周りでのTaylor展開を書いてみると：

v(\xi + \Delta \xi) \approx v(\xi) + \nabla v(\xi) \Delta \xi

ここから，以下のように$\nabla v(\xi)$を得ることで1階の微分を評価します：

\begin{align}
v(\xi + \Delta \xi) &= v(\xi) + \nabla v(\xi) \Delta \xi \\
\nabla v(\xi) \Delta \xi &= v(\xi + \Delta \xi) - v(\xi) \\
\therefore \nabla v(\xi) &= \frac{v(\xi + \Delta \xi) - v(\xi)}{\Delta \xi}
\end{align}

次に，2階の微分を評価します．$(\xi \pm \Delta \xi)$まわりでの$v$のTaylor展開は：

\begin{align}
v(\xi + \Delta \xi) &\approx v(\xi) + \nabla v(\xi) \Delta \xi + \frac{1}{2} \nabla^2 v(\xi) \Delta \xi^2 \\
v(\xi - \Delta \xi) &\approx v(\xi) - \nabla v(\xi) \Delta \xi + \frac{1}{2} \nabla^2 v(\xi) \Delta \xi^2
\end{align}

辺々の和を取ると，

\begin{align}
v(\xi + \Delta \xi) + v(\xi - \Delta \xi) &= 2 v(\xi) + \nabla^2 v(\xi) \Delta \xi^2 \\
\nabla^2 v(\xi) \Delta \xi^2 &= v(\xi + \Delta \xi) - 2 v(\xi) + v(\xi - \Delta \xi) \\
\therefore \nabla^2 v(\xi) &= \frac{v(\xi + \Delta \xi) - 2 v(\xi) + v(\xi - \Delta \xi)}{\Delta \xi^2}
\end{align}

以上より，$v(\xi)$の1階，および2階の導関数を近傍での$v (= v(\xi \pm \Delta \xi))$との差分で表現します．これを上記の熱方程式に適用すると，以下の有限差分近似を得ます：

\frac{u^{n+1}_{i} - u^{n}_{i}}{\Delta t} = D \frac{u^{n}_{i+1} - 2u^{n}_{i} + u^{n}_{i-1}}{\Delta x^2}

ここで，$u^{n}_{i}$はタイムステップ$n$，節点$i$での$u$（温度）であり，$\Delta t$は時間増分，$\Delta x$は空間解像度です．このとき，時間については前進差分（Foward-Time），空間については中心差分（Centered-Space）で離散化していることからFTCS法と呼ばれているそうです．
この子をじっと見つめてあげると，$n+1$ステップ目での物理量$u$は$n$ステップ目での$u$から計算できそうです．少しだけ式を操作して，

u^{n+1}_{i} = u^{n}_{i} + D \frac{\Delta t}{\Delta x^2} (u^{n}_{i+1} - 2u^{n}_{i} + u^{n}_{i-1})

ここで，$i$は任意の節点を示しているので，全ての節点$i$に対して漸化式的に時間積分を実行できます（境界条件の取り扱いは省きます）．

CFL (Courant–Friedrichs–Lewy) 条件拡散数の条件

陽解法は時間増分に敏感です． ~~時間増分に関する制約条件はCFL条件と呼ばれており，これが満足されていないと計算が発散します．~~ ← 誤解です．申し訳ありません．「正しい時間増分を選択しなければ計算が発散する」は正しいですが，今回は拡散方程式なので拡散数の条件が重要です（CFL条件は移流や非線形性をシミュレートするときに重要です）．陽解法の時間積分スキームを再掲します：

u^{n+1}_{i} = u^{n}_{i} + D \frac{\Delta t}{\Delta x^2} (u^{n}_{i+1} - 2u^{n}_{i} + u^{n}_{i-1})

ここで，$ d = D \frac{\Delta t}{\Delta x^2} $を「拡散数」と呼びます．熱伝導方程式や拡散方程式が安定である条件は，この$ d $が$ 1 / 2 $以下である必要があります．

\begin{align}
d = D \frac{\Delta t}{\Delta x^2} &\le \frac{1}{2} \\
\therefore \Delta t &\le \frac{\Delta x^2}{2D}
\end{align}

上記を満足する中で，できるだけ大きい$ \Delta t $を採用します．ここで，空間解像度が極めて高い（$\Delta x$が非常に小さい）ときや拡散係数$ D $が非常に大きいとき（あるいは流体の粘性が高いときなど），$\Delta t$に関する制約はとても厳しくなります．こうした状況では，次に紹介する陰解法が有効です．

陰解法 (Implicit method)

陰解法は，次のような形の時間積分を実行します：

x_{n+1} = x_n + h f(x_{n+1}, t_0 + (n+1)h)

右辺に$x_{n+1}$が登場する点が陽解法とは異なり，少し工夫した解き方をしなければうまくいきません．一方で，陰解法には時間増分に関する制約が無いので，陽解法より大きい$\Delta t$でも安定して計算できます（それが精度の良い解であるかは別の話）．

また，有限差分法や有限要素法では，節点$i$での物理量$x_i$は，近傍節点$j$での物理量$x_j$から計算することができます．これを連立方程式化すると，係数行列$A$，ソース項$b$，未知変数ベクトル$x$の関係は，$Ax=b$というシンプルな式で結ばれます．ここから，何らかの方法で$x$を求めることが陰解法のミソです．以下では，定常法（stationary method）と呼ばれる方法を紹介します．実用上は非定常な反復法（non-stationary iterative method）として，CG法（Conjugate gradient method）などが広く用いられますが，前段階として簡単な方法から紹介します．

Jacobi method

ヤコビの反復法は，シンプルかつパワフルな反復法です．

Ax=b

という，$n$次元の連立方程式（$A \in \mathbb{R}^{n \times n}, x, b \in \mathbb{R}^{n}$）に対して，$A$を狭義下三角行列$L$，対角行列$D$，狭義上三角行列$U$に分解して以下のような操作を行います．

\begin{align}
Ax &= b \\
(L + D + U)x &= b \\
D x &= b - (L + U)x \\
x &= D^{-1} (b - (L + U)x) \\
\therefore x^{(k+1)} &= D^{-1} (b - (L + U)x^{(k)}) 
\end{align}

Weighted Jacobi method

上記のJacobi methodに，新たなパラメータ（重み$w$）を導入したヤコビの重み付き反復法を紹介します．ヤコビの方法における更新式を，以下のように書き改めます．

\begin{align}
Ax &= b \\
x^{(k+1)} &= w D^{-1} (b - (L + U)x^{(k)}) + (1-w) x^{(k)}
\end{align}

ここで，重み$w$を$w=1.0$とすると先ほどのヤコビ法に帰着します．重み$w$には$w=2/3 \approx 0.67$程度の値がよく選ばれるようです．これを勉強するのに参考にできるものが少なかったのですが，直観的には$x^{(k+1)}$と$x^{(k)}$がある程度近いことを期待して項（$+ (1-w) x^{(k)}$）を追加しているものと解釈しています．

Gauss-Seidel method

ガウス-ザイデルの方法は，Jacobi法や重み付きJacobi法より高速です．先ほどと同様に，$n$次元の連立方程式（$A \in \mathbb{R}^{n \times n}, x, b \in \mathbb{R}^{n}$）を考え，やはり$A$を狭義下三角行列$L$，対角行列$D$，狭義上三角行列$U$に分解して以下の操作を実施します．

\begin{align}
Ax &= b \\
(L + D + U)x &= b \\
(L + D) x &= b - Ux \\
\therefore (L + D) x^{(k+1)} &= b - Ux^{(k)} \\
x^{(k+1)} &= (L + D)^{-1} (b - Ux^{(k)})
\end{align}

要素ごとに書き表すと「解いた変数を即，次の未知変数の計算に使う」という大胆な方法に書き換えることができます（$x^{(k+1)} = D^{-1} (b - Lx^{(k+1)} - Ux^{(k)})$）．この方法に書き換えると，$k$番目の解と$k+1$番目の解を混ぜて使うことができ，メモリ効率の良さでもJacobi法を上回ることになるため，数値解析ではこちらが好まれるようです．ただし，今回はPythonで実装することを念頭においています．Pythonはループで要素にアクセスするのが非常に遅く，その代わりスライス操作が高速です．スライス操作との相性が良い式ということで，上記の式展開を記しました．

数値実験

陽解法 vs. 陰解法

まずは，金属棒の熱伝導シミュレーションです．代表的な金属か分かりませんが，取り敢えず高価な金属（金，銀，銅）の材料パラメータは以下の通りです．

Material	Density (kg / m3)	Conductivity (W / m / K)	Capacity (J / kg / K)
Gold	19320	318	126
Silver	10490	429	233
Copper	8960	398	386

熱伝導効率が最も高いのは銀ですが，熱伝導効率に対するコストパフォーマンスは銅が最高なので，本稿では銅棒のシミュレーションを実施します．1mの銅棒に対して，実時間60分のシミュレーションを行います．初期条件は領域全域で20度の温度分布で，左端にはDirichlet条件（100度）を，右端にはNeumann条件を課します．また，$\Delta x = 10^{-2}$とします．

今回は時間増分$\Delta t$に対する感度を調べてみます．陽解法（FTCS法）と陰解法（ヤコビ法）で銅棒の熱伝導を計算します．$\Delta x = 10^{-2}$で空間を離散化しているので，CFL条件は$\Delta t \le 4.35 \times 10^{-1}$を要求しています．本実験では，$\Delta t$を$10^{-3}$から$10^{2}$までlogスケールで変化させながら計算してみます．

以下の表には，時間増分と陽解法（FTCS法）・陰解法（ヤコビ法）で計算終了までにかかった時間をまとめます：

Case	時間増分$\Delta t$	陽解法（FTCS法）	陰解法（ヤコビ法）
1	$10^{-3}$	15.32 s	43.68 s
2	$10^{-2}$	1.51 s	4.43 s
3	$10^{-1}$	0.15 s	0.59 s
4	$10^{0}$	(diverged)	0.11 s
5	$10^{1}$	(diverged)	0.04 s

チューニングには自信がないですが，同じ時間増分なら陽解法の方が速そうですね！ただし，陽解法では時間増分に応じて解けたり解けなかったりしました．一方で，陰解法では$\Delta t = 10^{1}$などの大胆な時間増分でも安定して解くことができました．

具体的な解を見てみるために，最も細かい時間増分を採用したCase 1と最も大きい時間増分を採用したCase 5での解を比較します：

Case	$\Delta t$	陽解法（FTCS法）	陰解法（ヤコビ法）
1	$10^{-3}$
5	$10^{1}$

以上より，Dirichlet条件が与えられた左端境界では$u$がべったり100度になっており，Neumann条件が課された右端境界では$u$の勾配がゼロになっています．また，Case 1では陽解法と陰解法でほぼ同一の解を得ることができました．Case 5の陽解法ではトンデモな解が出てしまいましたが，陰解法では解けました．ただし，Case 1と比較すると温度の伝播が遅いようで，安定に解けることと精度が良いことは別の話ということが分かります．

ヤコビ法 vs. 重み付きヤコビ法 vs. ガウス-ザイデル法

こちらは簡単な検証で，連立方程式：

\begin{bmatrix}
   3 & 2 & 1 \\
   1 & 4 & 1 \\
   2 & 2 & 5
\end{bmatrix}
\begin{bmatrix}
   x_1 \\
   x_2 \\
   x_3
\end{bmatrix}
=
\begin{bmatrix}
   10 \\
   12 \\
   21
\end{bmatrix}

をヤコビ法，重み付きヤコビ法，ガウス-ザイデル法の3つの方法で解きます．解は$x = [x_1, x_2, x_3]^\top = [1, 2, 3]^\top$です．初期値$x^0 = [0, 0, 0]^\top$の下で解くと，次のような動きになりました．

以上より，各手法できちんと解に収束している様子が確認できます．なお，ここでは，以下に定義する残差$R$が$10^{-6}$以下になると反復を打ち切ることにしました：

R = \frac{\sqrt{\sum_j (x^{(k+1)}_j - x^{(k)}_j)^2}}{\sqrt{\sum_j (x^{(k)}_j)^2}}

残差$R$と反復回数との関係を以下のグラフに示します．

どの手法も着実に収束していますが，Jacobiよりweighted Jacobi，weighted JacobiよりGauss-Seidelが速そうです（Jacobi法の振動のひどさがすごいですね，今回だけなのでしょうか）取り敢えず．一般的に認知されている性能と同様の結果を得られました．

終わりに

1次元の熱方程式を対象に陽解法と陰解法を紹介し，数値シミュレーションを実施しました． ~~CFL条件~~ 拡散数の条件を満足する$\Delta t$と，これを満足しない$\Delta t$を用いたシミュレーションから，陽解法が発散する様子と陰解法で安定して解ける様子を確認しました．

おまけ的に，紹介した反復法で連立方程式を解きました．どの方法もきちんと収束することを確認して，紹介した3つの方法の中ではGauss-Seidel法が最速であることを確認しました．

参考文献

[1] こちらの記事
[2] FTCS法
[3] ヤコビ法
[4] ガウス-ザイデル法

付録（Pythonコード）

数値実験のPythonコードです．参考まで．
（Pythonはelement-wiseな操作よりもslice操作が高速なので，これを利用しています．）

陽解法 vs. 陰解法

"""
Explicit (FTCS) vs. Implicit (Jacobi)
"""

import numpy as np
import matplotlib.pyplot as plt

def main():
    # conputational domain
    xmin, xmax = 0., 1.               # 1 m
    tmin, tmax = 0., 60. * 60. * 1.   # 1 hour
    n_max = int(1e3)   # max iteration (for implicit)
    c_tol = 1e-6       # convergence tolerance (for implicit)

    # material parameter (copper)
    rho = 8960.   # density (kg / m3)
    lam = 398.    # conductivity (W / m / K)
    cap = 386.    # capacity
    D   = lam / (rho * cap)

    # discretization
    dx = 1e-2
    dt = 1e-1
    dt_star = dx ** 2 / (2. * D)   # diffusion condition
    beta = D * dt / (dx ** 2.)    # coef
    nx = 1 + int((xmax - xmin) / dx)
    nt = 1 + int((tmax - tmin) / dt)
    x = np.linspace(xmin, xmax, nx)

    # unknown vector
    u = np.ones(shape=(nx))

    # initial condition
    u_ic = 20.
    u *= u_ic

    # boundary condition
    u_bc = 100.
    u[0]  = u_bc    # Dirichlet
    u[-1] = u[-2]   # Neumann

    # conputational setting
    mode = 1   # 0 - explicit (FTCS), 1 - implicit (Jacobi)

    if mode == 0:
        print(">>>>> explicit")
        if dt < dt_star:
            print(">>>>> dt: %.3e, dt_star: %.3e" % (dt, dt_star))
            print(">>>>> CFL met")
        else:
            print(">>>>> dt: %.3e, dt_star: %.3e" % (dt, dt_star))
            print(">>>>> CFL NOT met, program terminating now")
            exit()
    elif mode == 1:
        print(">>>>> implicit")
        print(">>>>> dt: %.3e, dt_star: %.3e" % (dt, dt_star))

    # FDM: Finite Difference Method
    if mode == 0:   # explicit
        for n in range(0, nt-1):
            # copy of u
            v = np.copy(u)

            # slow: element-wise operation
            # for i in range(1, nx-1):
            #     u[i] = v[i] + beta * (v[i+1] - 2. * v[i] + v[i-1])

            # fast: slice operation
            u[1:-1] = v[1:-1] + beta * (v[2:] - 2. * v[1:-1] + v[:-2])

            # boundary condition
            u[0] = u_bc
            u[-1] = u[-2]

            # damp-out
            if n % int(nt / 10) == 0:
                print("step: %d / %d" % (n, nt))
                plt.figure(figsize=(8, 4))
                plt.plot(x, u)
                plt.xlim(-.1, 1.1)
                plt.ylim(0, 120)
                plt.title("t: %.1f min (step: %d / %d)" % (n * dt / 60., n, nt))
                plt.xlabel("x")
                plt.ylabel("u")
                plt.grid(alpha=.3)
                plt.savefig("./res/" + str(n) + ".png")
                plt.clf()
                plt.close()

    elif mode == 1:   # implicit
        for n in range(0, nt-1):
            v = np.copy(u)
            for n_ in range(n_max):
                w = np.copy(u)

                # slow
                # for i in range(1, nx-1):
                #     u[i] = 1. / (1. + 2. * beta) \
                #             * (v[i] + beta * (w[i+1] + w[i-1]))

                # fast
                u[1:-1] = 1. / (1. + 2. * beta) \
                        * (v[1:-1] + beta * (w[2:] + w[:-2]))

                # residual
                if n_ % 10 == 0:
                    r_ = np.sqrt(np.sum(u - w) ** 2) / np.sum(w ** 2)
                    if r_ < c_tol:
                        break

            # boundary condition
            u[0]  = u_bc
            u[-1] = u[-2]

            # damp-out
            if n % int(nt / 10) == 0:
                print("step: %d / %d" % (n, nt))
                plt.figure(figsize=(8, 4))
                plt.plot(x, u)
                plt.xlim(-.1, 1.1)
                plt.ylim(0, 120)
                plt.title("t: %.1f min (step: %d / %d)" % (n * dt / 60., n, nt))
                plt.xlabel("x")
                plt.ylabel("u")
                plt.grid(alpha=.3)
                plt.savefig("./res/" + str(n) + ".png")
                plt.clf()
                plt.close()

if __name__ == "__main__":
    main()

ヤコビ法 vs. 重み付きヤコビ法 vs. ガウス-ザイデル法

"""
Jacobi vs. weighted Jacobi vs. Gauss-Seidel
"""

import numpy as np

def Jacobi(A, b, x, n_max=int(1e3), c_tol=1e-6):
    """
    Jacobi method

    params:
        A: coef
        b: source
        x: unknown
    """

    D = np.diag(A)         # diag elements
    D = np.diagflat(D)     # diag matrix
    C = A - D              # C = L + U
    E = np.linalg.inv(D)

    for n in range(n_max):
        y = np.copy(x)
        x = np.matmul(E, (b - np.matmul(C, y)))
        res = np.sqrt(np.sum((x - y) ** 2)) / np.sqrt(np.sum(y ** 2))

        if n % 10 == 0:
            if res < c_tol:
                print(">>>>> converged, loop terminating now")
                break
    return x

def w_Jacobi(A, b, x, w=2./3., n_max=int(1e3), c_tol=1e-6):
    """
    weighted Jacobi method

    params:
        A: coef
        b: source
        x: unknown
    """

    D = np.diag(A)
    D = np.diagflat(D)
    C = A - D
    E = np.linalg.inv(D)

    for n in range(n_max):
        y = np.copy(x)
        x = w * np.matmul(E, (b - np.matmul(C, y))) + (1. - w) * y
        res = np.sqrt(np.sum((x - y) ** 2)) / np.sqrt(np.sum(y ** 2))

        if n % 10 == 0:
            if res < c_tol:
                print(">>>>> converged, loop terminating now")
                break
    return x

def Gauss_Seidel(A, b, x, n_max=int(1e3), c_tol=1e-6):
    """
    Gauss-Seidel method

    params:
        A: coef
        b: source
        x: unknown
    """

    D = np.diag(A)         # diag elements
    D = np.diagflat(D)     # diag matrix
    L = np.tril(A, k=-1)   # strictly lower triangular
    U = np.triu(A, k= 1)   # strictly upper triangular
    P = L + D
    Q = np.linalg.inv(P)

    for n in range(n_max):
        y = np.copy(x)
        x = np.matmul(Q, (b - np.matmul(U, y)))
        res = np.sqrt(np.sum((x - y) ** 2)) / np.sqrt(np.sum(y ** 2))

        if n % 10 == 0:
            if res < c_tol:
                print(">>>>> converged, loop terminating now")
                break
    return x

def main():
    A = np.array([[3., 2., 1.],
                  [1., 4., 1.], 
                  [2., 2., 5.]])
    b = np.array([10., 12., 21.])
    x = np.zeros(shape=(3))

    x1 = Jacobi(A, b, x)
    x2 = w_Jacobi(A, b, x)
    x3 = Gauss_Seidel(A, b, x)

    print("Jacobi         :", x1)
    print("weighted Jacobi:", x2)
    print("Gauss-Seidel   :", x3)

if __name__ == "__main__":
    main()

$h$の直観的理解は，$\frac{f(x+h) - f(x)}{h} \xrightarrow{h \rightarrow 0} \frac{df}{dx}$，あるいは$\sum_i f(x_i)h \xrightarrow{h \rightarrow 0} \int f(x)dx$より，離散的な表現を厳密な表現に近づけるパラメータと解釈しています（ただし計算上は，$h$が大きすぎると打ち切り誤差が，$h$を小さくしすぎると丸め誤差が効いてきます）． ↩