【コード付き】二次元放物形の偏微分方程式の数値解法【Python】

本記事では、二次元放物形偏微分方程式の数値解法について、分かりやすい具体例とともに掘り下げていきます。Pythonを活用したアプローチ方法を学びます。

本記事を通して偏微分方程式の数値解法の1つを会得しましょう！

注）差分法の一部の話だけにとどめています。誤差や境界条件などの詳細な議論は冗長化を避けるためにご紹介していません。

動作検証済み環境

macOS Monterey(12.7.1), python3.11.4

Pythonで学ぶ偏微分方程式の数値シミュレーションの技術書を販売中

¥2,500 → ¥1,000 今なら60%OFF!!

計算方法、グラフやアニメーションの作成方法、計算ログの残し方を惜しみなく解説しています！
理論と実装の溝を埋めてくれる良書です！

技術書ページへ

偏微分方程式の数値解法とは

偏微分方程式の数値解法は、偏微分方程式（PDE: Partial Differential Equations）の解を近似的に求めるための手法のことを指します。これらの方程式は、多くの場合、解析的な解が見つけられないため、数値的な手法が必要となります。以下に、主な数値解法をいくつか紹介します。

有限差分法（Finite Difference Method）: 空間や時間を離散的なグリッドに分割し、微分を差分に置き換えることにより近似します。この方法は直感的で実装が比較的簡単ですが、グリッドの選択が解の精度に大きく影響します。
有限要素法（Finite Element Method）: 問題の領域を小さな「要素」に分割し、各要素内で方程式を近似します。この方法は複雑な形状や境界条件を持つ問題に適しています。
有限体積法（Finite Volume Method）: 保存則を基にした方程式に特に適しており、領域を小さな「体積」に分割し、各体積内の平均値を計算します。
スペクトル法（Spectral Method）: 方程式を異なる周波数の波形の和として表し、これらの波形の係数を計算することで解を近似します。高い精度を達成できますが、滑らかな解が必要です。
差分時間領域法（Finite-Difference Time-Domain, FDTD）: 電磁気学で広く使われている、時間領域における有限差分法の一種です。

これらの方法は、物理学、工学、気象学など様々な分野で偏微分方程式の解を求めるために用いられます。それぞれの方法には特徴や適用できる問題の種類が異なり、問題の性質に応じて最適な方法を選択する必要があります。また、高度なコンピュータシミュレーションにおいては、これらの方法が組み合わされることもあります。

本記事では有限差分法（Finite Difference Method）を用いた数値解法についてご紹介します。

差分近似

偏微分方程式の数値解法においては差分法は非常に多く用いられます。

差分法では差分近似を用います。偏微分係数での差分近似の考え方の詳細は別の記事に譲り、ここでは本記事のメイントピックである放物形方程式で使用する差分近似のみを手短にご紹介します。

$u = u(x,y)$ が領域 $0 \leqq x \leqq a, 0 \leqq y \leqq b$ で与えられているものとします。この領域を以下の図のように格子（差分格子）に分割します。格子間隔は等間隔でもよいですが、一般的には図のように不等間隔格子を用いてもよいです。

$x$ 方向への格子番号を $i=0 \sim m$ 、$y$ 方向格子番号を $j=0 \sim n$ とします。点 $(x,y)$ の格子点の番号を図に示したように $(i,j)$ とし、そこでの $u$ の値 $u(x,y)$ を $u_{i,j}$ と表します。格子点 $(i,j)$ の隣の格子点は図に示したように $\Delta x_i,\Delta y_i$ などで表しておきます。

格子間隔が等間隔の場合は

$$
\begin{aligned}
\Delta x_{i-1} = \Delta x_{i} =h \\
\Delta y_{i-1} = \Delta y_{i} =k
\end{aligned}
$$

となり、1階偏微分係数は前進差分近似を用いて次のようになります。

$$
\begin{aligned}
\frac{∂u}{∂x} = \frac{1}{h}(u_{i+1,j} – u_{i,j}) \\
\frac{∂u}{∂y} = \frac{1}{k}(u_{i,j+1} – u_{i,j})
\end{aligned}
$$

