Add Pleiades as a second-order ODE benchmark

benchmarks/pleiades.jmd

@@ -0,0 +1,131 @@
+# Pleiades benchmark
+
+```julia
+using LinearAlgebra, Statistics, InteractiveUtils
+using DiffEqDevTools, ParameterizedFunctions, SciMLBase, OrdinaryDiffEq, Sundials, Plots
+using ModelingToolkit
+using ProbNumDiffEq
+
+# Plotting theme
+theme(:dao;
+    markerstrokewidth=0.5,
+    legend=:outertopright,
+    bottom_margin=5Plots.mm,
+    size = (1000, 400),
+)
+```
+
+### Pleiades problem definition
+
+```julia
+# first-order ODE
+@fastmath function pleiades(du, u, p, t)
+    v = view(u, 1:7)   # x
+    w = view(u, 8:14)  # y
+    x = view(u, 15:21) # x′
+    y = view(u, 22:28) # y′
+    du[15:21] .= v
+    du[22:28] .= w
+    @inbounds @simd ivdep for i in 1:14
+        du[i] = zero(eltype(u))
+    end
+    @inbounds @simd ivdep for i in 1:7
+        @inbounds @simd ivdep for j in 1:7
+            if i != j
+                r = ((x[i] - x[j])^2 + (y[i] - y[j])^2)^(3 / 2)
+                du[i] += j * (x[j] - x[i]) / r
+                du[7+i] += j * (y[j] - y[i]) / r
+            end
+        end
+    end
+end
+x0 = [3.0, 3.0, -1.0, -3.0, 2.0, -2.0, 2.0]
+y0 = [3.0, -3.0, 2.0, 0, 0, -4.0, 4.0]
+dx0 = [0, 0, 0, 0, 0, 1.75, -1.5]
+dy0 = [0, 0, 0, -1.25, 1, 0, 0]
+u0 = [dx0; dy0; x0; y0]
+tspan = (0.0, 3.0)
+prob1 = ODEProblem(pleiades, u0, tspan)
+
+# second-order ODE
+function pleiades2(ddu, du, u, p, t)
+    x = view(u, 1:7)
+    y = view(u, 8:14)
+    for i in 1:14
+        ddu[i] = zero(eltype(u))
+    end
+    for i in 1:7, j in 1:7
+        if i != j
+            r = ((x[i] - x[j])^2 + (y[i] - y[j])^2)^(3 / 2)
+            ddu[i] += j * (x[j] - x[i]) / r
+            ddu[7+i] += j * (y[j] - y[i]) / r
+        end
+    end
+end
+u0 = [x0; y0]
+du0 = [dx0; dy0]
+prob2 = SecondOrderODEProblem(pleiades2, du0, u0, tspan)
+probs = [prob1, prob2]
+
+ref_sol1 = solve(prob1, Vern9(), abstol=1/10^14, reltol=1/10^14, dense=false)
+ref_sol2 = solve(prob2, Vern9(), abstol=1/10^14, reltol=1/10^14, dense=false)
+ref_sols = [ref_sol1, ref_sol2]
+
+plot(ref_sol1, idxs=[(14+i,21+i) for i in 1:7], title="Pleiades Solution", legend=false)
+scatter!(ref_sol1.u[end][15:21], ref_sol1.u[end][22:end], color=1:7)
+```
+
+## First-order ODE vs. second-order ODE
+```julia
+DENSE = false;
+SAVE_EVERYSTEP = false;
+
+_setups = [
+  "EK0(3) (1st order ODE)" => Dict(:alg => EK0(order=3, smooth=DENSE), :prob_choice => 1)
+  "EK0(5) (1st order ODE)" => Dict(:alg => EK0(order=5, smooth=DENSE), :prob_choice => 1)
+  "EK1(3) (1st order ODE)" => Dict(:alg => EK1(order=3, smooth=DENSE), :prob_choice => 1)
+  "EK1(5) (1st order ODE)" => Dict(:alg => EK1(order=5, smooth=DENSE), :prob_choice => 1)
+  "EK0(4) (2nd order ODE)" => Dict(:alg => EK0(order=4, smooth=DENSE), :prob_choice => 2)
+  "EK0(6) (2nd order ODE)" => Dict(:alg => EK0(order=6, smooth=DENSE), :prob_choice => 2)
+  "EK1(4) (2nd order ODE)" => Dict(:alg => EK1(order=4, smooth=DENSE), :prob_choice => 2)
+  "EK1(6) (2nd order ODE)" => Dict(:alg => EK1(order=6, smooth=DENSE), :prob_choice => 2)
+]
+
+labels = first.(_setups)
+setups = last.(_setups)
+
+abstols = 1.0 ./ 10.0 .^ (6:11)
+reltols = 1.0 ./ 10.0 .^ (3:8)
+
+wp = WorkPrecisionSet(
+    probs, abstols, reltols, setups;
+    names = labels,
+    #print_names = true,
+    appxsol = ref_sols,
+    dense = DENSE,
+    save_everystep = SAVE_EVERYSTEP,
+    numruns = 10,
+    maxiters = Int(1e7),
+    timeseries_errors = false,
+    verbose = false,
+)
+
+plot(wp, color=[1 1 2 2 3 3 4 4],
+     # xticks = 10.0 .^ (-16:1:5)
+)
+```
+
+
+## Conclusion
+
+- If the problem is a second-order ODE, _implement it as a second-order ODE_!
+  Just use `SecondOrderODEProblem`.
+
+
+## Appendix
+
+Computer information:
+
+```julia
+InteractiveUtils.versioninfo()
+```

