Juliaで常微分方程式を解きます。4段4次ルンゲクッタ編
はじめに
前回記事ではSciMLパッケージを用いて常微分方程式 を解くことを試みた。
これまでも何度か引用してきた1から始めるJuliaプログラミングをふたたび読み直していると、第4章の高速化の部分で4次ルンゲクッタ (Runge-Kutta) 法のコードが掲載されていたので、これを使って常微分方程式を解くことを考える。ただし、この関数は次の形の常微分方程式を対象としている：
function solve_ode(f, t0, y0, n, h)
t, y = t0, y0
for _ in 1:n
k1 = f(t , y )
k2 = f(t + h / 2, y + k1 / 2 * h)
k3 = f(t + h / 2, y + k2 / 2 * h)
k4 = f(t + h / 2, y + k3 * h)
t += h
y += h * (k1 + 2 * k2 + 2 * k3 + k4) / 6
end
return t, y
end
進藤 裕之、佐藤 建太著『1から始めるJuliaプログラミング』（コロナ社）, p.175より引用
この関数は関数fを引数に取るようになっていて、これが前述の常微分方程式の右辺に相当する。
しかしふと気づくとJuliaで関数を引数に取るやり方がわからなかった。
そこで何をやりたいかを頭の中で整理するのに、とりあえずFortranに書き起こしてみることにした。
Fortranで書いてみる
いにしえの記憶を頼りにしたらなんとなくできた。 特に文法を調べることなく書けたのでちょっとだけ嬉しい。
program main
implicit none
integer :: n
double precision :: s
n = 10
s = solve_ode(f, 0d0, 1d0, n, 1d0 / dble(n))
write(*,*) n, s, abs(s - exp(0.98d0))
n = 100
s = solve_ode(f, 0d0, 1d0, n, 1d0 / dble(n))
write(*,*) n, s, abs(s - exp(0.98d0))
n = 1000
s = solve_ode(f, 0d0, 1d0, n, 1d0 / dble(n))
write(*,*) n, s, abs(s - exp(0.98d0))
contains
double precision function f(t, y)
implicit none
double precision, intent(in) :: t, y
f = 0.98d0 * y
end function
double precision function solve_ode(f, t0, y0, n, h)
implicit none
double precision :: f
double precision :: t0, y0
integer :: n
double precision :: h
double precision :: t, y
integer :: i
double precision :: k1, k2, k3, k4
t = t0
y = y0
do i = 1, n
k1 = f(t, y)
k2 = f(t + h / 2d0, y + k1 / 2d0 * h)
k3 = f(t + h / 2d0, y + k2 / 2d0 * h)
k4 = f(t + h / 2d0, y + k3 * h)
t = t + h
y = y + h * (k1 + 2d0 * k2 + 2d0 * k3 + k4) / 6d0
end do
solve_ode = y
end function
end program
ここで関数fは解くべき微分方程式の右辺であり、double precisionで宣言している。関数solve_ode内で同じくdouble precisionとして宣言した仮引数fに渡すことができる。
すると宣言の通り関数fの型はdouble precisionなのだと思う。
このコードを実行して以下の出力を得た：
$ ./a.exe
10 2.6644543920994486 1.8498299683145092E-006
100 2.6644562417303468 1.9907009374264817E-010
1000 2.6644562419293938 2.3092638912203256E-014
先人の研究によれば、4段4次ルンゲクッタ法は5次精度ということで、確かに誤差はすごい勢いで小さくなる。
Juliaで書いてみる
Fortranでコードを書くことで何をやれば良いのかが理解できたので、改めてJuliaでコードを書いてみる。最終的に得られたコードは以下：
using Printf
function makef(c)
return f(t, y) = c * y
end
function solve_ode(f, t0, y0, n, h)
t, y = t0, y0
for _ in 1:n
k1 = f(t, y)
k2 = f(t + h / 2, y + k1 / 2 * h)
k3 = f(t + h / 2, y + k2 / 2 * h)
k4 = f(t + h / 2, y + k3 * h)
t += h
y += h * (k1 + 2 * k2 + 2 * k3 + k4) / 6
end
return t, y
end
function main()
f = makef(0.98)
n = 10
t, s = solve_ode(f, 0.0, 1.0, n, 1.0 / n)
@printf "%5i %23.15f %23.15e\n" n s (exp(0.98) - s)
n = 100
t, s = solve_ode(f, 0.0, 1.0, n, 1.0 / n)
@printf "%5i %23.15f %23.15e\n" n s (exp(0.98) - s)
n = 1000
t, s = solve_ode(f, 0.0, 1.0, n, 1.0 / n)
@printf "%5i %23.15f %23.15e\n" n s (exp(0.98) - s)
end
main()
- 関数
makefは係数cを引数に取って、関数（解くべき微分方程式の右辺）を返す関数。- 関数を返す関数はFortranでも頑張れば作れたように思うが、記憶があやしい。INTERFACE文とかを使うのだっただろうか…
- 関数
solve_odeは関数を引数に取る関数。解くべき微分方程式の右辺を渡してやる。- for文で添字にアンダースコアを使う流儀はこのあたりに由来するだろうか？ 変数として参照しようとしてもエラーになるので安全であり、単に反復回数を得たいだけの場合に便利。
typeof(solve_ode)はtypeof(solve_ode)を返し、supertype(typeof(solve_ode))はFunctionを返す。そういう型になる（Fortranからの類推でしばらく嵌ったポイント）solve_ode(f, t0, y0, n, h)::Float64というように関数名の後ろに型を指定できるが、これは関数の戻り値に対する型指定であって、関数そのものの型を指定するものではないことに注意。
- このように関数を引数に取る、または関数を返す関数を高階関数と呼んでいいのかな。高階関数の代表的なものに
mapがあるらしく、関数には匿名関数を使うことができる：たとえばmap(x -> x * x, array)など。 - 初期値と解析解はハードコーディングだが、解析解を用いても良いかもしれない：ふたたび匿名関数を用いて
t -> exp(0.98 * t)など。 - 出力のフォーマッティングには前回同様に
@printfを用いた。23.15とはかつてよく唱えた呪文のひとつ。
このコードを実行して以下の出力を得た。Fortranと同じ解が得られたので一安心。
$ julia solve_ode.jl 10 2.664454392099449 1.849829968314509e-06 100 2.664456241730347 1.990700937426482e-10 1000 2.664456241929394 2.309263891220326e-14
おわりに
常微分方程式を解くJuliaのサンプルコードを例にして、関数を引数に取る関数と関数を返す関数をそれぞれ作れるようになった。
今回はいにしえのFortranスキルを取っかかりにして問題を整理することができた。 Fortranは10年前ほど前にこのサイト（まだ残っていてめちゃくちゃ懐かしい。すごく参考になるのでおすすめです）を読みつつコード書きを繰り返し、文字通り "骨の髄まで" ある程度の練度を得たという感覚が残っている。 ある一つのテーマで練度を上げて、それを別のテーマで応用する、例えばプロゲーマーでいえばウメハラ氏にとってのゲームと麻雀であるとか、あるいはときど氏にとってのゲームと研究1の関係といったものがそうだけど、そうしたものは学びの理想形であると常々感じている。 それらと比較するのはおこがましいけれども、今回のFortranとJuliaの関係性、両者を行き来する体験はその一端を感じさせてくれるものであった。■
-
彼の著書である『東大卒プロゲーマー 論理は結局、情熱にかなわない (PHP新書)』を読むと、彼が一流のプロゲーマーであるだけでなく、一流の研究プレイヤーであることがよくわかる↩
