【Fortran】双曲型偏微分方程式の計算【数値計算】

はじめに

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.0d0

namelistファイルに書かれた変数の値から配列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ループで実装しましたが、大規模な計算になると効率の良い計算が必要になります。このあたりは、先人たちの計算ライブラリが公開されていてそれらを利用することで計算速度を向上させることができます。

次回は、行列の計算を例に、数値計算ライブラリの利用方法を紹介したいと思います。