benchmarks/runall.jl

@@ -4,6 +4,7 @@ FILES = [
     "lotkavolterra.jmd",
     "vanderpol.jmd",
     "rober.jmd",
+    "pleiades.jmd",
     "multi-language-wrappers.jmd",
 ]
 

docs/make.jl

@@ -50,6 +50,7 @@ makedocs(
             "Multi-Language Wrapper Benchmark" => "benchmarks/multi-language-wrappers.md",
             "Non-stiff ODE: Lotka-Volterra" => "benchmarks/lotkavolterra.md",
             "Stiff ODE: Van der Pol" => "benchmarks/vanderpol.md",
+            "Second order ODE: Pleiades" => "benchmarks/pleiades.md",
             "DAE: ROBER" => "benchmarks/rober.md",
         ],
         "Internals" => [

docs/src/benchmarks/pleiades.md

@@ -0,0 +1,161 @@
+# Pleiades benchmark
+
+```julia
+using LinearAlgebra, Statistics, InteractiveUtils
+using DiffEqDevTools, ParameterizedFunctions, SciMLBase, OrdinaryDiffEq, Sundials, Plots
+using ModelingToolkit
+using ProbNumDiffEq
+
+# Plotting theme
+theme(:dao;
+    markerstrokewidth=0.5,
+    legend=:outertopright,
+    bottom_margin=5Plots.mm,
+    size = (1000, 400),
+)
+```
+
+
+
+
+### Pleiades problem definition
+
+```julia
+# first-order ODE
+@fastmath function pleiades(du, u, p, t)
+    v = view(u, 1:7)   # x
+    w = view(u, 8:14)  # y
+    x = view(u, 15:21) # x′
+    y = view(u, 22:28) # y′
+    du[15:21] .= v
+    du[22:28] .= w
+    @inbounds @simd ivdep for i in 1:14
+        du[i] = zero(eltype(u))
+    end
+    @inbounds @simd ivdep for i in 1:7
+        @inbounds @simd ivdep for j in 1:7
+            if i != j
+                r = ((x[i] - x[j])^2 + (y[i] - y[j])^2)^(3 / 2)
+                du[i] += j * (x[j] - x[i]) / r
+                du[7+i] += j * (y[j] - y[i]) / r
+            end
+        end
+    end
+end
+x0 = [3.0, 3.0, -1.0, -3.0, 2.0, -2.0, 2.0]
+y0 = [3.0, -3.0, 2.0, 0, 0, -4.0, 4.0]
+dx0 = [0, 0, 0, 0, 0, 1.75, -1.5]
+dy0 = [0, 0, 0, -1.25, 1, 0, 0]
+u0 = [dx0; dy0; x0; y0]
+tspan = (0.0, 3.0)
+prob1 = ODEProblem(pleiades, u0, tspan)
+
+# second-order ODE
+function pleiades2(ddu, du, u, p, t)
+    x = view(u, 1:7)
+    y = view(u, 8:14)
+    for i in 1:14
+        ddu[i] = zero(eltype(u))
+    end
+    for i in 1:7, j in 1:7
+        if i != j
+            r = ((x[i] - x[j])^2 + (y[i] - y[j])^2)^(3 / 2)
+            ddu[i] += j * (x[j] - x[i]) / r
+            ddu[7+i] += j * (y[j] - y[i]) / r
+        end
+    end
+end
+u0 = [x0; y0]
+du0 = [dx0; dy0]
+prob2 = SecondOrderODEProblem(pleiades2, du0, u0, tspan)
+probs = [prob1, prob2]
+
+ref_sol1 = solve(prob1, Vern9(), abstol=1/10^14, reltol=1/10^14, dense=false)
+ref_sol2 = solve(prob2, Vern9(), abstol=1/10^14, reltol=1/10^14, dense=false)
+ref_sols = [ref_sol1, ref_sol2]
+
+plot(ref_sol1, idxs=[(14+i,21+i) for i in 1:7], title="Pleiades Solution", legend=false)
+scatter!(ref_sol1.u[end][15:21], ref_sol1.u[end][22:end], color=1:7)
+```
+
+![](figures/pleiades_2_1.svg)
+
+
+
+## First-order ODE vs. second-order ODE
+```julia
+DENSE = false;
+SAVE_EVERYSTEP = false;
+
+_setups = [
+  "EK0(3) (1st order ODE)" => Dict(:alg => EK0(order=3, smooth=DENSE), :prob_choice => 1)
+  "EK0(5) (1st order ODE)" => Dict(:alg => EK0(order=5, smooth=DENSE), :prob_choice => 1)
+  "EK1(3) (1st order ODE)" => Dict(:alg => EK1(order=3, smooth=DENSE), :prob_choice => 1)
+  "EK1(5) (1st order ODE)" => Dict(:alg => EK1(order=5, smooth=DENSE), :prob_choice => 1)
+  "EK0(4) (2nd order ODE)" => Dict(:alg => EK0(order=4, smooth=DENSE), :prob_choice => 2)
+  "EK0(6) (2nd order ODE)" => Dict(:alg => EK0(order=6, smooth=DENSE), :prob_choice => 2)
+  "EK1(4) (2nd order ODE)" => Dict(:alg => EK1(order=4, smooth=DENSE), :prob_choice => 2)
+  "EK1(6) (2nd order ODE)" => Dict(:alg => EK1(order=6, smooth=DENSE), :prob_choice => 2)
+]
+
+labels = first.(_setups)
+setups = last.(_setups)
+
+abstols = 1.0 ./ 10.0 .^ (6:11)
+reltols = 1.0 ./ 10.0 .^ (3:8)
+
+wp = WorkPrecisionSet(
+    probs, abstols, reltols, setups;
+    names = labels,
+    #print_names = true,
+    appxsol = ref_sols,
+    dense = DENSE,
+    save_everystep = SAVE_EVERYSTEP,
+    numruns = 10,
+    maxiters = Int(1e7),
+    timeseries_errors = false,
+    verbose = false,
+)
+
+plot(wp, color=[1 1 2 2 3 3 4 4],
+     # xticks = 10.0 .^ (-16:1:5)
+)
+```
+
+![](figures/pleiades_3_1.svg)
+
+
+
+
+## Conclusion
+
+- If the problem is a second-order ODE, _implement it as a second-order ODE_!
+  Just use `SecondOrderODEProblem`.
+
+
+## Appendix
+
+Computer information:
+
+```julia
+InteractiveUtils.versioninfo()
+```
+
+```
+Julia Version 1.9.3
+Commit bed2cd540a1 (2023-08-24 14:43 UTC)
+Build Info:
+  Official https://julialang.org/ release
+Platform Info:
+  OS: Linux (x86_64-linux-gnu)
+  CPU: 12 × Intel(R) Core(TM) i7-6800K CPU @ 3.40GHz
+  WORD_SIZE: 64
+  LIBM: libopenlibm
+  LLVM: libLLVM-14.0.6 (ORCJIT, broadwell)
+  Threads: 12 on 12 virtual cores
+Environment:
+  JULIA_NUM_THREADS = auto
+  JULIA_STACKTRACE_MINIMAL = true
+```
+
+