In [1]:
using DifferentialEquations, ODE, ODEInterfaceDiffEq, LSODA

f = (du,u,p,t) -> begin
x = view(u,1:7)   # x
y = view(u,8:14)  # y
v = view(u,15:21) # x′
w = view(u,22:28) # y′
du[1:7] .= v
du[8:14].= w
for i in 14:28
du[i] = zero(u[1])
end
for i=1:7,j=1:7
if i != j
r = ((x[i]-x[j])^2 + (y[i] - y[j])^2)^(3/2)
du[14+i] += j*(x[j] - x[i])/r
du[21+i] += j*(y[j] - y[i])/r
end
end
end

prob = ODEProblem(f,[3.0,3.0,-1.0,-3.0,2.0,-2.0,2.0,3.0,-3.0,2.0,0,0,-4.0,4.0,0,0,0,0,0,1.75,-1.5,0,0,0,-1.25,1,0,0],(0.0,3.0))

abstols = 1./10.^(6:9)
reltols = 1./10.^(3:6);

using Plots; gr()

Out[1]:
Plots.GRBackend()
In [3]:
sol = solve(prob,Vern8(),abstol=1/10^12,reltol=1/10^10,maxiters=1000000)
test_sol = TestSolution(sol);
plot(sol)

Out[3]:

## Low Order¶

ODE.jl had to be discarded. The error estimate is off since it throws errors and aborts and so that artificially lowers the error the the time is serverly diminished.

In [ ]:
setups = [Dict(:alg=>ode45())]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;appxsol=test_sol,save_everystep=false,numruns=100,maxiters=10000)
plot(wp)

In [4]:
setups = [Dict(:alg=>DP5())
Dict(:alg=>dopri5())
Dict(:alg=>Tsit5())
Dict(:alg=>Vern6())
]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;appxsol=test_sol,save_everystep=false,numruns=1000,maxiters=10000)
plot(wp)

Out[4]:

### Interpolation¶

In [5]:
setups = [Dict(:alg=>DP5())
Dict(:alg=>Tsit5())
Dict(:alg=>Vern6())
]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;appxsol=test_sol,numruns=1000,maxiters=10000,error_estimate=:L2,dense_errors=true)
plot(wp)

Out[5]:

## Higher Order¶

Once again ODE.jl had to be discarded since it errors.

In [ ]:
setups = [Dict(:alg=>ode78())]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;appxsol=test_sol,save_everystep=false,numruns=100,maxiters=1000)
plot(wp)

In [7]:
setups = [Dict(:alg=>DP8())
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=100,maxiters=1000)
plot(wp)

Out[7]:
In [4]:
setups = [Dict(:alg=>odex())
Dict(:alg=>Vern7())
Dict(:alg=>lsoda())
Dict(:alg=>Vern6())
Dict(:alg=>Tsit5())
Dict(:alg=>ddeabm())
Dict(:alg=>ARKODE(Sundials.Explicit(),order=6))
]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;appxsol=test_sol,save_everystep=false,numruns=20)
plot(wp)

Out[4]:

### Interpolations¶

In [8]:
setups = [Dict(:alg=>DP8())
Dict(:alg=>Vern7())
Dict(:alg=>Vern8())
Dict(:alg=>Vern6())
]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;appxsol=test_sol,numruns=1000,maxiters=1000,error_estimate=:L2,dense_errors=true)
plot(wp)

Out[8]:

## Conclusion¶

One big conclusion is that, once again, the ODE.jl algorithms fail to run on difficult problems. Its minimum timestep is essentially machine epsilon, and so this shows some fatal flaws in its timestepping algorithm. The OrdinaryDiffEq.jl algorithms come out as faster in each case than the ODEInterface algorithms. Overall, the Verner methods have a really good showing once again. The CVODE_Adams method does really well here when the tolerances are higher.