In [1]:
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()

Out[1]:
Plots.GRBackend()
In [2]:
sol = solve(prob,Vern7())
plot(sol)

Out[2]:
In [3]:
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)

Out[3]:

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.

## Higher Order¶

In [4]:
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)

Out[4]:
In [2]:
setups = [Dict(:alg=>Vern7())
Dict(:alg=>Vern8())
Dict(:alg=>odex())