merge issues

Nicholaswogan · Feb 24, 2022 · 1fe807f · 1fe807f
2 parents edb2192 + 6b37966
commit 1fe807f
Show file tree

Hide file tree

Showing 2 changed files with 33 additions and 15 deletions.
diff --git a/benchmark/README.md b/benchmark/README.md
@@ -14,8 +14,8 @@ Also, Scipy might be faster than numbalsoda for stiff problems with > ~500 ODEs
 | Scipy LSODA  | 595.029 ms  | 174x                   |
 | Scipy RK45   | 1744.073 ms | 510x                   |
 | Scipy DOP853 | 1052.442 ms | 308x                   |
-| Julia Tsit5  | 2.612 ms    | 0.76x                  |
-| Julia Vern8  | 1.803 ms    | 0.53x                  |
+| Julia Tsit5  | 1.340 ms    | 0.39x                  |
+| Julia Vern8  | 1.803 ms    | 0.21x                  |
 
 ## Stiff (Rober)
 

diff --git a/benchmark/benchmark.jl b/benchmark/benchmark.jl
@@ -1,5 +1,5 @@
 
-using OrdinaryDiffEq, StaticArrays, BenchmarkTools, Sundials, ModelingToolkit
+using OrdinaryDiffEq, StaticArrays, BenchmarkTools, LSODA, Sundials, ModelingToolkit
 
 println("Lorenz")
 # Lorenz
@@ -15,26 +15,27 @@ u0 = SA[1.0,0.0,0.0]
 p  = SA[10.0,28.0,8/3]
 tspan = (0.0,100.0)
 prob = ODEProblem(lorenz_static,u0,tspan,p)
-@btime solve(prob,Tsit5(),saveat=0.1,reltol=1.0e-8,abstol=1.0e-8) # 2.612 ms (30 allocations: 59.22 KiB)
-@btime solve(prob,Vern8(),saveat=0.1,reltol=1.0e-8,abstol=1.0e-8); # 1.803 ms (31 allocations: 65.70 KiB)
+@btime solve(prob,Tsit5(),saveat=0.1,reltol=1.0e-8,abstol=1.0e-8) # 1.340 ms (27 allocations: 64.55 KiB)
+@btime solve(prob,Vern8(),saveat=0.1,reltol=1.0e-8,abstol=1.0e-8); # 718.000 μs (28 allocations: 71.52 KiB)
 
 println("\nRober")
 # rober
 function rober(du,u,p,t)
     y₁,y₂,y₃ = u
     k₁,k₂,k₃ = p
     @inbounds begin
-    du[1] = -k₁*y₁+k₃*y₂*y₃
-    du[2] =  k₁*y₁-k₂*y₂^2-k₃*y₂*y₃
-    du[3] =  k₂*y₂^2
+        du[1] = -k₁*y₁+k₃*y₂*y₃
+        du[2] =  k₁*y₁-k₂*y₂^2-k₃*y₂*y₃
+        du[3] =  k₂*y₂^2
     end
     nothing
 end
-prob = ODEProblem(rober,[1.0,0.0,0.0],(0.0,1e5),[0.04,3e7,1e4]) 
+prob = ODEProblem(rober,[1.0,0.0,0.0],(0.0,1e5),[0.04,3e7,1e4])
 
-@btime solve(prob,Rodas5(),reltol=1.0e-8,abstol=1.0e-8, saveat = 1000) # 761.578 μs (633 allocations: 53.05 KiB)
-@btime solve(prob,TRBDF2(),reltol=1.0e-8,abstol=1.0e-8, saveat = 1000) # 5.765 ms (11415 allocations: 524.14 KiB)
-@btime solve(prob,CVODE_BDF(),reltol=1.0e-8,abstol=1.0e-8, saveat = 1000) # 1.115 ms (6954 allocations: 295.83 KiB)
+@btime solve(prob,Rodas5(),reltol=1.0e-8,abstol=1.0e-8, saveat = 1000) # 227.200 μs (1288 allocations: 52.56 KiB)
+@btime solve(prob,TRBDF2(),reltol=1.0e-8,abstol=1.0e-8, saveat = 1000) # 2.202 ms (3424 allocations: 194.41 KiB)
+@btime solve(prob,CVODE_BDF(),reltol=1.0e-8,abstol=1.0e-8, saveat = 1000) # 627.500 μs (7183 allocations: 214.33 KiB)
+@btime solve(prob,lsoda(),reltol=1.0e-8,abstol=1.0e-8, saveat = 1000) # 254.400 μs (2370 allocations: 127.09 KiB)
 
 # rober
 function rober_static(u,p,t)
@@ -46,8 +47,25 @@ function rober_static(u,p,t)
     SA[du1,du2,du3]
 end
 prob = ODEProblem{false}(rober_static,SA[1.0,0.0,0.0],(0.0,1e5),SA[0.04,3e7,1e4])
-@btime solve(prob,Rodas5(),reltol=1.0e-8,abstol=1.0e-8, saveat = 1000) # 172.953 μs (580 allocations: 53.36 KiB)
+@btime solve(prob,Rodas5(),reltol=1.0e-8,abstol=1.0e-8, saveat = 1000) # 42.300 μs (29 allocations: 9.75 KiB)
 
-prob2 = ODEProblem{false}(modelingtoolkitize(prob),SA[1.0,0.0,0.0],(0.0,1e5),SA[0.04,3e7,1e4],jac=true)
-@btime solve(prob2,Rodas5(), reltol=1.0e-8, abstol=1.0e-8, saveat = 1000); # 113.475 μs (400 allocations: 51.19 KiB)
+function rober_jac(u,p,t)
+  y₁,y₂,y₃ = u
+  k₁,k₂,k₃ = p
+  J11 = k₁ * -1
+  J21 = k₁
+  J31 = 0
+  J12 = y₃ * k₃
+  J22 = y₂ * k₂ * -2 + y₃ * k₃ * -1
+  J32 = y₂ * 2 * k₂
+  J13 = k₃ * y₂
+  J23 = k₃ * y₂ * -1
+  J33 = 0
+  SA[J11 J12 J13
+     J21 J22 J23
+     J31 J32 J33]
+end
 
+ff = ODEFunction(rober_static, jac=rober_jac)
+prob2 = ODEProblem{false}(ff,SA[1.0,0.0,0.0],(0.0,1e5),SA[0.04,3e7,1e4])
+@btime solve(prob2,Rodas5(), reltol=1.0e-8, abstol=1.0e-8, saveat = 1000); # 42.000 μs (30 allocations: 9.70 KiB)