また2階偏微分係数は中心差分近似を用いて次のようになります。

$$
\begin{align*}
\frac{∂^2u}{∂x∂y} &= \frac{1}{4hk}(u_{i+1,j+1} – u_{i+1,j-1} – u_{i-1,j+1} + u_{i-1,j-1}) \\
\frac{∂^2u}{∂x^2} &= \frac{1}{h^2}(u_{i+1,j} – 2u_{i,j} + u_{i-1,j}) \\
\frac{∂^2u}{∂y^2} &= \frac{1}{k^2}(u_{i,j+1} – 2u_{i,j} + u_{i,j-1})
\end{align*}
$$

これらの差分近似を用いて、偏微分方程式を差分方程式に変換します。

有限差分法（Finite Difference Method）を用いた数値解法

有限差分法を用いた偏微分方程式を数値的に解く手順としては、

偏微分方程式→差分方程式へ変換
差分方程式を使って、各要素の値を計算していく

となります。以下でそれぞれの手順について説明します。

偏微分方程式→差分方程式へ変換

放物形方程式

$$
\frac{∂u}{∂t} = \chi \Big( \frac{∂^2u}{∂x^2} + \frac{∂^2u}{∂y^2} \Big); \ \ \ \ \ \ 0 \leqq x \leqq a, 0 \leqq y \leqq b, 0 \leqq t < \infty \ \ \ \ (1)
$$

初期条件： $u(x,y,0)= f(x,y)$
境界条件：
$u(0,y,t) = g_1(y,t), \ \ \ u(a,y,t) = g_2(y,t)$
$u(x,0,t) = g_3(x,t), \ \ \ u(x,b,t) = g_4(x,t)$

のもとで解く場合を例に考えます。

ここで $t$ は時間、 $\chi$ は温度伝導率といわれる物性地で定数とします。（$\chi, a, b, g_1(y,t), g_2(y,t), g_3(x,t), g_4(x,t)$ は既知とします）

まずは式 $(1)$ を差分方程式で表していきましょう。そこで $x, y, t$ は次のようになります。

$$
t = k \Delta t, \ \ \ x = i \Delta x, \ \ \ y = j \Delta y
$$

差分方程式への変換は前節の差分近似した偏微分係数を流用します。すなわち

$$
\frac{(u_{i,j,k+1} – u_{i,j,k}) }{\Delta t} = \chi \Big\{ \frac{(u_{i+1,j,k} – 2u_{i,j,k} + u_{i-1,j,k}) }{(\Delta x)^2} + \frac{(u_{i,j+1,k} – 2u_{i,j,k} + u_{i,j-1,k}) }{(\Delta y)^2} \Big\}
$$

となります。この式を $u_{i,j,k+1}$ について整理すると

$$
u_{i,j,k+1} = u_{i,j,k} + \Delta t \chi \Big\{ \frac{u_{i+1,j,k} – 2u_{i,j,k} + u_{i-1,j,k} }{(\Delta x)^2} + \frac{u_{i,j+1,k} – 2u_{i,j,k} + u_{i,j-1,k}}{(\Delta y)^2} \Big\} \ \ \ \ \ (2)
$$

となります。

差分方程式を使って、各要素の値を計算していく

差分方程式 $(2)$ をよく眺めてみましょう。

この式で $k=0,1,2, … , l-1$ とおくと、 $k=0$ での $u_{i,j}$ は初期条件と境界条件より既知なので、$k=0$ 以降の時刻における $u_{i,j}$ が次々に求まります。

このように既知の時刻行上の $u$ 値より次の時刻行上の $u$ 値を直接求める方法を陽解法といいます。

また、差分方程式 $(2)$ を以下のように整理してみます。

$$
\begin{equation}
\begin{split}
u_{i,j,k+1} = \Delta t \chi \Big\{ \frac{u_{i+1,j,k} + u_{i-1,j,k} }{(\Delta x)^2} + \frac{u_{i,j+1,k} + u_{i,j-1,k}}{(\Delta y)^2} \Big\} \\ + \Big\{ 1 – 2 \Delta t \chi \Big( \frac{1}{(\Delta x)^2} + \frac{1}{(\Delta y)^2} \Big) \Big\} u_{i,j,k} \ \ \ \ \ (3)
\end{split}
\end{equation}
$$

