Rigid Body Work-Precision Diagrams (including ESERK4)

using Pkg
Pkg.update()
using OrdinaryDiffEq, ParameterizedFunctions, Plots, DiffEqDevTools

k(t) = 0.25*sin(t)^2

g = @ode_def RigidBody begin
  dy1  = I₁*y2*y3
  dy2  = I₂*y1*y3
  dy3  = I₃*y1*y2 + k(t)
end I₁ I₂ I₃

p = [-2.0,1.25,-0.5]
prob = ODEProblem(g,[1.0;0.0;0.9],(0.0,10.0),p)

abstols = 1.0 ./ 10.0 .^ (6:13)
reltols = 1.0 ./ 10.0 .^ (3:10);
sol = solve(prob,Vern7(),abstol=1/10^14,reltol=1/10^14)
test_sol = TestSolution(sol)
using Plots; gr()
plot(sol,dpi=200,linewidth=1)

469.1s

Julia

Setups

setups = [
          Dict(:alg=>DP5()),
          Dict(:alg=>BS5()),
          Dict(:alg=>Tsit5()),
          Dict(:alg=>Vern6()),
          #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 = ["DP5" "BS5" "Tsit5" "Vern6" "ROCK2" "ROCK4" "SERK2 PI" "SERK2 Predictive" "ESERK4" "ESERK5"]

3.2s

Julia

1×10 Array{String,2}: "DP5" "BS5" "Tsit5" "Vern6" "ROCK2" "ROCK4" … "ESERK4" "ESERK5"

Speed Only Tests

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

897.2s

Julia

Interpolation

wp = WorkPrecisionSet(prob,abstols,reltols,setups;names=Names,appxsol=test_sol,save_everystep=false,numruns=100,maxiters=Int(1e8),error_estimate=:L2,dense_errors=true)
plot(wp,dpi=200,linewidth=1,legend=:topleft,legendfontsize=6)

769.5s

Julia

As expected, the algorithms are all pretty matched on time for this problem. However, you can clearly see the OrdinaryDiffEq.jl algorithms solving to a much higher accuracy and still faster, especially when the interpolations are involved.