@@ -113,16 +113,17 @@ julia> @time MATLABDiffEq.eval_string("[t,u] = $(algstr)(f,tspan,u0,options);")
113
113
The following benchmarks demonstrate a ** 100x performance advantage for the
114
114
pure-Julia methods over the MATLAB ODE solvers** across a range of stiff and
115
115
non-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
117
120
118
121
``` julia
119
122
using ParameterizedFunctions, MATLABDiffEq, OrdinaryDiffEq, ODEInterface,
120
- ODEInterfaceDiffEq, Plots, Sundials
123
+ ODEInterfaceDiffEq, Plots, Sundials, SciPyDiffEq, deSolveDiffEq
121
124
using DiffEqDevTools
122
125
using LinearAlgebra
123
126
124
- # # Non-Stiff Problem 1: Lotka-Volterra
125
-
126
127
f = @ode_def_bare LotkaVolterra begin
127
128
dx = a* x - b* x* y
128
129
dy = - c* y + d* x* y
@@ -136,27 +137,48 @@ test_sol = TestSolution(sol)
136
137
137
138
setups = [Dict (:alg => DP5 ())
138
139
Dict (:alg => dopri5 ())
139
- Dict (:alg => BS5 ())
140
140
Dict (:alg => Tsit5 ())
141
- Dict (:alg => Vern6 ())
142
141
Dict (:alg => Vern7 ())
143
142
Dict (:alg => MATLABDiffEq. ode45 ())
144
143
Dict (:alg => MATLABDiffEq. ode113 ())
144
+ Dict (:alg => SciPyDiffEq. RK45 ())
145
+ Dict (:alg => SciPyDiffEq. LSODA ())
146
+ Dict (:alg => deSolveDiffEq. lsoda ())
147
+ Dict (:alg => deSolveDiffEq. ode45 ())
145
148
Dict (:alg => CVODE_Adams ())
146
149
]
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
+
147
165
abstols = 1.0 ./ 10.0 .^ (6 : 13 )
148
166
reltols = 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 ,
150
170
save_everystep= false ,numruns= 100 ,maxiters= 10000000 ,
151
171
timeseries_errors= false ,verbose= false )
152
172
plot (wp,title= " Non-stiff 1: Lotka-Volterra" )
153
173
savefig (" benchmark1.png" )
154
174
```
155
175
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
157
179
158
180
``` julia
159
- # # Non-Stiff Problem 2: Rigid Body
181
+
160
182
161
183
f = @ode_def_bare RigidBodyBench begin
162
184
dy1 = - 2 * y2* y3
@@ -169,61 +191,151 @@ test_sol = TestSolution(sol)
169
191
170
192
setups = [Dict (:alg => DP5 ())
171
193
Dict (:alg => dopri5 ())
172
- Dict (:alg => BS5 ())
173
194
Dict (:alg => Tsit5 ())
174
- Dict (:alg => Vern6 ())
175
195
Dict (:alg => Vern7 ())
176
196
Dict (:alg => MATLABDiffEq. ode45 ())
177
197
Dict (:alg => MATLABDiffEq. ode113 ())
198
+ Dict (:alg => SciPyDiffEq. RK45 ())
199
+ Dict (:alg => SciPyDiffEq. LSODA ())
200
+ Dict (:alg => deSolveDiffEq. lsoda ())
201
+ Dict (:alg => deSolveDiffEq. ode45 ())
178
202
Dict (:alg => CVODE_Adams ())
179
203
]
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
+
180
219
abstols = 1.0 ./ 10.0 .^ (6 : 13 )
181
220
reltols = 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 ,
183
224
save_everystep= false ,numruns= 100 ,maxiters= 10000000 ,
184
225
timeseries_errors= false ,verbose= false )
185
226
plot (wp,title= " Non-stiff 2: Rigid-Body" )
186
227
savefig (" benchmark2.png" )
187
228
```
188
229
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 )
190
231
232
+ #### Stiff Problem 1: ROBER Shorter and Simpler for SciPy
191
233
192
234
``` julia
193
- # # Stiff Problem 1: ROBER
235
+
194
236
195
237
rober = @ode_def begin
196
238
dy₁ = - k₁* y₁+ k₃* y₂* y₃
197
239
dy₂ = k₁* y₁- k₂* y₂^ 2 - k₃* y₂* y₃
198
240
dy₃ = k₂* y₂^ 2
199
241
end 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
+
200
293
prob = ODEProblem (rober,[1.0 ,0.0 ,0.0 ],(0.0 ,1e5 ),[0.04 ,3e7 ,1e4 ])
201
294
sol = solve (prob,CVODE_BDF (),abstol= 1 / 10 ^ 14 ,reltol= 1 / 10 ^ 14 )
202
295
test_sol = TestSolution (sol)
203
296
204
297
abstols = 1.0 ./ 10.0 .^ (5 : 8 )
205
298
reltols = 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 ())
212
306
Dict (:alg => MATLABDiffEq. ode15s ())
307
+ # Dict(:alg=>SciPyDiffEq.LSODA())
308
+ # Dict(:alg=>SciPyDiffEq.BDF())
309
+ Dict (:alg => deSolveDiffEq. lsoda ())
213
310
]
214
311
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
+
215
325
wp = 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" )
219
332
```
220
333
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 )
222
335
336
+ #### Stiff Problem 2: HIRES
223
337
224
338
``` julia
225
- # # Stiff Problem 2: HIRES
226
-
227
339
f = @ode_def Hires begin
228
340
dy1 = - 1.71 * y1 + 0.43 * y2 + 8.32 * y3 + 0.0007
229
341
dy2 = 1.71 * y1 - 8.75 * y2
@@ -246,26 +358,37 @@ test_sol = TestSolution(sol)
246
358
247
359
abstols = 1.0 ./ 10.0 .^ (5 : 8 )
248
360
reltols = 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 ())
255
367
Dict (:alg => MATLABDiffEq. ode15s ())
368
+ Dict (:alg => SciPyDiffEq. LSODA ())
369
+ Dict (:alg => SciPyDiffEq. BDF ())
370
+ Dict (:alg => deSolveDiffEq. lsoda ())
256
371
]
257
372
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
+
258
386
wp = 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 )
260
390
plot (wp,title= " Stiff 2: Hires" )
261
391
savefig (" benchmark4.png" )
262
392
```
263
393
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