本記事では、非線形の偏微分方程式の数値解法について、分かりやすい具体例とともに掘り下げていきます。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）を用いた数値解法
有限差分法を用いた偏微分方程式を数値的に解く手順としては、
- 偏微分方程式→差分方程式へ変換
- 差分方程式を使って、各要素の値を計算していく
となります。以下でそれぞれの手順について説明します。
偏微分方程式→差分方程式へ変換
※少し前置きが長くなります。
非線形偏微分方程式の例として、非定常熱伝導方程式で伝導率 $\lambda$ が温度 $u(x,t)$ の関数 $\lambda=\lambda (u)$ である場合を示します。 非定常熱伝導方程式は
$$
c \rho \frac{\partial u}{\partial t} = \frac{\partial}{\partial x} \Big( \lambda \frac{\partial u}{\partial x} \Big) + \frac{\partial}{\partial y} \Big( \lambda \frac{\partial u}{\partial y} \Big) + \frac{\partial}{\partial z} \Big( \lambda \frac{\partial u}{\partial z} \Big)
$$
のように表され、一次元の場合は次のようになります。
$$
\begin{align*} c \rho \frac{\partial u}{\partial t} &= \frac{\partial}{\partial x} \Big( \lambda \frac{\partial u}{\partial x} \Big) \\ \therefore \frac{\partial u}{\partial t} &= \frac{1}{c \rho} \Big( \frac{\partial \lambda}{\partial x} \frac{\partial u}{\partial x} + \lambda \frac{\partial^2 u}{\partial x^2} \Big) \end{align*}
$$
$c, \rho$ は比熱と密度です。ここで、 $\lambda = au + b$ （ $a, b$ は定数 ）とすると、
$$
\frac{\partial u}{\partial t} = \frac{1}{c \rho} \Big\{ a \Big( \frac{\partial u}{\partial x}\Big)^2 + (au+b) \frac{\partial^2 u}{\partial x^2} \Big\} \ \ \ \ \ (1)
$$
初期条件： $u(x, 0) = f(x)$
境界条件： $u(0,t) = g_1(t), u(l,t) = g_2(t)$
$0 \leqq x \leqq l, 0 \leqq t < \infty$ とします。
$(1)$ の非線形偏微分方程式を差分方程式に直していきましょう。 差分方程式への変換は前節の差分近似した偏微分係数を流用します。すなわち
$$
\begin{align*}
\frac{u_{i,j+1} – u_{i,j}}{k} &= \frac{1}{c \rho} \Big\{ \Big( \frac{u_{i+1,j} – u_{i-1,j} }{2h}\Big)^2 + (a u_{i,j}+b) \frac{(u_{i+1,j} – 2u_{i,j} + u_{i-1,j}) }{h^2} \Big\} \\ u_{i,j+1} &= u_{i,j} + \frac{k}{c \rho} \Big\{ \frac{a}{4h^2} ( u_{i+1,j} – u_{i-1,j})^2 + \frac{1}{h^2} (a u_{i,j}+b) (u_{i+1,j} – 2u_{i,j} + u_{i-1,j}) \Big\} \ \ \ \ (2)
\end{align*}
$$
となります。
差分方程式を使って、各要素の値を計算していく
差分方程式 $(2)$ をよく眺めてみましょう。
この式で $i=1,2, … , l-1$ とおくと、 $j=0$ での $u_{i,j}$ は初期条件と境界条件より既知なので、$j=0$ 以降の時刻における $u_{i,j}$ が次々に求まります。
このように既知の時刻行上の $u$ 値より次の時刻行上の $u$ 値を直接求める方法を陽解法といいます。
また、差分方程式 $(2)$ を以下のように整理してみます。
$$
\begin{align*} u_{i,j+1} &= u_{i,j} + \frac{ak}{4h^2} \frac{1}{c \rho} ( u_{i+1,j} – u_{i-1,j})^2 + \frac{k}{c \rho h^2} \lambda(u) (u_{i+1,j} – 2u_{i,j} + u_{i-1,j}) \\ &= \Big( 1 – \frac{2k}{h^2}\frac{1}{c \rho} \lambda(u) \Big)u_{i,j} + \frac{ak}{4h^2}\frac{1}{c \rho}(u_{i+1,j} – u_{i-1,j})^2 + \frac{k}{h^2}\frac{1}{c \rho}\lambda(u)(u_{i+1,j} + u_{i-1,j}) \ \ \ \ (3) \end{align*}
$$
いま、温度はすべて正の値をもつものとすれば、式 $(3)$ の右辺第ニ項、三項は正となり、右辺第一項は以下の条件をもちます。
$$
\begin{align*} 1 – \frac{2k}{h^2}\frac{1}{c \rho} \lambda(u) &\geqq 0 \\ \therefore \frac{k}{h^2}\frac{\lambda(u)}{c \rho} &\leqq \frac{1}{2} \ \ \ \ \ \ (4) \end{align*}
$$
陽解法は簡単でわかりやすいのですが、式 $(4)$ の条件より $h$ に対して $k$ を非常に小さくとる必要があり、計算時間が大幅に長くなる不便が生じます。
式 $(4)$ に満たすように $h,k$ を定めていても、 $\lambda(u)$ は時間変化することからこの条件を満たさなくなる可能性もあります。なので、時間ステップを増すごとに条件を満たすかを確認し、 $k$ の値を変えるか、 あるいは条件を満たさなくなれば計算を中止するようなプログラムにしておくのもよいかもしれません。しかし、もし $\lambda$ の最大値 $\lambda _{max}$ が推定できるならば、 $h$ を定めた時 $k$ は式 $(4)$ より
$$ k \leqq \frac{h^2 c \rho}{2 \lambda _{max}} $$
となるように定めておけばすべての $t$ に対して $k$ の値が有効となります。
例題
前提
熱伝導率 $\lambda$ が温度の関数であるとき、一次元非定常熱電動方程式は次で表されます。
$$
\frac{\partial u}{\partial t} = \frac{1}{c \rho} \Big( \frac{\partial \lambda}{\partial x} \frac{\partial u}{\partial x} + \lambda \frac{\partial^2 u}{\partial x^2} \Big)
$$
いまこの式を適当な基準温度 $u_0$ 、基準長さ $l_0$ 、基準熱伝導率 $\lambda_0$ を用いて無次元化すると次のようになります。
$$
\frac{\partial U}{\partial \tau} = \frac{\partial R}{\partial X} \frac{\partial U}{\partial X} + R \frac{\partial^2 U}{\partial X^2} \ \ \ \ (5) $$
ただし、
$$
U=u/u_0, \ \ \ X=x/l_0, \ \ \ R=\lambda/\lambda_0, \ \ \ \tau=\lambda_0 t/ c \rho l_0^2 \ \ \ \ (6)
$$
であり、これらはすべて無次元量であり、無次元温度、無次元座標、無次元熱伝導率、無次元時間を表します。
問題
ここで式 $(5)$ の問題として次の場合を考えます。細い棒の温度 $U(x, \tau)$ が最初 $0$ でした。この棒の両端の温度を $U=1+\sin2\pi \tau$ に保つとき、棒の中心の温度の時間的変化を調べましょう。ただし、棒の長さは $1$ 、棒の周囲からの熱移動はないものとします。無次元熱伝導率 $R$ は $U$ の一次式として
$$ R=U+1 \ \ \ \ (7) $$
と表されるものとします。
ヒント
温度分布 $U(X,\tau)$ に対する基礎方程式と条件は次式のようになります。
$$
\frac{\partial U}{\partial \tau} = \Big( \frac{\partial U}{\partial X}\Big)^2 + (U + 1)\frac{\partial^2 U}{\partial X^2}; \ \ \ \ \ \ \ 0 \leqq X \leqq 1, 0 \leqq \tau < \infty
$$
初期条件： $U(X, 0) = 0$
境界条件： $U(0,\tau) = U(1,\tau) = 1 + \sin 2 \pi \tau$
格子番号を $i=0 \sim m, j=0 \sim n$ とし、任意格子点の $U$ の値を $U_{i,j}$ とします。
上でもご紹介したように、まずはこの偏微分方程式を差分方程式に変換し、次に差分方程式を解くためのプログラムを実装していきましょう。
実装方法
コードは以下のように実装できます。
グラフの作成、結果のテキストファイル出力もしているので少しコードが長いですが、重要な部分は calculate_equation 関数になります。
以下の実装を main.py という名前でどこかのフォルダ（ここでは~/Desktop/labcode/python/numerical_calculation/non_linear_equation とします）に保存しましょう。
import matplotlib.pyplot as plt
from matplotlib.animation import ArtistAnimation
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
def _set_boundary_condition(
*,
array_2d: list,
grid_counts_x: int,
grid_counts_t: int,
time_delta: float,
) -> None:
"""境界条件の設定"""
for j in range(grid_counts_t + 1):
array_2d[0][j] = 1 + np.sin(2 * np.pi * j * time_delta)
array_2d[grid_counts_x][j] = 1 + np.sin(2 * np.pi * j * time_delta)
def calculate_equation(*, M: int, N: int, H: float, K: 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,
time_delta=K,
)
# 計算
KH = K / (H**2)
K4 = KH / 4
for j in range(N):
for i in range(1, M):
A = U[i][j] + K4 * (U[i + 1][j] - U[i - 1][j]) ** 2
B = KH * (U[i][j] + 1) * (U[i + 1][j] - 2 * U[i][j] + U[i - 1][j])
U[i][j + 1] = A + B
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 (a.u.)")
plt.ylabel("Temperature (a.u.)")
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,
1.05,
f"tau = {time:.3f} (a.u.)",
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,
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_t: {grid_counts_t}\n")
file.write(f"grid_space: {grid_space}\n")
file.write(f"time_delta: {time_delta}\n\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(
grid_counts_t: int, grid_space: float, time_delta: float
) -> None:
"""パラメータの妥当性をチェックする"""
def get_max_u_value() -> float:
max_u = 0.0
for j in range(grid_counts_t + 1):
u = 1 + np.sin(2 * np.pi * j * time_delta)
if max_u < u:
max_u = u
return max_u
max_r_value = get_max_u_value() + 1
r = (grid_space**2) / (2 * max_r_value)
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 = 1000 # 格子点番号nの値
grid_space: float = 0.1 # 刻み幅 (a.u.)
time_delta: float = 0.001 # 刻み幅 (a.u.)
validate_parameters(
grid_counts_t=grid_counts_t,
grid_space=grid_space,
time_delta=time_delta,
)
array_2d = calculate_equation(
M=grid_counts_x,
N=grid_counts_t,
H=grid_space,
K=time_delta,
)
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,
time_delta=time_delta,
)プログラムを実行する
ターミナルを開き、
$ cd Desktop/labcode/python/numerical_calculation/non_linear_equationと入力し、ディレクトリを移動します。あとは以下のコマンドを実行するだけです。（ $ マークは無視してください）
$ python main.py実行結果
~/Desktop/LabCode/python/numerical_calculation/non_linear_equation にanimation.gifというファイルができて，以下のようなgif画像になっていれば成功です！
両端の温度がじわじわ変化するのにつられて、内側の領域の温度が変化している様子が可視化されていますね。
また、 calculated_result.txt というテキストファイルも以下のように作成されているはずです。今回は簡単な情報しか出力していませんが、このように計算結果と計算に用いたパラメータを出力しておくことで、第三者でも計算結果を評価することが可能になります。自分のためにも、他の人のためにも、ぜひ結果は出力する癖をつけましょう。
# This file shows calculated result as below.
# Calculation Parameters.
grid_counts_x: 10
grid_counts_t: 1000
grid_space: 0.1
time_delta: 0.001
# Calculated Matrix U.
1.0 1.0062831439655588 1.0125660398833527 1.0188484397154083 1.0251300954433376 1.0314107590781283 1.0376901826699345 1.0439681183178648 1.0502443181797696 1.0565185344820245 1.0627905195293135 1.0690600257144058 1.0753268055279328 1.0815906115681575 1.
...コードの解説
このPythonコードは、非線形方程式（non-linear equation）の数値解法に関するもので、特定の条件下で方程式の近似解を計算し、その結果をアニメーションとテキストファイルに出力するものです。コードは以下に記すいくつかの主要な部分に分かれています：
1. Pythonにおけるメイン実行ブロック
if __name__ == "__main__":
"""
非線形方程式の数値解法（陽解法）
"""
grid_counts_x: int = 10 # 格子点番号mの値
grid_counts_t: int = 1000 # 格子点番号nの値
grid_space: float = 0.1 # 刻み幅 (a.u.)
time_delta: float = 0.001 # 刻み幅 (a.u.)
validate_parameters(
grid_counts_t=grid_counts_t,
grid_space=grid_space,
time_delta=time_delta,
)
array_2d = calculate_equation(
M=grid_counts_x,
N=grid_counts_t,
H=grid_space,
K=time_delta,
)
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,
time_delta=time_delta,
)1 パラメータの設定:
grid_counts_x: int = 10 # 格子点番号mの値
grid_counts_t: int = 1000 # 格子点番号nの値
grid_space: float = 0.1 # 刻み幅 (a.u.)
time_delta: float = 0.001 # 刻み幅 (a.u.)ここでは、数値シミュレーションに必要なパラメータが設定されています。これには格子点数（grid_counts, grid_counts_t）、格子の刻み幅（grid_space）、時間の刻み幅（time_delta）が含まれます。
2 安定性パラメータの計算と検証:
validate_parameters(
grid_counts_t=grid_counts_t,
grid_space=grid_space,
time_delta=time_delta,
)この部分では、数値解法の安定性に影響を与えるパラメータ（grid_counts_t , grid_space , time_delta）をvalidate_parameters 関数に渡して、その値を検証しています。
3 方程式の計算:
array_2d = calculate_equation(
M=grid_counts_x,
N=grid_counts_t,
H=grid_space,
K=time_delta,
)calculate_equation 関数は、設定されたパラメータを用いて非線形方程式の数値解を計算します。結果は2次元配列 array_2d に格納され、これは時間ごとの温度分布を表します。
4 結果のアニメーション表示:
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,
)この関数は、計算された温度分布の時間発展をアニメーションとして表示します。
5 結果のファイル出力:
output_result_file(
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 関数は計算結果をファイルに出力します。これにより、後で結果を分析したり、他のソフトウェアで表示したりできます。
2. パラメータの妥当性チェック
def validate_parameters(
grid_counts_t: int, grid_space: float, time_delta: float
) -> None:
"""パラメータの妥当性をチェックする"""
def get_max_u_value() -> float:
max_u = 0.0
for j in range(grid_counts_t + 1):
u = 1 + np.sin(2 * np.pi * j * time_delta)
if max_u < u:
max_u = u
return max_u
max_r_value = get_max_u_value() + 1
r = (grid_space**2) / (2 * max_r_value)
if r > 0.5:
raise ValueError("R must be less than 0.5.")この関数 validate_parameters は、指定されたパラメータ（grid_counts_t, grid_space, time_delta）が特定の条件を満たしているかどうかを確認するために使用されます。この関数は数値計算やシミュレーションで重要な役割を果たします。具体的な動作は以下の通りです。
- パラメータの説明:
grid_counts_t(整数): 時間に関する格子点の数を表します。grid_space(浮動小数点数): 空間に関する格子の間隔を表します。time_delta(浮動小数点数): 時間に関する格子の間隔を表します。
- 内部関数
get_max_u_value: この関数は、与えられた時間に関する格子点数 (grid_counts_t) に基づいて、特定の計算を行い、その結果として最大のu値を返します。具体的には、1 + np.sin(2 * np.pi * j * time_delta)の最大値を求めます。これは、時間ステップごとに正弦波の値を計算し、その中で最大値を探すことを意味します。 - 変数
max_r_valueとrの計算:max_r_valueは、get_max_u_valueから得られる最大値に1を加えたものです。rは、空間の格子間隔の二乗 (grid_space ** 2) を2 * max_r_valueで割った値です。
- 条件チェック: 最後に、計算された
rの値が0.5を超える場合、ValueErrorが発生します。これは、数値計算の安定性に関連する重要なチェックです。rが大きすぎると、計算が不安定になるか、間違った結果を生じる可能性があります。
3. 差分方程式の計算
def calculate_equation(*, M: int, N: int, H: float, K: 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,
time_delta=K,
)
# 計算
KH = K / (H**2)
K4 = KH / 4
for j in range(N):
for i in range(1, M):
A = U[i][j] + K4 * (U[i + 1][j] - U[i - 1][j]) ** 2
B = KH * (U[i][j] + 1) * (U[i + 1][j] - 2 * U[i][j] + U[i - 1][j])
U[i][j + 1] = A + B
return Uこの calculate_equation 関数は、特定の差分方程式を計算するための関数です。関数の詳細を説明します。
- 関数の引数:
M(整数): X方向の格子点の数を表します。N(整数): tau（時間）方向の格子点の数を表します。H(浮動小数点数): X方向の格子の間隔を表します。K(浮動小数点数): tau（時間）方向の格子の間隔を表します。
- U値のリストの初期化:
- 二次元リスト
Uは、各格子点での数値解を格納するために使われます。最初に、すべての値を0.0に設定しています。
- 二次元リスト
- 初期条件と境界条件の設定:
_set_initial_condition関数と_set_boundary_condition関数を使用して、Uの初期条件と境界条件を設定します。これらの関数は、Uの特定の要素に初期値を割り当てることで、計算の出発点を定義します。
- 差分方程式の計算:
- この部分は、与えられた差分方程式を具体的に計算する部分です。
KH = K / (H**2)とK4 = KH / 4は、計算に必要な係数を定義しています。これらは差分方程式の特性を表し、空間と時間の格子間隔に基づいています。- 二重の
forループを使用して、各格子点での数値解を計算しています。 AとBは、差分方程式の特定の項を計算するために使用されます。これらの項は、隣接する格子点の値に基づいており、差分方程式の数学的表現をコードに変換しています。- 最終的に、
U[i][j + 1] = A + Bで、新しい時間ステップにおける各格子点の数値解を更新しています。
- 関数の戻り値:
- 関数は、計算された数値解を含む二次元リスト
Uを返します。
- 関数は、計算された数値解を含む二次元リスト
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 (a.u.)")
plt.ylabel("Temperature (a.u.)")
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,
1.05,
f"tau = {time:.3f} (a.u.)",
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 関数は、matplotlibライブラリを使って2次元配列のデータからアニメーションを作成するための関数です。関数の動作は以下の通りです。
- 引数の説明:
array_2d: 数値データを含む2次元リスト。grid_counts_x: X軸方向の格子点数。grid_counts_t: 時間軸方向の格子点数。grid_space: X軸方向の格子の間隔。time_delta: 時間軸方向の格子の間隔。
- 配列の準備:
array_2dをNumPy配列に変換し、転置（.T）することで、x軸と時間軸が適切に整列されるようにします。- x軸の値は
np.arange(grid_counts_x + 1) * grid_spaceを使って生成され、格子点ごとのx座標を表します。
- アニメーションの設定:
- 出力ファイルのパス（
path_gif）と、フレームを保存するための空のリスト（frames）を用意します。 plt.subplotsを使って、図と軸のオブジェクトを作成します。x軸とy軸にラベルを付け、グリッドを表示します。
- 出力ファイルのパス（
- 各フレームの作成:
- 時間ステップ
jごとに、array_2d[j]のデータをプロットします。このプロットは点線スタイルで、青色のマーカーを持つ線グラフです。 - 各フレームには、現在の時間を表すテキスト（
tau = {time:.3f} (a.u.)）も追加されます。
- 時間ステップ
- アニメーションの生成と保存:
ArtistAnimationを使用して、フレームのリストからアニメーションを作成します。各フレーム間のインターバルは100ミリ秒に設定されています。- 作成されたアニメーションはGIFファイルとして保存されます。
5. 結果のファイル出力
def output_result_file(
array_2d: list,
grid_counts_x: int,
grid_counts_t: int,
grid_space: 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_t: {grid_counts_t}\n")
file.write(f"grid_space: {grid_space}\n")
file.write(f"time_delta: {time_delta}\n\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次元の数値配列データをテキストファイルに出力するための関数です。関数の動作は以下の通りです。
- 引数の説明:
array_2d: 数値データを含む2次元リスト。grid_counts_x: X軸方向の格子点数。grid_counts_t: 時間軸方向の格子点数。grid_space: X軸方向の格子の間隔。time_delta: 時間軸方向の格子の間隔。
- テキストファイルの作成と書き込み:
with open("./calculated_result.txt", "w") as fileを使用して、書き込みモードでcalculated_result.txtファイルを開きます。- ファイルには、まず計算結果に関するヘッダ情報として、計算パラメータ（
grid_counts_x,grid_counts_t,grid_space,time_delta）が書き込まれます。
- 2次元配列データの書き込み:
- ファイルにはさらに、計算された2次元配列
array_2dの各行が、スペース区切りの文字列として書き込まれます。これはfor row in array_2dループを使って行われ、各行はstr.joinメソッドを用いて文字列に変換されます。
- ファイルにはさらに、計算された2次元配列
この関数の主な目的は、数値計算の結果を保存し、後で分析やレビュー、または他のプログラムでの使用ができるようにすることです。テキストファイル形式での保存は、データの普遍性とアクセスの容易さを提供し、さまざまなソフトウェアやプラットフォームでの互換性を保証します。
最後に
有限差分法を用いた非線形偏微分方程式の数値解法についてご紹介しました。
Pythonで書いた計算プログラムも惜しみなく公開していますので、ぜひ動かしながら理解を深めてください。
今後は非線形方程式以外の偏微分方程式についても数値的に解けるように、新たな記事を執筆していきますので、そちらもご覧いただけたら幸いです。
Pythonで学ぶ偏微分方程式の数値シミュレーションの技術書を販売中
¥2,500 → ¥1,000 今なら60%OFF!!
計算方法、グラフやアニメーションの作成方法、計算ログの残し方を惜しみなく解説しています！
理論と実装の溝を埋めてくれる良書です！
