Skip to content

Commit 616b04f

Browse files
nathanaelboschnathanaelbosch
authored
Add Pleiades as a second-order ODE benchmark (#250)
1 parent 4fa8fd6 commit 616b04f

36 files changed

+8471
-5085
lines changed

‎benchmarks/lotkavolterra.jmd‎

Lines changed: 14 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -4,7 +4,7 @@ Adapted from
44
[SciMLBenchmarks.jl](https://docs.sciml.ai/SciMLBenchmarksOutput/stable/NonStiffODE/LotkaVolterra_wpd/).
55

66
```julia
7-
using LinearAlgebra, Statistics, InteractiveUtils
7+
using LinearAlgebra, Statistics
88
using DiffEqDevTools, ParameterizedFunctions, SciMLBase, OrdinaryDiffEq, Plots
99
using ProbNumDiffEq
1010

@@ -215,6 +215,7 @@ wp = WorkPrecisionSet(
215215
plot(wp, color=[2 2 2 4 4 4 5 5 5], xticks = 10.0 .^ (-16:1:5))
216216
```
217217

218+
218219
## Conclusion
219220

220221
- **Use the EK1!** It seems to be strictly better than the EK0 here.
@@ -224,11 +225,21 @@ plot(wp, color=[2 2 2 4 4 4 5 5 5], xticks = 10.0 .^ (-16:1:5))
224225
- Most likely, the default choice of `diffusionmodel=DynamicDiffusion` and `initialization=TaylorModeInit` are fine.
225226

226227

227-
228228
## Appendix
229229

230230
Computer information:
231-
232231
```julia
232+
using InteractiveUtils
233233
InteractiveUtils.versioninfo()
234234
```
235+
236+
Package Information:
237+
```julia
238+
using Pkg
239+
Pkg.status()
240+
```
241+
242+
And the full manifest:
243+
```julia
244+
Pkg.status(mode=Pkg.PKGMODE_MANIFEST)
245+
```

‎benchmarks/multi-language-wrappers.jmd‎

Lines changed: 13 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -5,7 +5,7 @@ Adapted from
55

66
```julia
77
# Imports
8-
using LinearAlgebra, Statistics, InteractiveUtils
8+
using LinearAlgebra, Statistics
99
using StaticArrays, DiffEqDevTools, ParameterizedFunctions, Plots, SciMLBase, OrdinaryDiffEq
1010
using ODEInterface, ODEInterfaceDiffEq, Sundials, SciPyDiffEq, deSolveDiffEq, MATLABDiffEq, LSODA
1111
using LoggingExtras
@@ -330,7 +330,18 @@ plot(
330330
## Appendix
331331

332332
Computer information:
333-
334333
```julia
334+
using InteractiveUtils
335335
InteractiveUtils.versioninfo()
336336
```
337+
338+
Package Information:
339+
```julia
340+
using Pkg
341+
Pkg.status()
342+
```
343+
344+
And the full manifest:
345+
```julia
346+
Pkg.status(mode=Pkg.PKGMODE_MANIFEST)
347+
```

‎benchmarks/pleiades.jmd‎

Lines changed: 141 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,141 @@
1+
# Pleiades benchmark
2+
3+
```julia
4+
using LinearAlgebra, Statistics
5+
using DiffEqDevTools, ParameterizedFunctions, SciMLBase, OrdinaryDiffEq, Sundials, Plots
6+
using ModelingToolkit
7+
using ProbNumDiffEq
8+
9+
# Plotting theme
10+
theme(:dao;
11+
markerstrokewidth=0.5,
12+
legend=:outertopright,
13+
bottom_margin=5Plots.mm,
14+
size = (1000, 400),
15+
)
16+
```
17+
18+
### Pleiades problem definition
19+
20+
```julia
21+
# first-order ODE
22+
@fastmath function pleiades(du, u, p, t)
23+
v = view(u, 1:7) # x
24+
w = view(u, 8:14) # y
25+
x = view(u, 15:21) # x′
26+
y = view(u, 22:28) # y′
27+
du[15:21] .= v
28+
du[22:28] .= w
29+
@inbounds @simd ivdep for i in 1:14
30+
du[i] = zero(eltype(u))
31+
end
32+
@inbounds @simd ivdep for i in 1:7
33+
@inbounds @simd ivdep for j in 1:7
34+
if i != j
35+
r = ((x[i] - x[j])^2 + (y[i] - y[j])^2)^(3 / 2)
36+
du[i] += j * (x[j] - x[i]) / r
37+
du[7+i] += j * (y[j] - y[i]) / r
38+
end
39+
end
40+
end
41+
end
42+
x0 = [3.0, 3.0, -1.0, -3.0, 2.0, -2.0, 2.0]
43+
y0 = [3.0, -3.0, 2.0, 0, 0, -4.0, 4.0]
44+
dx0 = [0, 0, 0, 0, 0, 1.75, -1.5]
45+
dy0 = [0, 0, 0, -1.25, 1, 0, 0]
46+
u0 = [dx0; dy0; x0; y0]
47+
tspan = (0.0, 3.0)
48+
prob1 = ODEProblem(pleiades, u0, tspan)
49+
50+
# second-order ODE
51+
function pleiades2(ddu, du, u, p, t)
52+
x = view(u, 1:7)
53+
y = view(u, 8:14)
54+
for i in 1:14
55+
ddu[i] = zero(eltype(u))
56+
end
57+
for i in 1:7, j in 1:7
58+
if i != j
59+
r = ((x[i] - x[j])^2 + (y[i] - y[j])^2)^(3 / 2)
60+
ddu[i] += j * (x[j] - x[i]) / r
61+
ddu[7+i] += j * (y[j] - y[i]) / r
62+
end
63+
end
64+
end
65+
u0 = [x0; y0]
66+
du0 = [dx0; dy0]
67+
prob2 = SecondOrderODEProblem(pleiades2, du0, u0, tspan)
68+
probs = [prob1, prob2]
69+
70+
ref_sol1 = solve(prob1, Vern9(), abstol=1/10^14, reltol=1/10^14, dense=false)
71+
ref_sol2 = solve(prob2, Vern9(), abstol=1/10^14, reltol=1/10^14, dense=false)
72+
ref_sols = [ref_sol1, ref_sol2]
73+
74+
plot(ref_sol1, idxs=[(14+i,21+i) for i in 1:7], title="Pleiades Solution", legend=false)
75+
scatter!(ref_sol1.u[end][15:21], ref_sol1.u[end][22:end], color=1:7)
76+
```
77+
78+
## First-order ODE vs. second-order ODE
79+
```julia
80+
DENSE = false;
81+
SAVE_EVERYSTEP = false;
82+
83+
_setups = [
84+
"EK0(3) (1st order ODE)" => Dict(:alg => EK0(order=3, smooth=DENSE), :prob_choice => 1)
85+
"EK0(5) (1st order ODE)" => Dict(:alg => EK0(order=5, smooth=DENSE), :prob_choice => 1)
86+
"EK1(3) (1st order ODE)" => Dict(:alg => EK1(order=3, smooth=DENSE), :prob_choice => 1)
87+
"EK1(5) (1st order ODE)" => Dict(:alg => EK1(order=5, smooth=DENSE), :prob_choice => 1)
88+
"EK0(4) (2nd order ODE)" => Dict(:alg => EK0(order=4, smooth=DENSE), :prob_choice => 2)
89+
"EK0(6) (2nd order ODE)" => Dict(:alg => EK0(order=6, smooth=DENSE), :prob_choice => 2)
90+
"EK1(4) (2nd order ODE)" => Dict(:alg => EK1(order=4, smooth=DENSE), :prob_choice => 2)
91+
"EK1(6) (2nd order ODE)" => Dict(:alg => EK1(order=6, smooth=DENSE), :prob_choice => 2)
92+
]
93+
94+
labels = first.(_setups)
95+
setups = last.(_setups)
96+
97+
abstols = 1.0 ./ 10.0 .^ (6:11)
98+
reltols = 1.0 ./ 10.0 .^ (3:8)
99+
100+
wp = WorkPrecisionSet(
101+
probs, abstols, reltols, setups;
102+
names = labels,
103+
#print_names = true,
104+
appxsol = ref_sols,
105+
dense = DENSE,
106+
save_everystep = SAVE_EVERYSTEP,
107+
numruns = 10,
108+
maxiters = Int(1e7),
109+
timeseries_errors = false,
110+
verbose = false,
111+
)
112+
113+
plot(wp, color=[1 1 2 2 3 3 4 4],
114+
# xticks = 10.0 .^ (-16:1:5)
115+
)
116+
```
117+
118+
119+
## Conclusion
120+
121+
- If the problem is a second-order ODE, _implement it as a second-order ODE_!
122+
123+
124+
## Appendix
125+
126+
Computer information:
127+
```julia
128+
using InteractiveUtils
129+
InteractiveUtils.versioninfo()
130+
```
131+
132+
Package Information:
133+
```julia
134+
using Pkg
135+
Pkg.status()
136+
```
137+
138+
And the full manifest:
139+
```julia
140+
Pkg.status(mode=Pkg.PKGMODE_MANIFEST)
141+
```

‎benchmarks/rober.jmd‎

Lines changed: 17 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -4,7 +4,7 @@ Adapted from
44
[SciMLBenchmarks.jl](https://docs.sciml.ai/SciMLBenchmarksOutput/stable/DAE/ROBERDAE/).
55

66
```julia
7-
using LinearAlgebra, Statistics, InteractiveUtils
7+
using LinearAlgebra, Statistics
88
using DiffEqDevTools, ParameterizedFunctions, SciMLBase, OrdinaryDiffEq, Sundials, Plots
99
using ModelingToolkit
1010
using ProbNumDiffEq
@@ -35,7 +35,7 @@ daeprob = DAEProblem(sys,[D(y₁)=>-0.04, D(y₂)=>0.04, D(y₃)=>0.0],[],(0.0,1
3535
odaeprob = ODAEProblem(structural_simplify(sys),[],(0.0,1e5)) # can't handle this yet
3636

3737
ref_sol = solve(daeprob,IDA(),abstol=1/10^14,reltol=1/10^14,dense=false)
38-
plot(ref_sol, vars=[y₁,y₂,y₃], title="ROBER Solution", legend=false, ylims=(0, 1))
38+
plot(ref_sol, idxs=[y₁,y₂,y₃], title="ROBER Solution", legend=false, ylims=(0, 1))
3939
```
4040

4141
## EK1 accross orders
@@ -46,16 +46,12 @@ SAVE_EVERYSTEP = false;
4646

4747
_setups = [
4848
"EK1($order)" => Dict(:alg => EK1(order=order, smooth=DENSE))
49-
for order in 2:6
49+
for order in 1:4
5050
]
5151

5252
labels = first.(_setups)
5353
setups = last.(_setups)
5454

55-
# works:
56-
# abstols = 1.0 ./ 10.0 .^ (4:10)
57-
# reltols = 1.0 ./ 10.0 .^ (1:7)
58-
# test:
5955
abstols = 1.0 ./ 10.0 .^ (4:9)
6056
reltols = 1.0 ./ 10.0 .^ (1:6)
6157

@@ -72,8 +68,7 @@ wp = WorkPrecisionSet(
7268
verbose = false,
7369
)
7470

75-
plot(wp, palette=Plots.palette([:blue, :red], length(_setups)), xticks = 10.0 .^ (-16:1:5),
76-
xlims = (2e-15, 3e-7), ylims = (1e-2, 6e-1))
71+
plot(wp, palette=Plots.palette([:blue, :red], length(_setups)), xticks = 10.0 .^ (-16:1:5))
7772
```
7873

7974

@@ -82,10 +77,22 @@ plot(wp, palette=Plots.palette([:blue, :red], length(_setups)), xticks = 10.0 .^
8277
- The `EK1` can solve mass-matrix DAEs! But it only really works well for low errors.
8378
- Order 3 seems to work well here. But the order-to-error-tolerance heuristic should in principle still hold: lower tolerance level ``\rightarrow`` higher order.
8479

80+
8581
## Appendix
8682

8783
Computer information:
88-
8984
```julia
85+
using InteractiveUtils
9086
InteractiveUtils.versioninfo()
9187
```
88+
89+
Package Information:
90+
```julia
91+
using Pkg
92+
Pkg.status()
93+
```
94+
95+
And the full manifest:
96+
```julia
97+
Pkg.status(mode=Pkg.PKGMODE_MANIFEST)
98+
```

‎benchmarks/runall.jl‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -4,15 +4,17 @@ FILES = [
44
"lotkavolterra.jmd",
55
"vanderpol.jmd",
66
"rober.jmd",
7+
"pleiades.jmd",
78
"multi-language-wrappers.jmd",
89
]
910

11+
filedir = @__DIR__
1012
for file in FILES
1113
@info "Weave file" file
1214
weave(
1315
file;
1416
doctype="github",
15-
out_path="../docs/src/benchmarks/",
17+
out_path=joinpath(filedir, "../docs/src/benchmarks/"),
1618
fig_ext=".svg",
1719
)
1820
end

‎benchmarks/vanderpol.jmd‎

Lines changed: 15 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -2,7 +2,7 @@
22

33

44
```julia
5-
using LinearAlgebra, Statistics, InteractiveUtils
5+
using LinearAlgebra, Statistics
66
using DiffEqDevTools, ParameterizedFunctions, SciMLBase, OrdinaryDiffEq, Plots
77
using ProbNumDiffEq
88

@@ -117,15 +117,28 @@ wp = WorkPrecisionSet(
117117
plot(wp, color=[1 1 1 1 2 2 2 2], xticks = 10.0 .^ (-16:1:5))
118118
```
119119

120+
120121
## Conclusion
121122

122123
- Use the `EK1` to solve stiff problems, with orders $\leq 6$ depending on the error tolerance.
123124
- When the problem is actually a second-order ODE, as is the case for the Van der Pol system here, _solve it as a second-order ODE_.
124125

126+
125127
## Appendix
126128

127129
Computer information:
128-
129130
```julia
131+
using InteractiveUtils
130132
InteractiveUtils.versioninfo()
131133
```
134+
135+
Package Information:
136+
```julia
137+
using Pkg
138+
Pkg.status()
139+
```
140+
141+
And the full manifest:
142+
```julia
143+
Pkg.status(mode=Pkg.PKGMODE_MANIFEST)
144+
```

‎docs/make.jl‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -50,6 +50,7 @@ makedocs(
5050
"Multi-Language Wrapper Benchmark" => "benchmarks/multi-language-wrappers.md",
5151
"Non-stiff ODE: Lotka-Volterra" => "benchmarks/lotkavolterra.md",
5252
"Stiff ODE: Van der Pol" => "benchmarks/vanderpol.md",
53+
"Second order ODE: Pleiades" => "benchmarks/pleiades.md",
5354
"DAE: ROBER" => "benchmarks/rober.md",
5455
],
5556
"Internals" => [

0 commit comments

Comments
 (0)