はじめに
Fortranは今日でも現役で数値計算に使用されています。多くの場合、Pythonなどのインタプリタ言語よりもコンパイラ言語であるFortranのほうが実行速度が速いです（あくまでなんの工夫もなく生の実装をした場合の比較ですが）。今回は双曲型偏微分方程式を解く問題を題材にして、Fortranで数値計算を体験してみましょう。
Pythonに計算条件を指定する部分 (namelistファイルの作成) と計算結果を描画する部分を担当させ、Fortranで計算した結果を読み込んで結果を描画してみましょう。
今回使用するプログラムの計算（偏微分方程式数値解法）の中身についてはこちらを参照してください。
動作検証済み環境
macOS Sonoma(14.7), gfortran (gcc version 14.1.0),
Processor: 2.4 GHz 8-Core Intel Core i9,
Memory: 32 GB 2400 MHz DDR4
Pythonで学ぶ偏微分方程式の数値シミュレーションの技術書を販売中
計算方法、グラフやアニメーションの作成方法、計算ログの残し方を惜しみなく解説しています！
理論と実装の溝を埋めてくれる良書です！
解きたい問題
今回解きたい方程式は次のようなものです。
次の $u(x, t)$ にかんする偏微分方程式を与えられた条件のもとで解いてください。
$$ \frac{\partial^2 u}{\partial t^2} = (x + 1)\frac{\partial^2 u}{\partial x^2} + xe^{-t};\quad 0\le x\le 1, \quad 0\le t< \infty $$
初期条件： $u(x, 0) = \sin \pi x, \frac{\partial u}{\partial t}(x, 0) = 0$
境界条件： $u(0, t)= 0, u(1, t) = 0$
方程式を解くプログラムをFortranで実装する
こちらのページで解説されているように方程式を差分化して数値計算で答えを求めるプログラムを作成します。
Fortranで実装すると次のようになります。適当な場所（~/Desktop/LabCode/fortran/04_hyperbolic-PDE とします）に適当な名前（ここではcalc-hyperbolic-PDE.f90 とします）で保存してください。
program main implicit none ! declare variables ! constants real(8), parameter :: pi=4.0d0*datan(1.0d0) ! parameters for calculation integer(4) :: M, N real(8) :: H, K ! variables for I/O file character(len=256) :: outfile integer(4) :: lnml=10, lout=61 integer(4) :: iost ! variables for calculation result real(8), allocatable :: U(:, :) ! working variables integer(4) :: i, j real(8) :: R, S, Q ! setting for namelist file namelist /input/ M, N, H, K ! read the namelist file open (lnml, file="input.nmlst", action="read", iostat=iost) if (iost /= 0) then write (*, *) "ERROR: cannot open namelist file!" stop end if read (lnml, nml=input) close (lnml) ! allocate and initialize allocate(U(M + 1, N + 1)) U = 0.0d0 ! varidate the parameter R R = 2.0d0 * (K / H) if (R > 1.0d0) then write(*, *) "Error: R must be less than 1.0." stop end if ! set the initial condition do i = 2, M U(i, 1) = sin(pi*(i-1)*H) end do ! set the boundary conditions do j = 1, N+1 U( 1, j) = 0.0d0 U(M+1, j) = 0.0d0 end do ! compute the equation ! calculation of U[i, 2] do i = 2, M R = (K / H) ** 2 * ((i - 1) * H + 1) S = 2.0d0 * (1.0d0 - R) Q = K ** 2 * (i - 1) * H * exp(0.0d0) U(i, 2) = (R * (U(i + 1, 1) + U(i - 1, 1)) + S * U(i, 1) + Q) / 2.0d0 end do ! calculation of U[i, j + 1] do j = 2, N do i = 2, M R = (K / H) ** 2 * ((i - 1) * H + 1) S = 2.0d0 * (1.0d0 - R) Q = K ** 2 * (i - 1) * H * exp(-dble(j - 1) * K) U(i, j + 1) = R * (U(i + 1, j) + U(i - 1, j)) & & + S * U(i, j) - U(i, j - 1) + Q end do end do ! save result outfile = "calculated_result.txt" open(unit=lout, file=outfile, status="replace") write(lout, "(a)") "# This file shows calculated result as below." write(lout, "(a)") "" write(lout, "(a)") "# Calculation Parameters." write(lout, "(a,i4)") "grid_counts_x: ", M write(lout, "(a,i4)") "grid_counts_t: ", N write(lout, "(a,es12.3)") "grid_space: ", H write(lout, "(a,es12.3)") "time_delta: ", K write(lout, "(a)") "" write(lout, "(a)") "# Calculated Matrix U." do j = 1, N + 1 do i = 1, M + 1 write(lout, "(f12.6, x)", advance="no") U(i, j) end do write(lout, *) end do close(lout) ! deallocate memory deallocate(U) ! end of program stop
end program mainコンパイル
保存したら、コンパイルしましょう。まず、プログラムを保存した場所に移動します：
$ cd ~/Desktop/LabCode/fortran/04_hyperbolic-PDE続いて、ソースファイルをコンパイルします。次のように打ち込んでください：
$ gfortran -o calc-hyperbolic-PDE calc-hyperbolic-PDE.f90コンパイルが正常に終了すると、実行ファイル calc-hyperbolic-PDE ができているはずです！
namelistファイルの用意と可視化をPythonにやらせる
特定の形式でパラメータを設定するnamelistファイルを作成し、計算を実行し、可視化をmatplotlibなりにやらせるという一連の操作をいちいちやるのは大変です。
そこで、前回やったように、これらをPythonスクリプトを書いて実行します。以下のようなスクリプトを書いて、main.pyという名前で~/Desktop/LabCode/fortran/04_hyperbolic-PDEに保存します。結果のアニメーションを作る部分は、こちらの記事のスクリプトからコピペしました。
import os
import matplotlib.pyplot as plt
from matplotlib.animation import ArtistAnimation
import numpy as np
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")
if __name__ == "__main__": """ 双曲形方程式の数値解法（陽解法） """ grid_counts_x: int = 10 # 格子点番号mの値 grid_counts_t: int = 40 # 格子点番号nの値 grid_space: float = 0.1 # 刻み幅 (meter) time_delta: float = 0.05 # 刻み幅 (hour) # namelistファイルの作成 nmlst_path = f"input.nmlst" with open(nmlst_path, "w") as fnmlst: fnmlst.write("&input\n") fnmlst.write(f"M={grid_counts_x},\n") fnmlst.write(f"N={grid_counts_t},\n") fnmlst.write(f"H={grid_space:.6e},\n".replace("e", "d")) fnmlst.write(f"K={time_delta:.6e},\n".replace("e", "d")) fnmlst.write(f"/\n") # 計算の実行 cmds = f"./calc-hyperbolic-PDE" print(cmds) os.system(cmds) # 計算結果を読み込む array_2d = np.genfromtxt("calculated_result.txt", skip_header=9) # アニメーションを作成する animate_time_evolution( array_2d=array_2d.T, grid_counts_x=grid_counts_x, grid_counts_t=grid_counts_t, grid_space=grid_space, time_delta=time_delta, )プログラムの実行
ターミナルを開き、
$ cd ~/Desktop/LabCode/fortran/04_hyperbolic-PDEと入力し、ディレクトリを移動します。あとは以下のコマンドを実行するだけです。（ $マークは無視してください）
$ python3 main.py実行結果
~/Desktop/LabCode/fortran/04_hyperbolic-PDEにinput.nmlst、calculated_result.txt、animation.gifというファイルができて、以下のようなgif画像になっていれば成功です！
実行時間の比較
計算をPythonで実行する場合と、Fortranでする場合で速度を比較してみましょう。計算を実行して、結果をファイルに出力するまでの時間を計測します。
実行時間を計測するには、time モジュールを使って
start_time = time.time()
""" 計測したい処理
"""
end_time = time.time()
print(f"{end_time - start_time}s")としました。
比較するプログラムは次の二つです。計算して計算結果をファイルに出力するまでの時間を計測します。
Fortranで計算する場合
import os
import time
if __name__ == "__main__": """ 双曲形方程式の数値解法（陽解法） """ grid_counts_x: int = 1000 # 格子点番号mの値 grid_counts_t: int = 4000 # 格子点番号nの値 grid_space: float = 0.001 # 刻み幅 (meter) time_delta: float = 0.0005 # 刻み幅 (hour) start_time = time.time() # namelistファイルの作成 nmlst_path = f"input.nmlst" with open(nmlst_path, "w") as fnmlst: fnmlst.write("&input\n") fnmlst.write(f"M={grid_counts_x},\n") fnmlst.write(f"N={grid_counts_t},\n") fnmlst.write(f"H={grid_space:.6e},\n".replace("e", "d")) fnmlst.write(f"K={time_delta:.6e},\n".replace("e", "d")) fnmlst.write(f"/\\n") # 計算の実行 cmds = f"./calc-hyperbolic-PDE" print(cmds) os.system(cmds) end_time = time.time() print(f"{end_time - start_time}s")Pythonで計算する場合
import time
import numpy as np
def _set_initial_condition( *, array_2d: list, grid_counts_x: int, grid_space: float,
) -> None: """ 初期条件の設定 境界以外のt=0のU値を設定する """ for i in range(1, grid_counts_x): array_2d[i][0] = np.sin(np.pi * i * grid_space)
def _set_boundary_condition( *, array_2d: list, grid_counts_x: int, grid_counts_t: int,
) -> None: """境界条件の設定""" for j in range(grid_counts_t + 1): array_2d[0][j] = 0.0 array_2d[grid_counts_x][j] = 0.0
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, grid_space=H, ) # 境界条件設定 _set_boundary_condition( array_2d=U, grid_counts_x=M, grid_counts_t=N, ) # U[i][1]の計算 for i in range(1, M): R = (K / H) ** 2 * (i * H + 1) S = 2 * (1.0 - R) Q = K**2 * i * H * np.exp(0) U[i][1] = (R * (U[i + 1][0] + U[i - 1][0]) + S * U[i][0] + Q) / 2.0 # U[i][j + 1]の計算 for j in range(1, N): for i in range(1, M): R = (K / H) ** 2 * (i * H + 1) S = 2 * (1.0 - R) Q = K**2 * i * H * np.exp(-j * K) U[i][j + 1] = ( R * (U[i + 1][j] + U[i - 1][j]) + S * U[i][j] - U[i][j - 1] + Q ) return U
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(R: float) -> None: """パラメータの妥当性をチェックする""" if R > 1.0: raise ValueError("R must be less than 1.0.")
if __name__ == "__main__": """ 双曲形方程式の数値解法（陽解法） """ grid_counts_x: int = 1000 # 格子点番号mの値 grid_counts_t: int = 4000 # 格子点番号nの値 grid_space: float = 0.001 # 刻み幅 (meter) time_delta: float = 0.0005 # 刻み幅 (hour) start_time = time.time() validate_parameters(R=2 * (time_delta / grid_space)) array_2d = calculate_equation( M=grid_counts_x, N=grid_counts_t, H=grid_space, K=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, ) end_time = time.time() print(f"{end_time - start_time}s")筆者の環境での（あくまで筆者の環境）実行時間は、
- Fortranを使った場合：1.701秒
- Pythonだけでやった場合：7.902秒
となりました。５回ずつ実行し、平均を取りました。5倍ぐらいFortranの方が速いことがわかります。
プログラムの解説
上に書いたソースプログラム calc-hyperbolic-PDE.f90 の内容を解説します。
program main implicit none ! declare variables ! constants real(8), parameter :: pi=4.0d0*datan(1.0d0) ! parameters for calculation integer(4) :: M, N real(8) :: H, K ! variables for I/O file character(len=256) :: outfile integer(4) :: lnml=10, lout=61 integer(4) :: iost ! variables for calculation result real(8), allocatable :: U(:, :) ! working variables integer(4) :: i, j real(8) :: R, S, Q ! setting for namelist file namelist /input/ M, N, H, K変数の宣言と、namelistファイルの宣言をしています。配列Uはこの後でサイズが決まりますので、allocatableをつけておきます。
! read the namelist file open (lnml, file="input.nmlst", action="read", iostat=iost) if (iost /= 0) then write (*, *) "ERROR: cannot open namelist file!" stop end if read (lnml, nml=input) close (lnml)namelistファイルに書かれた変数の値を読み込みます。
! allocate and initialize allocate(U(M + 1, N + 1)) U = 0.0d0namelistファイルに書かれた変数の値から配列Uのサイズが決まりましたので、メモリを確保し、配列の値を初期化します。
! varidate the parameter R R = 2.0d0 * (K / H) if (R > 1.0d0) then write(*, *) "Error: R must be less than 1.0." stop end if ! set the initial condition do i = 2, M U(i, 1) = sin(pi*(i-1)*H) end do ! set the boundary conditions do j = 1, N+1 U( 1, j) = 0.0d0 U(M+1, j) = 0.0d0 end do ! compute the equation ! calculation of U[i, 2] do i = 2, M R = (K / H) ** 2 * ((i - 1) * H + 1) S = 2.0d0 * (1.0d0 - R) Q = K ** 2 * (i - 1) * H * exp(0.0d0) U(i, 2) = (R * (U(i + 1, 1) + U(i - 1, 1)) + S * U(i, 1) + Q) / 2.0d0 end do ! calculation of U[i, j + 1] do j = 2, N do i = 2, M R = (K / H) ** 2 * ((i - 1) * H + 1) S = 2.0d0 * (1.0d0 - R) Q = K ** 2 * (i - 1) * H * exp(-dble(j - 1) * K) U(i, j + 1) = R * (U(i + 1, j) + U(i - 1, j)) & & + S * U(i, j) - U(i, j - 1) + Q end do end doパラメータの値のチェックをし、OKであったら初期条件と境界条件の設定をし、時間発展の計算をしていきます。計算式については、こちらをご覧ください。単純なdoループで構成しました。
! save result outfile = "calculated_result.txt" open(unit=lout, file=outfile, status="replace") write(lout, "(a)") "# This file shows calculated result as below." write(lout, "(a)") "" write(lout, "(a)") "# Calculation Parameters." write(lout, "(a,i4)") "grid_counts_x: ", M write(lout, "(a,i4)") "grid_counts_t: ", N write(lout, "(a,es12.3)") "grid_space: ", H write(lout, "(a,es12.3)") "time_delta: ", K write(lout, "(a)") "" write(lout, "(a)") "# Calculated Matrix U." do j = 1, N + 1 do i = 1, M + 1 write(lout, "(f12.6, x)", advance="no") U(i, j) end do write(lout, *) end do close(lout)計算結果を計算パラメータと合わせてテキストファイルに出力します。
! deallocate memory deallocate(U) ! end of program stop
end program main最後に、配列Uのために確保したメモリを開放し、プログラムを終了します。
最後に
今回はFortranで双曲型の偏微分方程式の数値計算をしてみました。数値計算といっても、これまでに紹介したファイルの入出力と、簡単なdoループで実装することができます。
筆者の環境で比較すると、Fortranで計算したほうが5倍ほど計算が速いことがわかりました。
今回は単純なdoループで実装しましたが、大規模な計算になると効率の良い計算が必要になります。このあたりは、先人たちの計算ライブラリが公開されていてそれらを利用することで計算速度を向上させることができます。
次回は、行列の計算を例に、数値計算ライブラリの利用方法を紹介したいと思います。