いま、温度はすべて正の値をもつものとすれば、式 $(3)$ の右辺第一項は正となり、右辺第二項は以下の条件をもちます。

$$
1 – 2 \Delta t \chi \Big( \frac{1}{(\Delta x)^2} + \frac{1}{(\Delta y)^2} \Big) \geqq 0
$$

$$
\therefore \Delta t \chi \Big( \frac{1}{(\Delta x)^2} + \frac{1}{(\Delta y)^2} \Big) \leqq \frac{1}{2} \ \ \ \ (4)
$$

陽解法は簡単でわかりやすいのですが、式 $(4)$ の条件より $\Delta x, \Delta y$ に対して $\Delta t$ を非常に小さくとる必要があり、計算時間が大幅に長くなる不便が生じます。この不便をなくすために有効な方法としてクランク-ニコルソンの陰解法があります。陰解法については、機会があれば別記事でご紹介します。

例題

温度 $u(x,y,t)$ が最初 $0$ ℃であった正方形板（辺の長さ $1$ $\mathrm{m}$ ）の各辺を、3辺で $u=0$ 、$y=1$ の辺で $u=4x(1-x)$ に保つとき、時間経過とともに、板内の温度分布はどのように変化するか調べよ。板の温度伝導率は $\chi = 0.07$ $\mathrm{m^2/h}$ とする。

温度分布 $u(x, y,t)$ に対する基礎方程式と条件は次式のようになります。

$$
\frac{∂u}{∂t} = \chi \Big( \frac{∂^2u}{∂x^2} + \frac{∂^2u}{∂y^2} \Big); \ \ \ \ \ \ 0 \leqq x \leqq 1, 0 \leqq y \leqq 1, 0 \leqq t < \infty \ \ \ \ (1)
$$

初期条件： $u(x,y,0)= 0$

境界条件：
$u(0,y,t) = 0, \ \ \ u(1,y,t) = 0$
$u(x,0,t) = 0, \ \ \ u(x,1,t) = 4x(1-x)$

格子番号を $i=0 \sim m, j=0 \sim n, k=0 \sim l$ とし、任意格子点の $u$ の値を $u_{i,j,k}$ とします。

上でもご紹介したように、まずはこの偏微分方程式を差分方程式に変換し、次に差分方程式を解くためのプログラムを実装していきましょう。

実装方法

コードは以下のように実装できます。

グラフの作成、結果のテキストファイル出力もしているので少しコードが長いですが、重要な部分は calculate_equation 関数になります。

以下の実装を main.py という名前でどこかのフォルダ（ここでは~/Desktop/labcode/python/numerical_calculation/2d_parabolic_equation とします）に保存しましょう。

import matplotlib.pyplot as plt
from matplotlib.animation import ArtistAnimation
import numpy as np


def _set_initial_condition(
    *,
    array_3d: list,
    grid_counts_x: int,
    grid_counts_y: int,
) -> None:
    """
    初期条件の設定

    境界以外のU値を全て0.0とおく
    """
    for i in range(1, grid_counts_x):
        for j in range(1, grid_counts_y):
            array_3d[0][i][j] = 0.0


def _set_boundary_condition(
    *,
    array_3d: list,
    grid_counts_x: int,
    grid_counts_y: int,
    grid_counts_t: int,
    grid_space_x: float,
) -> None:
    """境界条件の設定"""
    for k in range(grid_counts_t + 1):
        for i in range(grid_counts_x + 1):
            array_3d[k][i][0] = 0.0
            array_3d[k][i][grid_counts_y] = (
                4 * i * grid_counts_x * (1.0 - i * grid_space_x)
            )

        for j in range(grid_counts_y + 1):
            array_3d[k][0][j] = 0.0
            array_3d[k][grid_counts_x][j] = 0.0


