数値計算入門【2. 微分方程式の数値解法】

  1. 事前準備
  2. 微分方程式の数値解法

※2020 年度は 1 章, 2 章ともに一通り行います. 節の冒頭, もしくはコマンドの実行を指示する箇所には赤字で「※2020 年度 : 」から始まる注釈をつけているので, よく確認し, 実行可能なコマンドは joho18 にて行ってください.
例年はこの資料に沿って情報実験機を使った実習をしています

[2.1] Euler 法

ここでは,Euler 法を用いた微分方程式の数値解法を Fortran プログラムで実装することを考えます.

[2.1.1] Euler 法で常微分方程式を解く

※2020 年度 : [2.1.1] 節は一通りコマンドを実行してください.

サンプルプログラム test1.f90 をダウンロードします.

$ wget www.ep.sci.hokudai.ac.jp/~inex/y2020/0717/practical/prog/test1.f90

test1.f90 は $\frac{df(z)}{dz}=-(z-a)f(z),\ a=3,\ f(0)=0.01$ という微分方程式を,$0\leq z\leq 10$の範囲で解くためのプログラムになっています.

test1.f90 の中身を見てみましょう.

$ less test1.f90

1   program main
2 
3       implicit none
4 
5   ! 定数・変数の宣言 ---------------------------------------------
6       double precision:: f_0 = 1.0d-2     ! f(z=0)の値
7       double precision:: z_max = 1.0d1    ! zの最大値
8       double precision:: dz = 1.0d-2      ! zの刻み幅
9       double precision:: a = 3.0d0        ! パラメータ
10 
11      double precision:: z, f             ! 変数
12  ! --------------------------------------------------------------
13
14      open(10, file='data_z.dat')         ! 出力ファイルを開く
15
16  ! 変数の初期値 -------------------------------------------------
17      f = f_0
18      z = 0.0d0
19      write(10,*) z, f                    ! 初期値の出力
20  ! --------------------------------------------------------------
21
22  ! ループ計算 (Euler法) -----------------------------------------
23      do while (z <= z_max)
24          z = z+dz
25          f = f-(z-a)*f*dz
26          write(10,*) z, f
27      end do
28  ! --------------------------------------------------------------
29
30      close(10)                           ! 出力ファイルを閉じる
31
32  end program main

このプログラムには,どこで何を行っているかを分かりやすくするためにコメント文を挿入しています.コメント文は! (エクスクラメーションマーク) から始まります.前回の hello.f90 と異なる点を大まかに記します:

次に,このプログラムの計算の詳細に迫ってみましょう.25 行目に Euler 法による計算が記述されています.いま解くべき方程式は, \[\frac{df(z)}{dz}=-(z-a)f(z)\] ですから,これを微分の定義式を思い出して,小さな値 $dz$ で $z$ に離散化を施せば, \[\frac{f(z+dz)-f(z)}{dz}=-(z-a)f(z)\iff f(z+dz)=f(z)-(z-a)f(z)dz\] となります.この右の式こそが Euler 法による表式であり,プログラムの 25 行目の式です.

それでは,このプログラムを GFortran でコンパイルし,実行してみましょう.

$ gfortran test1.f90
$ ./a.out

すると,データファイル (data_z.dat) が出力されます.
$ ls

a.out  data_z.dat  test1.f90

このデータファイルの中身を見てみましょう.
$ less data_z.dat

  0.0000000000000000        1.0000000000000000E-002
  1.0000000000000000E-002   1.0299000000000001E-002
  2.0000000000000000E-002   1.0605910200000001E-002
  2.9999999999999999E-002   1.0920905732940000E-002
  4.0000000000000001E-002   1.1244164542635025E-002
  5.0000000000000003E-002   1.1575867396642759E-002
  6.0000000000000005E-002   1.1916197898104056E-002
  7.0000000000000007E-002   1.2265342496518505E-002
  8.0000000000000002E-002   1.2623490497416845E-002
  8.9999999999999997E-002   1.2990834070891676E-002
  9.9999999999999992E-002   1.3367568258947535E-002
  0.10999999999999999       1.3753890981631120E-002
  0.11999999999999998       1.4150003041902096E-002
  …… 略

1 列目に$z$,2 列目に$f(z)$が出力されていることがわかります.ここで,"E-002" が$10^{-2}$を表すことに注意してください.

