HIRES Work-Precision Diagrams (including ESERK4)

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

f = @ode_def Hires begin
  dy1 = -1.71*y1 + 0.43*y2 + 8.32*y3 + 0.0007
  dy2 = 1.71*y1 - 8.75*y2
  dy3 = -10.03*y3 + 0.43*y4 + 0.035*y5
  dy4 = 8.32*y2 + 1.71*y3 - 1.12*y4
  dy5 = -1.745*y5 + 0.43*y6 + 0.43*y7
  dy6 = -280.0*y6*y8 + 0.69*y4 + 1.71*y5 -
           0.43*y6 + 0.69*y7
  dy7 = 280.0*y6*y8 - 1.81*y7
  dy8 = -280.0*y6*y8 + 1.81*y7
end

u0 = zeros(8)
u0[1] = 1
u0[8] = 0.0057
prob = ODEProblem(f,u0,(0.0,321.8122))

sol = solve(prob,Rodas5(),abstol=1/10^14,reltol=1/10^14)
test_sol = TestSolution(sol)
abstols = 1.0 ./ 10.0 .^ (4:11)
reltols = 1.0 ./ 10.0 .^ (1:8);

439.2s

Julia

plot(sol,dpi=200,linewidth=1,legend=:topright,legendfontsize=6)

40.2s

Julia

plot(sol,tspan=(0.0,5.0),dpi=200,linewidth=1,legend=:topright,legendfontsize=6)

2.9s

Julia

High Tolerances

This is the speed when you just want the answer.

abstols = 1.0 ./ 10.0 .^ (5:8)
reltols = 1.0 ./ 10.0 .^ (2:5);
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" "RKC" "ROCK2" "ROCK4" "SERK2 PI" "SERK2 Predictive" "ESERK5"]
Names = ["Rosen23" "Rodas4" "ROCK2" "ROCK4" "SERK2 PI" "ESERK4" "ESERK5"]

2.7s

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)

120.1s

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)

40.3s

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)

39.8s

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)

387.7s

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)

431.7s

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)

437.0s

Julia

Conclusion

We see a decent overall performance as compared to the implicit methods. RKC and SERK2 with Predictive control are omitted because the former takes can't complete within the given iterations and the latter goes unstable.