def calculate_equation(
    *,
    M: int,
    N: int,
    L: int,
    DX: float,
    DY: float,
    DT: float,
    CHI: float,
) -> list:
    """
    差分方程式を計算する

    本プログラムの肝となる関数
    """
    # U値のリスト
    U: list = [
        [[0.0 for _ in range(N + 1)] for _ in range(M + 1)] for _ in range(L + 1)
    ]

    # 初期値設定
    _set_initial_condition(
        array_3d=U,
        grid_counts_x=M,
        grid_counts_y=N,
    )
    _set_boundary_condition(
        array_3d=U,
        grid_counts_x=M,
        grid_counts_y=N,
        grid_counts_t=L,
        grid_space_x=DX,
    )

    # 計算
    for k in range(L):
        for j in range(1, N):
            for i in range(1, M):
                A = (U[k][i + 1][j] - 2 * U[k][i][j] + U[k][i - 1][j]) / (DX**2)
                B = (U[k][i][j + 1] - 2 * U[k][i][j] + U[k][i][j - 1]) / (DY**2)
                U[k + 1][i][j] = U[k][i][j] + CHI * DT * (A + B)
    return U


def animate_time_evolution(
    *,
    array_3d: list,
    grid_counts_x: int,
    grid_counts_y: int,
    grid_counts_t: int,
    grid_space_x: float,
    grid_space_y: float,
    time_delta: float,
) -> None:
    """
    アニメーションを作成
    array_3dのデータをプロットし、時間発展をアニメーションで表示する
    """

    # アニメーションの作成
    fig, ax = plt.subplots(figsize=(5, 5))
    extent = [
        0 - grid_space_x * 0.5,
        grid_counts_x * grid_space_x + grid_space_x * 0.5,
        0 - grid_space_y * 0.5,
        grid_counts_y * grid_space_y + grid_space_y * 0.5,
    ]
    plt.xlabel("X (m)")
    plt.ylabel("Y (m)")
    path_gif = "./animation.gif"

    frames = []
    for k in range(grid_counts_t + 1):
        U_2dim = [
            [0.0 for _ in range(grid_counts_x + 1)] for _ in range(grid_counts_y + 1)
        ]
        for j in range(grid_counts_y + 1):
            for i in range(grid_counts_x + 1):
                U_2dim[i][j] = array_3d[k][i][j]

        frame = plt.imshow(
            np.array(U_2dim).T,
            cmap="viridis",
            interpolation="none",
            aspect="auto",
            origin="lower",
            extent=extent,
        )
        time = k * time_delta
        text = ax.text(
            0.05,
            1.05,
            f"t = {time:.2f} (h)",
            transform=ax.transAxes,
            fontsize=14,
            verticalalignment="top",
        )
        frames.append([frame] + [text])

    anim = ArtistAnimation(fig, frames, interval=100)
    anim.save(path_gif, writer="pillow")


def output_result_file(
    *,
    array_3d: list,
    grid_counts_x: int,
    grid_counts_y: int,
    grid_counts_t: int,
    grid_space_x: float,
    grid_space_y: float,
    time_delta: float,
) -> None:
    """2次元配列をテキストファイルに書き出す"""
    with open("./calculated_result.txt", "w") as file:
        file.write(f"# This file shows calculated result as below.\n\n")
        file.write(f"# Calculation Parameters.\n")
        file.write(f"grid_counts_x: {grid_counts_x}\n")
        file.write(f"grid_counts_y: {grid_counts_y}\n")
        file.write(f"grid_counts_t: {grid_counts_t}\n")
        file.write(f"grid_space_x: {grid_space_x}\n")
        file.write(f"grid_space_y: {grid_space_y}\n")
        file.write(f"time_delta: {time_delta}\n\n")

        file.write(f"# Calculated Matrix U at the each time.\n")
        for k in range(grid_counts_t + 1):
            file.write(f"time = {time_delta * k} h\n")

            for row in array_3d[k]:
                line = " ".join(map(str, row))
                file.write(line + "\n")
            file.write("\n")


