@@ -113,16 +113,17 @@ julia> @time MATLABDiffEq.eval_string("[t,u] = $(algstr)(f,tspan,u0,options);")
113113The following benchmarks demonstrate a ** 100x performance advantage for the
114114pure-Julia methods over the MATLAB ODE solvers** across a range of stiff and
115115non-stiff ODEs. These were ran with Julia 1.2, MATLAB 2019B, deSolve 1.2.5,
116- and SciPy 1.3.1 after verifying negligible overhead on interop.
116+ and SciPy 1.3.1 after verifying negligible overhead on interop. Note that the
117+ MATLAB solvers do outperform that of Python and R.
118+
119+ #### Non-Stiff Problem 1: Lotka-Volterra
117120
118121``` julia
119122using ParameterizedFunctions, MATLABDiffEq, OrdinaryDiffEq, ODEInterface,
120- ODEInterfaceDiffEq, Plots, Sundials
123+ ODEInterfaceDiffEq, Plots, Sundials, SciPyDiffEq, deSolveDiffEq
121124using DiffEqDevTools
122125using LinearAlgebra
123126
124- # # Non-Stiff Problem 1: Lotka-Volterra
125-
126127f = @ode_def_bare LotkaVolterra begin
127128 dx = a* x - b* x* y
128129 dy = - c* y + d* x* y
@@ -136,27 +137,48 @@ test_sol = TestSolution(sol)
136137
137138setups = [Dict (:alg => DP5 ())
138139 Dict (:alg => dopri5 ())
139- Dict (:alg => BS5 ())
140140 Dict (:alg => Tsit5 ())
141- Dict (:alg => Vern6 ())
142141 Dict (:alg => Vern7 ())
143142 Dict (:alg => MATLABDiffEq. ode45 ())
144143 Dict (:alg => MATLABDiffEq. ode113 ())
144+ Dict (:alg => SciPyDiffEq. RK45 ())
145+ Dict (:alg => SciPyDiffEq. LSODA ())
146+ Dict (:alg => deSolveDiffEq. lsoda ())
147+ Dict (:alg => deSolveDiffEq. ode45 ())
145148 Dict (:alg => CVODE_Adams ())
146149 ]
150+
151+ names = [
152+ " Julia: DP5"
153+ " Hairer: dopri5"
154+ " Julia: Tsit5"
155+ " Julia: Vern7"
156+ " MATLAB: ode45"
157+ " MATLAB: ode113"
158+ " SciPy: RK45"
159+ " SciPy: LSODA"
160+ " deSolve: lsoda"
161+ " deSolve: ode45"
162+ " Sundials: Adams"
163+ ]
164+
147165abstols = 1.0 ./ 10.0 .^ (6 : 13 )
148166reltols = 1.0 ./ 10.0 .^ (3 : 10 )
149- wp = WorkPrecisionSet (prob,abstols,reltols,setups;appxsol= test_sol,dense= false ,
167+ wp = WorkPrecisionSet (prob,abstols,reltols,setups;
168+ names = names,
169+ appxsol= test_sol,dense= false ,
150170 save_everystep= false ,numruns= 100 ,maxiters= 10000000 ,
151171 timeseries_errors= false ,verbose= false )
152172plot (wp,title= " Non-stiff 1: Lotka-Volterra" )
153173savefig (" benchmark1.png" )
154174```
155175
156- ![ benchmark1] ( https://user-images.githubusercontent.com/1814174/69478405-f9dac480-0dbf-11ea-8f8e-5572b86d2a97.png )
176+ ![ benchmark1] ( https://user-images.githubusercontent.com/1814174/69487806-157cb400-0e2e-11ea-876f-c519aed013c0.png )
177+
178+ #### Non-Stiff Problem 2: Rigid Body
157179
158180``` julia
159- # # Non-Stiff Problem 2: Rigid Body
181+
160182
161183f = @ode_def_bare RigidBodyBench begin
162184 dy1 = - 2 * y2* y3
@@ -169,61 +191,151 @@ test_sol = TestSolution(sol)
169191
170192setups = [Dict (:alg => DP5 ())
171193 Dict (:alg => dopri5 ())
172- Dict (:alg => BS5 ())
173194 Dict (:alg => Tsit5 ())
174- Dict (:alg => Vern6 ())
175195 Dict (:alg => Vern7 ())
176196 Dict (:alg => MATLABDiffEq. ode45 ())
177197 Dict (:alg => MATLABDiffEq. ode113 ())
198+ Dict (:alg => SciPyDiffEq. RK45 ())
199+ Dict (:alg => SciPyDiffEq. LSODA ())
200+ Dict (:alg => deSolveDiffEq. lsoda ())
201+ Dict (:alg => deSolveDiffEq. ode45 ())
178202 Dict (:alg => CVODE_Adams ())
179203 ]
204+
205+ names = [
206+ " Julia: DP5"
207+ " Hairer: dopri5"
208+ " Julia: Tsit5"
209+ " Julia: Vern7"
210+ " MATLAB: ode45"
211+ " MATLAB: ode113"
212+ " SciPy: RK45"
213+ " SciPy: LSODA"
214+ " deSolve: lsoda"
215+ " deSolve: ode45"
216+ " Sundials: Adams"
217+ ]
218+
180219abstols = 1.0 ./ 10.0 .^ (6 : 13 )
181220reltols = 1.0 ./ 10.0 .^ (3 : 10 )
182- wp = WorkPrecisionSet (prob,abstols,reltols,setups;appxsol= test_sol,dense= false ,
221+ wp = WorkPrecisionSet (prob,abstols,reltols,setups;
222+ names = names,
223+ appxsol= test_sol,dense= false ,
183224 save_everystep= false ,numruns= 100 ,maxiters= 10000000 ,
184225 timeseries_errors= false ,verbose= false )
185226plot (wp,title= " Non-stiff 2: Rigid-Body" )
186227savefig (" benchmark2.png" )
187228```
188229
189- ![ benchmark2] ( https://user-images.githubusercontent.com/1814174/69478406-fe9f7880-0dbf -11ea-8f15-a0a0f3e8717f .png )
230+ ![ benchmark2] ( https://user-images.githubusercontent.com/1814174/69487808-17467780-0e2e -11ea-9db2-324d4e319d07 .png )
190231
232+ #### Stiff Problem 1: ROBER Shorter and Simpler for SciPy
191233
192234``` julia
193- # # Stiff Problem 1: ROBER
235+
194236
195237rober = @ode_def begin
196238 dy₁ = - k₁* y₁+ k₃* y₂* y₃
197239 dy₂ = k₁* y₁- k₂* y₂^ 2 - k₃* y₂* y₃
198240 dy₃ = k₂* y₂^ 2
199241end k₁ k₂ k₃
242+ prob = ODEProblem (rober,[1.0 ,0.0 ,0.0 ],(0.0 ,1e3 ),[0.04 ,3e7 ,1e4 ])
243+ sol = solve (prob,CVODE_BDF (),abstol= 1 / 10 ^ 14 ,reltol= 1 / 10 ^ 14 )
244+ test_sol = TestSolution (sol)
245+
246+ abstols = 1.0 ./ 10.0 .^ (7 : 8 )
247+ reltols = 1.0 ./ 10.0 .^ (3 : 4 );
248+
249+ setups = [Dict (:alg => Rosenbrock23 ())
250+ Dict (:alg => TRBDF2 ())
251+ Dict (:alg => RadauIIA5 ())
252+ Dict (:alg => rodas ())
253+ Dict (:alg => radau ())
254+ Dict (:alg => MATLABDiffEq. ode23s ())
255+ Dict (:alg => MATLABDiffEq. ode15s ())
256+ Dict (:alg => SciPyDiffEq. LSODA ())
257+ Dict (:alg => SciPyDiffEq. BDF ())
258+ Dict (:alg => deSolveDiffEq. lsoda ())
259+ ]
260+
261+ names = [
262+ " Julia: Rosenbrock23"
263+ " Julia: TRBDF2"
264+ " Julia: radau"
265+ " Hairer: rodas"
266+ " Hairer: radau"
267+ " MATLAB: ode23s"
268+ " MATLAB: ode15s"
269+ " SciPy: LSODA"
270+ " SciPy: BDF"
271+ " deSolve: lsoda"
272+ ]
273+
274+ wp = WorkPrecisionSet (prob,abstols,reltols,setups;
275+ names = names,print_names = true ,
276+ dense= false ,verbose = false ,
277+ save_everystep= false ,appxsol= test_sol,
278+ maxiters= Int (1e5 ))
279+ plot (wp,title= " Stiff 1: ROBER Short" )
280+ savefig (" benchmark31.png" )
281+ ```
282+
283+ ![ benchmark31] ( https://user-images.githubusercontent.com/1814174/69501071-4b25a980-0ecf-11ea-99d1-b7274491498e.png )
284+
285+ #### Rober Standard length
286+
287+ ** SciPy Omitted due to failures at higher tolerances and because it's too slow to finish in a day!**
288+ See note below.
289+
290+ ``` julia
291+
292+
200293prob = ODEProblem (rober,[1.0 ,0.0 ,0.0 ],(0.0 ,1e5 ),[0.04 ,3e7 ,1e4 ])
201294sol = solve (prob,CVODE_BDF (),abstol= 1 / 10 ^ 14 ,reltol= 1 / 10 ^ 14 )
202295test_sol = TestSolution (sol)
203296
204297abstols = 1.0 ./ 10.0 .^ (5 : 8 )
205298reltols = 1.0 ./ 10.0 .^ (1 : 4 );
206- setups = [Dict (:alg => Rosenbrock23 ()),
207- Dict (:alg => TRBDF2 ()),
208- Dict (:alg => rodas ()),
209- Dict (:alg => radau ()),
210- Dict (:alg => RadauIIA5 ()),
211- Dict (:alg => MATLABDiffEq. ode23s ()),
299+
300+ setups = [Dict (:alg => Rosenbrock23 ())
301+ Dict (:alg => TRBDF2 ())
302+ Dict (:alg => RadauIIA5 ())
303+ Dict (:alg => rodas ())
304+ Dict (:alg => radau ())
305+ Dict (:alg => MATLABDiffEq. ode23s ())
212306 Dict (:alg => MATLABDiffEq. ode15s ())
307+ # Dict(:alg=>SciPyDiffEq.LSODA())
308+ # Dict(:alg=>SciPyDiffEq.BDF())
309+ Dict (:alg => deSolveDiffEq. lsoda ())
213310 ]
214311
312+ names = [
313+ " Julia: Rosenbrock23"
314+ " Julia: TRBDF2"
315+ " Julia: radau"
316+ " Hairer: rodas"
317+ " Hairer: radau"
318+ " MATLAB: ode23s"
319+ " MATLAB: ode15s"
320+ # "SciPy: LSODA"
321+ # "SciPy: BDF"
322+ " deSolve: lsoda"
323+ ]
324+
215325wp = WorkPrecisionSet (prob,abstols,reltols,setups;
216- save_everystep= false ,appxsol= test_sol,maxiters= Int (1e5 ),numruns= 100 )
217- plot (wp,title= " Stiff 1: ROBER" )
218- savefig (" benchmark3.png" )
326+ names = names,print_names = true ,
327+ dense= false ,verbose = false ,
328+ save_everystep= false ,appxsol= test_sol,
329+ maxiters= Int (1e5 ))
330+ plot (wp,title= " Stiff 1: ROBER Standard Length" )
331+ savefig (" benchmark32.png" )
219332```
220333
221- ![ benchmark3 ] ( https://user-images.githubusercontent.com/1814174/69478407-ffd0a580-0dbf -11ea-9311-909bd3938b42 .png )
334+ ![ benchmark32 ] ( https://user-images.githubusercontent.com/1814174/69501072-4b25a980-0ecf -11ea-82dd-47aa566ecbc2 .png )
222335
336+ #### Stiff Problem 2: HIRES
223337
224338``` julia
225- # # Stiff Problem 2: HIRES
226-
227339f = @ode_def Hires begin
228340 dy1 = - 1.71 * y1 + 0.43 * y2 + 8.32 * y3 + 0.0007
229341 dy2 = 1.71 * y1 - 8.75 * y2
@@ -246,26 +358,37 @@ test_sol = TestSolution(sol)
246358
247359abstols = 1.0 ./ 10.0 .^ (5 : 8 )
248360reltols = 1.0 ./ 10.0 .^ (1 : 4 );
249- setups = [Dict (:alg => Rosenbrock23 ()),
250- Dict (:alg => TRBDF2 ()),
251- Dict (:alg => rodas ()),
252- Dict (:alg => radau ()),
253- Dict (:alg => RadauIIA5 ()),
254- Dict (:alg => MATLABDiffEq. ode23s ()),
361+ setups = [Dict (:alg => Rosenbrock23 ())
362+ Dict (:alg => TRBDF2 ())
363+ Dict (:alg => RadauIIA5 ())
364+ Dict (:alg => rodas ())
365+ Dict (:alg => radau ())
366+ Dict (:alg => MATLABDiffEq. ode23s ())
255367 Dict (:alg => MATLABDiffEq. ode15s ())
368+ Dict (:alg => SciPyDiffEq. LSODA ())
369+ Dict (:alg => SciPyDiffEq. BDF ())
370+ Dict (:alg => deSolveDiffEq. lsoda ())
256371 ]
257372
373+ names = [
374+ " Julia: Rosenbrock23"
375+ " Julia: TRBDF2"
376+ " Julia: radau"
377+ " Hairer: rodas"
378+ " Hairer: radau"
379+ " MATLAB: ode23s"
380+ " MATLAB: ode15s"
381+ " SciPy: LSODA"
382+ " SciPy: BDF"
383+ " deSolve: lsoda"
384+ ]
385+
258386wp = WorkPrecisionSet (prob,abstols,reltols,setups;
259- save_everystep= false ,appxsol= test_sol,maxiters= Int (1e5 ),numruns= 100 )
387+ names = names,print_names = true ,
388+ save_everystep= false ,appxsol= test_sol,
389+ maxiters= Int (1e5 ),numruns= 100 )
260390plot (wp,title= " Stiff 2: Hires" )
261391savefig (" benchmark4.png" )
262392```
263393
264- ![ benchmark4] ( https://user-images.githubusercontent.com/1814174/69478409-019a6900-0dc0-11ea-9fce-c11a5e4de9a4.png )
265-
266-
267- Together this demonstrates that the algorithms like ` ode45 ` and ` ode15s ` are
268- not competitive performance-wise against the Fortran and Julia methods.
269- This shows that being able to run MATLAB ODE algorithms with MATLAB functions
270- is cute, but does not really have a practical use due to MATLAB's lack of
271- performance (and its [ pass by copy] ( https://www.mathworks.com/matlabcentral/answers/152-can-matlab-pass-by-reference ) for functions).
394+ ![ benchmark4] ( https://user-images.githubusercontent.com/1814174/69501114-bec7b680-0ecf-11ea-9095-7b7f2e98d514.png )
0 commit comments