julia(1.4.1)~数値積分の精度確認~

概要

前回の記事（以下）で常微分方程式df/dt = cos(t)の数値積分を行った．

Julia(1.4.1)~常微分方程式の計算~ - Tagebuch von Spargel

本記事ではsolverによる精度の比較を行う．

結果

解析解との誤差の平均，最大誤差，実行時間は下表のようになる． 計算過程は後述．

solver	reltol	abstol	max error	mean error	Benchmark time [us]	対応するmatlab solver
default	-	-	4.38e-04	9.77e-06	116.5	-
Tsit5()	-	-	4.38e-04	9.77e-06	25.7	ode45
BS3()	-	-	2.40e-04	1.17e-04	29.8	ode23
Vern7()	-	-	1.63e-06	-2.68e-08	41.9	-
Vern9()	-	-	1.87e-07	5.91e-10	67.3	-
Tsit5()	1e-12	1e-12	1.89e-12	-1.41e-13	137.4	pde45
BS3()	1e-12	1e-12	1.69e-12	8.05e-13	7331.4	ode23
Vern7()	1e-12	1e-12	2.42e-14	-5.16e-15	72.1	-
Vern9()	1e-12	1e-12	3.75e-13	-2.89e-15	85.5	-

誤差の検証

環境

• windows10
• "julia" => v"1.4.1"
  • "StatsBase" => v"0.33.0"
  • "DifferentialEquations" => v"6.13.0"
  • "IJulia" => v"1.21.2"
  • "PyPlot" => v"2.9.0"
  • "BenchmarkTools" => v"0.5.0"

モジュールのインポート

下記のコードでモジュールをインポートする． インストールされてないモジュールはPkg.add("hoge")でインストールしておく．

using Pkg
# Pkg.add("StatsBase")
# Pkg.add("DifferentialEquations")
# Pkg.add("PyPlot")
# Pkg.add("BenchmarkTools")

#モジュールの導入
using DifferentialEquations
using PyPlot
using Statistics
using BenchmarkTools

数値積分

数値積分の定義を行う．詳細は前回の記事参照

#関数の定義
f(u,p,t) = cos(t)
#初期値の定義
u0 = 0.0;
#積分区間の指定
tspan = (0.0,4*pi);
prob = ODEProblem(f,u0,tspan);

許容誤差を指定せずに積分実行

solverだけ指定して数値積分を実行する， solverは独断でTsit5() BS3() Vern7() Vern9()を選択． 角solverの詳細は公式ドキュメント参照

ODE Solvers · DifferentialEquations.jl

#積分の実行,出力間隔は0.1
sol = solve(prob,saveat=0.1);#default
sol_Tsit5 = solve(prob,Tsit5(),saveat=0.1);
sol_BS3 = solve(prob,BS3(),saveat=0.1);
sol_Vern7 = solve(prob,Vern7(),saveat=0.1);
sol_Vern9 = solve(prob,Vern9(),saveat=0.1);

解析解はy=sin(t)となるため，解析解と数値積分の差をプロットする．

#解析解
sol_true = sin.(sol.t);
#誤差のplot
plot(sol.t,sol_true-sol.u,label="default",color = "k")
plot(sol_Tsit5.t,sol_true-sol_Tsit5.u,label="Tsit5()")
plot(sol_BS3.t,sol_true-sol_BS3.u,label="BS3()")
plot(sol_Vern7.t,sol_true-sol_Vern7.u,label="Vern7()")
plot(sol_Vern9.t,sol_true-sol_Vern9.u,label="Vern9()")
legend()
ylabel("ERROR")
ylim([-1e-6,1e-6])

defaultとTsit5は完全に重なっている（solverが指定されないときはTsit5で計算される模様）． 誤差が約4e-4あるのは大きい感じがする． また，Vern7とVern9は他と比べて小さすぎるため拡大したのが右図． Vern9のほうが精度はいいっぽい．

許容誤差を指定

reltolとabstolを指定する．とりあえず，1e-12を指定する．

#許容誤差を指定
sol_BS3_tol = solve(prob,BS3(),reltol = 1e-12,abstol=1e-12,saveat=0.1);
sol_Tsit5_tol = solve(prob,Tsit5(),reltol = 1e-12,abstol=1e-12,saveat=0.1);
sol_Vern7_tol = solve(prob,Vern7(),reltol = 1e-12,abstol=1e-12,saveat=0.1);
sol_Vern9_tol = solve(prob,Vern9(),reltol = 1e-12,abstol=1e-12,saveat=0.1);
#結果のプロット
plot(sol_Tsit5_tol.t,sol_true-sol_Tsit5_tol.u,label="Tsit5()")
plot(sol_BS3_tol.t,sol_true-sol_BS3_tol.u,label="BS3()")
plot(sol_Vern7_tol.t,sol_true-sol_Vern7_tol.u,label="Vern7()")
plot(sol_Vern9_tol.t,sol_true-sol_Vern9_tol.u,label="Vern9()")
legend()
ylabel("ERROR")

縦軸の数値を見るとわかるようにかなり誤差が減っていることが確認できる．

やっぱりreltolとabstolは必要．

実行時間の計測

@benchmarkを使用することで実行時間を取得できる．何回か関数を実行して，平均値，最大値，使用メモリ量などを出してくれる． 以下のコードでrecord.timesとすると時間が取得できる．

#benchmarkの実行 時間はかかる
record = @benchmark solve(prob,saveat=0.1);#default
record_Tsit5 = @benchmark solve(prob,Tsit5(),saveat=0.1);
record_BS3 = @benchmark solve(prob,BS3(),saveat=0.1);
record_Vern7 = @benchmark solve(prob,Vern7(),saveat=0.1);
record_Vern9 = @benchmark solve(prob,Vern9(),saveat=0.1);
record_BS3_tol = @benchmark solve(prob,BS3(),reltol = 1e-12,abstol=1e-12,saveat=0.1);
record_Tsit5_tol = @benchmark solve(prob,Tsit5(),reltol = 1e-12,abstol=1e-12,saveat=0.1);
record_Vern7_tol = @benchmark solve(prob,Vern7(),reltol = 1e-12,abstol=1e-12,saveat=0.1);
record_Vern9_tol = @benchmark solve(prob,Vern9(),reltol = 1e-12,abstol=1e-12,saveat=0.1);

結果は本記事の最初の表を参照． 許容誤差を指定すると実行時間がかかるみたい． また，defaultはTsit5と同じsolverなのに実行時間がかかっているので，solverは明示したほうが良い．

結論

solverはVern7()でabstol=1e-12,reltol=1e-12にするのが良いのではないかと思う．

最後に検証コードを記載する

検証コード