def validate_parameters(R: float) -> None:
    """パラメータの妥当性をチェックする"""
    if R > 0.5:
        raise ValueError("R must be less than 0.5.")


if __name__ == "__main__":
    """
    ２次元放物形方程式の数値解法（陽解法）
    """
    grid_counts_x: int = 10  # 格子点数
    grid_counts_y: int = 10  # 格子点数
    grid_counts_t: int = 500  # 格子点数
    grid_space_x: float = 0.1  # 刻み幅 (m)
    grid_space_y: float = 0.1  # 刻み幅 (m)
    time_delta: float = 0.02  # 刻み幅 (h)
    CHI: float = 0.07  # 温伝導率 (m^2/h)
    r: float = CHI * (1 / grid_space_x**2 + 1 / grid_space_y**2) * time_delta

    validate_parameters(R=r)

    array_3d = calculate_equation(
        M=grid_counts_x,
        N=grid_counts_y,
        L=grid_counts_t,
        DX=grid_space_x,
        DY=grid_space_y,
        DT=time_delta,
        CHI=CHI,
    )
    animate_time_evolution(
        array_3d=array_3d,
        grid_counts_x=grid_counts_x,
        grid_counts_y=grid_counts_y,
        grid_counts_t=grid_counts_t,
        grid_space_x=grid_space_x,
        grid_space_y=grid_space_y,
        time_delta=time_delta,
    )
    output_result_file(
        array_3d=array_3d,
        grid_counts_x=grid_counts_x,
        grid_counts_y=grid_counts_y,
        grid_counts_t=grid_counts_t,
        grid_space_x=grid_space_x,
        grid_space_y=grid_space_y,
        time_delta=time_delta,
    )

プログラムを実行する

ターミナルを開き、

$ cd Desktop/labcode/python/numerical_calculation/2d_parabolic_equation

と入力し、ディレクトリを移動します。あとは以下のコマンドを実行するだけです。（ $ マークは無視してください）

$ python main.py

実行結果

~/Desktop/LabCode/python/numerical_calculation/2d_parabolic_equation にanimation.gifというファイルができて，以下のようなgif画像になっていれば成功です！

上端 $y=1.0$ の境界の温度がじわじわ下方へ広がっている様子が可視化されていますね。
$t=5.5$ $\mathrm{h}$ にはほぼ熱平衡に達していることもわかります。

また、 calculated_result.txt というテキストファイルも以下のように作成されているはずです。今回は簡単な情報しか出力していませんが、このように計算結果と計算に用いたパラメータを出力しておくことで、第三者でも計算結果を評価することが可能になります。自分のためにも、他の人のためにも、ぜひ結果は出力する癖をつけましょう。

# This file shows calculated result as below.

# Calculation Parameters.
grid_counts_x: 10
grid_counts_y: 10
grid_counts_t: 500
grid_space_x: 0.1
grid_space_y: 0.1
time_delta: 0.02

# Calculated Matrix U at the each time.
time = 0.0 h
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 36.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 64.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 84.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 96.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 100.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 95.99999999999997
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 83.99999999999999
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 63.999999999999986
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 35.99999999999999
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0

...

コードの解説

このPythonコードは、二次元放物形方程式（2d_parabolic equation）の数値解法に関するもので、特定の条件下で方程式の近似解を計算し、その結果をアニメーションとテキストファイルに出力するものです。コードは以下に記すいくつかの主要な部分に分かれています：

1. Pythonにおけるメイン実行ブロック

