301
301
# basic VariableRateJump test
302
302
let
303
303
N = 1000 # number of simulations for testing solve accuracy
304
- Random. seed! (rng, 1111 )
304
+ Random. seed! (rng, 1111 )
305
305
@variables A (t) B (t) C (t)
306
306
@parameters k
307
307
vrj = VariableRateJump (k * (sin (t) + 1 ), [A ~ A + 1 , C ~ C + 2 ])
@@ -425,17 +425,19 @@ end
425
425
426
426
# PDMP test
427
427
let
428
- Random. seed! (rng, 1111 )
428
+ Random. seed! (rng, 1111 )
429
429
@variables X (t) Y (t)
430
430
@parameters k1 k2
431
431
vrj1 = VariableRateJump (k1 * X, [X ~ X - 1 ]; save_positions = (false , false ))
432
432
vrj2 = VariableRateJump (k1, [Y ~ Y + 1 ]; save_positions = (false , false ))
433
- eqs = [D (X) ~ k2, D (Y) ~ - k2/ 10 * Y]
433
+ eqs = [D (X) ~ k2, D (Y) ~ - k2 / 10 * Y]
434
434
@named jsys = JumpSystem ([vrj1, vrj2, eqs[1 ], eqs[2 ]], t, [X, Y], [k1, k2])
435
435
jsys = complete (jsys)
436
- X0 = 0.0 ; Y0 = 0.0
436
+ X0 = 0.0
437
+ Y0 = 0.0
437
438
u0 = [X => X0, Y => Y0]
438
- k1val = 1.0 ; k2val = 20.0
439
+ k1val = 1.0
440
+ k2val = 20.0
439
441
p = [k1 => k1val, k2 => k2val]
440
442
tspan = (0.0 , 10.0 )
441
443
oprob = ODEProblem (jsys, u0, tspan, p)
@@ -447,27 +449,30 @@ let
447
449
Yv = zeros (length (times))
448
450
for n in 1 : Nsims
449
451
sol = solve (jprob, Tsit5 (); saveat = times)
450
- Xv .+ = sol[1 ,:]
451
- Yv .+ = sol[2 ,:]
452
+ Xv .+ = sol[1 , :]
453
+ Yv .+ = sol[2 , :]
452
454
end
453
- Xv ./= Nsims; Yv ./= Nsims;
455
+ Xv ./= Nsims
456
+ Yv ./= Nsims
454
457
455
458
Xact (t) = X0 * exp (- k1val * t) + (k2val / k1val) * (1 - exp (- k1val * t))
456
- Yact (t) = Y0 * exp (- k2val/ 10 * t) + (k1val / (k2val/ 10 )) * (1 - exp (- k2val/ 10 * t))
459
+ function Yact (t)
460
+ Y0 * exp (- k2val / 10 * t) + (k1val / (k2val / 10 )) * (1 - exp (- k2val / 10 * t))
461
+ end
457
462
@test all (abs .(Xv .- Xact .(times)) .<= 0.05 .* Xv)
458
- @test all (abs .(Yv .- Yact .(times)) .<= 0.1 .* Yv)
463
+ @test all (abs .(Yv .- Yact .(times)) .<= 0.1 .* Yv)
459
464
end
460
465
461
466
# that mixes ODEs and jump types, and then contin events
462
467
let
463
- Random. seed! (rng, 1111 )
468
+ Random. seed! (rng, 1111 )
464
469
@variables X (t) Y (t)
465
470
@parameters α β
466
471
vrj = VariableRateJump (β * X, [X ~ X - 1 ]; save_positions = (false , false ))
467
- crj = ConstantRateJump (β* Y, [Y ~ Y - 1 ])
472
+ crj = ConstantRateJump (β * Y, [Y ~ Y - 1 ])
468
473
maj = MassActionJump (α, [0 => 1 ], [Y => 1 ])
469
- eqs = [D (X) ~ α* (1 + Y)]
470
- @named jsys = JumpSystem ([maj, crj, vrj, eqs[1 ]], t, [X, Y], [α, β])
474
+ eqs = [D (X) ~ α * (1 + Y)]
475
+ @named jsys = JumpSystem ([maj, crj, vrj, eqs[1 ]], t, [X, Y], [α, β])
471
476
jsys = complete (jsys)
472
477
p = (α = 6.0 , β = 2.0 , X₀ = 2.0 , Y₀ = 1.0 )
473
478
u0map = [X => p. X₀, Y => p. Y₀]
@@ -481,32 +486,34 @@ let
481
486
Yv = zeros (length (times))
482
487
for n in 1 : Nsims
483
488
sol = solve (jprob, Tsit5 (); saveat = times)
484
- Xv .+ = sol[1 ,:]
485
- Yv .+ = sol[2 ,:]
489
+ Xv .+ = sol[1 , :]
490
+ Yv .+ = sol[2 , :]
486
491
end
487
- Xv ./= Nsims; Yv ./= Nsims;
492
+ Xv ./= Nsims
493
+ Yv ./= Nsims
488
494
489
495
function Yf (t, p)
490
496
local α, β, X₀, Y₀ = p
491
497
return (α / β) + (Y₀ - α / β) * exp (- β * t)
492
498
end
493
499
function Xf (t, p)
494
500
local α, β, X₀, Y₀ = p
495
- return (α / β) + (α^ 2 / β^ 2 ) + α * (Y₀ - α / β) * t * exp (- β * t) + (X₀ - α / β - α^ 2 / β^ 2 ) * exp (- β * t)
501
+ return (α / β) + (α^ 2 / β^ 2 ) + α * (Y₀ - α / β) * t * exp (- β * t) +
502
+ (X₀ - α / β - α^ 2 / β^ 2 ) * exp (- β * t)
496
503
end
497
- Xact = [Xf (t,p) for t in times]
498
- Yact = [Yf (t,p) for t in times]
504
+ Xact = [Xf (t, p) for t in times]
505
+ Yact = [Yf (t, p) for t in times]
499
506
@test all (abs .(Xv .- Xact) .<= 0.05 .* Xv)
500
- @test all (abs .(Yv .- Yact) .<= 0.05 .* Yv)
501
-
507
+ @test all (abs .(Yv .- Yact) .<= 0.05 .* Yv)
508
+
502
509
function affect! (integ, u, p, ctx)
503
510
savevalues! (integ, true )
504
511
terminate! (integ)
505
512
nothing
506
513
end
507
- cevents = [t ~ .2 ] => (affect!, [], [], [], nothing )
508
- @named jsys = JumpSystem ([maj, crj, vrj, eqs[1 ]], t, [X, Y], [α, β];
509
- continuous_events = cevents)
514
+ cevents = [t ~ 0 .2 ] => (affect!, [], [], [], nothing )
515
+ @named jsys = JumpSystem ([maj, crj, vrj, eqs[1 ]], t, [X, Y], [α, β];
516
+ continuous_events = cevents)
510
517
jsys = complete (jsys)
511
518
tspan = (0.0 , 200.0 )
512
519
oprob = ODEProblem (jsys, u0map, tspan, pmap)
516
523
for n in 1 : Nsims
517
524
sol = solve (jprob, Tsit5 (), saveat = tspan[2 ])
518
525
@test sol. retcode == ReturnCode. Terminated
519
- Xsamp += sol[1 ,end ]
526
+ Xsamp += sol[1 , end ]
520
527
end
521
528
Xsamp /= Nsims
522
- @test abs (Xsamp - Xf (.2 ,p) < .05 * Xf (.2 ,p))
523
- end
529
+ @test abs (Xsamp - Xf (0 .2 , p) < 0 .05 * Xf (0 .2 , p))
530
+ end
0 commit comments