# 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]


## 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)