update benchmarks

Author: ChrisRackauckas

README.md

@@ -109,31 +109,155 @@ julia> @time MATLABDiffEq.eval_string("[t,u] = $(algstr)(f,tspan,u0,options);")
 
 ## Benchmark
 
-MATLABDiffEq.jl will be included into [DiffEqBenchmarks.jl](https://github.com/JuliaDiffEq/DiffEqBenchmarks.jl). However, the benchmarks are not encouraging to MATLAB at all.  Instead, they show that the cost of evaluating functions in MATLAB are too large (note that the functions are re-defined in MATLAB, meaning that they are true MATLAB functions and not interop functions).
-
-Running benchmarks at various tolerance using the same Lotka-Volterra problem from before, we get the following:
+The following benchmarks demonstrate a **100x performance advantage for the
+pure-Julia methods over the MATLAB ODE solvers** across a range of stiff and
+non-stiff ODEs. These were ran with Julia 1.2 and MATLAB 2019B after verifying
+negligible overhead on interop.
 
 ```julia
-using OrdinaryDiffEq, ODEInterfaceDiffEq, Plots, ODE
+using ParameterizedFunctions, MATLABDiffEq, OrdinaryDiffEq, ODEInterface,
+      ODEInterfaceDiffEq, Plots, Sundials
 using DiffEqDevTools
-abstols = 1./10.^(6:13)
-reltols = 1./10.^(3:10)
+using LinearAlgebra
+
+## Non-Stiff Problem 1: Lotka-Volterra
+
+f = @ode_def_bare LotkaVolterra begin
+  dx = a*x - b*x*y
+  dy = -c*y + d*x*y
+end a b c d
+p = [1.5,1,3,1]
+tspan = (0.0,10.0)
+u0 = [1.0,1.0]
+prob = ODEProblem(f,u0,tspan,p)
+sol = solve(prob,Vern7(),abstol=1/10^14,reltol=1/10^14)
+test_sol = TestSolution(sol)
+
+setups = [Dict(:alg=>DP5())
+          Dict(:alg=>dopri5())
+          Dict(:alg=>BS5())
+          Dict(:alg=>Tsit5())
+          Dict(:alg=>Vern6())
+          Dict(:alg=>Vern7())
+          Dict(:alg=>MATLABDiffEq.ode45())
+          Dict(:alg=>MATLABDiffEq.ode113())
+          Dict(:alg=>CVODE_Adams())
+  ]
+abstols = 1.0 ./ 10.0 .^ (6:13)
+reltols = 1.0 ./ 10.0 .^ (3:10)
+wp = WorkPrecisionSet(prob,abstols,reltols,setups;appxsol=test_sol,dense=false,save_everystep=false,numruns=100,maxiters=10000000,ttimeseries_errors=false,verbose=false)
+plot(wp,title="Non-stiff 1: Lotka-Volterra")
+savefig("benchmark1.png")
+```
+
+![benchmark1](https://user-images.githubusercontent.com/1814174/69478405-f9dac480-0dbf-11ea-8f8e-5572b86d2a97.png)
+
+```julia
+## Non-Stiff Problem 2: Rigid Body
+
+f = @ode_def_bare RigidBodyBench begin
+  dy1  = -2*y2*y3
+  dy2  = 1.25*y1*y3
+  dy3  = -0.5*y1*y2 + 0.25*sin(t)^2
+end
+prob = ODEProblem(f,[1.0;0.0;0.9],(0.0,100.0))
 sol = solve(prob,Vern7(),abstol=1/10^14,reltol=1/10^14)
 test_sol = TestSolution(sol)
-plotly()
+
 setups = [Dict(:alg=>DP5())
           Dict(:alg=>dopri5())
           Dict(:alg=>BS5())
           Dict(:alg=>Tsit5())
           Dict(:alg=>Vern6())
           Dict(:alg=>Vern7())
           Dict(:alg=>MATLABDiffEq.ode45())
+          Dict(:alg=>MATLABDiffEq.ode113())
+          Dict(:alg=>CVODE_Adams())
   ]
-wp = WorkPrecisionSet(prob,abstols,reltols,setups;appxsol=test_sol,dense=false,save_everystep=false,numruns=1000,maxiters=10000000,ttimeseries_errors=false,verbose=false)
-plot(wp)
+abstols = 1.0 ./ 10.0 .^ (6:13)
+reltols = 1.0 ./ 10.0 .^ (3:10)
+wp = WorkPrecisionSet(prob,abstols,reltols,setups;appxsol=test_sol,dense=false,save_everystep=false,numruns=100,maxiters=10000000,ttimeseries_errors=false,verbose=false)
+plot(wp,title="Non-stiff 2: Rigid-Body")
+savefig("benchmark2.png")
+```
+
+![benchmark2](https://user-images.githubusercontent.com/1814174/69478406-fe9f7880-0dbf-11ea-8f15-a0a0f3e8717f.png)
+
+
+```julia
+## Stiff Problem 1: ROBER
+
+rober = @ode_def begin
+  dy₁ = -k₁*y₁+k₃*y₂*y₃
+  dy₂ =  k₁*y₁-k₂*y₂^2-k₃*y₂*y₃
+  dy₃ =  k₂*y₂^2
+end k₁ k₂ k₃
+prob = ODEProblem(rober,[1.0,0.0,0.0],(0.0,1e5),[0.04,3e7,1e4])
+sol = solve(prob,CVODE_BDF(),abstol=1/10^14,reltol=1/10^14)
+test_sol = TestSolution(sol)
+
+abstols = 1.0 ./ 10.0 .^ (5:8)
+reltols = 1.0 ./ 10.0 .^ (1:4);
+setups = [Dict(:alg=>Rosenbrock23()),
+          Dict(:alg=>TRBDF2()),
+          Dict(:alg=>rodas()),
+          Dict(:alg=>radau()),
+          Dict(:alg=>RadauIIA5()),
+          Dict(:alg=>MATLABDiffEq.ode23s()),
+          Dict(:alg=>MATLABDiffEq.ode15s())
+          ]
+
+wp = WorkPrecisionSet(prob,abstols,reltols,setups;
+                      save_everystep=false,appxsol=test_sol,maxiters=Int(1e5),numruns=100)
+plot(wp,title="Stiff 1: ROBER")
+savefig("benchmark3.png")
 ```
 
-![Benchmark](assets/matlab_bench.png)
+![benchmark3](https://user-images.githubusercontent.com/1814174/69478407-ffd0a580-0dbf-11ea-9311-909bd3938b42.png)
+
+
+```julia
+## Stiff Problem 2: HIRES
+
+f = @ode_def Hires begin
+  dy1 = -1.71*y1 + 0.43*y2 + 8.32*y3 + 0.0007
+  dy2 = 1.71*y1 - 8.75*y2
+  dy3 = -10.03*y3 + 0.43*y4 + 0.035*y5
+  dy4 = 8.32*y2 + 1.71*y3 - 1.12*y4
+  dy5 = -1.745*y5 + 0.43*y6 + 0.43*y7
+  dy6 = -280.0*y6*y8 + 0.69*y4 + 1.71*y5 -
+           0.43*y6 + 0.69*y7
+  dy7 = 280.0*y6*y8 - 1.81*y7
+  dy8 = -280.0*y6*y8 + 1.81*y7
+end
+
+u0 = zeros(8)
+u0[1] = 1
+u0[8] = 0.0057
+prob = ODEProblem(f,u0,(0.0,321.8122))
+
+sol = solve(prob,Rodas5(),abstol=1/10^14,reltol=1/10^14)
+test_sol = TestSolution(sol)
+
+abstols = 1.0 ./ 10.0 .^ (5:8)
+reltols = 1.0 ./ 10.0 .^ (1:4);
+setups = [Dict(:alg=>Rosenbrock23()),
+          Dict(:alg=>TRBDF2()),
+          Dict(:alg=>rodas()),
+          Dict(:alg=>radau()),
+          Dict(:alg=>RadauIIA5()),
+          Dict(:alg=>MATLABDiffEq.ode23s()),
+          Dict(:alg=>MATLABDiffEq.ode15s())
+          ]
+
+wp = WorkPrecisionSet(prob,abstols,reltols,setups;
+                      save_everystep=false,appxsol=test_sol,maxiters=Int(1e5),numruns=100)
+plot(wp,title="Stiff 2: Hires")
+savefig("benchmark4.png")
+```
+
+![benchmark4](https://user-images.githubusercontent.com/1814174/69478409-019a6900-0dc0-11ea-9fce-c11a5e4de9a4.png)
+
 
 This shows that being able to run MATLAB ODE algorithms with MATLAB functions
 is cute, but does not really have a practical use due to MATLAB's lack of