Fitzhugh-Nagumo Work-precision Diagrams (including ESERK4)

The purpose of this is to see how the errors scale on a standard nonlinear problem.

Problem

using Pkg
Pkg.update()
using Pkg, OrdinaryDiffEq, ParameterizedFunctions, DiffEqDevTools, Plots

f = @ode_def FitzhughNagumo begin
  dv = v - v^3/3 -w + l
  dw = τinv*(v +  a - b*w)
end a b τinv l

p = [0.7,0.8,1/12.5,0.5]
prob = ODEProblem(f,[1.0;1.0],(0.0,10.0),p)

abstols = 1.0 ./ 10.0 .^ (6:13)
reltols = 1.0 ./ 10.0 .^ (3:10);

sol = solve(prob,Vern7(),abstol=1/10^14,reltol=1/10^14)
test_sol = TestSolution(sol)
using Plots; gr()
GRBackend()
plot(sol,dpi=200,linewidth=1,legend=:bottomleft,legendfontsize=6)

Setups

setups = [
          Dict(:alg=>DP5()),
          Dict(:alg=>BS5()),
          Dict(:alg=>Tsit5()),
          Dict(:alg=>Vern6()),
          #Dict(:alg=>RKC()),
          Dict(:alg=>ROCK2()),
          Dict(:alg=>ROCK4()),
          Dict(:alg=>SERK2(controller=:PI)),
          Dict(:alg=>SERK2(controller=:Predictive)),
  				Dict(:alg=>ESERK4()),
          Dict(:alg=>ESERK5())
          ]
Names = ["DP5" "BS5" "Tsit5" "Vern6" "ROCK2" "ROCK4" "SERK2 PI" "SERK2 Predictive" "ESERK4" "ESERK5"]
1×10 Array{String,2}: "DP5" "BS5" "Tsit5" "Vern6" "ROCK2" "ROCK4" … "ESERK4" "ESERK5"

Speed only Tests

wp = WorkPrecisionSet(prob,abstols,reltols,setups;names=Names,appxsol=test_sol,save_everystep=false,numruns=100,maxiters=Int(1e6))
plot(wp,dpi=200,linewidth=1,legend=:topright,legendfontsize=6)
wp = WorkPrecisionSet(prob,abstols,reltols,setups;names=Names,appxsol=test_sol,save_everystep=false,numruns=100,maxiters=Int(1e6),error_estimate=:L2,dense_errors=true)
plot(wp,dpi=200,linewidth=1,legend=:topright,legendfontsize=6)

WP diagram after removing ROCK4 and BS5. Just to see clearly whats going on.

setups = [
          Dict(:alg=>DP5()),
          #Dict(:alg=>BS5()),
          Dict(:alg=>Tsit5()),
          Dict(:alg=>Vern6()),
          #Dict(:alg=>RKC()),
          Dict(:alg=>ROCK2()),
          #Dict(:alg=>ROCK4()),
          Dict(:alg=>SERK2(controller=:PI)),
          Dict(:alg=>SERK2(controller=:Predictive)),
          Dict(:alg=>ESERK4()),
  				Dict(:alg=>ESERK5())
          ]
Names = ["DP5" "Tsit5" "Vern6" "ROCK2" "SERK2 PI" "SERK2 Predictive" "ESERK4" "ESERK5"]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;names=Names,appxsol=test_sol,save_everystep=false,numruns=100,maxiters=Int(1e6),error_estimate=:L2,dense_errors=true)
plot(wp,dpi=200,linewidth=1,legend=:topleft,legendfontsize=6)

Thus we see that Stabilized methods are showing respectable preformance as compared to top performers of OrdinaryDiffEq.jl. RKC is commented out because it requires more iterations.