@@ -37,24 +37,24 @@ but are not exported. Thus to use `ode45`, you would specify the algorithm as
37
37
using MATLABDiffEq, ParameterizedFunctions
38
38
39
39
f = @ode_def LotkaVolterra begin
40
- dx = 1.5 x - x* y
41
- dy = - 3 y + x* y
40
+ dx = 1.5 x - x* y
41
+ dy = - 3 y + x* y
42
42
end
43
43
44
- tspan = (0.0 ,10.0 )
45
- u0 = [1.0 ,1.0 ]
46
- prob = ODEProblem (f,u0,tspan)
47
- sol = solve (prob,MATLABDiffEq. ode45 ())
44
+ tspan = (0.0 , 10.0 )
45
+ u0 = [1.0 , 1.0 ]
46
+ prob = ODEProblem (f, u0, tspan)
47
+ sol = solve (prob, MATLABDiffEq. ode45 ())
48
48
49
- function lorenz (du,u,p, t)
50
- du[1 ] = 10.0 (u[2 ]- u[1 ])
51
- du[2 ] = u[1 ]* (28.0 - u[3 ]) - u[2 ]
52
- du[3 ] = u[1 ]* u[2 ] - (8 / 3 )* u[3 ]
49
+ function lorenz (du, u, p, t)
50
+ du[1 ] = 10.0 (u[2 ]- u[1 ])
51
+ du[2 ] = u[1 ]* (28.0 - u[3 ]) - u[2 ]
52
+ du[3 ] = u[1 ]* u[2 ] - (8 / 3 )* u[3 ]
53
53
end
54
- u0 = [1.0 ;0.0 ;0.0 ]
55
- tspan = (0.0 ,100.0 )
56
- prob = ODEProblem (lorenz,u0,tspan)
57
- sol = solve (prob,MATLABDiffEq. ode45 ())
54
+ u0 = [1.0 ; 0.0 ; 0.0 ]
55
+ tspan = (0.0 , 100.0 )
56
+ prob = ODEProblem (lorenz, u0, tspan)
57
+ sol = solve (prob, MATLABDiffEq. ode45 ())
58
58
```
59
59
60
60
## Measuring Overhead
@@ -82,28 +82,28 @@ Generally, for long enough problems the overhead is minimal. Example:
82
82
``` julia
83
83
using DiffEqBase, ParameterizedFunctions, MATLABDiffEq
84
84
f = @ode_def_bare RigidBodyBench begin
85
- dy1 = - 2 * y2* y3
86
- dy2 = 1.25 * y1* y3
87
- dy3 = - 0.5 * y1* y2 + 0.25 * sin (t)^ 2
85
+ dy1 = - 2 * y2* y3
86
+ dy2 = 1.25 * y1* y3
87
+ dy3 = - 0.5 * y1* y2 + 0.25 * sin (t)^ 2
88
88
end
89
- prob = ODEProblem (f,[1.0 ;0.0 ;0.9 ],(0.0 ,100.0 ))
89
+ prob = ODEProblem (f, [1.0 ; 0.0 ; 0.9 ], (0.0 , 100.0 ))
90
90
alg = MATLABDiffEq. ode45 ()
91
91
algstr = string (typeof (alg). name. name)
92
92
```
93
93
94
94
For this, we get the following:
95
95
96
96
``` julia
97
- julia> @time sol = solve (prob,alg);
97
+ julia> @time sol = solve (prob, alg);
98
98
0.063918 seconds (38.84 k allocations: 1.556 MB)
99
99
100
- julia> @time sol = solve (prob,alg);
100
+ julia> @time sol = solve (prob, alg);
101
101
0.062600 seconds (38.84 k allocations: 1.556 MB)
102
102
103
- julia> @time sol = solve (prob,alg);
103
+ julia> @time sol = solve (prob, alg);
104
104
0.061438 seconds (38.84 k allocations: 1.556 MB)
105
105
106
- julia> @time sol = solve (prob,alg);
106
+ julia> @time sol = solve (prob, alg);
107
107
0.065460 seconds (38.84 k allocations: 1.556 MB)
108
108
109
109
julia> @time MATLABDiffEq. eval_string (" [t,u] = $(algstr) (diffeqf,tspan,u0,options);" )
@@ -131,14 +131,14 @@ MATLAB solvers do outperform that of Python and R.
131
131
132
132
``` julia
133
133
f = @ode_def_bare LotkaVolterra begin
134
- dx = a* x - b* x* y
135
- dy = - c* y + d* x* y
134
+ dx = a* x - b* x* y
135
+ dy = - c* y + d* x* y
136
136
end a b c d
137
- p = [1.5 ,1 , 3 , 1 ]
138
- tspan = (0.0 ,10.0 )
139
- u0 = [1.0 ,1.0 ]
140
- prob = ODEProblem (f,u0,tspan,p)
141
- sol = solve (prob,Vern7 (),abstol= 1 / 10 ^ 14 ,reltol= 1 / 10 ^ 14 )
137
+ p = [1.5 , 1 , 3 , 1 ]
138
+ tspan = (0.0 , 10.0 )
139
+ u0 = [1.0 , 1.0 ]
140
+ prob = ODEProblem (f, u0, tspan, p)
141
+ sol = solve (prob, Vern7 (), abstol = 1 / 10 ^ 14 , reltol = 1 / 10 ^ 14 )
142
142
test_sol = TestSolution (sol)
143
143
144
144
setups = [Dict (:alg => DP5 ())
@@ -152,32 +152,29 @@ setups = [Dict(:alg=>DP5())
152
152
Dict (:alg => SciPyDiffEq. odeint ())
153
153
Dict (:alg => deSolveDiffEq. lsoda ())
154
154
Dict (:alg => deSolveDiffEq. ode45 ())
155
- Dict (:alg => CVODE_Adams ())
156
- ]
157
-
158
- names = [
159
- " Julia: DP5"
160
- " Hairer: dopri5"
161
- " Julia: Tsit5"
162
- " Julia: Vern7"
163
- " MATLAB: ode45"
164
- " MATLAB: ode113"
165
- " SciPy: RK45"
166
- " SciPy: LSODA"
167
- " SciPy: odeint"
168
- " deSolve: lsoda"
169
- " deSolve: ode45"
170
- " Sundials: Adams"
171
- ]
155
+ Dict (:alg => CVODE_Adams ())]
156
+
157
+ names = [" Julia: DP5"
158
+ " Hairer: dopri5"
159
+ " Julia: Tsit5"
160
+ " Julia: Vern7"
161
+ " MATLAB: ode45"
162
+ " MATLAB: ode113"
163
+ " SciPy: RK45"
164
+ " SciPy: LSODA"
165
+ " SciPy: odeint"
166
+ " deSolve: lsoda"
167
+ " deSolve: ode45"
168
+ " Sundials: Adams" ]
172
169
173
170
abstols = 1.0 ./ 10.0 .^ (6 : 13 )
174
171
reltols = 1.0 ./ 10.0 .^ (3 : 10 )
175
- wp = WorkPrecisionSet (prob,abstols,reltols,setups;
176
- names = names,
177
- appxsol= test_sol,dense= false ,
178
- save_everystep= false ,numruns= 100 ,maxiters= 10000000 ,
179
- timeseries_errors= false ,verbose= false )
180
- plot (wp,title= " Non-stiff 1: Lotka-Volterra" )
172
+ wp = WorkPrecisionSet (prob, abstols, reltols, setups;
173
+ names = names,
174
+ appxsol = test_sol, dense = false ,
175
+ save_everystep = false , numruns = 100 , maxiters = 10000000 ,
176
+ timeseries_errors = false , verbose = false )
177
+ plot (wp, title = " Non-stiff 1: Lotka-Volterra" )
181
178
```
182
179
183
180
![ ] ( https://user-images.githubusercontent.com/1814174/71537082-ef42ac00-28e4-11ea-9acc-67dfd9990221.png )
@@ -186,12 +183,12 @@ plot(wp,title="Non-stiff 1: Lotka-Volterra")
186
183
187
184
``` julia
188
185
f = @ode_def_bare RigidBodyBench begin
189
- dy1 = - 2 * y2* y3
190
- dy2 = 1.25 * y1* y3
191
- dy3 = - 0.5 * y1* y2 + 0.25 * sin (t)^ 2
186
+ dy1 = - 2 * y2* y3
187
+ dy2 = 1.25 * y1* y3
188
+ dy3 = - 0.5 * y1* y2 + 0.25 * sin (t)^ 2
192
189
end
193
- prob = ODEProblem (f,[1.0 ;0.0 ;0.9 ],(0.0 ,100.0 ))
194
- sol = solve (prob,Vern7 (),abstol= 1 / 10 ^ 14 ,reltol= 1 / 10 ^ 14 )
190
+ prob = ODEProblem (f, [1.0 ; 0.0 ; 0.9 ], (0.0 , 100.0 ))
191
+ sol = solve (prob, Vern7 (), abstol = 1 / 10 ^ 14 , reltol = 1 / 10 ^ 14 )
195
192
test_sol = TestSolution (sol)
196
193
197
194
setups = [Dict (:alg => DP5 ())
@@ -205,32 +202,29 @@ setups = [Dict(:alg=>DP5())
205
202
Dict (:alg => SciPyDiffEq. odeint ())
206
203
Dict (:alg => deSolveDiffEq. lsoda ())
207
204
Dict (:alg => deSolveDiffEq. ode45 ())
208
- Dict (:alg => CVODE_Adams ())
209
- ]
210
-
211
- names = [
212
- " Julia: DP5"
213
- " Hairer: dopri5"
214
- " Julia: Tsit5"
215
- " Julia: Vern7"
216
- " MATLAB: ode45"
217
- " MATLAB: ode113"
218
- " SciPy: RK45"
219
- " SciPy: LSODA"
220
- " SciPy: odeint"
221
- " deSolve: lsoda"
222
- " deSolve: ode45"
223
- " Sundials: Adams"
224
- ]
205
+ Dict (:alg => CVODE_Adams ())]
206
+
207
+ names = [" Julia: DP5"
208
+ " Hairer: dopri5"
209
+ " Julia: Tsit5"
210
+ " Julia: Vern7"
211
+ " MATLAB: ode45"
212
+ " MATLAB: ode113"
213
+ " SciPy: RK45"
214
+ " SciPy: LSODA"
215
+ " SciPy: odeint"
216
+ " deSolve: lsoda"
217
+ " deSolve: ode45"
218
+ " Sundials: Adams" ]
225
219
226
220
abstols = 1.0 ./ 10.0 .^ (6 : 13 )
227
221
reltols = 1.0 ./ 10.0 .^ (3 : 10 )
228
- wp = WorkPrecisionSet (prob,abstols,reltols,setups;
229
- names = names,
230
- appxsol= test_sol,dense= false ,
231
- save_everystep= false ,numruns= 100 ,maxiters= 10000000 ,
232
- timeseries_errors= false ,verbose= false )
233
- plot (wp,title= " Non-stiff 2: Rigid-Body" )
222
+ wp = WorkPrecisionSet (prob, abstols, reltols, setups;
223
+ names = names,
224
+ appxsol = test_sol, dense = false ,
225
+ save_everystep = false , numruns = 100 , maxiters = 10000000 ,
226
+ timeseries_errors = false , verbose = false )
227
+ plot (wp, title = " Non-stiff 2: Rigid-Body" )
234
228
```
235
229
236
230
![ ] ( https://user-images.githubusercontent.com/1814174/71537083-ef42ac00-28e4-11ea-8dc7-a5dca0518baf.png )
@@ -239,12 +233,12 @@ plot(wp,title="Non-stiff 2: Rigid-Body")
239
233
240
234
``` julia
241
235
rober = @ode_def begin
242
- dy₁ = - k₁* y₁+ k₃* y₂* y₃
243
- dy₂ = k₁* y₁- k₂* y₂^ 2 - k₃* y₂* y₃
244
- dy₃ = k₂* y₂^ 2
236
+ dy₁ = - k₁* y₁+ k₃* y₂* y₃
237
+ dy₂ = k₁* y₁- k₂* y₂^ 2 - k₃* y₂* y₃
238
+ dy₃ = k₂* y₂^ 2
245
239
end k₁ k₂ k₃
246
- prob = ODEProblem (rober,[1.0 ,0.0 ,0.0 ],(0.0 ,1e5 ),[0.04 ,3e7 ,1e4 ])
247
- sol = solve (prob,CVODE_BDF (),abstol= 1 / 10 ^ 14 ,reltol= 1 / 10 ^ 14 )
240
+ prob = ODEProblem (rober, [1.0 , 0.0 , 0.0 ], (0.0 , 1e5 ), [0.04 , 3e7 , 1e4 ])
241
+ sol = solve (prob, CVODE_BDF (), abstol = 1 / 10 ^ 14 , reltol = 1 / 10 ^ 14 )
248
242
test_sol = TestSolution (sol)
249
243
250
244
abstols = 1.0 ./ 10.0 .^ (7 : 8 )
@@ -261,30 +255,27 @@ setups = [Dict(:alg=>Rosenbrock23())
261
255
Dict (:alg => SciPyDiffEq. BDF ())
262
256
Dict (:alg => SciPyDiffEq. odeint ())
263
257
Dict (:alg => deSolveDiffEq. lsoda ())
264
- Dict (:alg => CVODE_BDF ())
265
- ]
266
-
267
- names = [
268
- " Julia: Rosenbrock23"
269
- " Julia: TRBDF2"
270
- " Julia: radau"
271
- " Hairer: rodas"
272
- " Hairer: radau"
273
- " MATLAB: ode23s"
274
- " MATLAB: ode15s"
275
- " SciPy: LSODA"
276
- " SciPy: BDF"
277
- " SciPy: odeint"
278
- " deSolve: lsoda"
279
- " Sundials: CVODE"
280
- ]
281
-
282
- wp = WorkPrecisionSet (prob,abstols,reltols,setups;
283
- names = names,print_names = true ,
284
- dense= false ,verbose = false ,
285
- save_everystep= false ,appxsol= test_sol,
286
- maxiters= Int (1e5 ))
287
- plot (wp,title= " Stiff 1: ROBER" , legend= :topleft )
258
+ Dict (:alg => CVODE_BDF ())]
259
+
260
+ names = [" Julia: Rosenbrock23"
261
+ " Julia: TRBDF2"
262
+ " Julia: radau"
263
+ " Hairer: rodas"
264
+ " Hairer: radau"
265
+ " MATLAB: ode23s"
266
+ " MATLAB: ode15s"
267
+ " SciPy: LSODA"
268
+ " SciPy: BDF"
269
+ " SciPy: odeint"
270
+ " deSolve: lsoda"
271
+ " Sundials: CVODE" ]
272
+
273
+ wp = WorkPrecisionSet (prob, abstols, reltols, setups;
274
+ names = names, print_names = true ,
275
+ dense = false , verbose = false ,
276
+ save_everystep = false , appxsol = test_sol,
277
+ maxiters = Int (1e5 ))
278
+ plot (wp, title = " Stiff 1: ROBER" , legend = :topleft )
288
279
```
289
280
290
281
![ ] ( https://user-images.githubusercontent.com/1814174/71537080-ef42ac00-28e4-11ea-9abd-37631cd18ad9.png )
@@ -293,23 +284,23 @@ plot(wp,title="Stiff 1: ROBER", legend=:topleft)
293
284
294
285
``` julia
295
286
f = @ode_def Hires begin
296
- dy1 = - 1.71 * y1 + 0.43 * y2 + 8.32 * y3 + 0.0007
297
- dy2 = 1.71 * y1 - 8.75 * y2
298
- dy3 = - 10.03 * y3 + 0.43 * y4 + 0.035 * y5
299
- dy4 = 8.32 * y2 + 1.71 * y3 - 1.12 * y4
300
- dy5 = - 1.745 * y5 + 0.43 * y6 + 0.43 * y7
301
- dy6 = - 280.0 * y6* y8 + 0.69 * y4 + 1.71 * y5 -
302
- 0.43 * y6 + 0.69 * y7
303
- dy7 = 280.0 * y6* y8 - 1.81 * y7
304
- dy8 = - 280.0 * y6* y8 + 1.81 * y7
287
+ dy1 = - 1.71 * y1 + 0.43 * y2 + 8.32 * y3 + 0.0007
288
+ dy2 = 1.71 * y1 - 8.75 * y2
289
+ dy3 = - 10.03 * y3 + 0.43 * y4 + 0.035 * y5
290
+ dy4 = 8.32 * y2 + 1.71 * y3 - 1.12 * y4
291
+ dy5 = - 1.745 * y5 + 0.43 * y6 + 0.43 * y7
292
+ dy6 = - 280.0 * y6* y8 + 0.69 * y4 + 1.71 * y5 -
293
+ 0.43 * y6 + 0.69 * y7
294
+ dy7 = 280.0 * y6* y8 - 1.81 * y7
295
+ dy8 = - 280.0 * y6* y8 + 1.81 * y7
305
296
end
306
297
307
298
u0 = zeros (8 )
308
299
u0[1 ] = 1
309
300
u0[8 ] = 0.0057
310
- prob = ODEProblem (f,u0,(0.0 ,321.8122 ))
301
+ prob = ODEProblem (f, u0, (0.0 , 321.8122 ))
311
302
312
- sol = solve (prob,Rodas5 (),abstol= 1 / 10 ^ 14 ,reltol= 1 / 10 ^ 14 )
303
+ sol = solve (prob, Rodas5 (), abstol = 1 / 10 ^ 14 , reltol = 1 / 10 ^ 14 )
313
304
test_sol = TestSolution (sol)
314
305
315
306
abstols = 1.0 ./ 10.0 .^ (5 : 8 )
@@ -326,29 +317,26 @@ setups = [Dict(:alg=>Rosenbrock23())
326
317
Dict (:alg => SciPyDiffEq. BDF ())
327
318
Dict (:alg => SciPyDiffEq. odeint ())
328
319
Dict (:alg => deSolveDiffEq. lsoda ())
329
- Dict (:alg => CVODE_BDF ())
330
- ]
331
-
332
- names = [
333
- " Julia: Rosenbrock23"
334
- " Julia: TRBDF2"
335
- " Julia: radau"
336
- " Hairer: rodas"
337
- " Hairer: radau"
338
- " MATLAB: ode23s"
339
- " MATLAB: ode15s"
340
- " SciPy: LSODA"
341
- " SciPy: BDF"
342
- " SciPy: odeint"
343
- " deSolve: lsoda"
344
- " Sundials: CVODE"
345
- ]
346
-
347
- wp = WorkPrecisionSet (prob,abstols,reltols,setups;
348
- names = names,print_names = true ,
349
- save_everystep= false ,appxsol= test_sol,
350
- maxiters= Int (1e5 ),numruns= 100 )
351
- plot (wp,title= " Stiff 2: Hires" ,legend= :topleft )
320
+ Dict (:alg => CVODE_BDF ())]
321
+
322
+ names = [" Julia: Rosenbrock23"
323
+ " Julia: TRBDF2"
324
+ " Julia: radau"
325
+ " Hairer: rodas"
326
+ " Hairer: radau"
327
+ " MATLAB: ode23s"
328
+ " MATLAB: ode15s"
329
+ " SciPy: LSODA"
330
+ " SciPy: BDF"
331
+ " SciPy: odeint"
332
+ " deSolve: lsoda"
333
+ " Sundials: CVODE" ]
334
+
335
+ wp = WorkPrecisionSet (prob, abstols, reltols, setups;
336
+ names = names, print_names = true ,
337
+ save_everystep = false , appxsol = test_sol,
338
+ maxiters = Int (1e5 ), numruns = 100 )
339
+ plot (wp, title = " Stiff 2: Hires" , legend = :topleft )
352
340
```
353
341
354
342
![ ] ( https://user-images.githubusercontent.com/1814174/71537081-ef42ac00-28e4-11ea-950f-59c762ce9a69.png )
0 commit comments