diff --git a/benchmarks/hodgkinhuxley.jmd b/benchmarks/hodgkinhuxley.jmd
index dc612d9e7..d99f78d29 100644
--- a/benchmarks/hodgkinhuxley.jmd
+++ b/benchmarks/hodgkinhuxley.jmd
@@ -9,7 +9,7 @@
- [Results are similar for fixed time steps.](@ref hh_fixed_steps)
-```julia
+```julia, results="hidden"
using LinearAlgebra, Statistics, Distributions
using DiffEqDevTools, SciMLBase, OrdinaryDiffEq, Plots, SimpleUnPack
using ProbNumDiffEq
@@ -111,6 +111,33 @@ wp = WorkPrecisionSet(
)
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)
```
### With smoothing
@@ -146,6 +173,7 @@ wp = WorkPrecisionSet(
)
plot(wp, title="Adaptive steps - with smoothing", color=colors)
+plot!(ref_wp_dense, x=:final, color=:gray, alpha=0.7, linestyle=:dash)
```
```@raw html
@@ -153,6 +181,7 @@ plot(wp, title="Adaptive steps - with smoothing", color=colors)
```
```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
diff --git a/benchmarks/lotkavolterra.jmd b/benchmarks/lotkavolterra.jmd
index 738332518..0f39dc556 100644
--- a/benchmarks/lotkavolterra.jmd
+++ b/benchmarks/lotkavolterra.jmd
@@ -20,7 +20,7 @@
Benchmark adapted from
[SciMLBenchmarks.jl](https://docs.sciml.ai/SciMLBenchmarksOutput/stable/NonStiffODE/LotkaVolterra_wpd/).
-```julia
+```julia, results="hidden"
using LinearAlgebra, Statistics, Distributions
using DiffEqDevTools, ParameterizedFunctions, SciMLBase, OrdinaryDiffEq, Plots
using ProbNumDiffEq
@@ -32,6 +32,14 @@ Plots.theme(
margin=5Plots.mm,
xticks=10.0 .^ (-16:1:16),
)
+
+function plot_chisq_interval!(df, q=0.01)
+ dist = Chisq(df)
+ low, high, mid = quantile(dist, [q, 1-q])..., mean(dist)
+ hline!([low, high], linestyle=:dash, color=:black, label="",
+ fill_between=[high nothing], fillcolor=:green, fillalpha=0.15)
+ hline!([mid], linestyle=:solid, color=:black, label="")
+end
```
```julia
@@ -75,6 +83,32 @@ wp = WorkPrecisionSet(
)
plot(wp, palette=Plots.palette([:blue, :red], length(_setups)))
+
+_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, reltols, ref_setups;
+ names = ref_labels,
+ appxsol = test_sol,
+ dense = false,
+ save_everystep = false,
+ maxiters = Int(1e7),
+)
+ref_wp_dense = WorkPrecisionSet(
+ prob, abstols, reltols, 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)
```
```@raw html
@@ -105,6 +139,7 @@ wp = WorkPrecisionSet(
)
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
@@ -115,6 +150,7 @@ plot(wp, x=:l2, palette=Plots.palette([:blue, :red], length(_setups)))
```
```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
@@ -147,6 +183,7 @@ wp = WorkPrecisionSet(
)
plot(wp, palette=Plots.palette([:blue, :red], length(_setups)))
+plot!(ref_wp_final, x=:final, color=:gray, alpha=0.7, linestyle=:dash)
```
```@raw html
@@ -177,6 +214,7 @@ wp = WorkPrecisionSet(
)
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
@@ -187,6 +225,7 @@ plot(wp, x=:l2, palette=Plots.palette([:blue, :red], length(_setups)))
```
```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
@@ -225,6 +264,7 @@ wp = WorkPrecisionSet(
)
plot(wp, color=[1 1 1 1 2 2 2 2])
+plot!(ref_wp_final, x=:final, color=:gray, alpha=0.7, linestyle=:dash)
```
```@raw html
@@ -261,6 +301,7 @@ wp = WorkPrecisionSet(
)
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
@@ -269,15 +310,6 @@ plot(wp, x=:L2, color=[1 1 1 1 2 2 2 2])
## [`EK0`](@ref) vs. [`EK1`](@ref): Calibration
```julia
plot(wp, x=:final, y=:chi2_final, color=[1 1 1 1 2 2 2 2], yguide="Chi-squared (final)")
-
-# Should be distributed according to a Chi-squared distribution:
-function plot_chisq_interval!(df, q=0.01)
- dist = Chisq(df)
- low, high, mid = quantile(dist, [q, 1-q])..., mean(dist)
- hline!([low, high], linestyle=:dash, color=:black, label="",
- fill_between=[high nothing], fillcolor=:green, fillalpha=0.15)
- hline!([mid], linestyle=:solid, color=:black, label="")
-end
plot_chisq_interval!(2)
```
diff --git a/benchmarks/orego.jmd b/benchmarks/orego.jmd
index e7553f835..11c90dcd7 100644
--- a/benchmarks/orego.jmd
+++ b/benchmarks/orego.jmd
@@ -2,13 +2,14 @@
!!! note "Summary"
- TODO
+ - [**The `EK1` is able to solve mass-matrix DAEs.** To achieve low error, use order 4 or higher.](@ref orego_results)
+ - The order-to-error-tolerance heuristic holds: lower tolerance level ``\rightarrow`` higher order.
Adapted from
[SciMLBenchmarks.jl](https://docs.sciml.ai/SciMLBenchmarksOutput/stable/DAE/OregoDAE/).
-```julia
+```julia, results="hidden"
using LinearAlgebra, Statistics, Distributions
using DiffEqDevTools, ParameterizedFunctions, SciMLBase, OrdinaryDiffEq, Sundials, Plots
using ModelingToolkit
@@ -21,9 +22,15 @@ Plots.theme(
margin=5Plots.mm,
xticks=10.0 .^ (-16:1:16),
)
-```
-### OREGO problem definition
+function plot_chisq_interval!(df, q=0.01)
+ dist = Chisq(df)
+ low, high, mid = quantile(dist, [q, 1-q])..., mean(dist)
+ hline!([low, high], linestyle=:dash, color=:black, label="",
+ fill_between=[high nothing], fillcolor=:green, fillalpha=0.15)
+ hline!([mid], linestyle=:solid, color=:black, label="")
+end
+```
```julia
@variables t y1(t)=1.0 y2(t)=2.0 y3(t)=3.0
@@ -47,7 +54,7 @@ ref_sol = solve(daeprob,IDA(),abstol=1/10^14,reltol=1/10^14)
plot(ref_sol, title="OREGO Solution", legend=false, xticks=:auto)
```
-## EK1 accross orders
+## [`EK1` across orders](@id orego_results)
```julia
DENSE = false;
@@ -55,14 +62,14 @@ SAVE_EVERYSTEP = false;
_setups = [
"EK1($order)" => Dict(:alg => EK1(order=order, smooth=DENSE))
- for order in 2:4
+ for order in 2:6
]
labels = first.(_setups)
setups = last.(_setups)
-abstols = 1.0 ./ 10.0 .^ (4:8)
-reltols = 1.0 ./ 10.0 .^ (1:5)
+abstols = 1.0 ./ 10.0 .^ (6:10)
+reltols = 1.0 ./ 10.0 .^ (3:7)
wp = WorkPrecisionSet(
mmprob, abstols, reltols, setups;
@@ -75,6 +82,24 @@ wp = WorkPrecisionSet(
)
plot(wp, palette=Plots.palette([:blue, :red], length(_setups)))
+
+
+_ref_setups = [
+ "Rosenbrock23" => Dict(:alg => Rosenbrock23())
+ "Rodas4P" => Dict(:alg => Rodas4P())
+ "RadauIIA" => Dict(:alg => RadauIIA5())
+]
+ref_labels = first.(_ref_setups)
+ref_setups = last.(_ref_setups)
+ref_wp = WorkPrecisionSet(
+ mmprob, abstols, reltols, ref_setups;
+ names = ref_labels,
+ appxsol = ref_sol,
+ dense = DENSE,
+ save_everystep = SAVE_EVERYSTEP,
+ maxiters = Int(1e7),
+)
+plot!(ref_wp, x=:final, color=:gray, alpha=0.7, linestyle=:dash)
```
### Calibration
@@ -82,14 +107,6 @@ plot(wp, palette=Plots.palette([:blue, :red], length(_setups)))
```julia
plot(wp; x=:final, y=:chi2_final, yguide="Chi-squared (final)",
palette=Plots.palette([:blue, :red], length(_setups)))
-
-function plot_chisq_interval!(df, q=0.01)
- dist = Chisq(df)
- low, high, mid = quantile(dist, [q, 1-q])..., mean(dist)
- hline!([low, high], linestyle=:dash, color=:black, label="",
- fill_between=[high nothing], fillcolor=:green, fillalpha=0.15)
- hline!([mid], linestyle=:solid, color=:black, label="")
-end
plot_chisq_interval!(3)
```
diff --git a/benchmarks/pleiades.jmd b/benchmarks/pleiades.jmd
index 64b5e7b58..c13dcab9e 100644
--- a/benchmarks/pleiades.jmd
+++ b/benchmarks/pleiades.jmd
@@ -3,11 +3,11 @@
!!! note "Summary"
Pleiades is a medium-dimensional, non-stiff, second-order ODE. We see that:
- - The `EK0` is _much_ faster than the `EK1` as it scales linearly with the ODE dimension.
- - If the problem is a second-order ODE, _implement it as a second-order ODE_!
+ - [**The `EK0` is _much_ faster than the `EK1` as it scales linearly with the ODE dimension.**](@ref pleiades_results)
+ - [**If the problem is a second-order ODE, _implement it as a second-order ODE_!**](@ref pleiades_results)
-```julia
+```julia, results="hidden"
using LinearAlgebra, Statistics, Distributions
using DiffEqDevTools, ParameterizedFunctions, SciMLBase, OrdinaryDiffEq, Sundials, Plots, ODEInterfaceDiffEq
using ModelingToolkit
@@ -20,9 +20,15 @@ Plots.theme(
margin=5Plots.mm,
xticks=10.0 .^ (-16:1:16),
)
-```
-### Pleiades problem definition
+function plot_chisq_interval!(df, q=0.01)
+ dist = Chisq(df)
+ low, high, mid = quantile(dist, [q, 1-q])..., mean(dist)
+ hline!([low, high], linestyle=:dash, color=:black, label="",
+ fill_between=[high nothing], fillcolor=:green, fillalpha=0.15)
+ hline!([mid], linestyle=:solid, color=:black, label="")
+end
+```
```julia
# first-order ODE
@@ -83,7 +89,7 @@ plot(ref_sol1, idxs=[(14+i,21+i) for i in 1:7], title="Pleiades Solution", legen
scatter!(ref_sol1.u[end][15:21], ref_sol1.u[end][22:end], color=1:7)
```
-## EK0 vs EK1 & first-order vs. second-order
+## [`EK0` vs `EK1` & first-order vs. second-order](@id pleiades_results)
```julia
DENSE = false;
SAVE_EVERYSTEP = false;
@@ -97,9 +103,6 @@ _setups = [
"EK1(5) (1st order ODE)" => Dict(:alg => EK1(order=5, smooth=DENSE), :prob_choice => 1)
"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)
- "Classic: Tsit5" => Dict(:alg => Tsit5(), :prob_choice => 1)
- "Classic: RadauIIA5" => Dict(:alg => RadauIIA5(), :prob_choice => 1)
- "Classic: DPRKN6" => Dict(:alg => DPRKN6(), :prob_choice => 2)
]
labels = first.(_setups)
@@ -118,30 +121,32 @@ wp = WorkPrecisionSet(
maxiters = Int(1e7),
)
-plot(wp,
- color=[1 1 1 1 2 2 2 2 3 3 3],
- linestyle=[:solid :solid :dash :dash :solid :solid :dash :dash :solid :solid :dash],
- alpha=[1 1 1 1 1 1 1 1 0.5 0.5 0.5],
+color = [1 1 1 1 2 2 2 2]
+linestyle = [:solid :solid :dash :dash :solid :solid :dash :dash]
+plot(wp; color, linestyle)
+
+_ref_setups = [
+ "Classic: Tsit5" => Dict(:alg => Tsit5(), :prob_choice => 1)
+ "Classic: RadauIIA5" => Dict(:alg => RadauIIA5(), :prob_choice => 1)
+ "Classic: DPRKN6" => Dict(:alg => DPRKN6(), :prob_choice => 2)
+]
+ref_labels = first.(_ref_setups)
+ref_setups = last.(_ref_setups)
+ref_wp = WorkPrecisionSet(
+ probs, abstols ./ 1000, reltols ./ 1000, ref_setups;
+ names = ref_labels,
+ appxsol = ref_sols,
+ dense = false,
+ save_everystep = false,
+ maxiters = Int(1e7),
)
+plot!(ref_wp, x=:final, color=:gray, alpha=0.7, linestyle=[:solid :solid :dash])
```
### Calibration
```julia
-# We can only evaluate calibration for the PN solvers
-_wp = WorkPrecisionSet(wp.wps[1:end-3], wp.N-3, wp.abstols, wp.reltols, wp.prob, wp.setups[1:end-3],
- wp.names[1:end-3], wp.error_estimate, wp.numruns)
-plot(_wp; x=:final, y=:chi2_final,
- color=[1 1 1 1 2 2 2 2],
- linestyle=[:solid :solid :dash :dash :solid :solid :dash :dash],
-)
-function plot_chisq_interval!(df, q=0.01)
- dist = Chisq(df)
- low, high, mid = quantile(dist, [q, 1-q])..., mean(dist)
- hline!([low, high], linestyle=:dash, color=:black, label="",
- fill_between=[high nothing], fillcolor=:green, fillalpha=0.15)
- hline!([mid], linestyle=:solid, color=:black, label="")
-end
+plot(wp; x=:final, y=:chi2_final, color, linestyle)
plot_chisq_interval!(length(u0)*2)
```
diff --git a/benchmarks/rober.jmd b/benchmarks/rober.jmd
index df379859e..f24a9accf 100644
--- a/benchmarks/rober.jmd
+++ b/benchmarks/rober.jmd
@@ -2,14 +2,14 @@
!!! note "Summary"
- - The `EK1` can solve mass-matrix DAEs! But it only really works well for low tolerances.
- - 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.
+ - [**The `EK1` can solve mass-matrix DAEs.** But for this problem, it only works well for low tolerances.](@ref rober_results)
+ - For this problem it only works well for low tolerances, but the order-to-error-tolerance heuristic should in principle still hold: lower tolerance level ``\rightarrow`` higher order.
Adapted from
[SciMLBenchmarks.jl](https://docs.sciml.ai/SciMLBenchmarksOutput/stable/DAE/ROBERDAE/).
-```julia
+```julia, results="hidden"
using LinearAlgebra, Statistics, Distributions
using DiffEqDevTools, ParameterizedFunctions, SciMLBase, OrdinaryDiffEq, Sundials, Plots
using ModelingToolkit
@@ -22,9 +22,15 @@ Plots.theme(
margin=5Plots.mm,
xticks=10.0 .^ (-16:1:16),
)
-```
-### ROBER problem definition
+function plot_chisq_interval!(df, q=0.01)
+ dist = Chisq(df)
+ low, high, mid = quantile(dist, [q, 1-q])..., mean(dist)
+ hline!([low, high], linestyle=:dash, color=:black, label="",
+ fill_between=[high nothing], fillcolor=:green, fillalpha=0.15)
+ hline!([mid], linestyle=:solid, color=:black, label="")
+end
+```
```julia
@variables t y₁(t)=1.0 y₂(t)=0.0 y₃(t)=0.0
@@ -44,7 +50,7 @@ ref_sol = solve(daeprob,IDA(),abstol=1/10^14,reltol=1/10^14)
plot(ref_sol, idxs=[y₁,y₂,y₃], title="ROBER Solution", legend=false, ylims=(0, 1), xticks=:auto)
```
-## EK1 across orders
+## [`EK1` across orders](@id rober_results)
```julia
DENSE = false;
@@ -58,8 +64,8 @@ _setups = [
labels = first.(_setups)
setups = last.(_setups)
-abstols = 1.0 ./ 10.0 .^ (4:8)
-reltols = 1.0 ./ 10.0 .^ (1:5)
+abstols = 1.0 ./ 10.0 .^ (4:7)
+reltols = 1.0 ./ 10.0 .^ (1:4)
wp = WorkPrecisionSet(
mmprob, abstols, reltols, setups;
@@ -72,6 +78,24 @@ wp = WorkPrecisionSet(
)
plot(wp, palette=Plots.palette([:blue, :red], length(_setups)))
+
+
+_ref_setups = [
+ "Rosenbrock23" => Dict(:alg => Rosenbrock23())
+ "Rodas4P" => Dict(:alg => Rodas4P())
+ "RadauIIA" => Dict(:alg => RadauIIA5())
+]
+ref_labels = first.(_ref_setups)
+ref_setups = last.(_ref_setups)
+ref_wp = WorkPrecisionSet(
+ mmprob, abstols ./ 10000, reltols ./ 10000, ref_setups;
+ names = ref_labels,
+ appxsol = ref_sol,
+ dense = DENSE,
+ save_everystep = SAVE_EVERYSTEP,
+ maxiters = Int(1e7),
+)
+plot!(ref_wp, x=:final, color=:gray, alpha=0.7, linestyle=:dash)
```
### Calibration
@@ -79,14 +103,6 @@ plot(wp, palette=Plots.palette([:blue, :red], length(_setups)))
```julia
plot(wp; x=:final, y=:chi2_final, yguide="Chi-squared (final)",
palette=Plots.palette([:blue, :red], length(_setups)))
-
-function plot_chisq_interval!(df, q=0.01)
- dist = Chisq(df)
- low, high, mid = quantile(dist, [q, 1-q])..., mean(dist)
- hline!([low, high], linestyle=:dash, color=:black, label="",
- fill_between=[high nothing], fillcolor=:green, fillalpha=0.15)
- hline!([mid], linestyle=:solid, color=:black, label="")
-end
plot_chisq_interval!(3)
```
diff --git a/benchmarks/vanderpol.jmd b/benchmarks/vanderpol.jmd
index 6cf0eb5f4..ed7e37683 100644
--- a/benchmarks/vanderpol.jmd
+++ b/benchmarks/vanderpol.jmd
@@ -3,11 +3,13 @@
!!! note "Summary"
Van der Pol is a low-dimensional, stiff, second-order ODE. We see that:
- - The `EK1` is very well able to solve stiff problems.
- - Since Van der Pol is actually a second-order ODE, _do solve it as a second-order ODE_.
+ - [**The `EK1` is very well able to solve stiff problems.**](@ref vdp_main_results)
+ - [**Since Van der Pol is actually a second-order ODE, _do solve it as a second-order ODE_.**](@ref vdp_second_order)
+ - [**Use the `TaylorInit` or `ForwardDiffInit` initialization.**](@ref vdp_initialization)
+ While `SimpleInit` works well for lower orders, it fails for higher orders. And since Taylor-mode initialization is fast and works well, there is no reason not to use it.
-```julia
+```julia, results="hidden"
using LinearAlgebra, Statistics, Distributions
using DiffEqDevTools, ParameterizedFunctions, SciMLBase, OrdinaryDiffEq, Plots
using ProbNumDiffEq
@@ -19,9 +21,15 @@ Plots.theme(
margin=5Plots.mm,
xticks = 10.0 .^ (-16:1:16)
)
-```
-### Van der Pol problem definition
+function plot_chisq_interval!(df, q=0.01)
+ dist = Chisq(df)
+ low, high, mid = quantile(dist, [q, 1-q])..., mean(dist)
+ hline!([low, high], linestyle=:dash, color=:black, label="",
+ fill_between=[high nothing], fillcolor=:green, fillalpha=0.15)
+ hline!([mid], linestyle=:solid, color=:black, label="")
+end
+```
```julia
function vanderpol!(du, u, p, t)
@@ -37,11 +45,8 @@ test_sol = solve(prob, RadauIIA5(), abstol=1/10^14, reltol=1/10^14)
plot(test_sol, title="Van der Pol Solution", legend=false, ylims=(-5, 5), xticks=:auto)
```
-## EK1 across orders
-
-### Final value only
+## [`EK1` across orders](@id vdp_main_results)
-This does not require smoothing or saving intermediate values.
```julia
DENSE = false;
SAVE_EVERYSTEP = false;
@@ -68,11 +73,37 @@ wp = WorkPrecisionSet(
)
plot(wp, palette=Plots.palette([:blue, :red], length(_setups)))
-```
-### Full trajectory
+_ref_setups = [
+ "Rosenbrock23" => Dict(:alg => Rosenbrock23())
+ "Rodas4P" => Dict(:alg => Rodas4P())
+ "RadauIIA5" => Dict(:alg => RadauIIA5())
+]
+ref_labels = first.(_ref_setups)
+ref_setups = last.(_ref_setups)
+ref_wp_final = WorkPrecisionSet(
+ prob, abstols, reltols, ref_setups;
+ names = ref_labels,
+ appxsol = test_sol,
+ dense = false,
+ save_everystep = false,
+ maxiters = Int(1e7),
+)
+ref_wp_dense = WorkPrecisionSet(
+ prob, abstols, reltols, 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)
+```
-This does require smoothing or saving intermediate values.
+```@raw html
+Discrete time-series errors (l2):
+```
```julia
DENSE = true;
SAVE_EVERYSTEP = true;
@@ -99,30 +130,31 @@ wp = WorkPrecisionSet(
)
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
plot(wp; x=:L2, color)
+plot!(ref_wp_dense, x=:L2, color=:gray, alpha=0.7, linestyle=:dash)
+```
+```@raw html
+