得られた解をグラフに表してみましょう.今年度はグラフ描画のために gnuplot というオープンソースのグラフ描画ソフトウェアを使います. gnuplot は数式からのグラフ描画とデータファイルからのグラフ描画に対応し, 様々な出力形式 (PNG, JPEG, GIF, PDF, etc...) をサポートしています (gnuplot について[2019 年度資料]). 今回はグラフ描画に必要な諸々の設定 (データファイルの指定, 軸の設定等) を, バッチファイルと呼ばれるスクリプトを実行することで指示します. gnuplot バッチファイルをダウンロードし,実行すると,微分方程式の解のグラフ (out1.pdf) が生成されます.

$ wget www.ep.sci.hokudai.ac.jp/~inex/y2020/0717/practical/prog/plot1.gpl
$ gnuplot plot1.gpl  

生成されたグラフ (out1.pdf) を見てみましょう.
※2020 年度 : local として各自の計算機を使ってください.
local$ ftp hoge@joho18.ep.sci.hokudai.ac.jp
  > get out1.pdf
  > exit

グラフの横軸は$z$,縦軸は$f(z)$です.この微分方程式の解析解 ($f(z)=f(0)\exp\left\{-\frac{z^2-2az}{2}\right\}$) のグラフは,$z=3$に頂上がある山のような形をしています(Fig. 1).out1.pdf に表示される数値解は,特徴が一致しているでしょうか?また一致しない場所があるとすれば,どこでしょうか?
Fig. 1:微分方程式$\frac{df(z)}{dz}=-(z-a)f(z),\ a=3,\ f(0)=0.01$の解析解.
計算が成功していれば,これによく似た数値解のグラフが表示されるはずです.

なお,plot1.gpl を vi で開き,コメントアウトのマーク (行頭の"#") を消して plot1.gpl を再度実行すると,数値解のグラフと合わせて解析解のグラフも表示され,数値解と解析解の比較ができます.余裕がある人はこれらがどの程度一致しているかを確かめてみましょう.このとき,ou1.pdf を開いたまま gnuplot を実行しても,ファイルが更新されません.out1.pdf を閉じた状態で gnuplot を実行するようにしてください.

[2.1.2] Euler 法で運動方程式を解く

※2020 年度 : [2.1.2] 節は一通りコマンドを実行してください.
ここまでは1階微分方程式を解く例を見てきましたが,自然界に現れる方程式の多くは運動方程式であり,これは基本的に2階微分方程式です.ここでは,運動方程式を解くサンプルプログラムを見てみます.

サンプルプログラム test2.f90 をダウンロードします.

$ wget www.ep.sci.hokudai.ac.jp/~inex/y2020/0717/practical/prog/test2.f90

このプログラムは,太陽周りの天体の運動を計算するためのプログラムで, \begin{align} &\frac{d^2 x}{dt^2}=-\frac{GM}{r^2}\frac{x}{r},\\ &\frac{d^2y}{dt^2}=-\frac{GM}{r^2}\frac{y}{r},\\ &r=\sqrt{x^2+y^2},\\ \end{align} という微分方程式を解いています.ここで,$x,y$は太陽を中心とした天体の座標位置,$r$は太陽からの距離,$G$は万有引力定数,$M$は太陽質量です.このプログラムでは,天体の初期位置は太陽から 1 AU 離れたところにあり,公転速度も地球とほぼ同じとしているため,軌道はほぼ地球と同じ円軌道となります.

ここで,test2.f90 の中身を見てみましょう.

$ less test2.f90

1   program main
2 
3       implicit none
4 
5   ! 定数・変数の宣言 ---------------------------------------------
6       double precision:: M = 1.9884D30          ! 太陽質量 [kg]
7       double precision:: G = 6.67408D-11        ! 万有引力定数
8 
9       double precision:: t_max = 2.0D0*3.1536D7 ! 時間の最大値 (2年)
10      double precision:: dt = 1.0D3             ! 時間の刻み幅 (1000sec)
11
12      double precision:: x, y, u, v, t, r       ! 変数
13      double precision:: uu, vv                 ! 変数
14  ! --------------------------------------------------------------
15
16      open(10, file='data_planet.dat', status='replace')  ! 出力ファイルを開く
17
18  ! 変数の初期値 -------------------------------------------------
19      t = 0.0D0
20      x = 1.496D11                              ! 1AU
21      y = 0.0D0
22      u = 0.0D0
23      v = 2.9788D4                              ! 29.8km/sec
24      write(10,*) t, x, y                       ! 初期値の出力
25  ! --------------------------------------------------------------
26
27  ! ループ計算 (Euler法) -----------------------------------------
28      do while (t <= t_max)
29          uu = u
30          vv = v
31
32          t = t+dt
33
34          r = sqrt(x*x+y*y)
35          u = u-G*M*x/(r*r*r)*dt
36          v = v-G*M*y/(r*r*r)*dt
37
38          x = x+uu*dt
39          y = y+vv*dt
40          write(10,*) t, x, y
41      end do
42  ! --------------------------------------------------------------
43
44      close(10)                                        ! 出力ファイルを閉じる
45
46  end program main

