Dynamical Systems

In the following notebook we are going to explore how to define and visualize dynamical systems in Julia. The ecosystem is wide and very flexible, so we'll first start with a simple example, and then consider something a tad more complex.

using

we load some package to make our life simpler.
we take a look at what Pluto.jl does for us, and why that makes reproducibility easy peasy.

using OrdinaryDiffEq, Plots

The story of 🐱, 😄 and their troubled love

We are going to talk about two characters: a cat, and a smiley-person.

Being the data people that we are, we'll measure their love (expressed as a real variable u) across time (t). When u is positive, there's love. When u is negative, well... not so much love.

# At the begin of the story, smiley-person really loves his cat.
# The cat, however, is not impressed.
🐱₀, 😄₀ = [-1.0, 1.0];

@show 🐱₀, 😄₀

(-1.0, 1.0)

a, b = [1.0, .5];

function cat_love(du,u,p,t)
    🐱, 😄 = u # the value of love
    a, b = p # the parameters of the love story
  # the cat being a cat, runs away from overly affectionate people
    du[1] = d🐱 = - a*😄
  # the smiley person love the cat more when the cat is lovely
    du[2] = d😄 = b*🐱
end;

begin
    p = [a, b]
    u₀ = [🐱₀,😄₀]
end

2-element Vector{Float64}:
 -1.0
  1.0

prob = ODEProblem(cat_love, # the equation system
    u₀, # initial state
    (0.0,30.0), # time interval
    p # parameters
)

ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 30.0)
u0: 2-element Vector{Float64}:
 -1.0
  1.0

solution = solve(prob,Tsit5(), saveat = 0.1)

	timestamp	value1	value2
1	0.0	-1.0	1.0
2	0.1	-1.09742	0.947543
3	0.2	-1.18935	0.89035
4	0.3	-1.27534	0.828707
5	0.4	-1.35495	0.762922
6	0.5	-1.4278	0.693324
7	0.6	-1.4935	0.620261
8	0.7	-1.55175	0.544099
9	0.8	-1.60223	0.465216
10	0.9	-1.64471	0.384008
...

trajectories = Array(solution)'

301×2 adjoint(::Matrix{Float64}) with eltype Float64:
 -1.0        1.0
 -1.09742    0.947543
 -1.18935    0.89035
 -1.27534    0.828707
 -1.35495    0.762922
 -1.4278     0.693324
 -1.4935     0.620261
  ⋮         
 -0.85362   -1.06643
 -0.744942  -1.10641
 -0.632529  -1.14085
 -0.516944  -1.16959
 -0.398768  -1.19248
 -0.27861   -1.20942

plot( solution.t, trajectories, label = ["🐱" "😄"])

plot(solution,idxs=(1,2))

And now for something different

we are going to consider some sources of noise, and casuality: stochastic differential equations

using DifferentialEquations

function sde_cat_love(du,u,p,t)
    🐱, 😄 = u
    a, b = p
    du[1] = d🐱 = - a*😄 - 1.0
    du[2] = d😄 = b*🐱
end

sde_cat_love (generic function with 1 method)

function σ_cat_love(du,u,p,t)
  du[1] = 0.4
  du[2] = 0.4
end

σ_cat_love (generic function with 1 method)

prob_sde_cat_lover = SDEProblem(sde_cat_love, # the equation system
    σ_cat_love, # the (diagonal) noise
    u₀, # initial state
    (0.0,30.0), # time interval
    p # parameters
)

SDEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 30.0)
u0: 2-element Vector{Float64}:
 -1.0
  1.0

sde_sol = solve(prob_sde_cat_lover)

	timestamp	value1	value2
1	0.0	-1.0	1.0
2	0.000570296	-0.996901	1.01415
3	0.000684356	-1.00193	1.01398
4	0.000812672	-1.00696	1.01311
5	0.000957029	-0.994853	1.01654
6	0.00111943	-0.99314	1.01157
7	0.00130213	-0.993825	1.01095
8	0.00150767	-0.994438	1.01363
9	0.0017389	-1.00304	1.00714
10	0.00199903	-1.00758	1.00445
...

plot(sde_sol.t, Array(sde_sol)')

plot(sde_sol,idxs=(1,2))

ensembleprob = EnsembleProblem(prob_sde_cat_lover)

EnsembleProblem with problem SDEProblem

sols = solve(ensembleprob,EnsembleThreads(),trajectories=100)

EnsembleSolution Solution of length 100 with uType:
RODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, NoiseProcess{Float64, 2, Float64, Vector{Float64}, Vector{Float64}, Vector{Vector{Float64}}, typeof(DiffEqNoiseProcess.INPLACE_WHITE_NOISE_DIST), typeof(DiffEqNoiseProcess.INPLACE_WHITE_NOISE_BRIDGE), true, ResettableStacks.ResettableStack{Tuple{Float64, Vector{Float64}, Vector{Float64}}, true}, ResettableStacks.ResettableStack{Tuple{Float64, Vector{Float64}, Vector{Float64}}, true}, RSWM{Float64}, Nothing, RandomNumbers.Xorshifts.Xoroshiro128Plus}, SDEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, Nothing, SDEFunction{true, SciMLBase.FullSpecialize, typeof(sde_cat_love), typeof(σ_cat_love), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, typeof(σ_cat_love), Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, Nothing}, SOSRI, StochasticDiffEq.LinearInterpolationData{Vector{Vector{Float64}}, Vector{Float64}}, DiffEqBase.DEStats, Nothing}

begin
    using DifferentialEquations.EnsembleAnalysis
    summ = EnsembleSummary(sols,0:0.1:30.)
    plot(summ,labels="Middle 95%")
    summ = EnsembleSummary(sols,0:0.1:30.;quantiles=[0.25,0.75])
    plot!(summ,labels="Middle 50%",legend=true)
end

Built with Julia 1.8.3 and

DifferentialEquations 7.6.0
OrdinaryDiffEq 6.33.3
Plots 1.36.6

To run this tutorial locally, download [this file](/tutorials/01DiffPluto.jl) and open it with Pluto.jl.

Julia 2022 AusDm

Dynamical Systems

using

The story of 🐱, 😄 and their troubled love

And now for something different