diff --git a/Project.toml b/Project.toml
index e37edd243..9bc3e4b46 100644
--- a/Project.toml
+++ b/Project.toml
@@ -32,6 +32,12 @@ TaylorIntegration = "92b13dbe-c966-51a2-8445-caca9f8a7d42"
TaylorSeries = "6aa5eb33-94cf-58f4-a9d0-e4b2c4fc25ea"
ToeplitzMatrices = "c751599d-da0a-543b-9d20-d0a503d91d24"
+[weakdeps]
+DiffEqDevTools = "f3b72e0c-5b89-59e1-b016-84e28bfd966d"
+
+[extensions]
+DiffEqDevToolsExt = "DiffEqDevTools"
+
[compat]
ArrayAllocators = "0.3"
DiffEqBase = "6.122"
@@ -59,6 +65,7 @@ Statistics = "1"
StructArrays = "0.4, 0.5, 0.6"
TaylorIntegration = "0.8, 0.9, 0.10, 0.11, 0.12, 0.13, 0.14"
TaylorSeries = "0.10, 0.11, 0.12, 0.13, 0.14, 0.15"
+Test = "1"
ToeplitzMatrices = "0.7, 0.8"
julia = "1.6"
diff --git a/benchmarks/Project.toml b/benchmarks/Project.toml
index c7d748b6c..4abd083bc 100644
--- a/benchmarks/Project.toml
+++ b/benchmarks/Project.toml
@@ -1,5 +1,6 @@
[deps]
DiffEqDevTools = "f3b72e0c-5b89-59e1-b016-84e28bfd966d"
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
IJulia = "7073ff75-c697-5162-941a-fcdaad2a7d2a"
LSODA = "7f56f5a3-f504-529b-bc02-0b1fe5e64312"
LoggingExtras = "e6f89c97-d47a-5376-807f-9c37f3926c36"
@@ -18,3 +19,6 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4"
Weave = "44d3d7a6-8a23-5bf8-98c5-b353f8df5ec9"
deSolveDiffEq = "0518478a-701b-4815-806c-24ad5cb92f09"
+
+[compat]
+DiffEqDevTools = "2.42.0"
diff --git a/benchmarks/hodgkinhuxley.jmd b/benchmarks/hodgkinhuxley.jmd
index 67951d590..d99f78d29 100644
--- a/benchmarks/hodgkinhuxley.jmd
+++ b/benchmarks/hodgkinhuxley.jmd
@@ -1,21 +1,29 @@
# Hodgkin-Huxley benchmark
-```julia
-using LinearAlgebra, Statistics
+
+!!! note "Summary"
+ Hodgkin-Huxley is a four-dimensional ODE, which can be stiff or non-stiff depending on the parameters;
+ here we consider a non-stiff version. We see that:
+ - [**`EK0` is the fastest solver**.](@ref hh_solver_comparison)
+ - [**`RosenbrockExpEK` is slowest; but suffers less from smoothing than `EK0` and `EK1`**.](@ref hh_solver_comparison)
+ - [Results are similar for fixed time steps.](@ref hh_fixed_steps)
+
+
+```julia, results="hidden"
+using LinearAlgebra, Statistics, Distributions
using DiffEqDevTools, SciMLBase, OrdinaryDiffEq, Plots, SimpleUnPack
using ProbNumDiffEq
-# Plotting theme
-theme(:dao;
+Plots.theme(
+ :dao;
markerstrokewidth=0.5,
legend=:outertopright,
- bottom_margin=5Plots.mm,
- size = (1000, 400),
+ margin=5Plots.mm,
+ xticks=10.0 .^ (-16:1:16),
+ yticks=10.0 .^ (-6:1:5),
)
```
-### Hodgkin-Huxley problem definition
-
```julia
αm(V, VT) = -0.32 * (V - VT - 13) / (exp(-(V - VT - 13) / 4) - 1)
βm(V, VT) = 0.28 * (V - VT - 40) / (exp((V - VT - 40) / 5) - 1)
@@ -26,14 +34,14 @@ theme(:dao;
αh(V, VT) = 0.128 * exp(-(V - VT - 17) / 18)
βh(V, VT) = 4 / (1 + exp(-(V - VT - 40) / 5))
-Inj(t) = (10 <= t <= 90) ? 500one(t) : zero(t)
+Inj(t) = (5 <= t <= 40) ? 500one(t) : zero(t)
function f(du, u, p, t)
@unpack gNa, gK, ENa, EK, area, C, Eleak, VT, gleak = p
V, m, n, h = u
- I_inj = Inj(t) * 1e-6 # uA
+ I_inj = Inj(t) * 1e-6
du[2] = dmdt = (αm(V, VT) * (1 - m) - βm(V, VT) * m)
du[3] = dndt = (αn(V, VT) * (1 - n) - βn(V, VT) * n)
@@ -53,9 +61,9 @@ n_inf(V, VT) = 1 / (1 + βn(V, VT) / αn(V, VT))
h_inf(V, VT) = 1 / (1 + βh(V, VT) / αh(V, VT))
u0 = [p.V0, m_inf(p.V0, p.VT), n_inf(p.V0, p.VT), h_inf(p.V0, p.VT)]
-prob = ODEProblem{true,SciMLBase.FullSpecialize()}(f, u0, (0.0, 100.0), p)
+prob = ODEProblem{true,SciMLBase.FullSpecialize()}(f, u0, (0.0, 50.0), p)
-test_sol = solve(prob, Vern7(), abstol=1/10^14, reltol=1/10^14, dense=false)
+test_sol = solve(prob, Vern7(), abstol=1/10^14, reltol=1/10^14)
plot(test_sol,
legend=false,
layout=(4,1),
@@ -64,10 +72,13 @@ plot(test_sol,
xlabel=["" "" "" "t"],
size = (1000, 600),
color=[1 2 3 4],
+ xticks=:auto, yticks=:auto
)
```
-## Adaptive steps - no smoothing
+## [Solver comparison: `EK0` vs. `EK1` vs `RosenbrockExpEK`](@id hh_solver_comparison)
+
+### Without smoothing
```julia
DENSE = SAVE_EVERYSTEP = false
@@ -92,26 +103,44 @@ reltols = 1.0 ./ 10.0 .^ (3:7)
wp = WorkPrecisionSet(
prob, abstols, reltols, setups;
names = labels,
- #print_names = true,
appxsol = test_sol,
dense = DENSE,
save_everystep = SAVE_EVERYSTEP,
- numruns = 10,
maxiters = Int(1e7),
- timeseries_errors = false,
- verbose = false,
+ numruns = 5,
)
-plot(
- wp,
- title = "Adaptive steps - no smoothing",
- color = colors,
- xticks = 10.0 .^ (-16:1:5),
- yticks = 10.0 .^ (-6:1:5),
+plot(wp, title="Adaptive steps - no smoothing", color=colors)
+
+
+_ref_setups = [
+ "Tsit5" => Dict(:alg => Tsit5())
+ "Vern7" => Dict(:alg => Vern7())
+ "RadauIIA5" => Dict(:alg => RadauIIA5())
+]
+ref_labels = first.(_ref_setups)
+ref_setups = last.(_ref_setups)
+ref_wp_final = WorkPrecisionSet(
+ prob, abstols ./ 1000, reltols ./ 1000, ref_setups;
+ names = ref_labels,
+ appxsol = test_sol,
+ dense = false,
+ save_everystep = false,
+ maxiters = Int(1e7),
)
+ref_wp_dense = WorkPrecisionSet(
+ prob, abstols ./ 1000, reltols ./ 1000, ref_setups;
+ names = ref_labels,
+ appxsol = test_sol,
+ dense = true,
+ save_everystep = true,
+ maxiters = Int(1e7),
+)
+
+plot!(ref_wp_final, x=:final, color=:gray, alpha=0.7, linestyle=:dash)
```
-## Adaptive steps - with smoothing
+### With smoothing
```julia
DENSE = SAVE_EVERYSTEP = true
@@ -136,28 +165,48 @@ reltols = 1.0 ./ 10.0 .^ (3:7)
wp = WorkPrecisionSet(
prob, abstols, reltols, setups;
names = labels,
- #print_names = true,
appxsol = test_sol,
dense = DENSE,
save_everystep = SAVE_EVERYSTEP,
- numruns = 10,
maxiters = Int(1e7),
- timeseries_errors = false,
- verbose = false,
+ numruns = 5,
)
-plot(
- wp,
- title = "Adaptive steps - with smoothing",
- color = colors,
- xticks = 10.0 .^ (-16:1:5),
- yticks = 10.0 .^ (-6:1:5),
-)
+plot(wp, title="Adaptive steps - with smoothing", color=colors)
+plot!(ref_wp_dense, x=:final, color=:gray, alpha=0.7, linestyle=:dash)
```
+```@raw html
+Interoplation errors (L2):
+```
+```julia
+plot(wp, x=:L2, title="Adaptive steps - with smoothing", color=colors)
+plot!(ref_wp_dense, x=:L2, color=:gray, alpha=0.7, linestyle=:dash)
+```
+```@raw html
+With smoothing:
+```
```julia
DENSE = SAVE_EVERYSTEP = true
@@ -222,41 +262,21 @@ wp = WorkPrecisionSet(
prob, abstols, reltols, setups;
adaptive = false,
names = labels,
- #print_names = true,
appxsol = test_sol,
dense = DENSE,
save_everystep = SAVE_EVERYSTEP,
- numruns = 10,
maxiters = Int(1e7),
- timeseries_errors = false,
- verbose = false,
-)
-
-plot(
- wp,
- title = "Fixed steps - with smoothing",
- color = colors,
- xticks = 10.0 .^ (-16:1:5),
- yticks = 10.0 .^ (-6:1:5),
+ numruns = 5,
)
-```
-
-
-## Appendix
-Computer information:
-```julia
-using InteractiveUtils
-InteractiveUtils.versioninfo()
+plot(wp, title="Fixed steps - with smoothing", color=colors)
```
-
-Package Information:
-```julia
-using Pkg
-Pkg.status()
+```@raw html
+Discrete time-series errors (l2):
+```
+```julia
+DENSE = true;
+SAVE_EVERYSTEP = true;
+
+_setups = [
+ "EK0(order=$order)" => Dict(:alg => EK0(order=order, smooth=DENSE))
+ for order in 2:7
+]
+
+labels = first.(_setups)
+setups = last.(_setups)
+
+abstols = 1.0 ./ 10.0 .^ (4:12)
+reltols = 1.0 ./ 10.0 .^ (1:9)
+
+wp = WorkPrecisionSet(
+ prob, abstols, reltols, setups;
+ names = labels,
+ appxsol = test_sol,
+ dense = DENSE,
+ save_everystep = SAVE_EVERYSTEP,
+ maxiters = Int(1e7),
+)
+
+plot(wp, x=:l2, palette=Plots.palette([:blue, :red], length(_setups)))
+plot!(ref_wp_dense, x=:l2, color=:gray, alpha=0.7, linestyle=:dash)
+```
+```@raw html
+Interoplation errors (L2):
+```
+```julia
+plot(wp, x=:L2, palette=Plots.palette([:blue, :red], length(_setups)))
+plot!(ref_wp_dense, x=:L2, color=:gray, alpha=0.7, linestyle=:dash)
+```
+```@raw html
+Discrete time-series errors (l2):
+```
+```julia
+DENSE = true;
+SAVE_EVERYSTEP = true;
+
+_setups = [
+ "EK1(order=$order)" => Dict(:alg => EK1(order=order, smooth=DENSE))
+ for order in 2:7
+]
+
+labels = first.(_setups)
+setups = last.(_setups)
+
+abstols = 1.0 ./ 10.0 .^ (4:12)
+reltols = 1.0 ./ 10.0 .^ (1:9)
wp = WorkPrecisionSet(
prob, abstols, reltols, setups;
names = labels,
- #print_names = true,
appxsol = test_sol,
dense = DENSE,
save_everystep = SAVE_EVERYSTEP,
- numruns = 10,
maxiters = Int(1e7),
- timeseries_errors = false,
- verbose = false,
)
-plot(wp, palette=Plots.palette([:blue, :red], length(_setups)), xticks = 10.0 .^ (-16:1:5))
+plot(wp, x=:l2, palette=Plots.palette([:blue, :red], length(_setups)))
+plot!(ref_wp_dense, x=:l2, color=:gray, alpha=0.7, linestyle=:dash)
+```
+```@raw html
+Interoplation errors (L2):
+```
+```julia
+plot(wp, x=:L2, palette=Plots.palette([:blue, :red], length(_setups)))
+plot!(ref_wp_dense, x=:L2, color=:gray, alpha=0.7, linestyle=:dash)
+```
+```@raw html
+Interoplation errors (L2):
+```
+```julia
+DENSE = true;
+SAVE_EVERYSTEP = true;
+
+_setups = [
+ "EK0(order=2)" => Dict(:alg => EK0(order=2, smooth=DENSE))
+ "EK0(order=3)" => Dict(:alg => EK0(order=3, smooth=DENSE))
+ "EK0(order=4)" => Dict(:alg => EK0(order=4, smooth=DENSE))
+ "EK0(order=5)" => Dict(:alg => EK0(order=5, smooth=DENSE))
+ "EK1(order=2)" => Dict(:alg => EK1(order=2, smooth=DENSE))
+ "EK1(order=3)" => Dict(:alg => EK1(order=3, smooth=DENSE))
+ "EK1(order=4)" => Dict(:alg => EK1(order=4, smooth=DENSE))
+ "EK1(order=5)" => Dict(:alg => EK1(order=5, smooth=DENSE))
+]
+
+labels = first.(_setups)
+setups = last.(_setups)
+
+abstols = 1.0 ./ 10.0 .^ (4:12)
+reltols = 1.0 ./ 10.0 .^ (1:9)
wp = WorkPrecisionSet(
prob, abstols, reltols, setups;
names = labels,
- #print_names = true,
appxsol = test_sol,
dense = DENSE,
save_everystep = SAVE_EVERYSTEP,
- numruns = 10,
maxiters = Int(1e7),
- timeseries_errors = false,
- verbose = false,
)
-plot(wp, color=[1 1 1 1 2 2 2 2], xticks = 10.0 .^ (-16:1:5))
+plot(wp, x=:L2, color=[1 1 1 1 2 2 2 2])
+plot!(ref_wp_dense, x=:L2, color=:gray, alpha=0.7, linestyle=:dash)
```
+```@raw html
+Code:
+ ```
+ ```julia""",
+ codeend="""
+ ```
+ ```@raw html
+ Computer information:
+ ```
+ ```julia
+ using InteractiveUtils
+ InteractiveUtils.versioninfo()
+ ```
+ ```
+ $(chomp(vinfo))
+ ```
+ ```@raw html
+ Package information:
+ ```
+ ```julia
+ using Pkg
+ Pkg.status()
+ ```
+ ```
+ $(chomp(proj))
+ ```
+ ```@raw html
+ Full manifest:
+ ```
+ ```julia
+ Pkg.status(mode=Pkg.PKGMODE_MANIFEST)
+ ```
+ ```
+ $(chomp(mani))
+ ```
+ ```@raw html
+ Discrete time-series errors (l2):
+```
+```julia
+DENSE = true;
+SAVE_EVERYSTEP = true;
+
+_setups = [
+ "EK1($order)" => Dict(:alg => EK1(order=order, smooth=DENSE))
+ for order in 3:7
+]
+
+labels = first.(_setups)
+setups = last.(_setups)
+
+abstols = 1.0 ./ 10.0 .^ (6:11)
+reltols = 1.0 ./ 10.0 .^ (3:8)
+
+wp = WorkPrecisionSet(
+ prob, abstols, reltols, setups;
+ names = labels,
+ appxsol = test_sol,
+ dense = DENSE,
+ save_everystep = SAVE_EVERYSTEP,
+ maxiters = Int(1e7),
+ numruns = 5,
+)
+
+plot(wp, x=:l2, palette=Plots.palette([:blue, :red], length(_setups)))
+plot!(ref_wp_dense, x=:l2, color=:gray, alpha=0.7, linestyle=:dash)
+```
+```@raw html
+Interpolation errors (L2):
+```
+```julia
+plot(wp, x=:L2, palette=Plots.palette([:blue, :red], length(_setups)))
+plot!(ref_wp_dense, x=:L2, color=:gray, alpha=0.7, linestyle=:dash)
+```
+```@raw html
+Interpolation errors (L2):
+```
```julia
-using InteractiveUtils
-InteractiveUtils.versioninfo()
+plot(wp; x=:L2, color)
+plot!(ref_wp_dense, x=:L2, color=:gray, alpha=0.7, linestyle=:dash)
+```
+```@raw html
+