[Fortranのはじめ方] Part24: 数値積分・常微分方程式の解法

科学技術計算において、解析的に解けない積分や微分方程式を扱うことは非常に多いです。Fortranはこのような数値計算を得意としており、効率的なプログラムを作成できます。このステップでは、数値積分と常微分方程式の基本的な解法について学びましょう。

数値積分 ∫

数値積分は、定積分の値を近似的に求める手法です。解析的に積分できない関数や、実験データのように離散的な値しか得られない場合に用いられます。

代表的な手法

  • 区分求積法: 積分区間を短冊状に分割し、長方形の面積の和として近似します。最も単純な方法です。
  • 台形公式: 積分区間を分割し、各区間を台形で近似します。区分求積法よりも精度が高いことが多いです。
  • シンプソン則: 積分区間を分割し、各区間を二次関数で近似します。より高精度な結果が期待できます。

Fortranによる台形公式の実装例

簡単な例として、関数 f(x) = x^2 を区間 [0, 1] で台形公式を用いて数値積分するプログラムを見てみましょう。

<!DOCTYPE html>
<html>
<head> <title>Fortran Trapezoidal Rule Example</title>
</head>
<body>
<pre><code class="language-fortran">
program trapezoidal_rule implicit none integer, parameter :: dp = kind(1.0d0) ! 倍精度実数型 real(dp) :: a, b, integral_value integer :: n integer :: i real(dp) :: h, x, total_sum ! 積分区間と分割数 a = 0.0_dp b = 1.0_dp n = 1000 ! 分割数を増やすと精度が向上 ! 刻み幅の計算 h = (b - a) / real(n, dp) ! 関数値の合計を計算 (台形公式) total_sum = (func(a) + func(b)) / 2.0_dp do i = 1, n - 1 x = a + real(i, dp) * h total_sum = total_sum + func(x) end do ! 積分値の計算 integral_value = h * total_sum ! 結果の出力 write(*, '(a, f10.6)') 'Numerical integral: ', integral_value write(*, '(a, f10.6)') 'Analytical integral: ', 1.0_dp / 3.0_dp ! 解析解 (x^3 / 3)
contains function func(x) result(y) real(dp), intent(in) :: x real(dp) :: y y = x**2 ! 被積分関数 f(x) = x^2 end function func
end program trapezoidal_rule
</code></pre>
</body>
</html>

ポイント:

  • implicit none で暗黙の型宣言を無効化しています。
  • 倍精度実数型 (real(dp)) を使用して計算精度を確保しています。
  • contains 以下に関数を定義しています。
  • 分割数 n を大きくすると、より正確な値に近づきますが、計算時間が増加します。

常微分方程式 (ODE) の解法 dy/dx

常微分方程式は、未知関数とその導関数の関係を表す方程式です。物理現象、化学反応、生物の個体数変化など、様々な分野でモデル化に用いられます。

数値解法では、初期値問題（ある点での関数の値が与えられている）を扱うことが多いです。時刻や位置などを少しずつ進めながら、関数の値を逐次的に計算していきます。

代表的な手法

  • オイラー法: 最も単純な手法で、現在の点の傾きを使って次の点の値を線形に近似します。刻み幅を小さくしないと誤差が大きくなりやすいです。
  • ルンゲ＝クッタ法: オイラー法を改良し、より高精度な結果を得る手法です。特に4次のルンゲ＝クッタ法 (RK4) が広く使われています。

Fortranによるオイラー法の実装例

簡単な例として、dy/dx = y, y(0) = 1 という常微分方程式をオイラー法で解いてみましょう（解析解は y = exp(x)）。

