Google ColabのJulia言語でDifferentialEquationsパッケージを使って35段14次ルンゲクッタFeagen法（BigFloat使用）、オイラー法、4段4次のルンゲクッタ法、Tsit5、Dormand＆Princeの8次(全部Float64)でローレンツ方程式をアダプティブは切って固定刻み幅で計算して比較。

Google Colabでjuliaが試せるようになったのでたまに使っている。DifferentialEquationsパッケージの公式サイトを見ていたら

35段14次のFeagen法がデフォルトで使えると分かった。これは昔、PARI/GPで試していたもの。

 35段14次のルンゲクッタ法をPARI/GPに実装、ローレンツ方程式を計算し、通常の4次、そして8次、オイラー法と比較してみる。

簡単に使えそうなので同様なことをやってみる。使うのはEuler, RK4, デフォルトのTsit5, Dormand&Princeの8次のDP8、そしてFeagen14。

Google ColabにはDifferentialEquationsパッケージはインストールされてなかったのでまずこれ。

using Pkg

Pkg.add("DifferentialEquations")

これが

185 dependencies successfully precompiled in 1221 seconds. 460 already precompiled.

と、とんでもなく時間がかかった。毎回これをするのは現実的じゃないな…

それはさておきプログラムは簡単で、今回は比較するのにアダプティブは使わず、刻み幅0.01に固定していて、

using DifferentialEquations

using Plots

function lorenz!(du, u, p, t)

    du[1] = 10.0 * (u[2] - u[1])

    du[2] = u[1] * (28.0 - u[3]) - u[2]

    du[3] = u[1] * u[2] - (8.0/3.0) * u[3]

    nothing

end

tspan = (0.0, 100.0)

u0 = [1.0;0.0;0.0]

bigu0 = [big(1.0);big(0.0);big(0.0)]

euler = solve(ODEProblem(lorenz!, u0, tspan), Euler(), adaptive = false, dt = 0.01);

rk4 = solve(ODEProblem(lorenz!, u0, tspan), RK4(), adaptive = false, dt = 0.01);

tsit5 = solve(ODEProblem(lorenz!, u0, tspan), Tsit5(), adaptive = false, dt = 0.01);

dp8 = solve(ODEProblem(lorenz!, u0, tspan), DP8(), adaptive = false, dt = 0.01);

feagin14 = solve(ODEProblem(lorenz!, bigu0, tspan), Feagin14(), adaptive = false, dt = 0.01);

plot(euler.t, euler[1,:], label = "Euler", xlabel="t", ylabel="X",size = (1200, 800))

plot!(rk4.t, rk4[1,:], label = "RK4")

plot!(tsit5.t, tsit5[1,:], label = "Tsit5")

plot!(dp8.t,dp8[1,:], label = "DP8")

plot!(feagin14.t, feagin14[1,:], label = "Feagin14")

とするだけ。Feagen14はBigFloatを使っているが、初期値だけそうしておくと勝手に全部BigFloatで計算してくれて便利。

結果はこちら。

まあオイラー法はほぼ最初から全然ダメ…オイラー法だけ除くと25くらいからルンゲクッタ4次もずれてきている。

拡大するとこんな感じ。

ルンゲクッタ4次は25くらいから、Tsit5は30くらいから、Dormand⁻Prince8次は40くらいからずれてる。

まあ実際に使うときはTsit5でアダプティブありにしとけばいいのでは。

念のため、Feagin14でアダプティブあり、abstol=reltol=1e-40にした結果。

うーんやっぱり50くらいからずれ始めるのか、難しい…