Towards a more modular state init by adding more options (#266)

benchmarks/lotkavolterra.jmd

@@ -175,44 +175,50 @@ wp = WorkPrecisionSet(
 plot(wp, color=[2 2 2 3 3 3], xticks = 10.0 .^ (-16:1:5))
 ```
 
-## TaylorModeInit vs ClassicSolverInit
+## Comparison of the different initialization schemes
 
 ```julia
 DENSE = false;
 SAVE_EVERYSTEP = false;
 
-_setups = [
-  "EK1(2) TaylorInit" => Dict(:alg => EK1(order=2, smooth=DENSE, initialization=TaylorModeInit()))
-  "EK1(3) TaylorInit" => Dict(:alg => EK1(order=3, smooth=DENSE, initialization=TaylorModeInit()))
-  "EK1(5) TaylorInit" => Dict(:alg => EK1(order=5, smooth=DENSE, initialization=TaylorModeInit()))
-  "EK1(2) Tsit5Init" => Dict(:alg => EK1(order=2, smooth=DENSE, initialization=ClassicSolverInit()))
-  "EK1(3) Tsit5Init" => Dict(:alg => EK1(order=3, smooth=DENSE, initialization=ClassicSolverInit()))
-  "EK1(5) Tsit5Init" => Dict(:alg => EK1(order=5, smooth=DENSE, initialization=ClassicSolverInit()))
-  "EK1(2) Tsit5Init+ddu" => Dict(:alg => EK1(order=2, smooth=DENSE, initialization=ClassicSolverInit(init_on_ddu=true)))
-  "EK1(3) Tsit5Init+ddu" => Dict(:alg => EK1(order=3, smooth=DENSE, initialization=ClassicSolverInit(init_on_ddu=true)))
-  "EK1(5) Tsit5Init+ddu" => Dict(:alg => EK1(order=5, smooth=DENSE, initialization=ClassicSolverInit(init_on_ddu=true)))
-]
-
-labels = first.(_setups)
-setups = last.(_setups)
-
 abstols = 1.0 ./ 10.0 .^ (4:14)
 reltols = 1.0 ./ 10.0 .^ (1:11)
 
-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,
+orders = (2, 3, 5, 8)
+ps = []
+for o in orders
+    _setups = [
+      "EK1($o) TaylorInit" => Dict(:alg => EK1(order=o, smooth=DENSE, initialization=TaylorModeInit(o)))
+      "EK1($o) ForwardDiffInit" => Dict(:alg => EK1(order=o, smooth=DENSE, initialization=ForwardDiffInit(o)))
+      "EK1($o) SimpleInit" => Dict(:alg => EK1(order=o, smooth=DENSE, initialization=SimpleInit()))
+      "EK1($o) ClassicSolverInit" => Dict(:alg => EK1(order=o, smooth=DENSE, initialization=ClassicSolverInit()))
+    ]
+
+    labels = first.(_setups)
+    setups = last.(_setups)
+
+    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,
+    )
+
+    p = plot(wp, color=[2 4 5 6], xticks = 10.0 .^ (-16:1:5))
+    push!(ps, p)
+end
+plot(
+    ps...,
+    layout=(length(orders), 1),
+    size = (1000, length(orders)*300),
+    xlabel=["" "" "" "Error"],
 )
-
-plot(wp, color=[2 2 2 4 4 4 5 5 5], xticks = 10.0 .^ (-16:1:5))
 ```
 
 

benchmarks/vanderpol.jmd

@@ -64,6 +64,90 @@ wp = WorkPrecisionSet(
 plot(wp, palette=Plots.palette([:blue, :red], length(_setups)), xticks = 10.0 .^ (-16:1:5))
 ```
 
+## Comparison of the different initialization schemes
+
+```julia
+DENSE = false;
+SAVE_EVERYSTEP = false;
+
+abstols = 1.0 ./ 10.0 .^ (6:11)
+reltols = 1.0 ./ 10.0 .^ (3:8)
+
+orders = (3, 5, 8)
+ps = []
+for o in orders
+    _setups = [
+      "EK1($o) TaylorInit" => Dict(:alg => EK1(order=o, smooth=DENSE, initialization=TaylorModeInit(o)))
+      "EK1($o) ForwardDiffInit" => Dict(:alg => EK1(order=o, smooth=DENSE, initialization=ForwardDiffInit(o)))
+      "EK1($o) SimpleInit" => Dict(:alg => EK1(order=o, smooth=DENSE, initialization=SimpleInit()))
+      # "EK1($o) ClassicSolverInit" => Dict(:alg => EK1(order=o, smooth=DENSE, initialization=ClassicSolverInit())) # unstable
+    ]
+
+    labels = first.(_setups)
+    setups = last.(_setups)
+
+    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,
+    )
+
+    p = plot(wp, color=[2 4 5 6], xticks = 10.0 .^ (-16:1:5))
+    push!(ps, p)
+end
+plot(
+    ps...,
+    layout=(length(orders), 1),
+    size = (1000, length(orders)*300),
+    xlabel=["" "" "" "Error"],
+)
+```
+
+```julia
+DENSE = false;
+SAVE_EVERYSTEP = false;
+
+abstols = 1.0 ./ 10.0 .^ (5:8)
+reltols = 1.0 ./ 10.0 .^ (2:5)
+
+_setups = [
+    "EK1(8) TaylorInit(8)" => Dict(:alg => EK1(order=8, smooth=DENSE, initialization=TaylorModeInit(8)))
+    "EK1(8) TaylorInit(7)" => Dict(:alg => EK1(order=8, smooth=DENSE, initialization=TaylorModeInit(7)))
+    "EK1(8) TaylorInit(6)" => Dict(:alg => EK1(order=8, smooth=DENSE, initialization=TaylorModeInit(6)))
+    "EK1(8) TaylorInit(5)" => Dict(:alg => EK1(order=8, smooth=DENSE, initialization=TaylorModeInit(5)))
+    "EK1(8) TaylorInit(4)" => Dict(:alg => EK1(order=8, smooth=DENSE, initialization=TaylorModeInit(4)))
+    "EK1(8) TaylorInit(3)" => Dict(:alg => EK1(order=8, smooth=DENSE, initialization=TaylorModeInit(3)))
+    "EK1(8) TaylorInit(2)" => Dict(:alg => EK1(order=8, smooth=DENSE, initialization=TaylorModeInit(2)))
+    # "EK1(8) TaylorInit(1)" => Dict(:alg => EK1(order=8, smooth=DENSE, initialization=TaylorModeInit(1))) # fails, see above
+]
+
+labels = first.(_setups)
+setups = last.(_setups)
+
+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))
+```
+
+
 ## Solving the first- vs second-order ODE
 
 ```julia

docs/Project.toml

@@ -12,3 +12,4 @@ OptimizationOptimJL = "36348300-93cb-4f02-beb5-3c3902f8871e"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 ProbNumDiffEq = "bf3e78b0-7d74-48a5-b855-9609533b56a5"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"