using DifferentialEquations, ParameterizedFunctions, ODE, ODEInterfaceDiffEq, LSODA
f(t) = 0.25*sin(t)^2
g = @ode_def RigidBody begin
dy1 = I₁*y2*y3
dy2 = I₂*y1*y3
dy3 = I₃*y1*y2 + f(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./10.^(6:13)
reltols = 1./10.^(3:10);
sol = solve(prob,Vern7(),abstol=1/10^14,reltol=1/10^14)
test_sol = TestSolution(sol)
using Plots; gr()
sol = solve(prob,Vern7())
plot(sol)
setups = [Dict(:alg=>DP5())
#Dict(:alg=>ode45()) # fails
Dict(:alg=>dopri5())
Dict(:alg=>Tsit5())
Dict(:alg=>Vern6())
]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;appxsol=test_sol,save_everystep=true,numruns=1000,maxiters=10000)
plot(wp)
The DifferentialEquations.jl algorithms once again pull ahead. This is the first benchmark we've ran where ode45
doesn't fail. However, it still doesn't do as well as Tsit5
. One reason why it does so well is that the maximum norm that ODE.jl uses (as opposed to the L2 norm of Sundials, DifferentialEquations, and ODEInterface) seems to do really well on this problem. dopri5
does surprisingly bad in this test.
setups = [Dict(:alg=>DP8())
#Dict(:alg=>ode78()) # fails
Dict(:alg=>Vern7())
Dict(:alg=>Vern8())
Dict(:alg=>dop853())
Dict(:alg=>Vern6())
]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;appxsol=test_sol,save_everystep=false,numruns=1000,maxiters=1000)
plot(wp)
setups = [Dict(:alg=>Vern7())
Dict(:alg=>Vern8())
Dict(:alg=>odex())
Dict(:alg=>CVODE_Adams())
Dict(:alg=>lsoda())
Dict(:alg=>ddeabm())
Dict(:alg=>ARKODE(Sundials.Explicit(),order=6))
]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;appxsol=test_sol,save_everystep=false,numruns=1000,maxiters=1000)
plot(wp)
Once again, the OrdinaryDiffEq.jl pull far ahead in terms of speed and accuracy.