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

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

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

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

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

動作検証済み環境

macOS Monterey(12.7.1), python3.11.4

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

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

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

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

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

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

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

放物形方程式

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

初期条件: $u(x,0)= f(x)$
境界条件: $u(0,t) = g_1(t), \ \ \ u(a,t) = g_2(t)$

のもとで解く場合を例に考えます。($\chi, a, g_1(t), g_2(t)$ は既知とします)

まずは式 $(1)$ を差分方程式で表していきましょう。差分方程式への変換は前節の差分近似した偏微分係数を流用します。すなわち

$$ \frac{(u_{i,j+1} – u_{i,j}) }{k} = \chi \frac{(u_{i+1,j} – 2u_{i,j} + u_{i-1,j}) }{h^2} $$

となります。この式を整理すると

$$ \begin{align*} u_{i,j+1} &= (1-2r) u_{i,j} + r(u_{i+1,j} + u_{i-1,j}) \ \ \ \ (2) \\ r &= \chi k/h^2 \end{align*} $$

となります。

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

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

たとえば、 $u_{1,1}$ を求めたいとき、式 $(2)$ を用いると

$$ u_{1,1} = (1-2r) u_{1,0} + r(u_{2,0} + u_{0,0}) $$

となります。下の図にも表しましたが $u_{0,0}, u_{1,0}, u_{2,0}$ が既知であれば $u_{1,1}$ は計算できることがわかります。実際に $u_{0,0}, u_{1,0}, u_{2,0}$ は初期条件と境界条件より既知です。

また $u_{2,1}$ も同様に $u_{1,0}, u_{2,0}, u_{3,0}$ から求めることができます。 1つ前の時刻行上の要素から次の時刻行上の各要素をすべて計算できるのです。

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

いま、温度はすべて正の値をもつものとすれば、式 $(2)$ で $r>0$ ( $\chi > 0$ ) なので右辺最後の項は正であり、求めようとする $u_{i,j+1}$ が正であるためには、少なくとも右辺第1項も正であれば良いです。すなわち

$$ 1-2r \geqq 0 \ \ \ \ \ \therefore r \leqq 1/2 \ \ \ \ (3) $$

の条件が成り立てば良いです。この条件は $t$ と $x$ の格子間隔 $k$ と $h$ は互いに依存していて、それぞれを独立に定めることができないことを示すものです。

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

例題

長さ $1$ $\mathrm{m}$の細い棒の温度が最初 $0$ ℃であったとする。この棒の一端を $0$ ℃、他端を $100$ ℃に保つものとすれば、その後の棒の温度分布はどのように変化するかを調べよ。ただし棒の周囲と外界との間の熱移動はないものとする。温度伝導率 $\chi = 0.07$ $\mathrm{m^2/h}$ とし、 $x$ 方向の格子間隔は $h=0.1$ として計算せよ。

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

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

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

境界条件: $u(0,t) = 0, \ \ \ u(1,t) = 100$

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

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

実装方法

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

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

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

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


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

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


def _set_boundary_condition(
    *,
    array_2d: list,
    grid_counts_x: int,
    grid_counts_t: int,
) -> None:
    """境界条件の設定"""
    # x=0上の境界条件
    for j in range(grid_counts_t + 1):
        array_2d[0][j] = 0.0

    # x=M上の境界条件
    for j in range(grid_counts_t + 1):
        array_2d[grid_counts_x][j] = 100.0


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

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

    # 初期値設定
    _set_initial_condition(
        array_2d=U,
        grid_counts_x=M,
    )
    _set_boundary_condition(
        array_2d=U,
        grid_counts_x=M,
        grid_counts_t=N,
    )

    # 計算
    for j in range(N):
        for i in range(1, M):
            U[i][j + 1] = (1.0 - 2 * R) * U[i][j] + R * (U[i + 1][j] + U[i - 1][j])
    return U


