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]:
0 1 2 3 -5 0 5 10 15 t u1(t) u2(t) u3(t) u4(t) u5(t) u6(t) u7(t) u8(t) u9(t) u10(t) u11(t) u12(t) u13(t) u14(t) u15(t) u16(t) u17(t) u18(t) u19(t) u20(t) u21(t) u22(t) u23(t) u24(t) u25(t) u26(t) u27(t) u28(t)

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]:
10 - 4 10 - 3 10 - 2 10 - 1 10 - 2.8 10 - 2.7 10 - 2.6 10 - 2.5 10 - 2.4 Error Time (s) OrdinaryDiffEq.DP5 ODEInterfaceDiffEq.dopri5 OrdinaryDiffEq.Tsit5 OrdinaryDiffEq.Vern6

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]:
10 - 4 10 - 3 10 - 2 10 - 1 10 - 2.8 10 - 2.7 10 - 2.6 10 - 2.5 10 - 2.4 Error Time (s) OrdinaryDiffEq.DP5 OrdinaryDiffEq.Tsit5 OrdinaryDiffEq.Vern6

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]:
10 - 4 10 - 3 10 - 2 10 - 1 10 0 10 - 2.7 10 - 2.6 10 - 2.5 10 - 2.4 Error Time (s) OrdinaryDiffEq.DP8 OrdinaryDiffEq.Vern7 OrdinaryDiffEq.Vern8 ODEInterfaceDiffEq.dop853 OrdinaryDiffEq.Vern6
In [4]:
setups = [Dict(:alg=>odex())
          Dict(:alg=>Vern7())
          Dict(:alg=>CVODE_Adams())
          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]:
10 - 4 10 - 3 10 - 2 10 - 1 10 0 10 1 10 - 2.7 10 - 2.6 10 - 2.5 10 - 2.4 10 - 2.3 Error Time (s) ODEInterfaceDiffEq.odex OrdinaryDiffEq.Vern7 Sundials.CVODE_Adams LSODA.lsoda OrdinaryDiffEq.Vern6 OrdinaryDiffEq.Tsit5 ODEInterfaceDiffEq.ddeabm Sundials.ARKODE

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]:
10 - 4 10 - 3 10 - 2 10 - 1 10 0 10 - 2.7 10 - 2.6 10 - 2.5 10 - 2.4 10 - 2.3 Error Time (s) OrdinaryDiffEq.DP8 OrdinaryDiffEq.Vern7 OrdinaryDiffEq.Vern8 OrdinaryDiffEq.Vern6

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.