if __name__ == "__main__":
    """
    ２次元放物形方程式の数値解法（陽解法）
    """
    grid_counts_x: int = 10  # 格子点数
    grid_counts_y: int = 10  # 格子点数
    grid_counts_t: int = 500  # 格子点数
    grid_space_x: float = 0.1  # 刻み幅 (m)
    grid_space_y: float = 0.1  # 刻み幅 (m)
    time_delta: float = 0.02  # 刻み幅 (h)
    CHI: float = 0.07  # 温伝導率 (m^2/h)
    r: float = CHI * (1 / grid_space_x**2 + 1 / grid_space_y**2) * time_delta

    validate_parameters(R=r)

    array_3d = calculate_equation(
        M=grid_counts_x,
        N=grid_counts_y,
        L=grid_counts_t,
        DX=grid_space_x,
        DY=grid_space_y,
        DT=time_delta,
        CHI=CHI,
    )
    animate_time_evolution(
        array_3d=array_3d,
        grid_counts_x=grid_counts_x,
        grid_counts_y=grid_counts_y,
        grid_counts_t=grid_counts_t,
        grid_space_x=grid_space_x,
        grid_space_y=grid_space_y,
        time_delta=time_delta,
    )
    output_result_file(
        array_3d=array_3d,
        grid_counts_x=grid_counts_x,
        grid_counts_y=grid_counts_y,
        grid_counts_t=grid_counts_t,
        grid_space_x=grid_space_x,
        grid_space_y=grid_space_y,
        time_delta=time_delta,
    )

パラメータの設定:
- ここでは、数値シミュレーションに必要なパラメータが設定されています。これには格子点数（grid_counts_x, grid_counts_y, grid_counts_t）、格子の刻み幅（grid_space_x, grid_space_y）、時間の刻み幅（time_delta）、そして温伝導率（CHI）が含まれます。
安定性パラメータの計算と検証:
- r: float = CHI * (1 / grid_space_x**2 + 1 / grid_space_y**2) * time_delta
  - validate_parameters(R=r) この部分では、数値解法の安定性に影響を与えるパラメータ r を計算し、validate_parameters 関数を使ってその値を検証しています。
方程式の計算:
- array_3d = calculate_equation(...)
  - calculate_equation 関数は、設定されたパラメータを用いて2次元放物形方程式の数値解を計算します。
  - 結果は3次元配列 array_3d に格納され、これは時間による空間の温度分布を表します。
結果のアニメーション表示:
- animate_time_evolution(...)
  - この関数は、計算された温度分布の時間発展をアニメーションとして表示します。
結果のファイル出力:
- output_result_file(...)
  - 最後に、output_result_file 関数は計算結果をファイルに出力します。これにより、後で結果を分析したり、他のソフトウェアで表示したりできます。

2. パラメータの妥当性チェック

def validate_parameters(R: float) -> None:
    """パラメータの妥当性をチェックする"""
    if R > 0.5:
        raise ValueError("R must be less than 0.5.")

このvalidate_parameters関数は、数値解法において重要なパラメータの妥当性をチェックするための関数です。具体的には、放物型偏微分方程式の数値解法において使用されるパラメータRの値が特定の条件を満たしているかを確認します。関数の動作は以下の通りです：

関数のパラメータ:
- R: 数値解法において使用されるパラメータです。この値は通常、空間と時間の刻み幅、および物理的な係数（例えば、熱伝導率）に基づいて計算されます。
妥当性チェック:
- この関数は、Rが0.5より大きい場合にエラーを発生させます。R > 0.5の場合、ValueErrorを投げて処理を中断します。
重要性:
- この妥当性チェックは、数値解法における安定性や正確性を確保するために重要です。特に、差分方程式を用いた数値解法では、Rの値によって計算の安定性が大きく左右されます。Rがある閾値より大きい場合、計算が不安定になるか、誤った結果を生む可能性があります。

3. 差分方程式の計算

def calculate_equation(
    *,
    M: int,
    N: int,
    L: int,
    DX: float,
    DY: float,
    DT: float,
    CHI: float,
) -> list:
    """
    差分方程式を計算する

    本プログラムの肝となる関数
    """
    # U値のリスト
    U: list = [
        [[0.0 for _ in range(N + 1)] for _ in range(M + 1)] for _ in range(L + 1)
    ]

    # 初期値設定
    _set_initial_condition(
        array_3d=U,
        grid_counts_x=M,
        grid_counts_y=N,
    )
    _set_boundary_condition(
        array_3d=U,
        grid_counts_x=M,
        grid_counts_y=N,
        grid_counts_t=L,
        grid_space_x=DX,
    )

    # 計算
    for k in range(L):
        for j in range(1, N):
            for i in range(1, M):
                A = (U[k][i + 1][j] - 2 * U[k][i][j] + U[k][i - 1][j]) / (DX**2)
                B = (U[k][i][j + 1] - 2 * U[k][i][j] + U[k][i][j - 1]) / (DY**2)
                U[k + 1][i][j] = U[k][i][j] + CHI * DT * (A + B)
    return U

