OREGO Work-Precision Diagrams (including ESERK4)

using Pkg
Pkg.update()
using OrdinaryDiffEq, DiffEqDevTools, ParameterizedFunctions, Plots
gr() #gr(fmt=:png)

f = @ode_def Orego begin
  dy1 = p1*(y2+y1*(1-p2*y1-y2))
  dy2 = (y3-(1+y1)*y2)/p1
  dy3 = p3*(y1-y3)
end p1 p2 p3

p = [77.27,8.375e-6,0.161]
prob = ODEProblem(f,[1.0,2.0,3.0],(0.0,30.0),p)
sol = solve(prob,Rodas5(),abstol=1/10^14,reltol=1/10^14)
test_sol = TestSolution(sol)

434.9s

Julia

retcode: Success Interpolation: 3rd order Hermite t: nothing u: nothing

High Tolerances

This is the speed when you just want the answer.

abstols = 1.0 ./ 10.0 .^ (5:8)
reltols = 1.0 ./ 10.0 .^ (3:6);
setups = [
          Dict(:alg=>Rosenbrock23()),
          Dict(:alg=>Rodas3()),
          #Dict(:alg=>RKC()),
          Dict(:alg=>ROCK2()),
          Dict(:alg=>ROCK4()),
          Dict(:alg=>SERK2(controller=:PI)),
          #Dict(:alg=>SERK2(controller=:Predictive)),
  				Dict(:alg=>ESERK4()),
          Dict(:alg=>ESERK5())
          ]
Names = ["Rosen23" "Rodas4" "ROCK2" "ROCK4" "SERK2 PI" "ESERK4" "ESERK5"]

3.2s

Julia

1×7 Array{String,2}: "Rosen23" "Rodas4" "ROCK2" "ROCK4" "SERK2 PI" "ESERK4" "ESERK5"

wp = WorkPrecisionSet(prob,abstols,reltols,setups; names=Names,
                      save_everystep=false,appxsol=test_sol,maxiters=Int(1e6),numruns=50)
plot(wp,dpi=200,linewidth=1,legend=:topright,legendfontsize=6)

300.3s

Julia

wp = WorkPrecisionSet(prob,abstols,reltols,setups;names=Names,dense = false,verbose=false,
                      appxsol=test_sol,maxiters=Int(1e6),error_estimate=:l2,numruns=50)
plot(wp,dpi=200,linewidth=1,legend=:topright,legendfontsize=6)

198.8s

Julia

wp = WorkPrecisionSet(prob,abstols,reltols,setups; names=Names,
                      appxsol=test_sol,maxiters=Int(1e6),error_estimate=:L2,numruns=50)
plot(wp,dpi=200,linewidth=1,legend=:topright,legendfontsize=6)

196.2s

Julia

Low Tolerances

This is the speed at lower tolerances, measuring what's good when accuracy is needed.

abstols = 1.0 ./ 10.0 .^ (7:13)
reltols = 1.0 ./ 10.0 .^ (4:10)
wp = WorkPrecisionSet(prob,abstols,reltols,setups; names=Names,
                      save_everystep=false,appxsol=test_sol,maxiters=Int(1e8),numruns=50)
plot(wp,dpi=200,linewidth=1,legend=:topright,legendfontsize=6)

842.5s

Julia

wp = WorkPrecisionSet(prob,abstols,reltols,setups;names=Names, verbose=false,
                      dense=false,appxsol=test_sol,maxiters=Int(1e8),error_estimate=:l2,numruns=50)
plot(wp,dpi=200,linewidth=1,legend=:topright,legendfontsize=6)

945.1s

Julia

wp = WorkPrecisionSet(prob,abstols,reltols,setups;names=Names, 
                      appxsol=test_sol,maxiters=Int(1e8),error_estimate=:L2,numruns=50)
plot(wp,dpi=200,linewidth=1,legend=:topright,legendfontsize=6)

913.3s

Julia

Conclusion

As we can see the stabilized methods are performing really well at low tolerances, considering that they are explicit methods. Also for high tolerances we have reasonable performance. RKC and SERK2 with Predictive control are omitted because the former takes can't complete within the given iterations and the latter goes unstable.