1
1
using ModelingToolkit
2
- using IfElse
3
- using Symbolics
2
+ using ModelingToolkit. Symbolics
4
3
using Symbolics: unwrap
5
4
using DiffEqBase, StaticArrays, LinearAlgebra
6
5
@@ -124,7 +123,7 @@ sd2sys = let
124
123
ODESystem (sd2eqs, t, name = :sd2 )
125
124
end
126
125
127
- sd2prob = ODEProblem {false} (sd2sys, [y[1 ] => 1.0 ; y[2 : 3 ] .=> 0.0 ], (0 , 40.0 ), dt = 1e-5 , cse = true )
126
+ sd2prob = ODEProblem {false} (sd2sys, [y[1 ] => 1.0 ; y[2 : 3 ] .=> 0.0 ], (0 , 40.0 ), dt = 1e-5 )
128
127
129
128
sd3sys = let
130
129
sd3eqs = [D (y[1 ]) ~ y[3 ] - 100 * (y[1 ] * y[2 ]),
@@ -135,7 +134,7 @@ sd3sys = let
135
134
ODESystem (sd3eqs, t, name = :sd3 )
136
135
end
137
136
138
- sd3prob = ODEProblem {false} (sd3sys, [y[1 : 2 ] .=> 1 ; y[3 : 4 ] .=> 0.0 ], (0 , 20.0 ), dt = 2.5e-5 , cse = true )
137
+ sd3prob = ODEProblem {false} (sd3sys, [y[1 : 2 ] .=> 1 ; y[3 : 4 ] .=> 0.0 ], (0 , 20.0 ), dt = 2.5e-5 )
139
138
140
139
sd4sys = let
141
140
sd4eqs = [D (y[1 ]) ~ - 0.013 y[1 ] - 1000 * (y[1 ] * y[3 ]),
@@ -145,7 +144,7 @@ sd4sys = let
145
144
ODESystem (sd4eqs, t, name = :sd4 )
146
145
end
147
146
148
- sd4prob = ODEProblem {false} (sd4sys, [y[1 : 2 ] .=> 1 ; y[3 ] => 0.0 ], (0 , 50.0 ), dt = 2.9e-4 , cse = true )
147
+ sd4prob = ODEProblem {false} (sd4sys, [y[1 : 2 ] .=> 1 ; y[3 ] => 0.0 ], (0 , 50.0 ), dt = 2.9e-4 )
149
148
150
149
sd5sys = let
151
150
sd5eqs = [D (y[1 ]) ~ 0.01 - (1 + (y[1 ] + 1000 ) * (y[1 ] + 1 )) * (0.01 + y[1 ] + y[2 ]),
@@ -154,7 +153,7 @@ sd5sys = let
154
153
ODESystem (sd5eqs, t, name = :sd5 )
155
154
end
156
155
157
- sd5prob = ODEProblem {false} (sd5sys, y[1 : 2 ] .=> 0.0 , (0 , 100.0 ), dt = 1e-4 , cse = true )
156
+ sd5prob = ODEProblem {false} (sd5sys, y[1 : 2 ] .=> 0.0 , (0 , 100.0 ), dt = 1e-4 )
158
157
159
158
sd6sys = let
160
159
sd6eqs = [D (y[1 ]) ~ - y[1 ] + 10 ^ 8 * y[3 ] * (1 - y[1 ]),
@@ -165,7 +164,7 @@ sd6sys = let
165
164
ODESystem (sd6eqs, t, name = :sd6 )
166
165
end
167
166
168
- sd6prob = ODEProblem {false} (sd6sys, [y[1 ] => 1.0 ; y[2 : 3 ] .=> 0.0 ], (0 , 1.0 ), dt = 3.3e-8 , cse = true )
167
+ sd6prob = ODEProblem {false} (sd6sys, [y[1 ] => 1.0 ; y[2 : 3 ] .=> 0.0 ], (0 , 1.0 ), dt = 3.3e-8 )
169
168
170
169
se1sys = let
171
170
Γ = 100
@@ -178,7 +177,7 @@ se1sys = let
178
177
ODESystem (se1eqs, t, name = :se1 )
179
178
end
180
179
181
- se1prob = ODEProblem {false} (se1sys, y .=> 0.0 , (0 , 1.0 ), dt = 6.8e-3 , cse = true )
180
+ se1prob = ODEProblem {false} (se1sys, y .=> 0.0 , (0 , 1.0 ), dt = 6.8e-3 )
182
181
183
182
se2sys = let
184
183
se2eqs = [D (y[1 ]) ~ y[2 ],
@@ -188,7 +187,7 @@ se2sys = let
188
187
ODESystem (se2eqs, t, name = :se2 )
189
188
end
190
189
191
- se2prob = ODEProblem {false} (se2sys, [y[1 ] => 2.0 , y[2 ] => 0.0 ], (0 , 1.0 ), dt = 1e-3 , cse = true )
190
+ se2prob = ODEProblem {false} (se2sys, [y[1 ] => 2.0 , y[2 ] => 0.0 ], (0 , 1.0 ), dt = 1e-3 )
192
191
193
192
se3sys = let
194
193
se3eqs = [D (y[1 ]) ~ - (55 + y[3 ]) * y[1 ] + 65 * y[2 ],
@@ -214,7 +213,7 @@ se4sys = let y = y
214
213
ODESystem (se4eqs, t, name = :se4 )
215
214
end
216
215
217
- se4prob = ODEProblem {false} (se4sys, [y[1 ] => 0.0 ; y[2 ] => - 2.0 ; y[3 : 4 ] .=> - 1.0 ], (0 , 1000.0 ), dt = 1e-3 , cse = true )
216
+ se4prob = ODEProblem {false} (se4sys, [y[1 ] => 0.0 ; y[2 ] => - 2.0 ; y[3 : 4 ] .=> - 1.0 ], (0 , 1000.0 ), dt = 1e-3 )
218
217
219
218
se5sys = let
220
219
se5eqs = [
@@ -227,7 +226,7 @@ se5sys = let
227
226
ODESystem (se5eqs, t, name = :se5 )
228
227
end
229
228
230
- se5prob = ODEProblem {false} (se5sys, [y[1 ] => 1.76e-3 ; y[2 : 4 ] .=> 0.0 ], (0 , 1000.0 ), dt = 1e-3 , cse = true )
229
+ se5prob = ODEProblem {false} (se5sys, [y[1 ] => 1.76e-3 ; y[2 : 4 ] .=> 0.0 ], (0 , 1000.0 ), dt = 1e-3 )
231
230
232
231
sf1sys = let
233
232
k = exp (20.7 - 1500 / y[1 ])
@@ -241,7 +240,7 @@ sf1sys = let
241
240
ODESystem (sf1eqs, t, name = :sf1 )
242
241
end
243
242
244
- sf1prob = ODEProblem {false} (sf1sys, [y[1 ] => 761.0 ; y[2 ] => 0.0 ; y[3 ] => 600.0 ; y[4 ] => 0.1 ], (0 , 1000.0 ), dt = 1e-4 , cse = true )
243
+ sf1prob = ODEProblem {false} (sf1sys, [y[1 ] => 761.0 ; y[2 ] => 0.0 ; y[3 ] => 600.0 ; y[4 ] => 0.1 ], (0 , 1000.0 ), dt = 1e-4 )
245
244
246
245
sf2sys = let
247
246
sf2eqs = [
@@ -252,7 +251,7 @@ sf2sys = let
252
251
ODESystem (sf2eqs, t, name = :sf2 )
253
252
end
254
253
255
- sf2prob = ODEProblem {false} (sf2sys, [y[1 ] => 1.0 ; y[2 ] => 0.0 ], (0 , 240.0 ), dt = 1e-2 , cse = true )
254
+ sf2prob = ODEProblem {false} (sf2sys, [y[1 ] => 1.0 ; y[2 ] => 0.0 ], (0 , 240.0 ), dt = 1e-2 )
256
255
257
256
sf3sys = let
258
257
sf3eqs = [
@@ -266,7 +265,7 @@ sf3sys = let
266
265
ODESystem (sf3eqs, t, name = :sf3 )
267
266
end
268
267
269
- sf3prob = ODEProblem {false} (sf3sys, [y[1 ] => 4e-6 ; y[2 ] => 1e-6 ; y[3 : 5 ] .=> 0.0 ], (0 , 100.0 ), dt = 1e-6 , cse = true )
268
+ sf3prob = ODEProblem {false} (sf3sys, [y[1 ] => 4e-6 ; y[2 ] => 1e-6 ; y[3 : 5 ] .=> 0.0 ], (0 , 100.0 ), dt = 1e-6 )
270
269
271
270
sf4sys = let
272
271
sf4eqs = [
@@ -278,7 +277,7 @@ sf4sys = let
278
277
ODESystem (sf4eqs, t, name = :sf4 )
279
278
end
280
279
281
- sf4prob = ODEProblem {false} (sf4sys, [y[1 ] => 4.0 ; y[2 ] => 1.1 ; y[3 ] => 4.0 ], (0 , 300.0 ), dt = 1e-3 , cse = true )
280
+ sf4prob = ODEProblem {false} (sf4sys, [y[1 ] => 4.0 ; y[2 ] => 1.1 ; y[3 ] => 4.0 ], (0 , 300.0 ), dt = 1e-3 )
282
281
283
282
sf5sys = let
284
283
sf5eqs = [
@@ -291,7 +290,7 @@ sf5sys = let
291
290
ODESystem (sf5eqs, t, name = :sf5 )
292
291
end
293
292
294
- sf5prob = ODEProblem {false} (sf5sys, [y[1 ] => 3.365e-7 ; y[2 ] => 8.261e-3 ; y[3 ] => 1.642e-3 ; y[4 ] => 9.38e-6 ], (0 , 100.0 ), dt = 1e-7 , cse = true )
293
+ sf5prob = ODEProblem {false} (sf5sys, [y[1 ] => 3.365e-7 ; y[2 ] => 8.261e-3 ; y[3 ] => 1.642e-3 ; y[4 ] => 9.38e-6 ], (0 , 100.0 ), dt = 1e-7 )
295
294
296
295
# Non-stiff
297
296
na1sys = let y = y[1 ]
@@ -300,39 +299,39 @@ na1sys = let y = y[1]
300
299
ODESystem (na1eqs, t, name = :na1 )
301
300
end
302
301
303
- na1prob = ODEProblem {false} (na1sys, [y[1 ] => 1 ], (0 , 20.0 ), cse = true )
302
+ na1prob = ODEProblem {false} (na1sys, [y[1 ] => 1 ], (0 , 20.0 ))
304
303
305
304
na2sys = let y = y[1 ]
306
305
na2eqs = D (y) ~ - y^ 3 / 2
307
306
308
307
ODESystem (na2eqs, t, name = :na2 )
309
308
end
310
309
311
- na2prob = ODEProblem {false} (na2sys, [y[1 ] => 1 ], (0 , 20.0 ), cse = true )
310
+ na2prob = ODEProblem {false} (na2sys, [y[1 ] => 1 ], (0 , 20.0 ))
312
311
313
312
na3sys = let y = y[1 ]
314
313
na3eqs = D (y) ~ y * cos (t)
315
314
316
315
ODESystem (na3eqs, t, name = :na3 )
317
316
end
318
317
319
- na3prob = ODEProblem {false} (na3sys, [y[1 ] => 1 ], (0 , 20.0 ), cse = true )
318
+ na3prob = ODEProblem {false} (na3sys, [y[1 ] => 1 ], (0 , 20.0 ))
320
319
321
320
na4sys = let y = y[1 ]
322
321
na4eqs = D (y) ~ y/ 4 * (1 - y/ 20 )
323
322
324
323
ODESystem (na4eqs, t, name = :na4 )
325
324
end
326
325
327
- na4prob = ODEProblem {false} (na4sys, [y[1 ] => 1 ], (0 , 20.0 ), cse = true )
326
+ na4prob = ODEProblem {false} (na4sys, [y[1 ] => 1 ], (0 , 20.0 ))
328
327
329
328
na5sys = let y = y[1 ]
330
329
na5eqs = D (y) ~ (y - t) / (y + t)
331
330
332
331
ODESystem (na5eqs, t, name = :na5 )
333
332
end
334
333
335
- na5prob = ODEProblem {false} (na5sys, [y[1 ] => 4 ], (0 , 20.0 ), cse = true )
334
+ na5prob = ODEProblem {false} (na5sys, [y[1 ] => 4 ], (0 , 20.0 ))
336
335
337
336
nb1sys = let
338
337
nb1eqs = [
@@ -343,7 +342,7 @@ nb1sys = let
343
342
ODESystem (nb1eqs, t, name = :nb1 )
344
343
end
345
344
346
- nb1prob = ODEProblem {false} (nb1sys, [y[1 ] => 1.0 , y[2 ] => 3 ], (0 , 20.0 ), cse = true )
345
+ nb1prob = ODEProblem {false} (nb1sys, [y[1 ] => 1.0 , y[2 ] => 3 ], (0 , 20.0 ))
347
346
348
347
nb2sys = let
349
348
nb2eqs = [
@@ -355,7 +354,7 @@ nb2sys = let
355
354
ODESystem (nb2eqs, t, name = :nb2 )
356
355
end
357
356
358
- nb2prob = ODEProblem {false} (nb2sys, [y[1 ] => 2.0 , y[2 ] => 0.0 , y[3 ] => 1.0 ], (0 , 20.0 ), cse = true )
357
+ nb2prob = ODEProblem {false} (nb2sys, [y[1 ] => 2.0 , y[2 ] => 0.0 , y[3 ] => 1.0 ], (0 , 20.0 ))
359
358
360
359
nb3sys = let
361
360
nb3eqs = [
@@ -367,7 +366,7 @@ nb3sys = let
367
366
ODESystem (nb3eqs, t, name = :nb3 )
368
367
end
369
368
370
- nb3prob = ODEProblem {false} (nb3sys, [y[1 ] => 1.0 ; y[2 : 3 ] .=> 0.0 ], (0 , 20.0 ), cse = true )
369
+ nb3prob = ODEProblem {false} (nb3sys, [y[1 ] => 1.0 ; y[2 : 3 ] .=> 0.0 ], (0 , 20.0 ))
371
370
372
371
nb4sys = let
373
372
r = sqrt (y[1 ]^ 2 + y[2 ]^ 2 )
@@ -380,7 +379,7 @@ nb4sys = let
380
379
ODESystem (nb4eqs, t, name = :nb4 )
381
380
end
382
381
383
- nb4prob = ODEProblem {false} (nb4sys, [y[1 ] => 3.0 ; y[2 : 3 ] .=> 0.0 ], (0 , 20.0 ), cse = true )
382
+ nb4prob = ODEProblem {false} (nb4sys, [y[1 ] => 3.0 ; y[2 : 3 ] .=> 0.0 ], (0 , 20.0 ))
384
383
385
384
nb5sys = let
386
385
nb5eqs = [
@@ -392,7 +391,7 @@ nb5sys = let
392
391
ODESystem (nb5eqs, t, name = :nb5 )
393
392
end
394
393
395
- nb5prob = ODEProblem {false} (nb5sys, [y[1 ] => 0.0 ; y[2 : 3 ] .=> 1.0 ], (0 , 20.0 ), cse = true )
394
+ nb5prob = ODEProblem {false} (nb5sys, [y[1 ] => 0.0 ; y[2 : 3 ] .=> 1.0 ], (0 , 20.0 ))
396
395
397
396
nc1sys = let y = y
398
397
n = 10
@@ -412,7 +411,7 @@ nc2sys = let y = y
412
411
ODESystem (nc2eqs, t, name = :nc2 )
413
412
end
414
413
415
- nc2prob = ODEProblem {false} (nc2sys, [y[1 ] => 1.0 ; y[2 : 10 ] .=> 0.0 ], (0 , 20.0 ), cse = true )
414
+ nc2prob = ODEProblem {false} (nc2sys, [y[1 ] => 1.0 ; y[2 : 10 ] .=> 0.0 ], (0 , 20.0 ))
416
415
417
416
nc3sys = let y = y
418
417
n = 10
@@ -422,7 +421,7 @@ nc3sys = let y = y
422
421
ODESystem (nc3eqs, t, name = :nc3 )
423
422
end
424
423
425
- nc3prob = ODEProblem {false} (nc3sys, [y[1 ] => 1.0 ; y[2 : 10 ] .=> 0.0 ], (0 , 20.0 ), cse = true )
424
+ nc3prob = ODEProblem {false} (nc3sys, [y[1 ] => 1.0 ; y[2 : 10 ] .=> 0.0 ], (0 , 20.0 ))
426
425
427
426
@variables y (t)[1 : 51 ]
428
427
y = collect (y)
@@ -434,7 +433,7 @@ nc4sys = let y = y
434
433
ODESystem (nc4eqs, t, name = :nc4 )
435
434
end
436
435
437
- nc4prob = ODEProblem {false} (nc4sys, [y[1 ] => 1.0 ; y[2 : 51 ] .=> 0.0 ], (0 , 20.0 ), cse = true )
436
+ nc4prob = ODEProblem {false} (nc4sys, [y[1 ] => 1.0 ; y[2 : 51 ] .=> 0.0 ], (0 , 20.0 ))
438
437
439
438
@variables y (t)[1 : 3 , 1 : 5 ]
440
439
y = collect (y)
@@ -470,7 +469,7 @@ ys′ = [-0.557160570446, 0.505696783289, 0.230578543901,
470
469
y0 = y .=> reshape (ys, 3 , 5 )
471
470
y0′ = D .(y) .=> reshape (ys′, 3 , 5 )
472
471
# The orginal paper has t_f = 20, but 1000 looks way better
473
- nc5prob = ODEProblem {false} (nc5sys, [y0; y0′], (0 , 20.0 ), cse = true )
472
+ nc5prob = ODEProblem {false} (nc5sys, [y0; y0′], (0 , 20.0 ))
474
473
475
474
@variables y (t)[1 : 4 ]
476
475
y = collect (y)
@@ -514,7 +513,7 @@ ne2sys = let
514
513
end
515
514
516
515
y0 = [y[1 ] => 2.0 ; y[2 ] => 0.0 ]
517
- ne2prob = ODEProblem {false} (ne2sys, y0, (0 , 20.0 ), cse = true )
516
+ ne2prob = ODEProblem {false} (ne2sys, y0, (0 , 20.0 ))
518
517
519
518
ne3sys = let
520
519
ne3eqs = [D (y[1 ]) ~ y[2 ],
0 commit comments