このcalculate_equation関数は、二次元放物形偏微分方程式の数値解法の核となる部分です。この関数は、熱伝導方程式のような方程式の差分近似を計算するために使用されます。具体的な動作は以下の通りです：

関数のパラメータ:
- def calculate_equation(*, M, N, L, DX, DY, DT, CHI) -> list:
  - M, N: 空間の格子点数（X軸およびY軸方向）
  - L: 時間の格子点数
  - DX, DY: 空間の刻み幅（X軸およびY軸方向）
  - DT: 時間の刻み幅
  - CHI: 温度伝導率
3次元リストの初期化:
- U = [[[0.0 for _ in range(N + 1)] for _ in range(M + 1)] for _ in range(L + 1)]
  - ここで3次元リスト U を初期化しています。このリストは各格子点における温度値を保持します。各次元はそれぞれ時間、X軸、Y軸に対応しています。
初期条件と境界条件の設定:
- _set_initial_condition関数を使用して、Uの初期条件を設定します。この例では、特定の境界以外の値を0.0に設定しています。
- _set_boundary_condition関数を使用して、Uの境界条件を設定します。この例では、一方の境界を0.0、もう一方を100.0に設定しています。
差分方程式の計算:
- 時間と空間の各格子点に対して、差分方程式を使用してUの値を更新します。ここでは、次の時間ステップにおけるある空間点の値は、現在の時間ステップにおけるその点と隣接する2点の値に依存しています。
結果の返却:
- 計算されたUの値の3次元リストを返します

4. アニメーションの作成

def animate_time_evolution(
    *,
    array_3d: list,
    grid_counts_x: int,
    grid_counts_y: int,
    grid_counts_t: int,
    grid_space_x: float,
    grid_space_y: float,
    time_delta: float,
) -> None:
    """
    アニメーションを作成
    array_3dのデータをプロットし、時間発展をアニメーションで表示する
    """

    # アニメーションの作成
    fig, ax = plt.subplots(figsize=(5, 5))
    extent = [
        0 - grid_space_x * 0.5,
        grid_counts_x * grid_space_x + grid_space_x * 0.5,
        0 - grid_space_y * 0.5,
        grid_counts_y * grid_space_y + grid_space_y * 0.5,
    ]
    plt.xlabel("X (m)")
    plt.ylabel("Y (m)")
    path_gif = "./animation.gif"

    frames = []
    for k in range(grid_counts_t + 1):
        U_2dim = [
            [0.0 for _ in range(grid_counts_x + 1)] for _ in range(grid_counts_y + 1)
        ]
        for j in range(grid_counts_y + 1):
            for i in range(grid_counts_x + 1):
                U_2dim[i][j] = array_3d[k][i][j]

        frame = plt.imshow(
            np.array(U_2dim).T,
            cmap="viridis",
            interpolation="none",
            aspect="auto",
            origin="lower",
            extent=extent,
        )
        time = k * time_delta
        text = ax.text(
            0.05,
            1.05,
            f"t = {time:.2f} (h)",
            transform=ax.transAxes,
            fontsize=14,
            verticalalignment="top",
        )
        frames.append([frame] + [text])

    anim = ArtistAnimation(fig, frames, interval=100)
    anim.save(path_gif, writer="pillow")

このanimate_time_evolution関数は、数値的に解かれた放物形偏微分方程式の結果をアニメーションとして視覚化するためのものです。具体的な機能とプロセスは以下の通りです：

関数の定義とパラメータ:
- 関数はいくつかのパラメータを受け取ります。これらは、アニメーションを作成するために必要な情報を提供します。
  - array_3d: 3次元のデータ配列。通常、時間による空間内の状態変化を表します。
  - grid_counts_x, grid_counts_y, grid_counts_t: X軸、Y軸、時間軸に沿った格子点の数。
  - grid_space_x, grid_space_y: X軸、Y軸に沿った格子の間隔。
  - time_delta: 時間ステップの大きさ。