Euler 法をはじめとする数値積分法では,2階の微分方程式を連立した2つの1階微分方程式に直して解きます.例えば,今回の運動方程式では, \[\frac{d^2x}{dt^2}=-\frac{GM}{r^2}\frac{x}{r}\longrightarrow \frac{dx}{dt}=u,\ \frac{du}{dt}=-\frac{GM}{r^2}\frac{x}{r}\] というふうに,速度変数$u$を用いれば2つの1階微分方程式に直すことができます.この test2.f90 のループ計算部分でも,u (速度の$x$成分), v (速度の$y$成分) を導入して2階微分方程式を連立1階微分方程式に直し,それを Euler 法で計算するという手順が見て取れます.

それでは,このプログラムをコンパイル・実行してみましょう.

$ gfortran test2.f90
$ ./a.out

すると,データファイル (data_planet.dat) が出力されます.得られたデータを可視化するために,gnuplot バッチファイル (plot_planet_orbit.gpl) をダウンロードし,gnuplot を実行します.
$ wget www.ep.sci.hokudai.ac.jp/~inex/y2020/0717/practical/prog/plot_planet_orbit.gpl
$ gnuplot plot_planet_orbit.gpl

これにより,グラフ (out_planet.pdf) が出力されますので見てみましょう.
local$ ftp hoge@joho18.ep.sci.hokudai.ac.jp
  > get out1.pdf
  > exit

このグラフでは中心に太陽位置が描かれ,横軸は位置座標$x$,縦軸は位置座標$y$で,惑星軌道を見下ろすような視点で見ることになります.太陽の周りを円を描いて回る軌道は得られたでしょうか? また,数値解は厳密には円軌道を描かないはずです.円軌道との違いは何でしょうか? なお,vi 等で plot_planet_orbit.pdf を開き,コメントマーク部を消去して再度 gnuplot を実行すれば,解析解もグラフに表示されるようになります.解析解は円軌道を表現するため,円軌道との違いが明瞭になります.

[2.1.3] 計算パラメータの変更

※2020 年度 : [2.1.3] 節は一通りコマンドを実行してください.
test2.f90 の初期条件や変数の値を書き換えると,様々な天体の軌道の計算ができます.例題1では,惑星探査機を例に,初期条件を変えて様々な軌道を計算してみましょう.また,Euler法をはじめとする数値積分では,計算で用いる変数の値を適切に設定しないと,思うような計算精度が得られないことがあります.例題2では,計算で用いる変数の値を変更し,数値解が正しく解を表現できない例を見てもらいます.

例題1:惑星探査機の軌道計算

惑星探査機の軌道を計算してみましょう.地球から打ち上げることを考え,初期位置は元のままに,初速度を変更することにします.この計算は,惑星たちからの重力を無視し,太陽の重力のみを考えて探査機の軌道を導くものになります(このような計算設定は,惑星探査機で要求される推進剤質量を概算する際によく用いられます).ここで探査機の初速度を地球の速度 ( = 約29.8 km/s,元の速度) から数 % 上下させるとき,探査機が火星軌道/金星軌道と交差するのは何 % 変えたときでしょうか?実際に 23 行目で v の値を書き換えて確かめてみましょう.

探査機の軌道解析用のプログラムを用意し,これを explorer.f90 とします.

$ cp test2.f90 explorer.f90

精度を上げるため,explorer.f90 の 10 行目で,時間の刻み幅 dt を,100 (=1.0D2) 秒と書き換えます.
double precision:: dt = 1.0D2

まず,速度を 3% 上昇させてみましょう.explorer.f90 の 23 行目の v の欄を以下のようにします.
v = 1.03D0*2.9788D4