def animate_time_evolution(
    *,
    array_2d: list,
    grid_counts_x: int,
    grid_counts_t: int,
    grid_space: float,
    time_delta: float,
) -> None:
    """
    点線プロットのアニメーションを作成
    array_1dのデータをプロットし、時間発展をアニメーションで表示する
    """
    # 2次元配列を転置
    array_2d = np.array(array_2d).T
    # x軸の値を生成
    x_values = np.arange(grid_counts_x + 1) * grid_space

    # アニメーションの作成
    path_gif = "./animation.gif"
    frames = []

    fig, ax = plt.subplots(figsize=(7, 5))
    plt.xlabel("X (m)")
    plt.ylabel("Temperature (deg.)")
    plt.grid(True)

    for j in range(grid_counts_t + 1):
        frame = plt.plot(
            x_values,
            array_2d[j],
            linestyle="--",
            marker="o",
            color="b",
        )

        time = j * time_delta
        text = ax.text(
            0.05,
            0.95,
            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_2d: list,
    grid_counts_x: int,
    grid_counts_t: int,
    grid_space: 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_t: {grid_counts_t}\n")
        file.write(f"grid_space: {grid_space}\n")

        file.write(f"# Calculated Matrix U.\n")
        for row in array_2d:
            line = " ".join(map(str, row))
            file.write(line + "\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  # 格子点番号mの値
    grid_counts_t: int = 150  # 格子点番号nの値
    grid_space: float = 0.1  # 刻み幅 (meter)
    time_delta: float = 0.07  # 刻み幅 (hour)
    CHI: float = 0.07  # 温度伝導率 (m^2/h)
    r: float = CHI * time_delta / grid_space**2

    validate_parameters(R=r)

    array_2d = calculate_equation(
        M=grid_counts_x,
        N=grid_counts_t,
        R=r,
    )
    animate_time_evolution(
        array_2d=array_2d,
        grid_counts_x=grid_counts_x,
        grid_counts_t=grid_counts_t,
        grid_space=grid_space,
        time_delta=time_delta,
    )
    output_result_file(
        array_2d=array_2d,
        grid_counts_x=grid_counts_x,
        grid_counts_t=grid_counts_t,
        grid_space=grid_space,
    )

プログラムを実行する

ターミナルを開き、

$ cd Desktop/labcode/python/numerical_calculation/parabolic_equation

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

$ python main.py

実行結果

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

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

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

# This file shows calculated result as below.

# Calculation Parameters.
grid_counts_x: 10
grid_counts_t: 150
grid_space: 0.1
# Calculated Matrix U.
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 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 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 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 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 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 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 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.1628413597910449 0.19215280455343298 0.5078696329163107 0.5768974339662957 0.9888295275754517 1.0944978820409426 1.551012363014935 1.6850160194332535 2.1500977769238085 2.3033870859117305 2.754869134748203 2.919486364331002 3.345056608216407 3.514655129233868 3.9084026904690203 4.078207216085826 4.438158013540345 4.604729318860016 4.931218010087263 5.092180311274705 5.3868172332113495 5.540609684963954 5.805642938710244 5.95131380404918 6.189245051403674 6.326289342213703 6.5396493340803055 6.667884291811582 6.859107522575876 6.978578771617452 7.14993879503712 7.260850421392517 7.414431656424816 7.517094533113063 7.654785487732547 7.749579327863821 7.873077918457643 7.960423572142136 8.071248831207127 8.151588191278279 8.2510949198846 8.32487646620043 8.414270798618634 8.481939317334676 8.562294039126904 8.624283432980405 8.696552430253067 8.753280817481818 8.818312359822977 8.870178867002531 8.928727619328772 8.976110426195419 9.028848195395613 9.072103502155386 9.119628784579485 9.15909045512667 9.201936880440936 9.237916576298625 9.276560354419628 9.30934801988395 9.34421449833513 9.374079091470692 9.40554852534445 9.432738915191198 9.461151543852182 9.485897513670679 9.511558029139374 9.534071340342898 9.557252822960534 9.577728305723989 9.598675693777135 9.617292339027777 9.636225490784298 9.653147524996902 9.670263924189564 9.685641853596362 9.70111900280575 9.715090617639314 9.729088158238858 9.741779490718843 9.754441082990708 9.76596731514476 9.777422307782452 9.787888627012553 9.798253541419555 9.80775594310794 9.817135794609015 9.825761832096076 9.834251307332098 9.842080790365985 9.849765297684728 9.856870940996789 9.863827548527153 9.870275572574643 9.876573846834459 9.882424533008871 9.888127289306016 9.893435492060838 9.898599466569156 9.903415084998082 9.908091537193874 9.91245994860796 9.91669520171351 9.920657659738284 9.924493585914668 9.928087585567274 9.931562041810226 9.934821653931884 9.937968873936242 9.940925051253663 9.94377599790988 9.94645685488684 9.949039537545568 9.951470606066904 9.953810366244928 9.956014829053355 9.958134597847254 9.96013350153121 9.962054031646955 9.963866480857336 9.965606555848199 9.9672498903045 9.968826513330644 9.97031646906416 9.971745033240532 9.973095889414207 9.974390331594787 9.975615044136827 9.976787983789322 9.977898306980993 9.978961171633896 9.97996776870085 9.980930907291535 9.981843450963119 9.982716236279243 9.983543500200788 9.984334421485624 9.98508436329505 9.985801109978839 9.9864809467905 9.987130484212924 9.987746761188347
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.33232930569601 0.38550199460737156 1.02862566699029 1.1566123291999375 1.9944726099920935 2.193308758141701 3.1206579701512576 3.375501575863173 4.319178482724783 4.613030878312764 5.528166108224426 5.84569179925722 6.707483430468953 7.036232647080694 7.832876709968047 8.163345229135603 8.891007896364549 9.21625746650859 9.875762089204208 10.190930512393795 10.78566046323644 11.087496612856585 11.622103561246869 11.90857335770403 12.388201582291204 12.658172328950265 13.088007239257207 13.341002663530563 13.726020074978866 13.962033920746805 14.306871876744431 14.526227847942398 14.835131934687686 15.038379387723607 15.315191014429155 15.503027792059532 15.751196595715033 15.924412273006087 16.14702114237609 16.306455540110484 16.50625133889599 16.652764424087223 16.83219034549924 16.96664061502511 17.127867862816757 17.251097045301773 17.39605461549686 17.50887707934032 17.63927906831294 17.742474734298103 17.85984498365045 17.954154844507844 18.0598489527994 18.14597252703566 18.241197376604852 18.319791590683835 18.405622594568168 18.477301711611847 18.55469800590542 18.62003431182767 18.6898521182397 18.74937714592651 18.812381517377624 18.866587642212874 18.923462786833383 18.972805066925787 19.024163426257065 19.069061591347165 19.115451828885053 19.15629234543832 19.198206382985013 19.23534454112701 19.27322376327294 19.306985745267784 19.341226476917605 19.37191137778076 19.402869725987394 19.430751505271836 19.45874764466545 19.484076995008298 19.509398965665984 19.532405088744788 19.55531216628481 19.576204450728373 19.59693014057001 19.61589973934602 19.63465444029971 19.65187574735489 19.668849123857502 19.684481151468045 19.699844248729402 19.71403190825121 19.72793904021881 19.740814329804287 19.753404766087687 19.765087869535066 19.77648734417518 19.787087645458605 19.79740970760804 19.807026725850406 19.816373949987575 19.825098199722618 19.833563270900715 19.841477052457172 19.8491437402447 19.85632186499792 19.863265898160982 19.869776353256487 19.876066205831847 19.881970762805995 19.8876683609894 19.89302313250743 19.89818449070986 19.903040439420394 19.907716232905287 19.912119636180524 19.916355716869766 19.92034859096789 19.924186452278228 19.927806939233463 19.931284135165413 19.934566855484768 19.93771737862154 19.940693752648055 19.94354837522296 19.946246915815095 19.94883349756429 19.951280076540677 19.953623842666502 19.955841933275934 19.95796572549852 19.95997662299627 19.96190112636228 19.963724148606712 19.96546809644894 19.967120766275926 19.968701125472666 19.97019933645987 19.971631474924365 19.97298964202229 19.974287480157376 19.975518676538957 19.9766948242176
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.678223072849 0.77317430304786 1.920659920001083 2.1262957992427545 3.5152739658757195 3.8178358434660695 5.290337400972351 5.66690704951038 7.125862027204491 7.553039712522395 8.943586897888785 9.400956744889763 10.69527289481059 11.166397265202294 12.353212895919171 12.825523840023768 13.903313619289829 14.367583209590299 15.340284720896463 15.78999572320062 16.664377608256746 17.095132757280105 17.87921058105342 18.28822967519746 18.990336888100217 19.3760500680007 20.004313709164084 20.36604232869532 20.92810581226307 21.26581676787664 21.768711249525133 22.0828302292142 22.532933404409487 22.82420416846433 23.227248806176576 23.496627858256222 23.857737042622233 24.10631524579289 24.43005044259536 24.658995034516966 24.949408777147525 25.159920782822336 25.420609278110287 25.613892529952885 25.848045631006375 26.025284535819278 26.235731827384974 26.398075724224455 26.58732823653496 26.735880719113787 26.90616822990047 27.041980204088972 27.195284331506418 27.319349875803905 27.457433286367937 27.570687606317552 27.695119710837847 27.79843864634371 27.910618163790055 28.00481883535354 28.105993588443937 28.191835864935875 28.283120143418117 28.36130868784629 28.44369848279663 28.51488518927639 28.589271568226497 28.654058247190594 28.72123910771391 28.780180310393057 28.84087072005222 28.894476619719494 28.94931792346859 28.998057191589215 29.04762504393458 29.091927675497807 29.13673913383447 29.17699918875484 29.217518986092234 29.254097223411684 29.290743322915873 29.32396971220252 29.357118232335647 29.387294332615795 29.417283919881722 29.444685121017283 29.471820837177972 29.496698462093885 29.521255243983646 29.543838512783307 29.56606425532105 29.58656211428019 29.606680420792387 29.625283240635675 29.643495879006732 29.66037702786069 29.676866126195087 29.692183423266847 29.707113434569226 29.721010490996804 29.73452995275853 29.747137406273072 29.759380517716867 29.77081716779848 29.78190520480803 29.79227905494136 29.802321640332316 29.811730853806736 29.820827098530017 29.82936087400558 29.83760040307046 29.84533977586397 29.852803651191568 29.859822225942132 29.86658377697868 29.872948397040076 29.87907396874723 29.88484532733541 29.890394954109624 29.895628151913478 29.90065616504844 29.905401218696664 29.90995679417597 29.914259099645754 29.91838675232242 29.92228750707999 29.926027536653322 29.92956412403353 29.932953017661717 29.936159356725625 29.93923015260456 29.94213701646044 29.944919632248467 29.947554937583778 29.950076467150925 29.95246553750004 29.954750519123095 29.956916324187826 29.958986982994308 29.960950356141282 29.962826823321013 29.964606659200943 29.966307170249973 29.96792060431825
0.0 0.0 0.0 0.0 0.0 0.0 1.3841287200999999 1.550224166512 3.5558266819368995 3.875482905071675 6.058615251154331 6.491408822123645 8.64630429634062 9.155875530615306 11.19061526702075 11.748013673072142 13.624752341026928 14.207550768487213 15.915208912005768 16.5063322972858 18.047383815684363 18.634093030685396 20.01774194817349 20.590770969240932 21.82926364793067 22.382089531007104 23.488682183842187 24.01691684498369 25.00476387535095 25.505657394442174 26.387227652214904 26.85927017339536 27.646069249962004 28.088676761890632 28.791145751945166 29.204414992860002 29.831928729444808 30.21644772052604 30.77736649024578 31.134069091600747 31.635816479204742 31.96587139213903 32.415022221197596 32.71974874863525 33.122117970977826 33.40292247324913 33.76365003406715 34.02197835229642 34.34560756407103 34.582909748138654 34.873458185944266 35.09116269762622 35.35218547732322 35.5516806846508 35.78632644991227 35.96894772552312 36.18000790213235 36.34702901976509 36.53698099068555 36.689608806867206 36.86065367751052 37.00002531090351 37.15412090647984 37.28130280141199 37.42019248882608 37.53618088446516 37.66141875274124 37.76714118443151 37.88011405751418 37.9764315989264 38.07837829698912 38.16608831619623 38.25811652835227 38.33795578169883 38.42105686511048 38.49370479281361 38.56876677079731 38.63484888990545 38.70266788451335 38.76275919993389 38.82404950224297 38.878677876403884 38.93408082986712 38.983730267276684 39.03382211547615 39.07893593082988 39.12423476032973 39.1652186085513 39.206190499827066 39.24341525404926 39.280479738265 39.31428420767908 39.34781911403739 39.378512598103626 39.408858365295735 39.43672304428518 39.46418655994941 39.489479724396794 39.51433774822322 39.537293871789245 39.55979609078939 39.580628752400926 39.6010005107243 39.619904172797476 39.6383489131797 39.65550056332796 39.6722020126748 39.68776267663511 39.7028868042806 39.717002937918494 39.7306997116545 39.743504479879974 39.75590944186446 39.76752389214513 39.778759574193884 39.789293712121015 39.799470907617845 39.809024681689024 39.81824358936629 39.82690779181563 39.83525904492246 39.843116135070424 39.85068172792416 39.85780658414696 39.86466070672802 39.87112131276943 39.877331102844565 39.88318917381898 39.888815395042826 39.894126948111634 39.89922460164276 39.90404047599182 39.90865935235198 39.91302568275728 39.91721085994805 39.92116950789277 39.92496180114992 39.928550747149266 39.93198711515302 39.93524081565436 39.93835472751373 39.94130443946884 39.94412620635329 39.94680028230645 39.94935735720193 39.95178151350216 39.95409876221479 39.95629632274204 39.95839626895682 39.960388386549305 39.962291433469346
0.0 0.0 0.0 0.0 0.0 2.82475249 3.1072277390000003 6.515291618185 6.990838699876501 10.28567802500483 10.874186888717569 13.865289544120499 14.514714043653305 17.173943734555085 17.85187125198294 20.200244233133237 20.885849515368776 22.956531258511266 23.635835334795306 25.46239479853359 26.12573611106649 27.73874646644774 28.379407038093508 29.805764398150682 30.41916017418213 31.682246693842405 32.265398565180696 33.385470429398914 33.93662666503654 34.931227182877315 35.44957771327409 36.33391585984692 36.819369068512344 37.60664999840194 38.05965578930693 38.76136445220144 39.18277332587264 39.80891637462845 40.199867243838526 40.75917903920832 41.12101024491308 41.62112832311897 41.95530745646071 42.402922132898425 42.71099110896133 43.11197319944987 43.39550568306482 43.75501570657717 44.01558451688631 44.33816621368867 44.57731876373556 44.866979314648106 45.086219495005885 45.3464984506135 45.54727365643856 45.78130226881585 45.96499450889119 46.175546893140854 46.34346711642837 46.53300444680619 46.6863893841008 46.85709814273771 46.997109094258754 47.15093423364754 47.27865734281678 47.41733109141798 47.53377873480295 47.658845664261264 47.764958661137825 47.877797540276966 47.97444794533013 48.076290827465286 48.164285119155835 48.25623404295005 48.33631655997552 48.419358188084125 48.492214698769196 48.56723317120338 48.63349448690014 48.70128272600956 48.761528290776134 48.82279796083897 48.87755936670492 48.93294966235046 48.982714053123445 49.03279946638353 49.078012803822396 49.12330999883115 49.16438017361936 49.2053540802848 49.242653857004164 49.27972307988115 49.31359287045797 49.34713449615603 49.37788496031277 49.408238835738196 49.43615331006934 49.463625854342794 49.48896261393366 49.513830218704776 49.536824576884996 49.559336642781595 49.580202895778065 49.60058454671053 49.61951777074207 49.63797228259045 49.65514999141464 49.67186096713515 49.687444638283836 49.70257795758137 49.71671443556251 49.73042000389985 49.74324278854402 49.75565610732276 49.767286535248644 49.77853011244023 49.78907843933456 49.79926305760812 49.80882944868358 49.81805530612709 49.826730741755775 49.835088478577966 49.84295558171192 49.85052720481414 49.85766099640923 49.8645207123989 49.87098930066221 49.87720426671982 49.88306947561105 49.888700476600604 49.89401841863785 49.89912047794726 49.90394207600214 49.908565006801595 49.91293646921915 49.917125372118015 49.921088625164714 49.92488433762001 49.92847741895037 49.931916921222225 49.93517433776023 49.9382911197142 49.941244173070125 49.94406856568468 49.94674564797162 49.949305123015385 49.95173198569084 49.95405142668302 49.95625142481981 49.958353372073226 49.96034768625909
0.0 0.0 0.0 0.0 5.764800999999999 6.22598508 11.7855591644 12.450863318208 17.150012605833098 17.896911642310922 21.794049659756173 22.564526387612872 25.810123160153292 26.575537484265997 29.305725157660017 30.051669246736246 32.372827798164685 33.09184659809459 35.08413413256094 35.77221508046077 37.49594627457362 38.150870456589196 39.652005515885314 40.27258594141033 41.58678429157701 42.172509716551765 43.32797589390331 43.878832657972566 44.89828642990738 45.414655196612 46.31669950203293 46.79928238030954 47.59936466234757 48.04912484900506 48.76022432237177 49.17833202107592 49.81146087894499 50.1992439454861 50.763820641043104 51.12272001886643 51.6268531498171 51.958383485408575 52.40909203139263 52.71480771320833 53.11819504386023 53.39966164731378 53.7610552575979 54.01982614312147 54.343891463175794 54.581489107905824 54.87232332298102 55.09022487117567 55.35143505474492 55.55106153458655 55.785830275979734 55.968538934711304 56.1796798590167 56.346759099756454 56.536764120493274 56.6894305710199 56.860510312399065 56.999907610329146 57.15402613793314 57.28122507235093 57.42012984701401 57.53612955017686 57.661377348412245 57.76710728057379 57.880086691576054 57.976409206225824 58.0783602102735 58.16607352585791 58.258104574845554 58.33794601241652 58.42104896526628 58.493698339854774 58.56876155007615 58.63484462739266 58.702664434408206 58.76275638426192 58.82404722229468 58.878676016429075 58.93407932322749 58.983729038593594 59.03382111987294 59.07893511915937 59.12423410243497 59.165218072352246 59.206190065096216 59.24341489982516 59.28047945100248 59.31428397366856 59.347818924221535 59.37851244350792 59.40885823987149 59.43672294215274 59.46418647707386 59.4894796569232 59.51433769346269 59.537293827212594 59.55979605460631 59.58062872295098 59.601000486816446 59.61990415334101 59.63834889738276 59.65550055047373 59.67220200223716 59.68776266814271 59.702886797384075 59.7170029323078 59.73069970709774 59.74350447617312 59.75590943885367 59.76752388969608 59.778759572204585 59.789293710502974 59.799470906303455 59.809024680620006 59.81824358849784 59.826907791109335 59.83525904434866 59.84311613460378 59.85068172754504 59.85780658383865 59.86466070647752 59.87112131256574 59.87733110267906 59.8831891736844 59.88881539493347 59.894126948022716 59.89922460157051 59.904040475933066 59.90865935230424 59.91302568271846 59.917210859916516 59.92116950786713 59.924961801129086 59.92855074713232 59.93198711513925 59.93524081564316 59.93835472750462 59.94130443946144 59.94412620634728 59.94680028230156 59.94935735719797 59.95178151349894 59.95409876221217 59.95629632273991 59.95839626895508 59.9603883865479 59.9622914334682
0.0 0.0 0.0 11.764899999999999 12.470794000000001 20.973287229999997 21.821654169000002 27.976535605221 28.83347025970362 33.46148815236243 34.28631391621661 37.8874301138408 38.66760232586288 41.54894281233206 42.28191393900287 44.640152528560556 45.327191386939994 47.29305486516026 47.93667976277521 49.59985412588356 50.20273640137784 51.62653538406966 52.191094658013895 53.42144493398974 53.949766419124444 55.02086697875407 55.51475070031828 56.452712788454 56.913760004379476 57.7389899561037 58.16868430768282 58.89746988452325 59.29723818813686 59.942823114642295 60.31407387320017 60.887399299887264 61.23154244536654 61.74176885421519 62.06022174809121 62.51510412291367 62.80928857633639 63.21545200684931 63.486786161050475 63.84993274382824 64.09982061846023 64.42488807534937 64.65470868013827 64.94599437569501 65.15709156197555 65.41835121656216 65.61202491935842 65.84655243498742 66.02405159375259 66.23474449895384 66.3972617092588 66.58667544612798 66.73534325146824 66.90573665481453 67.0416253093064 67.19499902457885 67.31911552218727 67.45724468527565 67.57053284368465 67.69499504301979 67.7983364392987 67.91053576043653 68.0047513337798 68.10593912342155 68.19179128252402 68.28308414574137 68.36127924171568 68.44367469157226 68.51486573986378 68.58925584484169 68.65404540032947 68.72122871656698 68.780171824478 68.84086385300144 68.8944710142708 68.94931338543827 68.9980534887814 69.04762204507719 69.0919252294768 69.13673715214152 69.17699757292127 69.21751767657928 69.25409615598075 69.29074245759574 69.32396900704032 69.35711766054374 69.3872938667683 69.41728354205388 69.44468481326345 69.47182058752016 69.4966982587797 69.52125507901854 69.54383837846478 69.5660641463189 69.58656202554255 69.60668034876882 69.62528318201063 69.64349583141724 69.66037698912943 69.6768660947506 69.69218339767849 69.70711341379257 69.7210104740914 69.73452993903062 69.74713739510416 69.75938050864636 69.77081716041947 69.78190519881487 69.7922790500662 69.80232163637247 69.81173085058582 69.82082709591364 69.82936087187755 69.83760040134176 69.84533977445801 69.85280365004938 69.85982222501322 69.866583776224 69.87294839642635 69.8790739682486 69.88484532692993 69.89039495378017 69.89562815164558 69.90065616483076 69.90540121851966 69.90995679403215 69.9142590995288 69.91838675222739 69.92228750700271 69.92602753659054 69.92956412398247 69.93295301762025 69.93615935669189 69.93923015257715 69.94213701643814 69.94491963223034 69.94755493756905 69.95007646713894 69.95246553749031 69.95475051911518 69.9569163241814 69.95898698298909 69.96095035613705 69.96282682331756 69.96460665919814 69.96630717024767 69.96792060431639
0.0 0.0 24.009999999999998 24.9704 36.528814 37.45195048 44.418731696500004 45.251054125871995 49.96186238838803 50.70938254783064 54.12779490252233 54.80252202786889 57.40536819865242 58.01820758258344 60.07083453576957 60.630755929526586 62.29352125779747 62.80778333144011 64.18345878708817 64.65786516211101 65.81523873775723 66.25436209510978 67.24110273070542 67.64851103908896 68.4986680561201 68.87715020333232 69.61573384986737 69.96750560449762 70.61337633747245 70.94025196371605 71.50798661708869 71.81151106958104 72.3126319721005 72.59417598958133 73.03797531522672 73.29879952865632 73.6929023971935 73.93419762799738 74.28495423232802 74.50786466251748 74.82062896308612 75.02626370012916 75.30559576576498 75.49503302879916 75.74484916630337 75.9191361458113 76.14282270585565 76.30297319251453 76.50347462931937 76.65046577990472 76.83035409865214 76.96512317995257 77.12665365357915 77.25009523945369 77.39525178920968 77.5082156419125 77.63874828340204 77.74203799454988 77.8594940803082 77.95386645284889 78.05961698317626 78.14578208350581 78.24104403787521 78.31966582242113 78.40552123798712 78.47721865098859 78.55463101229381 78.61997945423354 78.68980783922169 78.7493409137759 78.81235225245702 78.86656371094278 78.92344344565326 78.97278925989495 79.02415064404065 79.06905115024045 79.11544338158076 79.15628544854772 79.19820080059765 79.23533998527407 79.27322007423913 79.30698273576535 79.3412240391235 79.37190938972977 79.40286811506758 79.43075019196137 79.45874658016936 79.48407612742001 79.50939826225671 79.53240451559815 79.55531170148427 79.57620407209139 79.59692983344151 79.61589948920492 79.63465423735902 79.65187558210114 79.66884898976204 79.68448104229348 79.699844160125 79.71403183612468 79.72793898167335 79.7408142821533 79.75340472740396 79.76508783805384 79.7764873186152 79.78708762466003 79.79740969071956 79.80702671210942 79.81637393882876 79.82509819064431 79.83356326352772 79.84147704645935 79.84914373537316 79.85632186103528 79.86326589494222 79.86977635063843 79.87606620370512 79.88197076107629 79.88766835958423 79.89302313136464 79.89818448978143 79.90304043866536 79.90771623229186 79.91211963568168 79.91635571646447 79.92034859063831 79.92418645201045 79.9278069390157 79.93128413498847 79.93456685534089 79.93771737850463 79.94069375255299 79.94354837514574 79.94624691575228 79.94883349751326 79.95128007649917 79.95362384263278 79.95584193324851 79.95796572547624 79.95997662297815 79.96190112634754 79.96372414859474 79.9654680964392 79.96712076626802 79.96870112546623 79.97019933645466 79.97163147492012 79.97298964201885 79.97428748015456 79.97551867653668 79.97669482421571
0.0 49.0 49.980000000000004 61.7645 62.470786 68.14853458 68.71442642679999 72.139467059821 72.6158058628737 74.93362868756762 75.34627002218836 77.0295449026797 77.39382669170935 78.67650695117386 79.00245185448937 80.01475795961687 80.30936556466035 81.13001272761397 81.39841408695793 82.07786308741235 82.32391119118265 82.8959452053247 83.1225563307103 83.61059146465986 83.8199822384468 84.24074699226777 84.4346185394782 84.80040195722457 84.98008578534831 85.30015612106847 85.46672658464222 85.74824797406632 85.90260538357605 86.15124177400078 86.29417107037487 86.5144913258686 86.64670159555898 86.842456206536 86.96460596184944 87.13891969307771 87.25163207849512 87.4071408334821 87.51101202973292 87.64996216581949 87.74556542742798 87.8698874000372 87.95777445944827 88.06913861505824 88.14983963663327 88.24969936109916 88.3237222193753 88.41334795272707 88.4811773172313 88.56168383659842 88.62378034406427 88.69614898359403 88.752948644209 88.81804563175118 88.86995952996448 88.9285512899503 88.97596558769497 89.02873163351026 89.07200785358806 89.11955173563061 89.15902728769898 89.20188595236766 89.23787485803177 89.2765266931846 89.30932046643812 89.3441922505474 89.37406089276112 89.40553382155916 89.43272689479315 89.46114182626597 89.48588957387386 89.5115516070574 89.53406609575897 89.55724857888976 89.57772484136618 89.59867288912018 89.6172900505667 89.63622363738851 89.6531460132728 89.67026269943598 89.6856408549563 89.70111819348224 89.7150899579307 89.7290876234416 89.74177905490463 89.75444072960387 89.76596702723518 89.777422074272 89.78788843681022 89.79825338712254 89.80775581745287 89.81713569265497 89.82576174908266 89.83425123996506 89.84208073552311 89.84976525317171 89.85687090476452 89.86382751911522 89.87027554863741 89.87657382740069 89.8824245171944 89.88812727646533 89.89343548161273 89.89859945808485 89.9034150780953 89.908091531588 89.91245994404747 89.91669519800953 89.92065765672528 89.92449358346735 89.92808758357663 89.93156204019321 89.9348216526167 89.93796887286786 89.94092505038473 89.94377599720396 89.94645685431276 89.94903953707917 89.9514706056876 89.95381036593677 89.95601482880275 89.95813459764365 89.96013350136563 89.96205403151242 89.96386648074794 89.96560655575931 89.96724989023222 89.96882651327192 89.9703164690164 89.97174503320174 89.97309588938266 89.97439033156914 89.975615044116 89.9767879837724 89.97789830696722 89.9789611716227 89.97996776869175 89.98093090728413 89.9818434509571 89.98271623627434 89.98354350019682 89.98433442148239 89.98508436329242 89.9858011099767 89.98648094678876 89.98713048421149 89.98774676118721
100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0

コードの解説

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

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

if __name__ == "__main__":
    """
    放物形方程式の数値解法(陽解法)
    """
    grid_counts_x: int = 10  # 格子点番号mの値
    grid_counts_t: int = 150  # 格子点番号nの値
    grid_space: float = 0.1  # 刻み幅 (meter)
    time_delta: float = 0.07  # 刻み幅 (hour)
    CHI: float = 0.07  # 温度伝導率 (m^2/h)
    r: float = CHI * time_delta / grid_space**2

    validate_parameters(R=r)

    array_2d = calculate_equation(
        M=grid_counts_x,
        N=grid_counts_t,
        R=r,
    )
    animate_time_evolution(
        array_2d=array_2d,
        grid_counts_x=grid_counts_x,
        grid_counts_t=grid_counts_t,
        grid_space=grid_space,
        time_delta=time_delta,
    )
    output_result_file(
        array_2d=array_2d,
        grid_counts_x=grid_counts_x,
        grid_counts_t=grid_counts_t,
        grid_space=grid_space,
    )

このコードは、Pythonにおけるメイン実行ブロック(if __name__ == "__main__")の中に記述されており、放物形方程式の数値解法を実装するスクリプトの主要部分を構成しています。Pythonのプログラムが実装されると、まずこの実行ブロック内のコードが実行されます。

ここでは、陽解法を用いて数値解を計算し、その結果をアニメーションおよびファイルに出力するプロセスが定義されています。以下に各ステップを詳しく説明します。

処理の流れ

  1. グリッドと時間の設定:
    • grid_counts_x: 空間方向の格子点の数を定義します。
    • grid_counts_t: 時間方向の格子点の数を定義します。
    • grid_space: 空間方向の格子点間の距離(刻み幅)です。単位はメートル(m)。
    • time_delta: 時間方向の格子点間の時間間隔(刻み幅)です。単位は時間(h)。
    • CHI: 温度伝導率を定義します。これは熱がどのくらい早く伝わるかを示す値です。
    • r: 陽解法で使用される定数で、温度伝導率、時間の刻み幅、空間の刻み幅に依存します。
  2. パラメータの妥当性チェック:
    • validate_parameters関数を使って、計算に使用するパラメータが妥当かどうかをチェックします。ここでは、rが0.5未満であることを確認します。
  3. 差分方程式の計算:
    • calculate_equation関数を使用して、差分方程式を計算します。この関数は、空間と時間の格子点の数、およびrの値に基づいて数値解を計算します。
  4. アニメーションの作成:
    • animate_time_evolution関数を使用して、計算結果を基にアニメーションを作成します。これにより、時間が経過するにつれて解がどのように変化するかを視覚的に確認できます。
  5. 結果のファイル出力:
    • 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の値が特定の条件を満たしているかを確認します。関数の動作は以下の通りです:

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

3. 差分方程式の計算

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

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

    # 初期値設定
    _set_initial_condition(
        array_2d=U,
        grid_counts_x=M,
    )
    _set_boundary_condition(
        array_2d=U,
        grid_counts_x=M,
        grid_counts_t=N,
    )

    # 計算
    for j in range(N):
        for i in range(1, M):
            U[i][j + 1] = (1.0 - 2 * R) * U[i][j] + R * (U[i + 1][j] + U[i - 1][j])
    return U

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

  1. 関数のパラメータ:
    • M: 空間方向の格子点の数です。
    • N: 時間方向の格子点の数です。
    • R: 差分方程式の計算に使用する定数で、主に時間刻み幅と空間刻み幅、温度伝導率に依存します。
  2. U値の初期化:
    • U: 2次元リスト(配列)で、計算される各格子点の値を格納します。初期値は全て0.0に設定されています。
  3. 初期条件の設定:
    • _set_initial_condition関数を使用して、Uの初期条件を設定します。この例では、特定の境界以外の値を0.0に設定しています。
  4. 境界条件の設定:
    • _set_boundary_condition関数を使用して、Uの境界条件を設定します。この例では、一方の境界を0.0、もう一方を100.0に設定しています。
  5. 差分方程式の計算:
    • 時間と空間の各格子点に対して、差分方程式を使用してUの値を更新します。ここでは、次の時間ステップにおけるある空間点の値は、現在の時間ステップにおけるその点と隣接する2点の値に依存しています。
  6. 結果の返却:
    • 計算されたUの値の2次元リストを返します

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

def animate_time_evolution(
    *,
    array_2d: list,
    grid_counts_x: int,
    grid_counts_t: int,
    grid_space: float,
    time_delta: float,
) -> None:
    """
    点線プロットのアニメーションを作成
    array_1dのデータをプロットし、時間発展をアニメーションで表示する
    """
    # 2次元配列を転置
    array_2d = np.array(array_2d).T
    # x軸の値を生成
    x_values = np.arange(grid_counts_x + 1) * grid_space

    # アニメーションの作成
    path_gif = "./animation.gif"
    frames = []

    fig, ax = plt.subplots(figsize=(7, 5))
    plt.xlabel("X (m)")
    plt.ylabel("Temperature (deg.)")
    plt.grid(True)

    for j in range(grid_counts_t + 1):
        frame = plt.plot(
            x_values,
            array_2d[j],
            linestyle="--",
            marker="o",
            color="b",
        )

        time = j * time_delta
        text = ax.text(
            0.05,
            0.95,
            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関数は、数値的に解かれた放物形偏微分方程式の結果をアニメーションとして視覚化するためのものです。具体的な機能とプロセスは以下の通りです:

  1. パラメータ:
    • array_2d: 2次元の数値データを格納したリスト。これはcalculate_equation関数で計算された結果を含む。
    • grid_counts_x: 空間方向の格子点の数。
    • grid_counts_t: 時間方向の格子点の数。
    • grid_space: 空間方向の格子点間の距離(刻み幅)。
    • time_delta: 時間方向の格子点間の時間間隔(刻み幅)。
  2. 2次元配列の転置とx軸値の生成:
    • 数値データをnumpy配列に変換し、転置します。これにより、行が時間ステップ、列が空間ステップに対応するようになります。
    • x軸の値を生成します。これは、空間格子点に対応する実際の位置を示します。
  3. アニメーションの準備:
    • 出力するアニメーションファイルのパスを設定します(ここでは"./animation.gif")。
    • 各フレームを格納するための空のリストframesを作成します。
  4. グラフの設定:
    • matplotlibを使用してグラフを作成します。
    • x軸(位置)とy軸(温度)のラベルを設定します。
    • グリッドを表示するためにplt.grid(True)を呼び出します。
  5. 各時間ステップでのフレームの生成:
    • 時間の各ステップに対して、その時点での温度分布をプロットします。
    • 点線とマーカーでプロットを行い、色を青("b")に設定します。
    • 各フレームに時間を表示するテキストを追加します。
  6. アニメーションの作成と保存:
    • ArtistAnimationを使用して、フレームのリストからアニメーションを作成します。
    • 作成したアニメーションをGIF形式で保存します。

5. 結果のファイル出力

def output_result_file(
    array_2d: list,
    grid_counts_x: int,
    grid_counts_t: int,
    grid_space: 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_t: {grid_counts_t}\n")
        file.write(f"grid_space: {grid_space}\n")

        file.write(f"# Calculated Matrix U.\n")
        for row in array_2d:
            line = " ".join(map(str, row))
            file.write(line + "\n")

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

  1. 関数のパラメータ:
    • array_2d: 2次元配列(リスト)。これは計算結果を含みます。
    • grid_counts_x: 空間方向の格子点の数。
    • grid_counts_t: 時間方向の格子点の数。
    • grid_space: 空間方向の格子点間の距離(刻み幅)。
  2. テキストファイルの作成と書き込み:
    • open関数を使用して、"./calculated_result.txt"という名前の新しいテキストファイルを書き込みモードで開きます。
    • withステートメントを使うことで、ファイルの正しいクローズを保証しています。
  3. 計算パラメータの書き込み:
    • ファイルに計算に使用されたパラメータ(格子点の数、格子間の距離)を書き込みます。これにより、ファイルの内容を見た人が計算条件を容易に理解できます。
  4. 計算結果の書き込み:
    • 計算された2次元配列(array_2d)の各行をテキストファイルに書き込みます。
    • 各行の数値は空白で区切られ、str関数を使用して文字列に変換されます。

最後に

有限差分法を用いた放物形偏微分方程式の数値解法についてご紹介しました。

Pythonで書いた計算プログラムも惜しみなく公開していますので、ぜひ動かしながら理解を深めてください。

今後は放物形方程式以外の偏微分方程式についても数値的に解けるように、新たな記事を執筆していきますので、そちらもご覧いただけたら幸いです。

コメントを残す

メールアドレスが公開されることはありません。 が付いている欄は必須項目です