アニメーションの基本設定:
- matplotlib を使ってアニメーションの基本的なフレームワークを設定しています。ここで、軸ラベルと画像の大きさ、保存先のパスなどを設定します。
フレームの準備:
- frames = [] for k in range(grid_counts_t + 1): ...
  - ここでは時間ステップごとにアニメーションの各フレームを準備しています。各フレームは、特定の時間における状態（温度分布など）を表す2次元の配列です。
データの変換とプロット:
- U_2dim = ...
  - この部分では、3次元データ配列から特定の時間ステップの2次元データを抽出し、それをプロットしています。
フレームとテキストの追加:
- frame = plt.imshow(...) text = ax.text(...) frames.append([frame] + [text])
  - imshow を使って2次元データを画像として表示し、時間を示すテキストも追加しています。
アニメーションの作成と保存:
- anim = ArtistAnimation(fig, frames, interval=100) anim.save(path_gif, writer="pillow")
  - ArtistAnimation を使用して、準備したフレームからアニメーションを作成し、GIFファイルとして保存しています。

5. 結果のファイル出力

def output_result_file(
    *,
    array_3d: list,
    grid_counts_x: int,
    grid_counts_y: int,
    grid_counts_t: int,
    grid_space_x: float,
    grid_space_y: float,
    time_delta: float,
) -> None:
    """2次元配列をテキストファイルに書き出す"""
    with open("./calculated_result.txt", "w") as file:
        file.write(f"# This file shows calculated result as below.\n\n")
        file.write(f"# Calculation Parameters.\n")
        file.write(f"grid_counts_x: {grid_counts_x}\n")
        file.write(f"grid_counts_y: {grid_counts_y}\n")
        file.write(f"grid_counts_t: {grid_counts_t}\n")
        file.write(f"grid_space_x: {grid_space_x}\n")
        file.write(f"grid_space_y: {grid_space_y}\n")
        file.write(f"time_delta: {time_delta}\n\n")

        file.write(f"# Calculated Matrix U at the each time.\n")
        for k in range(grid_counts_t + 1):
            file.write(f"time = {time_delta * k} h\n")

            for row in array_3d[k]:
                line = " ".join(map(str, row))
                file.write(line + "\n")
            file.write("\n")

このoutput_result_file関数は、計算された2次元配列のデータをテキストファイルに書き出すためのものです。具体的には、放物形偏微分方程式の数値解法によって得られた結果を保存するために使用されます。関数の動作は以下の通りです： 1.

関数の定義とパラメータ:
- 関数はいくつかのパラメータを受け取ります。これらは、書き出すデータとそのデータの生成に使用されたパラメータに関する情報です。
  - array_3d: 3次元のデータ配列。通常、時間による空間内の状態変化を表します。
  - grid_counts_x, grid_counts_y, grid_counts_t: X軸、Y軸、時間軸に沿った格子点の数。
  - grid_space_x, grid_space_y: X軸、Y軸に沿った格子の間隔。
  - time_delta: 時間ステップの大きさ。
ファイルへの書き出し:
- with open("./calculated_result.txt", "w") as file: ...
  - ここでは、calculated_result.txt という名前のファイルを作成し、書き込みモードで開いています。with open(...) 構文は、ファイル操作を行う際に一般的に使用され、ファイルの適切な開閉を保証します。
シミュレーションパラメータの書き込み:
- file.write(f"# Calculation Parameters.\\\\n") ...
  - まず、シミュレーションに使用されたパラメータをファイルに書き込んでいます。これにより、後で結果を分析する際に参照できるようになります。
計算結果の書き出し:
- file.write( ...
  - 各時間ステップにおける計算結果をファイルに書き出しています。これは、時間ステップごとに3次元配列の各「スライス」（2次元配列）を取り出し、それを文字列に変換してファイルに書き込んでいます。各行は配列の一行を表し、各値はスペースで区切られています。