SciML Workshop Exercise Solutions

Chris Rackauckas
using DifferentialEquations
using Sundials
using BenchmarkTools
using Plots

Problem 1: Investigating Sources of Randomness and Uncertainty in a Biological System

Part 1: Simulating the Oregonator ODE model

using DifferentialEquations, Plots
function orego(du,u,p,t)
  s,q,w = p
  y1,y2,y3 = u
  du[1] = s*(y2+y1*(1-q*y1-y2))
  du[2] = (y3-(1+y1)*y2)/s
  du[3] = w*(y1-y3)
end
p = [77.27,8.375e-6,0.161]
prob = ODEProblem(orego,[1.0,2.0,3.0],(0.0,360.0),p)
sol = solve(prob)
plot(sol)
plot(sol,vars=(1,2,3))

Part 2: Investigating Stiffness

using BenchmarkTools
prob = ODEProblem(orego,[1.0,2.0,3.0],(0.0,50.0),p)
@btime sol = solve(prob,Tsit5())
945.686 ms (8723143 allocations: 920.67 MiB)
retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 872306-element Array{Float64,1}:
  0.0
  0.01618926718934831
  0.02355386004837834
  0.03818038870154586
  0.050503515877727514
  0.06810672932191658
  0.08676359998206734
  0.11145368602241688
  0.14105967462147356
  0.18104879156165962
  ⋮
 49.99977330536325
 49.99980456142745
 49.999835817515255
 49.999867073624586
 49.999898329755446
 49.99992958590576
 49.99996084207554
 49.999992098266844
 50.0
u: 872306-element Array{Array{Float64,1},1}:
 [1.0, 2.0, 3.0]
 [1.7128564042197614, 1.9996098373795999, 2.9959141611121862]
 [1.8376268914687968, 1.9993653073090198, 2.994474646468457]
 [1.9480445809808178, 1.9988333244430836, 2.991907642632475]
 [1.9807789479174538, 1.998364632682339, 2.989876120098015]
 [1.996520358969301, 1.9976843022063284, 2.9870473687154533]
 [2.0012471416469095, 1.9969587120867922, 2.9840850652644586]
 [2.003267094253373, 1.9959962346456372, 2.980190667568818]
 [2.0046071951018165, 1.9948405279663373, 2.9755485736940304]
 [2.0062040975915965, 1.9932773146432707, 2.969322732597494]
 ⋮
 [1.00114451241949, 1453.0173573419604, 414.83224206133156]
 [1.0011445128905938, 1453.0163492345089, 414.8301595725294]
 [1.001144513536549, 1453.0153411262695, 414.82807709263454]
 [1.001144514166616, 1453.014333017309, 414.82599462178484]
 [1.0011445147807905, 1453.013324907627, 414.8239121599803]
 [1.0011445151883325, 1453.0123167972909, 414.82182970735875]
 [1.0011445153892404, 1453.0113086863003, 414.8197472639202]
 [1.001144515574252, 1453.0103005745884, 414.8176648295267]
 [1.0008765717435082, 1453.0100456736175, 414.8171383809634]
@btime sol = solve(prob,Rodas5())
479.696 μs (1909 allocations: 130.16 KiB)
retcode: Success
Interpolation: 3rd order Hermite
t: 110-element Array{Float64,1}:
  0.0
  0.019615259849088615
  0.029598314714131158
  0.04705295553350644
  0.06489958093933189
  0.08933251171067431
  0.12069400166576917
  0.16655311655246774
  0.24089140897016648
  0.39558909491704786
  ⋮
 26.756905610888992
 27.982111658219903
 29.768997154114096
 32.21837697976615
 35.093850201346655
 38.49798110093118
 42.33811919585127
 46.60842194880463
 50.0
u: 110-element Array{Array{Float64,1},1}:
 [1.0, 2.0, 3.0]
 [1.7804115041903392, 1.9994992840408727, 2.995224421497252]
 [1.898773632635922, 1.9991507098568697, 2.9933805881501456]
 [1.9745775749460168, 1.9984968888022705, 2.9904382700551317]
 [1.9949959346655894, 1.9978087397951183, 2.9875591897227847]
 [2.0015958931121642, 1.9968586608477479, 2.98367866834122]
 [2.003748190575679, 1.9956356930368464, 2.9787387129953866]
 [2.0056429388535917, 1.9938442509772465, 2.9715736894920433]
 [2.0085949421229565, 1.9909335157971781, 2.960099467726684]
 [2.014815188384092, 1.9848502001186519, 2.936770263171178]
 ⋮
 [1.0009510454262696, 1052.1681949981978, 17454.97704553619]
 [1.000790082105047, 1266.4223517298105, 14330.342720311946]
 [1.0006713873660182, 1490.2781714142227, 10747.93771088393]
 [1.000598803847115, 1670.9447027478102, 7245.705166049239]
 [1.000568993307521, 1758.4723173221284, 4560.988616721742]
 [1.000569273504183, 1757.6100577789323, 2636.982996349979]
 [1.000594225030407, 1683.8471494545056, 1421.4818618119598]
 [1.0006409946157637, 1561.0560213127278, 715.2527024515273]
 [1.0006887475677544, 1452.8969192375328, 414.7220773988324]

(Optional) Part 3: Specifying Analytical Jacobians (I)

(Optional) Part 4: Automatic Symbolicification and Analytical Jacobian Calculations

Part 5: Adding stochasticity with stochastic differential equations

function orego(du,u,p,t)
  s,q,w = p
  y1,y2,y3 = u
  du[1] = s*(y2+y1*(1-q*y1-y2))
  du[2] = (y3-(1+y1)*y2)/s
  du[3] = w*(y1-y3)
end
function g(du,u,p,t)
  du[1] = 0.1u[1]
  du[2] = 0.1u[2]
  du[3] = 0.1u[3]
end
p = [77.27,8.375e-6,0.161]
prob = SDEProblem(orego,g,[1.0,2.0,3.0],(0.0,30.0),p)
sol = solve(prob,SOSRI())
plot(sol)