[Fortranのはじめ方] Part27: シミュレーション（熱伝導、粒子運動など）

Fortranを使って物理現象をコンピュータ上で再現してみよう！

これまでのステップでFortranの基本から応用まで学んできましたね。いよいよ集大成として、Fortranが得意とする科学技術計算の中でも代表的な「シミュレーション」に挑戦してみましょう。ここでは、比較的シンプルな例として「熱伝導」と「粒子運動」のシミュレーションを取り上げます。

シミュレーションとは？
物理法則や数学モデルに基づいて、現実の現象をコンピュータ上で模擬実験することです。Fortranはその数値計算能力の高さから、古くから様々なシミュレーションに用いられてきました。

1. 熱伝導シミュレーション
- 考え方（1次元・陽解法）
- Fortranコード例 (1次元熱伝導 – 陽解法)
2. 粒子運動シミュレーション
- 考え方（オイラー法）
- Fortranコード例 (鉛直投げ上げ – オイラー法)
3. シミュレーションを行う上での考慮事項
4. 次のステップと参考情報

1. 熱伝導シミュレーション

物体の温度が時間とともにどのように変化するかをシミュレーションします。ここでは、1次元の棒における熱の伝わり方を考えてみましょう。簡単な方法として有限差分法を用います。これは、空間と時間を小さな区間に区切り（離散化）、各点での温度変化を計算していく手法です。

考え方（1次元・陽解法）

棒をいくつかの点（格子点）に分割し、各点の温度を配列で管理します。ある時刻 `t` における格子点 `i` の温度を `T(i, t)` とすると、次の時刻 `t+dt` における温度 `T(i, t+dt)` は、隣接する点の温度を使って以下のように近似的に計算できます（陽解法）。

T(i, t+dt) = T(i, t) + alpha * dt / dx**2 * (T(i+1, t) - 2*T(i, t) + T(i-1, t))

ここで、`alpha` は熱拡散率、`dt` は時間ステップ幅、`dx` は格子点の間隔です。これをすべての格子点について、時間を進めながら繰り返し計算します。

Fortranコード例 (1次元熱伝導 – 陽解法)

簡単な一次元熱伝導のシミュレーションコードの例です。棒の両端の温度は固定されているとします。

program heat_conduction_1d implicit none integer, parameter :: N = 100 ! 格子点の数 real, parameter :: L = 1.0 ! 棒の長さ real, parameter :: D = 0.01 ! 熱拡散率 (alpha) real, parameter :: T_end = 100.0 ! シミュレーション終了時間 real, parameter :: dx = L / real(N-1) ! 格子間隔 real, parameter :: dt = 0.1 * dx**2 / D ! 時間刻み (安定性条件から決定) integer, parameter :: nsteps = int(T_end / dt) ! 時間ステップ数 integer :: i, n real, dimension(0:N+1) :: T, T_new ! 温度配列 (境界条件用に+2) real :: time ! 初期条件: 中央部のみ高温、両端は0度 T = 0.0 T(N/2) = 100.0 ! --- 時間発展ループ --- time = 0.0 do n = 1, nsteps time = time + dt ! 内部点の温度更新 (陽解法) do i = 1, N T_new(i) = T(i) + D * dt / dx**2 * (T(i+1) - 2.0 * T(i) + T(i-1)) end do ! 境界条件 (両端0度固定) T_new(0) = 0.0 T_new(N+1) = 0.0 ! 温度配列の更新 T = T_new ! --- 定期的に結果を出力 (例: 100ステップごと) --- if (mod(n, 100) == 0) then write(*, '(F8.3, 5F10.3)') time, T(1), T(N/4), T(N/2), T(3*N/4), T(N) end if end do write(*,*) 'Simulation finished.'
end program heat_conduction_1d

注意点

陽解法では、時間刻み `dt` が小さくないと計算が不安定になることがあります（CFL条件）。上記のコードでは、安定性を保つために `dt` を `0.1 * dx**2 / D` と比較的小さく設定しています。陰解法など、より安定な手法もありますが、ここでは扱いません。

2. 粒子運動シミュレーション

質点（粒子）が力（例えば重力）を受けてどのように運動するかをシミュレーションします。ここでは、簡単な例として、重力下での投げ上げ運動を考えます。基本的な方法としてオイラー法を用います。これは、微小時間 `dt` 後の速度と位置を、現在の速度と加速度（力から計算）を使って近似的に求める方法です。

考え方（オイラー法）

時刻 `t` における粒子の位置を `x(t)`、速度を `v(t)`、加速度を `a(t)` とします。加速度 `a(t)` は、粒子にはたらく力 `F(t)` と質量 `m` から `a(t) = F(t) / m` で求まります。オイラー法では、微小時間 `dt` 後の速度と位置を以下のように近似します。