これは,元の速度を 1.03 倍するという意味です.変更したら,コンパイルし,実行ファイルを実行します.
$ gfortran explorer.f90
$ ./a.out  

グラフ描画のバッチファイルは,新しいファイル (plot_planet_orbit2.gpl) をダウンロードして使いましょう.

$ wget www.ep.sci.hokudai.ac.jp/~inex/y2020/0717/practical/prog/plot_planet_orbit2.gpl
$ gnuplot plot_planet_orbit2.gpl

このバッチファイルを使ってグラフ (out_planet2.pdf) を出力すると,先ほどのグラフと同様の座標系で他の惑星の軌道も描かれるため,惑星軌道への到達の判定に便利です.

軌道交差が起こるためには±10 % 以下の速度変化で充分です.なので,全てを逐一確かめるのではなく,「アタリ」をつけて計算を実行しましょう.

例題2:Euler 法の計算誤差

Euler 法で運動方程式を計算した結果では,誤差が積み重なってエネルギーが厳密には保存しません.例えば今回の場合,エネルギーは単調に増加してしまい,計算上,地球がどんどん太陽から離れてしまいます.この影響は,時間刻みを大きく取ったときに顕著です.ここでは,時間刻みで誤差が積み重なる様を見てみましょう.

解析用のプログラムを用意し,これを err.f90 とします.

$ cp test2.f90 err.f90

誤差の影響を分かりやすくするため,err.f90 の9行目で,時間の最大値を10年に書き換えます.
double precision:: t_max = 1.0D1*3.1536D7

ここで,10行目を書き換え,時間の刻み幅を 1000秒 (1.0D3) から,1万秒 (1.0D4),10万秒 (1.0D5),100万秒 (1.0D6) とします.例えば,1万秒に書き換える場合は,変数 dt の定義行 (10行目) を以下のようにします.
double precision:: dt = 1.0D4

書き換えたら,計算実行・グラフ描画をしてみましょう.グラフ描画には,先ほどの plot_planet_orbit2.gpl を用いてください.
$ gfortran err.f90
$ ./a.out
$ gnuplot plot_planet_orbit2.gpl

例えば,地球が 10 年以内に火星軌道以遠まで弾き飛ばされてしまうような時間刻みはどれくらいでしょうか? Euler 法によるエネルギー増大を味わってみてください.また,すべて終わって時間が余った人は,付録にチャレンジしてみましょう.

[2.1.4] 応用課題 :静水圧平衡の式のプログラムを書いてみよう

※2020 年度 : この課題は必須ではありませんが, 取り組むことで評定が上がる可能性があります.

サンプルプログラムの「定数・変数の宣言」と「変数の初期値」,「ループ計算」の部分を適切に書き換えれば,色々な微分方程式を解くことができるようになります. test1.f90 を,大気圧の鉛直プロファイルが満たす微分方程式(静水圧平衡の式): \[\frac{dp(z)}{dz}=-\frac{mg}{kT}p(z)\] を解くためのプログラムへと書き換えてみましょう.ここで,$p(z)$は高度$z$における大気圧,$m$は大気分子の平均質量,$g$は重力加速度,$k$はBoltzmann定数,$T$は大気の温度[K]です.解く際には,高度範囲は$0\sim 10^4\,\mathrm{m}$,高度刻み ($dz$) は$1000\,\mathrm{m}$としましょう.なお,書き換えに当たっては以下の定数を用いてください.

静水圧平衡の式を解くためには,どのように書き換えればいいでしょうか?必要な変数や定数を宣言し,Euler法を思い出して微分方程式を離散化してみましょう.

書き換えが終わったら,プログラムを実行して,データファイル (data_z.dat) を生成してみましょう.

$ gfortran test1.f90
$ ./a.out

データファイルを生成したら,gnuplot バッチファイル (plot_atm_profile.gpl) をダウンロードしてグラフ (out_atm_profile.pdf) を描き,それを表示しましょう.
$ wget www.ep.sci.hokudai.ac.jp/~inex/y2020/0717/practical/prog/plot_atm_profile.gpl
$ gnuplot plot_atm_profile.gpl
local$ sftp hoge@joho18.ep.sci.hokudai.ac.jp
  > get out1.pdf
  > exit

提出内容

以下に示す suu のレポートタグに所定のものを投稿してください.

情報実験第 9 回のページへ戻る

最終更新日: 2020/07/17 関口 太郎 修正
Copyright © 2000-2020 inex