<!DOCTYPE html>
<html>
<head> <title>Fortran Euler Method Example</title>
</head>
<body>
<pre><code class="language-fortran">
program euler_method implicit none integer, parameter :: dp = kind(1.0d0) real(dp) :: x, y, x_end, h integer :: i, n_steps ! 初期条件と計算範囲、刻み幅 x = 0.0_dp y = 1.0_dp x_end = 2.0_dp h = 0.01_dp ! 刻み幅 n_steps = nint((x_end - x) / h) ! ステップ数 write(*, '(a, 2(a, f8.4))') 'Initial: ', 'x=', x, ' y=', y ! オイラー法による逐次計算 do i = 1, n_steps y = y + h * dydx(x, y) ! 次のyを計算 x = x + h ! xを更新 if (mod(i, nint(0.1_dp/h)) == 0) then ! 0.1ごとに途中結果を表示 write(*, '(a, i5, 2(a, f8.4), a, f8.4)') 'Step', i, ': x=', x, ' y=', y, & ' (Exact=', exp(x), ')' end if end do
contains function dydx(x, y) result(deriv) real(dp), intent(in) :: x, y real(dp) :: deriv deriv = y ! 微分方程式 dy/dx = y end function dydx
end program euler_method
</code></pre>
</body>
</html>

注意:

オイラー法は単純ですが、誤差が累積しやすいです。より精度が必要な場合は、ルンゲ＝クッタ法などの高次手法や、後述するライブラリの利用を検討しましょう。刻み幅 h を小さくすると精度は向上しますが、計算コストが増加します。

数値計算ライブラリの活用

複雑な積分や微分方程式を解く場合、自分でアルゴリズムを実装するよりも、最適化され、信頼性の高い数値計算ライブラリを利用するのが一般的です。

代表的なライブラリ

ライブラリ名種類主な機能ライセンス
NAG Library商用数値積分、ODE、PDE、最適化、統計など広範な数値計算ルーチンを提供プロプライエタリ
IMSL (by Rogue Wave)商用NAGと同様、広範な数値計算・統計解析ルーチンを提供プロプライエタリ
QUADPACKオープンソース1次元数値積分 (多くの数値計算ライブラリの基礎)Public Domain (一部改変版は異なる場合あり)
ODEPACKオープンソース常微分方程式の初期値問題ソルバー集 (LSODAなど)Public Domain / BSD-like
SUNDIALSオープンソース大規模な常微分方程式、微分代数方程式 (DAE) システムのソルバー (CVODE, IDAなど)BSD
SciPy (Python経由)オープンソースPythonの科学技術計算ライブラリ。Fortranで書かれたライブラリ (QUADPACK, ODEPACKなど) を内部で使用。scipy.integrateモジュールで数値積分やODEソルバーを提供。BSD
GNU Scientific Library (GSL)オープンソースC言語ライブラリだが、Fortranからのインターフェースも存在する。数値積分、ODE、乱数生成など多数の機能。GPL

ライブラリ利用のメリット:

  • 信頼性・精度: 専門家によって開発・検証されており、安定した高精度な計算が期待できます。
  • 効率: アルゴリズムが高度に最適化されており、自分で実装するよりも高速な場合が多いです。
  • 開発時間の短縮: 実績のあるルーチンを利用することで、複雑なアルゴリズムの実装にかかる時間を節約できます。

多くのライブラリでは、特定の問題設定（スティッフ方程式など）に特化した効率的なソルバーも提供されています。

まとめ

今回は、Fortranを用いた数値積分と常微分方程式の基本的な解法について学びました。

  • 数値積分では台形公式やシンプソン則、常微分方程式ではオイラー法やルンゲ＝クッタ法などの基本的なアルゴリズムがあります。
  • 簡単な問題であれば自作することも可能ですが、実用的な計算では信頼性と効率の高い数値計算ライブラリを利用することが推奨されます。
  • Fortranから利用できるライブラリには、商用のNAGやIMSL、オープンソースのQUADPACK, ODEPACK, SUNDIALS, GSLなどがあります。

これらの知識を活用して、次のステップではファイルを使ったデータの可視化など、より実践的な応用に進んでいきましょう！

参考情報

コメントを残す

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

前の記事

Netcat (nc) 徹底解説：ネットワークの万能ツールを使…

次の記事

[Fortranのはじめ方] Part25: ファイルによるデ…