v(t+dt) = v(t) + a(t) * dt
x(t+dt) = x(t) + v(t) * dt

これを時間を進めながら繰り返し計算します。

Fortranコード例 (鉛直投げ上げ – オイラー法)

地面から初速度 `v0` で真上にボールを投げ上げるシミュレーションの例です。空気抵抗は無視します。

program particle_motion_1d implicit none real, parameter :: g = 9.8 ! 重力加速度 (m/s^2) real, parameter :: v0 = 20.0 ! 初速度 (m/s) real, parameter :: t_end = 5.0 ! シミュレーション終了時間 (s) real, parameter :: dt = 0.01 ! 時間刻み (s) integer, parameter :: nsteps = int(t_end / dt) ! 時間ステップ数 integer :: n real :: time, y, v ! 時刻, 高さ, 速度 ! 初期条件 time = 0.0 y = 0.0 v = v0 ! ヘッダー出力 write(*, '(A8, A15, A15)') 'Time(s)', 'Height(m)', 'Velocity(m/s)' write(*, '(F8.3, F15.6, F15.6)') time, y, v ! --- 時間発展ループ --- do n = 1, nsteps ! 加速度 (一定: -g) real :: a = -g ! 速度と位置の更新 (オイラー法) v = v + a * dt y = y + v * dt ! 更新前の速度 v(t) を使うことに注意 (より正確には v(t+dt)*dt や平均を使う方法もある) ! 時間更新 time = time + dt ! 地面より下に行ったら停止 (簡単な衝突判定) if (y < 0.0) then y = 0.0 v = 0.0 write(*, '(F8.3, F15.6, F15.6)') time, y, v write(*,*) 'Hit the ground!' exit ! ループを抜ける end if ! --- 定期的に結果を出力 (例: 10ステップごと) --- if (mod(n, 10) == 0) then write(*, '(F8.3, F15.6, F15.6)') time, y, v end if end do write(*,*) 'Simulation finished.'
end program particle_motion_1d

改善のヒント

精度向上: オイラー法は最も単純ですが、誤差が蓄積しやすいです。より精度の高いルンゲ＝クッタ法などを使うと、より現実に近い結果が得られます。
多次元へ: 2次元（放物運動）や3次元に拡張することも可能です。位置、速度、力をベクトルとして扱います。
相互作用: 複数の粒子間の力（万有引力、クーロン力など）を考慮すると、より複雑で興味深いシミュレーション（N体シミュレーション）が可能です。

3. シミュレーションを行う上での考慮事項

数値計算手法の選択: 問題に応じて適切なアルゴリズム（陽解法/陰解法、オイラー法/ルンゲ＝クッタ法など）を選ぶ必要があります。精度、安定性、計算コストを考慮します。
パラメータ設定: 時間刻み `dt` や格子間隔 `dx` は、計算の精度や安定性に大きく影響します。小さすぎると計算時間が膨大になり、大きすぎると誤差が大きくなったり不安定になったりします。
境界条件: シミュレーション領域の端で物理量がどう振る舞うかを設定します（例：温度固定、断熱、周期的境界など）。
可視化: シミュレーション結果は数値の羅列であることが多いです。Step 7で学んだGnuplotなどのツールを使ってグラフやアニメーションにすることで、結果を理解しやすくなります。
検証: 解析解（もしあれば）と比較したり、物理的な保存則（エネルギー保存則など）が成り立っているかを確認したりして、シミュレーションの妥当性を検証することが重要です。

4. 次のステップと参考情報

今回紹介したのは非常に簡単な例ですが、Fortranによるシミュレーションの基本的な流れは掴めたでしょうか？

発展: 熱伝導を2次元や3次元に拡張する、粒子の相互作用（衝突、引力など）を導入するなど、モデルを複雑にしてみましょう。
ライブラリ活用: より高度な数値計算（連立一次方程式の解法、高速フーリエ変換など）には、BLAS/LAPACKやFFTWといった最適化されたライブラリの利用も検討しましょう（Step 7参照）。
学習リソース:
- Fortran-lang 公式サイト: Fortranに関する最新情報やチュートリアルがあります。
  https://fortran-lang.org/
  https://fortran-lang.org/learn/ (学習リソース)
- 数値計算の教科書: 有限差分法、有限要素法、数値積分法など、シミュレーションに使われる様々な手法について詳しく解説されています。大学の講義資料なども参考になります。
- オンラインチュートリアル: “Fortran simulation tutorial”, “Finite Difference Method Fortran”, “Particle Simulation Euler method Fortran”などで検索すると、様々な解説記事やコード例が見つかります。

Fortranを使ったシミュレーションは奥が深く、様々な分野で活用されています。ぜひ興味のある現象を対象に、自分だけのシミュレーションプログラムを作成してみてください！