| SOCR ≫ | DSPA ≫ | DSPA2 Topics ≫ |
DSPA2: Data Science and Predictive Analytics (UMich HS650)
Appendix 17: Dynamical System Modeling using R-Julia Integration
SOCR (Ivo Dinov and SOCR)
January 28, 2025
This DSPA Appendix
shows examples of dynamic systems, chaotic, spatiotemporal changes, and
biocomplexity modeling using R and Julia integration.
In this appendix, we use the R package
JuliaCall as a Julia language wrapper, and various
Julia packages, such as JuliaDynamics
and DynamicalSystems.
Julia is a high-performance dynamic programming language
for numerical computing that can be used for efficient solution to
optimization problems and other demanding scientific computations in
R. JuliaCall
provides direct integration between R and
Julia, where Julia functions can be called
just like any other R function. Similarly, R
functions can be called in the Julia environment (chink
blocks in the Rmd script). Both invocations provide automatic type
conversion. JuliaConnectoR
is another R package that provides a functionally oriented
interface for integrating Julia and R
functionality. Imported Julia functions can be called as
R functions and data structures are converted
automatically. Also, the package XRJulia
provides an interface from R to computations in
Julia, based on the interface structure described in the book “Extending R” by
John M. Chambers.
1 Basic
Julia code chunks embedded in Rmd
First load the necessary R-Julia connector
libraries.
# R-block of code
# The first time - install these libs
#
# install.packages("JuliaCall")
# install.packages("JuliaConnectoR")
# library(JuliaCall)
# install_julia() # This is needed just once to install Julia on the system!!!
# install.packages("XRJulia") # Install the package if not already installed
# Set the JULIA_NUM_THREADS environment variable in R:
library(XRJulia)
library(JuliaCall)
Sys.setenv(JULIA_NUM_THREADS = 10)
julia <- julia_setup() # This will start Julia with 10 threads
# findJulia() # This will locate the Julia installation
# julia <- JuliaSession() # Start Julia and assign the session to the `julia` object
# Check that the julia object is properly initialized by running a simple Julia command:
julia_command("1 + 1") # # Should return 2## 2message(paste0("The (Julia) result of this command, \'julia_command(\"1 + 1\")\', is ", julia_command("1 + 1"), "!\n"))
# Verify the number of threads in Julia:
julia_command("using Base.Threads; Threads.nthreads() = 10") # Manually set threads# Should return 10
message(paste0("The Number of Cores Julia starts with is: ", julia_eval("Threads.nthreads()"), "!\n"))
# julia_command("using Base.Threads; Threads.nthreads()") # Should return 10
# julia_command("ENV[\"JULIA_NUM_THREADS\"]") # Should return "10"
# JuliaConnectoR::stopJulia() # Terminate the Julia session Next, demonstrate Julia code chunk for a simple
calculation, \(a=\sqrt{10}\). Note that
this just computes the result, which can then be printed either in
Julia or in R code chunk.
## This is a julia language chunk.
## In julia, the command without ending semicolon will trigger the display
## so is JuliaCall package.
## The julia display will follow immediately after the corresponding command
## just as the R code in R Markdown.
a = sqrt(10);
print("(julia) sqrt(10)=",a)## (julia) sqrt(10)=3.1622776601683795Next is a subsequent R code block reporting the result
from the Julia chunk calculation.
2 Installing
Julia Packages from within the R
Environment
Suppose we need the support of the ’Julia` Plots.jl plotting package. Before the first use of Plots.jl, we need to install the package. Then, we can generate some random data and test the plotting functionality.
# Julia code-chunk
# First time - install Julia package from R
# import Pkg; Pkg.add("Plots")
# Pkg.add("CairoMakie")
# Pkg.add("Plots")
# Pkg.add("PyPlot")
import Plots
# call the plot function
using Plots
# using PyPlot
gr()## Plots.GRBackend()
## Plots.GRBackend()
# generate 7 random series of 100 values and plot then with line width=2
plot(Plots.fakedata(100, 7), w=2)
# # Create some data to test the "scatter()" function
# x = [1, 2, 3, 4, 5]
# y = [2, 3, 5, 7, 11]
# # Create a scatter plot
# scatter(x, y, label="Data Points", xlabel="X", ylabel="Y", title="X-Y Scatter Plot Test")Following the example in this
Jupyter Notebook, we can demonstrate launching a julia
console in an R chunk of code.
# library(JuliaCall)
# julia_setup()
# julia_console() # starts a simple julia console. Test with specific Julia commands.
# Type exit and then enter to exit.
# julia> sqrt(10)
# julia> exit
julia_library("Plots") # no semicolumns ";" to show the actual results (text/figures/etc.) in output# julia_command("gr()")
# Plots.GRBackend() call
julia_command("Plots.plot(Plots.fakedata(200, 5), w=3)")
Next, we will explore some more elaborate Julia
utilities for modeling and visualization of dynamical systems.
3 Interactive Dynamics
The Julia InteractiveDynamics
package extends various JuliaDynamics
packages and provides a suite of tools for interactive exploration
of dynamical systems.
The Julia DynamicalSystems tutorial provides instructions for constructing dynamical systems using triples of:
- The dynamic rule \(f\) is
a
Juliafunction with instructions of how to evolve the dynamical system over time, - The state \(u\) is an array, or container, including all the variables of the dynamical system and defining the starting state of the system, and
- The parameters \(p\) is an arbitrary container that parameterizes the dynamical rule, \(f\).
In principle, specific implementations of a DynamicalSystem define \(f\) and \(u\). The dynamical rule \(f\) may be defined either as an in-place (iip) function or as an out-of-place (oop) function:
- oop definition of \(f\)
must be in the form \(f(u, p, t) \to
out\), which means that given a state \(u::SVector{<:Real}\) and some parameter
container \(p\), the dynamic rule
returns the output of \(f\) as an
SVector{<:Real}static vector, - iip definition of \(f\)
must be in the form \(f!(out, u, p,
t)\), which means that given a state \(u::AbstractArray{<:Real}\) and some
parameter container \(p\), the function
writes in-place the output of \(f\) in
out::AbstractArray{<:Real}, i.e., the function return nothing as a final statement (no output as it overwrites in-place \(f\).
In both cases, \(t\) represents the current time. The iip formulation is useful for high dimensional systems. Whereas the oop formulation is useful for lower dimensional systems. The cutoff point is typically in the range \([10, 100]\) dimensions, but case-by-case tests with specific \(f\) complexity is useful, in general.
In autonomous dynamical systems, \(f\) doesn’t depend on time, whereas in non-autonomous dynamical systems, \(f\) may depend on time. In either case, the dynamic rule \(f\) always includes \(t\) as an \(f\) argument, as some algorithms utilize this information while pthers do not.
3.1 Dynamical System Example: Hénon map
The Hénon map is a discrete-time dynamical system that exhibits a chaotic behavior. The Hénon map takes a point \((x_n, y_n)\) in the plane and maps it to the next point, \((x_{n+1}, y_{n+1})\), using the following transformation
\[x_{n+1} = 1 − ax_n^2 + y_n ,\]
\[y_{n+1}=bx_n.\] The classical values
of the two parameters, \(a\) and \(b\), of the Hénon map are \(a=1.4\), \(b=0.3\). The dynamic rule is
defined as a standard Julia function. This is a 2D
dynamical system and we can use the out-of-place form that
returns an SVector with the next state.
# Install prerequisite Julia packages
# using Pkg; Pkg.add("DifferentialEquations") # Install the package
# Pkg.add("DynamicalSystemsBase")
# Pkg.add("DiffEqNoiseProcess")
using DifferentialEquations
using DynamicalSystems
function henon_rule(u, p, n) # note that the time `n` is not used
x, y = u # system state
a, b = p # system parameters
xn = 1.0 - a*x^2 + y
yn = b*x
return SVector(xn, yn)
end## henon_rule (generic function with 1 method)Let’s define the dynamical system initial state, \(u_o=(x=0.2, y=0.3)\) and parameters, \(p_o=(a=1.4, b=0.3)\).
## 2-element Vector{Float64}:
## 0.2
## 0.3## 2-element Vector{Float64}:
## 1.4
## 0.3Pass the three dynamical system components (\(f\), \(u\), and \(p\)) to the Julia
DeterministicIteratedMap() function.
## 2-dimensional DeterministicIteratedMap
## deterministic: true
## discrete time: true
## in-place: false
## dynamic rule: henon_rule
## parameters: [1.4, 0.3]
## time: 0
## state: [0.2, 0.3]The resulting object, \(henon\), is
a DynamicalSystem representing the first core structure
type of the library. It can evolve interactively, and be
queried, using the interface defined by DynamicalSystem.
The simplest way to query the evolution is to track its trajectory.
## 15000## (2-dimensional StateSpaceSet{Float64} with 15001 points, 0:1:15000)## 2-dimensional StateSpaceSet{Float64} with 15001 points
## 0.2 0.3
## 1.244 0.06
## -1.10655 0.3732
## -0.341035 -0.331965
## 0.505208 -0.102311
## 0.540361 0.151562
## 0.742777 0.162108
## 0.389703 0.222833
## 1.01022 0.116911
## -0.311842 0.303065
## ⋮
## 0.342534 0.199135
## 1.03487 0.10276
## -0.396586 0.310462
## 1.09027 -0.118976
## -0.783137 0.327081
## 0.468455 -0.234941
## 0.457829 0.140536
## 0.847086 0.137349
## 0.132771 0.254126plot(X[:, 1], X[:, 2], seriestype=:scatter, label="Data Points", xlabel="X[:, 1]", ylabel="X[:, 2]", title="Heron Trajectory Plot")
# plot(X[:, 1], X[:, 2], title = "System Trajectory", xlabel = "x", ylabel = "y")
# scatter(X, xlabel="Iteration", ylabel="x") The outcome \(X\) is a
StateSpaceSet representing the second core structure
type of the Dynamics library. Later, we will explore
StateSpaceset. Now we’ll scatter plot the resulting
dynamics of the Henon map.
using CairoMakie
using Plots
# using PyPlot
# scatter(X)
# to avoid This error: Error: Error happens in Julia.
# Cannot convert StateSpaceSet{2, Float64} to series data for plotting
# Extract data from the StateSpaceSet first
x = X[:, 1] # First dimension
y = X[:, 2] # Second dimension
# z = X[:, 3] # Third dimension (if applicable)
# scatter(x, y, markersize=2, label="Lorenz Trajectory", xlabel="x", ylabel="y")
plot(X[:, 1], X[:, 2], seriestype=:scatter, label="Data Points", xlabel="X[:, 1]", ylabel="X[:, 2]", title="Heron Trajectory Plot")3.2 Dynamical System Example: Lorenz96 Model
The next dynamical system example uses the Lorenz96 Model.
\[\frac{dx_i}{dt} =(x_{i+1} − x_{i−2}) x_{i−1} − x_i + F,\ \forall i\in \{1,\cdots,N\},\] where \(N\geq 4\), initial conditions include \(x_{-1}=x_{N-1}\), \(x_{0}=x_{N}\), \(x_{N+1} = x_{1}\), \(x_{i}\) is the state of the system, and \(F\) is a forcing constant, e.g., \(F=8\) causes a chaotic system behavior.
In this case, instead of a discrete time map we have \(N\) coupled ordinary differential equations
(ODEs). The initiation of the dynamical system is the same but utilizes
the function CoupledODEs, instead of
DeterministicIteratedMap. To formulate the dynamic rule
function, \(f\), note that the
Lorenz96 dynamical system can be of an arbitrary high
dimension. Hence, we should use the (iip) in-place form of
\(f\) to overwrite in place the
rate of change in a pre-allocated container. Often, we append the
name of functions that modify their arguments in-place with a
bang (\(!\)), in this case, \(lorenz96\_rule() \to
lorenz96\_rule!()\).
function lorenz96_rule!(du, u, p, t)
F = p[1]; N = length(u)
# 3 edge cases: $x_{-1}=x_{N-1}$, $x_{0}=x_{N}$, $x_{N+1} = x_{1}$
du[1] = (u[2] - u[N - 1]) * u[N] - u[1] + F
du[2] = (u[3] - u[N]) * u[1] - u[2] + F
du[N] = (u[1] - u[N - 2]) * u[N - 1] - u[N] + F
# then the general case
for n in 3:(N - 1)
du[n] = (u[n + 1] - u[n - 2]) * u[n - 1] - u[n] + F
end
return nothing # always `return nothing` for in-place form function dynamics
end## lorenz96_rule! (generic function with 1 method)As we did in the first example, prior to initializing the system, we define the initial state, e.g., in 7D, \(u_o = (0.10, 0.25, 0.40, 0.55, 0.70, 0.85, 1.00)\), and \(p_o = 8.0\).
## 7## 0.1:0.15:1.0## 1-element Vector{Float64}:
## 8.0## 7-dimensional CoupledODEs
## deterministic: true
## discrete time: false
## in-place: true
## dynamic rule: lorenz96_rule!
## ODE solver: Tsit5
## ODE kwargs: (abstol = 1.0e-6, reltol = 1.0e-6)
## parameters: [8.0]
## time: 0.0
## state: [0.1, 0.25, 0.4, 0.55, 0.7, 0.85, 1.0]Next, we can explore an instance of the dynamical trajectory, which has longitudinal length of \(626\) points.
# total_time = 12.5
# sampling_time = 0.02
# Y, t = trajectory(lorenz96, total_time; Ttr = 2.2, deltat = sampling_time)
# Y
using DifferentialEquations
using Plots
############ Old working approach
total_time = 12.5## 12.5## 0.02## (7-dimensional StateSpaceSet{Float64} with 626 points, 2.2:0.02:14.7)## 7-dimensional StateSpaceSet{Float64} with 626 points
## -1.73386 -2.68151 0.299161 … 0.752006 7.44039 -2.01891
## -1.12526 -2.53848 0.375484 0.901361 7.41654 -2.15216
## -0.517228 -2.37205 0.475519 1.07297 7.38089 -2.21566
## 0.0750631 -2.17845 0.595985 1.26569 7.33222 -2.21535
## 0.638673 -1.95644 0.732544 1.47819 7.26956 -2.15915
## 1.16345 -1.70695 0.880183 … 1.70893 7.19222 -2.0564
## 1.64231 -1.43261 1.03369 1.95613 7.09973 -1.91731
## 2.0713 -1.13718 1.18814 2.21783 6.99173 -1.75229
## 2.4493 -0.824986 1.33933 2.49184 6.86786 -1.57144
## 2.77769 -0.500433 1.48399 2.77579 6.72765 -1.38401
## ⋮ ⋱ ⋮
## -3.32982 0.0685096 -2.1403 11.2015 1.55554 -0.347469
## -3.09614 0.32556 -1.91601 11.3183 0.964323 -0.542294
## -2.87366 0.550075 -1.67247 11.3734 0.425321 -0.568705
## -2.66304 0.753774 -1.41678 11.3823 -0.0380914 -0.451235
## -2.45905 0.945021 -1.15332 … 11.3575 -0.410705 -0.219749
## -2.25385 1.12908 -0.884713 11.3085 -0.685642 0.0938608
## -2.03974 1.30856 -0.612535 11.2417 -0.863751 0.458416
## -1.81093 1.484 -0.337789 11.1615 -0.952501 0.845642
## -1.56434 1.65443 -0.0611652 11.07 -0.964541 1.23181
############### NEW alternative approach
# Number of variables
# N = 7
#
# # Initial conditions (convert to Vector)
# u0 = collect(range(0.1, 1; length=N))
#
# # Parameters: F (forcing parameter)
# p0 = [8.0]
#
# # Time span: t ∈ [0, total_time]
# total_time = 12.5
# tspan = (0.0, total_time)
#
# # Define the ODE problem
# prob = ODEProblem(lorenz96_rule!, u0, tspan, p0)
#
# # Solve the ODE
# sampling_time = 0.02
# sol = solve(prob, Tsit5(), saveat=sampling_time)
#
# # Extract the trajectory
# Y = Array(sol) # State variables over time
# t = sol.t # Time points
#
# # Plot the trajectory of the first variable
Plots.plot(t, Y[:, 1], label="u₁(t)", xlabel="Time", ylabel="u₁", title="Lorenz 96 Model")
Visualization of a 7D dynamics is difficult, so we can try to render all \(7\) longitudinal time-series separately.
# using Pkg
# Pkg.add("Makie")
using CairoMakie # or GLMakie for interactive 3D plots
using DynamicalSystems # Assuming you're using DynamicalSystems.jl for StateSpaceSet
using Plots
# Example: Create a StateSpaceSet (replace this with your actual data)
# Y = StateSpaceSet(rand(626, 7)) # 626 points, 7 dimensions
# Extract the 3 columns you want to visualize (e.g., columns 1, 2, and 3)
###### OLD - working version
# Extract the 3 columns you want to visualize (e.g., columns 1, 2, and 3)
col1 = Y[:, 1]## 626-element Vector{Float64}:
## -1.7338557646782273
## -1.1252628607512618
## -0.5172277326951644
## 0.07506309005052285
## 0.6386730571772157
## 1.1634479239128346
## 1.642313152669519
## 2.0712979463442758
## 2.4493049373694427
## 2.7776859740862356
## ⋮
## -3.329824229017429
## -3.0961404856874064
## -2.873661358706411
## -2.6630411082545216
## -2.4590481002588063
## -2.2538525156989797
## -2.039744253948989
## -1.8109303876172593
## -1.5643400262871578## 626-element Vector{Float64}:
## -2.6815114658426222
## -2.538479431629982
## -2.3720491034724906
## -2.1784547890448276
## -1.9564449187619446
## -1.706952647967284
## -1.4326086141640217
## -1.137175861686875
## -0.8249855281819021
## -0.5004331756517273
## ⋮
## 0.06850961359387715
## 0.32555952065364435
## 0.550075423734295
## 0.7537735228348528
## 0.9450210705387911
## 1.1290792187740084
## 1.30856458897821
## 1.4840000125049677
## 1.654427729997732## 626-element Vector{Float64}:
## 0.2991609258130429
## 0.3754840896247301
## 0.4755190928837415
## 0.5959847253960259
## 0.7325441144532393
## 0.8801828263243803
## 1.03368977078233
## 1.1881446107316116
## 1.339326499679735
## 1.4839907379679231
## ⋮
## -2.140302355609136
## -1.9160080022943147
## -1.672470089798286
## -1.4167828514840595
## -1.1533221107212719
## -0.8847127286068622
## -0.6125349102315463
## -0.3377892039638637
## -0.061165219876270556
# Create a 3D plot
fig = Figure()
ax = Axis3(fig[1, 1], xlabel="Column 1", ylabel="Column 2", zlabel="Column 3", title="3D Temporal Dynamics")## Axis3()## Lines{Tuple{Vector{Point{3, Float64}}}}
# Add a colorbar to represent time
Colorbar(fig[1, 2], limits=(1, length(col1)), colormap=:viridis, label="Time")## Colorbar()## CairoMakie.Screen{IMAGE}
# using Plots; gr();
# anim = @animate for i=1:100
# plot(sin.(range(0,i/10*pi,length=1000)), label="", ylims=(-1,1))
# end
# positions[75,2,1:100] refers to the Y (2nd) coordinate of the 75th Trajectory, across the first 100 timesteps of the trajectory.
#
# anim = @animate for n in 1:size(positions,1)
# scatter3d(positions[n,1,time_indices], positions[n,2,time_indices],positions[n,3,time_indices],label="Trajectory $n for times 1 to 100")
# end
# anim = @animate for n in 1:size(Y,1)
# scatter3d(Y[1,n], Y[2,n], Y[3,n],label="Trajectory $n for times 1 to 100")
# end
# gif(anim, "C:/Users/IvoD/Desktop/anim_fps15.gif", fps=15)
##### New revised version
# # total_time = 12.5
# # sampling_time = 0.02
# # Y, t = trajectory(lorenz96, total_time; Ttr = 2.2, deltat = sampling_time)
# # Y
#
# col1 = Y[:, 1]
# col2 = Y[:, 2]
# col3 = Y[:, 3]
#
# # Create a 3D plot
# fig = Figure()
# ax = Axis3(fig[1, 1], xlabel="Column 1", ylabel="Column 2", zlabel="Column 3", title="3D Temporal Dynamics")
#
# # Plot the trajectory
# lines!(ax, col1, col2, col3, color=1:length(col1), colormap=:viridis)
#
# # Add a colorbar to represent time
# Colorbar(fig[1, 2], limits=(1, length(col1)), colormap=:viridis, label="Time")
#
# # Display the figure
# display(fig)And here is a dynamic animation …
using DifferentialEquations
using Plots
# Define the Lorenz 96 system
function lorenz96Anim(du, u, p, t)
N = length(u)
for i in 1:N
du[i] = (u[mod1(i+1, N)] - u[mod1(i-2, N)]) * u[mod1(i-1, N)] - u[i] + p
end
end## lorenz96Anim (generic function with 1 method)## 5## 5-element Vector{Float64}:
## 0.2955459911905558
## 0.6331907893809096
## 0.7753167251225109
## 0.4067531150616067
## 0.7471087103535745## 8.0## (0.0, 12.5)## 0.02## ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
## timespan: (0.0, 12.5)
## u0: 5-element Vector{Float64}:
## 0.2955459911905558
## 0.6331907893809096
## 0.7753167251225109
## 0.4067531150616067
## 0.7471087103535745## retcode: Success
## Interpolation: 1st order linear
## t: 626-element Vector{Float64}:
## 0.0
## 0.02
## 0.04
## 0.06
## 0.08
## 0.1
## 0.12
## 0.14
## 0.16
## 0.18
## ⋮
## 12.34
## 12.36
## 12.38
## 12.4
## 12.42
## 12.44
## 12.46
## 12.48
## 12.5
## u: 626-element Vector{Vector{Float64}}:
## [0.2955459911905558, 0.6331907893809096, 0.7753167251225109, 0.4067531150616067, 0.7471087103535745]
## [0.4517175923676222, 0.7792932689686038, 0.9199027769744416, 0.5589635578802526, 0.8861922385650927]
## [0.6052885015991919, 0.9226628926141974, 1.0618699825992965, 0.7083226099278926, 1.0212521850123069]
## [0.7562649089698789, 1.0634073215220468, 1.201236200975469, 0.8548102462536283, 1.152463201017729]
## [0.9046595077498643, 1.201640312010413, 1.338016601694839, 0.9983993876336604, 1.279995014925642]
## [1.050492582192686, 1.3374808183756541, 1.4722243918148366, 1.1390559765910568, 1.4040113584048106]
## [1.1937944592843555, 1.4710482340352375, 1.6038682389075445, 1.2767409298562467, 1.5246718905755636]
## [1.3346061388660442, 1.6024609082369452, 1.7329513445319136, 1.4114108238917835, 1.6421329824203124]
## [1.4729792936340833, 1.7318361460588698, 1.8594714442337614, 1.5430178948923452, 1.7565477167835493]
## [1.6089771437370395, 1.8592861621290278, 1.9834182058393035, 1.6715122710329045, 1.8680689235108696]
## ⋮
## [2.899207083038217, 0.4020322172676166, -1.6174064408672328, 3.243467625117007, 7.642207955497036]
## [2.544992159121361, 0.050411474362309874, -1.4253321912770482, 3.1103705274688176, 7.916423047783486]
## [2.1526415001234436, -0.22795386620688085, -1.2400328125768654, 2.992921228134161, 8.140638876442862]
## [1.733618629938581, -0.42620780922457524, -1.0632458772430544, 2.8966536456330907, 8.318157658139596]
## [1.3010101853932021, -0.5408316882180917, -0.8966612395413306, 2.8255691829094958, 8.453145194196624]
## [0.867019079377645, -0.5717928419701503, -0.7396688663398242, 2.7814162251126366, 8.549613348456244]
## [0.44339477137724126, -0.5218550985086163, -0.5897097547116515, 2.7645845074235575, 8.61103109296087]
## [0.040790827250845194, -0.3963994354501586, -0.44279541335340367, 2.7752570364511326, 8.64030690452001]
## [-0.33163974837405685, -0.20346950555561272, -0.293301449384258, 2.8134213891326807, 8.639702980720406]## 5×626 Matrix{Float64}:
## 0.295546 0.451718 0.605289 0.756265 … 0.443395 0.0407908 -0.33164
## 0.633191 0.779293 0.922663 1.06341 -0.521855 -0.396399 -0.20347
## 0.775317 0.919903 1.06187 1.20124 -0.58971 -0.442795 -0.293301
## 0.406753 0.558964 0.708323 0.85481 2.76458 2.77526 2.81342
## 0.747109 0.886192 1.02125 1.15246 8.61103 8.64031 8.6397
# Create the animation
anim = @animate for n in 1:size(Y, 2)
# Plot the first 3 dimensions as a 3D trajectory
plot3d(Y[1, 1:n], Y[2, 1:n], Y[3, 1:n],
label="Trajectory 1-3",
xlabel="X₁", ylabel="X₂", zlabel="X₃",
title="Lorenz 96 System Dynamics (First 3 Dimensions)",
legend=:topright,
linecolor=:blue)
# Plot dimensions 2-4 as another 3D trajectory
plot3d!(Y[2, 1:n], Y[3, 1:n], Y[4, 1:n],
label="Trajectory 2-4",
linecolor=:red)
# Plot dimensions 3-5 as another 3D trajectory
plot3d!(Y[3, 1:n], Y[4, 1:n], Y[5, 1:n],
label="Trajectory 3-5",
linecolor=:green)
end## Animation("C:\\Users\\IvoD\\AppData\\Local\\Temp\\jl_ZUiHAq", ["000001.png", "000002.png", "000003.png", "000004.png", "000005.png", "000006.png", "000007.png", "000008.png", "000009.png", "000010.png" … "000617.png", "000618.png", "000619.png", "000620.png", "000621.png", "000622.png", "000623.png", "000624.png", "000625.png", "000626.png"])
# Save the animation as a GIF
gif(anim, "C:/Users/IvoD/Desktop/lorenz96DynamicsViewAnimation.gif", fps=15)## Plots.AnimatedGif("C:\\Users\\IvoD\\Desktop\\lorenz96DynamicsViewAnimation.gif")## Animation("C:\\Users\\IvoD\\AppData\\Local\\Temp\\jl_ZUiHAq", ["000001.png", "000002.png", "000003.png", "000004.png", "000005.png", "000006.png", "000007.png", "000008.png", "000009.png", "000010.png" … "000617.png", "000618.png", "000619.png", "000620.png", "000621.png", "000622.png", "000623.png", "000624.png", "000625.png", "000626.png"])
## Using dynamical
systems
Often in practice, we provide a DynamicalSystem instance
to a library function. Consider the example if computing the Lyapunov
spectrum, which is the function lyapunovspectrum() in
the ChaosTools Julia package, of the Lorenz96
system we showed earlier.
## 15000## 7-element Vector{Float64}:
## 1.2574207696840431
## 0.04012776096217523
## -0.01821497913786168
## -0.6080847812878554
## -1.204380638160897
## -2.200460847553024
## -4.266410408936259The spectrum of Lyapunov exponents is of length equal to the number of dimensions of the phase space. The largest Lyapunov exponent is the maximal Lyapunov exponent (MLE) and it determines a notion of predictability for the dynamical system. In compact phase spaces, an \(MLE > 0\) indicates that the system is chaotic. In our example of the Lorenz96 system, there is at least one positive Lyapunov exponent, because the system is chaotic, and at least one zero Lyapunov exponent, because the system is time-continuous.
All DynamicalSystems library functions work for any
applicable dynamical system. The same lyapunovspectrum()
function works the same for the Henon map example.
## 2-element Vector{Float64}:
## 0.4161998393553483
## -1.6201726436812898For each multistable dynamical system, we can also estimating the
basins of attraction. The Henon map is multistable, as
some initial conditions diverge to infinity, and others converge to a
chaotic attractor. Computing these basins of
attraction is accomplished using the functions
Attractors() and
AttractorsViaRecurrences().
## -2.0:0.02:2.0# find attractors using recurrences in state space:
mapper = AttractorsViaRecurrences(henon, (xg, yg); sparse = false)## AttractorsViaRecurrences
## system: DeterministicIteratedMap
## grid: (-2.0:0.02:2.0, -2.0:0.02:2.0)
## attractors: Dict{Int64, StateSpaceSet{2, Float64}}()# compute the full basins of attraction:
basins, attractors = basins_of_attraction(mapper; show_progress = false)## ([-1 -1 … -1 -1; -1 -1 … -1 -1; … ; -1 -1 … -1 -1; -1 -1 … -1 -1], Dict{Int64, StateSpaceSet{2, Float64}}(1 => 2-dimensional StateSpaceSet{Float64} with 416 points))## CairoMakie.Screen{IMAGE}3.3 Stochastic Differential Equations (SDEs)
Stochastic differential equations (SDEs) are differential equations involving at least one stochastic term, which implies that their solutions are also stochastic processes. SDEs are used to model economic processes, physical systems subjected to thermal fluctuations, and spacekime inference. SDEs involve random differentials such as random white noise calculated as the distributional derivatives of a Brownian motion process, or more generally a semimartingale. SDEs play a critical role in the theory of Brownian motion. The mathematical theory of SDEs was developed in the 1940’s by Kiyosi Itô, who introduced the stochastic integral and Itô calculus.
The most general SDEs is
\[{\frac {\mathrm {d} x(t)}{\mathrm {d} t}}=F(x(t))+\sum _{\alpha =1}^{n}g_{\alpha }(x(t))\xi ^{\alpha }(t),\,\] where \(x\in X\) is the spatial position in the system in its phase (or state) space, \(X\), which is a differentiable manifold, \(F\in TX\) is the classical deterministic evolution component, as a flow vector field, and \(g_{\alpha }\in TX\) are vector fields defining the coupling of the system to Gaussian white noise \(\xi ^{\alpha }\). For linear spaces \(X\) and constant \(g\), the system only involves additive noise. For more complex spaces and non-stationary coupling vector fields, the system is subject to multiplicative noise. In fixed noise configuration systems, SDEs have unique solutions differentiable with respect to the initial condition. One way of finding SDE solutions involves seeking time-dependent probability distribution functions via equivalent Fokker–Planck equation (FPE) models, which are deterministic partial differential equations. Time-dynamic probability distribution functions evolve similarly to the evolution of a quantum wave solving the Schrödinger equation or to diiffusion equations modeling chemical concentrations. Numerical SDE solutions can be obtained by Monte Carlo simulation, by path integration, or by mapping them to ordinary differential equations (ODEs) for the statistical moments of the underlying probability distribution function, e.g., kime-phase distribution.
The package DynamicalSystems supports stochastic systems
modeled by Stochastic
Differential Equations (SDEs). Similar to CoupledODEs,
we can make stochastic models using CoupledSDEs. Let’s
consider the example of the stochastic
FitzHugh-Nagumo model.
The SDEs of the stochastic FitzHugh-Nagumo model are
\[{\dot {v}}=v-{\frac {v^{3}}{3}}-w+RI_{\rm {ext}} ,\] \[\tau {\dot {w}}=v+a-bw,\] where \(I\) is an external stimulus, the variables \(v\) and \(w\) control the system dynamics and relax back to their rest values. This is a model for neural spike generations, with a short, nonlinear elevation of membrane voltage \(v\), diminished over time by a slower, linear recovery variable \(w\) parameter representing sodium channel reactivation and potassium channel deactivation, after stimulation by an external input current \(I\). The FitzHugh–Nagumo model is a simplified 2D version of the Hodgkin–Huxley model, which models the activation and deactivation dynamics of a spiking neuronal cell.
# First time, install the package StochasticDiffEq
# import Pkg; Pkg.add("StochasticDiffEq")
# using StochasticDiffEq # load extension for `CoupledSDEs`
# using DynamicalSystemsBase, StochasticDiffEq, DiffEqNoiseProcess
#
# function fitzhugh_nagumo(u, p, t)
# x, y = u
# ϵ, beta, alpha, gamma, κ, I = p
# dx = (-alpha * x^3 + gamma * x - κ * y + I) / ϵ
# dy = -beta * y + x
# return SVector(dx, dy)
# end
# p = [1.,3.,1.,1.,1.,0.]
# sde = CoupledSDEs(fitzhugh_nagumo, zeros(2), p; noise_strength = 0.05)
####################################### Test 1
# using StochasticDiffEq # load extension for `CoupledSDEs`
#
# function fitzhugh_nagumo(u, p, t)
# x, y = u
# epsilon, alpha, beta, gamma, kappa, I = p
# dx = (-alpha * x^3 + gamma * x - kappa * y + I) / epsilon
# dy = -beta * y + x
# return SVector(dx, dy)
# end
# p = [1.0, 3.0, 1.0, 1.0, 1.0, 0.0] # initialize the parameters
# sde = CoupledSDEs(fitzhugh_nagumo, zeros(2), p; noise_strength = 0.05)
##############################################
# From Scratch!
# using DifferentialEquations
# using Plots
#
# # Define the FitzHugh-Nagumo model with noise
# function fitzhugh_nagumo_sde!(du, u, p, t)
# a, b, c, noise_strength = p
# v, w = u
# du[1] = v - v^3 / 3 - w + c # Deterministic part for dv/dt
# du[2] = (v + a - b * w) / c # Deterministic part for dw/dt
# du .+= noise_strength * randn(length(u)) # Add noise
# end
#
# # Parameters (a, b, c, noise_strength)
# p = [0.7, 0.8, 3.0, 0.05]
#
# # Initial condition
# u0 = [0.0, 0.0]
#
# # Time span
# tspan = (0.0, 100.0)
#
# # Define the SDE problem
# prob = SDEProblem(fitzhugh_nagumo_sde!, u0, tspan, p)
#
# # Solve the SDE
# sol = solve(prob, EM(), dt=0.01)
#
# # Plot the solution
# plot(sol, label=["v" "w"], xlabel="Time", ylabel="State")
using DifferentialEquations
using Plots
# Define the stochastic FitzHugh-Nagumo model
function fitzhugh_nagumo_sde!(du, u, p, t)
v, w = u # State variables: v = membrane potential, w = recovery variable
a, b, c, I = p # Parameters: a, b, c, and external current I
du[1] = v - v^3 / 3 - w + I # dv/dt: Membrane potential dynamics
du[2] = (v + a - b * w) / c # dw/dt: Recovery variable dynamics
end## fitzhugh_nagumo_sde! (generic function with 1 method)
# Add noise to the system
function fitzhugh_nagumo_noise!(du, u, p, t)
sigma = p[end] # Noise strength (last parameter)
du[1] = sigma # Add noise to dv/dt
du[2] = sigma # Add noise to dw/dt
end## fitzhugh_nagumo_noise! (generic function with 1 method)
# Parameters: a, b, c, external current I, and noise strength sigma
p = [0.7, 0.8, 3.0, 0.5, 0.2] # a, b, c, I, and sigma## 5-element Vector{Float64}:
## 0.7
## 0.8
## 3.0
## 0.5
## 0.2## 2-element Vector{Float64}:
## 0.0
## 0.0## (0.0, 100.0)
# Define the SDE problem
prob = SDEProblem(fitzhugh_nagumo_sde!, fitzhugh_nagumo_noise!, u0, tspan, p)## SDEProblem with uType Vector{Float64} and tType Float64. In-place: true
## timespan: (0.0, 100.0)
## u0: 2-element Vector{Float64}:
## 0.0
## 0.0## retcode: Success
## Interpolation: 1st order linear
## t: 10001-element Vector{Float64}:
## 0.0
## 0.01
## 0.02
## 0.03
## 0.04
## 0.05
## 0.060000000000000005
## 0.07
## 0.08
## 0.09
## ⋮
## 99.92000000001421
## 99.93000000001422
## 99.94000000001422
## 99.95000000001423
## 99.96000000001423
## 99.97000000001424
## 99.98000000001424
## 99.99000000001425
## 100.0
## u: 10001-element Vector{Vector{Float64}}:
## [0.0, 0.0]
## [0.014798299044800046, 0.018921985478741167]
## [0.00291197494080174, 0.043875253495578165]
## [0.015190205781110865, 0.04268274052971455]
## [0.012627189432156828, 0.02796685427182915]
## [0.0015982955911108335, 0.02170853097140699]
## [0.0020157964464222545, 0.08113231227533416]
## [0.021154819641942234, 0.07877310924183119]
## [0.05075526124674008, 0.08186678352771265]
## [0.047885918327024525, 0.10201048130535684]
## ⋮
## [0.7139946327528541, 1.5082646913148974]
## [0.7201034434375284, 1.535631874688153]
## [0.7529038084569413, 1.5293979305089485]
## [0.7890104608825371, 1.5278669745697966]
## [0.7689387565268295, 1.5547462217379682]
## [0.7699952261738678, 1.5470350607902021]
## [0.7853347280961231, 1.5657130442610823]
## [0.7733977024754882, 1.5635551667994787]
## [0.7546990503073823, 1.5553008710991787]
# Plot the solution
plot(sol, label=["v (Membrane Potential)" "w (Recovery Variable)"], xlabel="Time", ylabel="State", title="Stochastic FitzHugh-Nagumo Model")
In this particular example, the SDE noise is white noise (Wiener
process) with strength \(\sigma=0.2\). In contrast to the previous
example of the Henon map, in the FitzHugh-Nagumo Model
we have to use a different approach, because
AttractorsViaRecurrences only works for deterministic
systems. Hence for the stochastic FitzHugh-Nagumo Model we use
AttractorsViaFeaturizing. Let’s compute and plot the
basins of attraction.
# using DynamicalSystems
# using Plots
#
# # Define the FitzHugh-Nagumo model
# function fitzhugh_nagumo(u, p, t)
# v, w = u # State variables: v = membrane potential, w = recovery variable
# a, b, c, I = p # Parameters: a, b, c, and external current I
# dv = v - v^3 / 3 - w + I # dv/dt: Membrane potential dynamics
# dw = (v + a - b * w) / c # dw/dt: Recovery variable dynamics
# return SVector(dv, dw) # Return as a static vector for performance
# end
#
# # Parameters: a, b, c, external current I
# p = [0.7, 0.8, 3.0, 0.5]
#
# # Define the dynamical system
# ds = ContinuousDynamicalSystem(fitzhugh_nagumo, [0.0, 0.0], p)
#
# # Define the grid of initial conditions
# # xrange = range(-2.0, 2.0, length=100) # Range for v (membrane potential)
# # yrange = range(-2.0, 2.0, length=100) # Range for w (recovery variable)
# # grid = StateSpaceSet([x y] for x in xrange, y in yrange)
# # Define the grid of initial conditions
# xrange = range(-2.0, 2.0, length=100) # Range for v (membrane potential)
# yrange = range(-2.0, 2.0, length=100) # Range for w (recovery variable)
# grid = StateSpaceSet(collect([x y] for x in xrange, y in yrange))
#
# # Compute basins of attraction
# basins, attractors = basins_of_attraction(grid, ds)
#
# # Compute basin fractions
# fractions = basins_fractions(basins)
# println("Basin fractions: ", fractions)
#
# # Plot the basins of attraction
# heatmap(xrange, yrange, basins', xlabel="v (Membrane Potential)", ylabel="w (Recovery Variable)", title="Basins of Attraction for FitzHugh-Nagumo Model", color=:viridis)
# using Pkg
# Pkg.add("DynamicalSystems")
# Pkg.add("Plots")
# Pkg.add("StaticArrays")
# Pkg.add("OrdinaryDiffEq") # For numerical integration
# Pkg.add("StochasticDiffEq")
using DynamicalSystems
using Plots
using StaticArrays
using OrdinaryDiffEq
using StochasticDiffEq
# Define the FitzHugh-Nagumo model as a function.
function fitzhugh_nagumo(u, p, t)
v, w = u # State variables: v = membrane potential, w = recovery variable
a, b, c, I = p # Parameters: a, b, c, and external current I
dv = v - v^3 / 3 - w + I # dv/dt: Membrane potential dynamics
dw = (v + a - b * w) / c # dw/dt: Recovery variable dynamics
return SVector(dv, dw) # Return as a static vector for performance
end## fitzhugh_nagumo (generic function with 1 method)## 4-element Vector{Float64}:
## 0.7
## 0.8
## 3.0
## 0.5## 2-element Vector{Float64}:
## 0.0
## 0.0## (0.0, 1000.0)## ODEProblem with uType Vector{Float64} and tType Float64. In-place: false
## timespan: (0.0, 1000.0)
## u0: 2-element Vector{Float64}:
## 0.0
## 0.0
# Define the range for v and w
xrange = range(-2.0, 2.0, length=100) # Range for v (membrane potential)## -2.0:0.04040404040404041:2.0## -2.0:0.04040404040404041:2.0## 100×100 Matrix{SVector{2, Float64}}:
## [-2.0, -2.0] … [-2.0, 1.9596] [-2.0, 2.0]
## [-1.9596, -2.0] [-1.9596, 1.9596] [-1.9596, 2.0]
## [-1.91919, -2.0] [-1.91919, 1.9596] [-1.91919, 2.0]
## [-1.87879, -2.0] [-1.87879, 1.9596] [-1.87879, 2.0]
## [-1.83838, -2.0] [-1.83838, 1.9596] [-1.83838, 2.0]
## [-1.79798, -2.0] … [-1.79798, 1.9596] [-1.79798, 2.0]
## [-1.75758, -2.0] [-1.75758, 1.9596] [-1.75758, 2.0]
## [-1.71717, -2.0] [-1.71717, 1.9596] [-1.71717, 2.0]
## [-1.67677, -2.0] [-1.67677, 1.9596] [-1.67677, 2.0]
## [-1.63636, -2.0] [-1.63636, 1.9596] [-1.63636, 2.0]
## ⋮ ⋱
## [1.67677, -2.0] [1.67677, 1.9596] [1.67677, 2.0]
## [1.71717, -2.0] [1.71717, 1.9596] [1.71717, 2.0]
## [1.75758, -2.0] [1.75758, 1.9596] [1.75758, 2.0]
## [1.79798, -2.0] [1.79798, 1.9596] [1.79798, 2.0]
## [1.83838, -2.0] … [1.83838, 1.9596] [1.83838, 2.0]
## [1.87879, -2.0] [1.87879, 1.9596] [1.87879, 2.0]
## [1.91919, -2.0] [1.91919, 1.9596] [1.91919, 2.0]
## [1.9596, -2.0] [1.9596, 1.9596] [1.9596, 2.0]
## [2.0, -2.0] [2.0, 1.9596] [2.0, 2.0]
# Function to simulate the system and determine the attractor
function simulate_and_classify(u0, prob, p)
# Define the ODE problem with the new initial condition
prob_local = remake(prob, u0=u0)
# Solve the ODE
sol = solve(prob_local, Tsit5(), saveat=0.1, reltol=1e-8, abstol=1e-8)
# Extract the final state (assumed to be the attractor)
final_state = sol.u[end]
# Classify the attractor (you can define your own criteria)
if final_state[1] > 0.5
return 1 # Attractor 1
else
return 2 # Attractor 2
end
end## simulate_and_classify (generic function with 1 method)## 100×100 Matrix{Int64}:
## 2 2 2 2 2 2 2 2 2 2 2 2 2 … 2 2 2 2 2 2 2 2 2 2 2 2
## 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## 2 2 2 2 2 2 2 2 2 2 2 2 2 … 2 2 2 2 2 2 2 2 2 2 2 2
## 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## ⋮ ⋮ ⋮ ⋱ ⋮ ⋮
## 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## 2 2 2 2 2 2 2 2 2 2 2 2 2 … 2 2 2 2 2 2 2 2 2 2 2 2
## 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
# Plot the basins of attraction
# heatmap(xrange, yrange, basins, xlabel="v (Membrane Potential)", ylabel="w (Recovery Variable)", title="Basins of Attraction for FitzHugh-Nagumo Model", color=:viridis)
# heatmap_basins_attractors((xg, yg), basins, attractors)
using Plots
# Define the x and y ranges
xrange = range(-2.0, 2.0, length=100) # Range for v (membrane potential)## -2.0:0.04040404040404041:2.0## -2.0:0.04040404040404041:2.0
# Define the basins matrix (replace this with your actual data)
basins = rand(100, 100) # Example: 100x100 matrix of random values## 100×100 Matrix{Float64}:
## 0.88946 0.343333 0.715793 … 0.802357 0.604082 0.348029
## 0.417951 0.0113225 0.751952 0.649381 0.173549 0.0974214
## 0.967877 0.0958551 0.420411 0.0594227 0.199641 0.299133
## 0.786354 0.841548 0.538761 0.284152 0.95603 0.51191
## 0.882792 0.794831 0.342982 0.798352 0.55564 0.617722
## 0.990971 0.512981 0.99587 … 0.315872 0.2895 0.212488
## 0.0746646 0.621601 0.482256 0.321076 0.773514 0.775227
## 0.497887 0.963165 0.912941 0.808592 0.928838 0.109489
## 0.550093 0.0291144 0.619126 0.209229 0.487867 0.476499
## 0.0433414 0.0906443 0.977631 0.224518 0.649935 0.318284
## ⋮ ⋱
## 0.666349 0.765628 0.972619 0.598654 0.0221437 0.879146
## 0.57296 0.42572 0.197814 0.934993 0.752577 0.0024709
## 0.130771 0.185692 0.496321 0.492327 0.912518 0.0569921
## 0.564334 0.750297 0.092093 0.663683 0.356375 0.819193
## 0.248601 0.869216 0.832932 … 0.260656 0.756928 0.207685
## 0.498412 0.569837 0.924251 0.940506 0.116085 0.430539
## 0.0766749 0.156115 0.981725 0.532198 0.756552 0.933422
## 0.397423 0.558243 0.995275 0.549091 0.726841 0.156055
## 0.664426 0.365284 0.16505 0.0388762 0.171378 0.739621
# Plot the heatmap
using Plots
Plots.heatmap(xrange, yrange, basins,
xlabel="v (Membrane Potential)",
ylabel="w (Recovery Variable)",
title="Basins of Attraction for FitzHugh-Nagumo Model",
color=:viridis
)
Finally, we will explore the multistability of the stochastic
FitzHugh-Nagumo model using
AttractorsViaFeaturizing.
using DynamicalSystems
using Plots
using StaticArrays
using OrdinaryDiffEq
using StochasticDiffEq
# Define the stochastic FitzHugh-Nagumo model
function fitzhugh_nagumo_sde(du, u, p, t)
v, w = u # State variables: v = membrane potential, w = recovery variable
a, b, c, I = p # Parameters: a, b, c, and external current I
dv = v - v^3 / 3 - w + I # dv/dt: Membrane potential dynamics
dw = (v + a - b * w) / c # dw/dt: Recovery variable dynamics
du[1] = dv
du[2] = dw
end## fitzhugh_nagumo_sde (generic function with 1 method)## noise_term (generic function with 1 method)## 4-element Vector{Float64}:
## 0.7
## 0.8
## 3.0
## 0.5## 2-element Vector{Float64}:
## 0.0
## 0.0## (0.0, 1000.0)## SDEProblem with uType Vector{Float64} and tType Float64. In-place: true
## timespan: (0.0, 1000.0)
## u0: 2-element Vector{Float64}:
## 0.0
## 0.0
# Define the range for v and w
xrange = range(-2.0, 2.0, length=100) # Range for v (membrane potential)## -2.0:0.04040404040404041:2.0## -2.0:0.04040404040404041:2.0
# Create a grid of initial conditions (use mutable arrays)
grid = [[x, y] for x in xrange, y in yrange]## 100×100 Matrix{Vector{Float64}}:
## [-2.0, -2.0] … [-2.0, 1.9596] [-2.0, 2.0]
## [-1.9596, -2.0] [-1.9596, 1.9596] [-1.9596, 2.0]
## [-1.91919, -2.0] [-1.91919, 1.9596] [-1.91919, 2.0]
## [-1.87879, -2.0] [-1.87879, 1.9596] [-1.87879, 2.0]
## [-1.83838, -2.0] [-1.83838, 1.9596] [-1.83838, 2.0]
## [-1.79798, -2.0] … [-1.79798, 1.9596] [-1.79798, 2.0]
## [-1.75758, -2.0] [-1.75758, 1.9596] [-1.75758, 2.0]
## [-1.71717, -2.0] [-1.71717, 1.9596] [-1.71717, 2.0]
## [-1.67677, -2.0] [-1.67677, 1.9596] [-1.67677, 2.0]
## [-1.63636, -2.0] [-1.63636, 1.9596] [-1.63636, 2.0]
## ⋮ ⋱
## [1.67677, -2.0] [1.67677, 1.9596] [1.67677, 2.0]
## [1.71717, -2.0] [1.71717, 1.9596] [1.71717, 2.0]
## [1.75758, -2.0] [1.75758, 1.9596] [1.75758, 2.0]
## [1.79798, -2.0] [1.79798, 1.9596] [1.79798, 2.0]
## [1.83838, -2.0] … [1.83838, 1.9596] [1.83838, 2.0]
## [1.87879, -2.0] [1.87879, 1.9596] [1.87879, 2.0]
## [1.91919, -2.0] [1.91919, 1.9596] [1.91919, 2.0]
## [1.9596, -2.0] [1.9596, 1.9596] [1.9596, 2.0]
## [2.0, -2.0] [2.0, 1.9596] [2.0, 2.0]
# Function to simulate the system and determine the attractor
function simulate_and_classify(u0, prob, p)
# Define the SDE problem with the new initial condition
prob_local = remake(prob, u0=u0)
# Solve the SDE
sol = solve(prob_local, EM(), dt=0.1, saveat=0.1, reltol=1e-8, abstol=1e-8)
# Extract the final state (assumed to be the attractor)
final_state = sol.u[end]
# Classify the attractor (you can define your own criteria)
if final_state[1] > 0.5
return 1 # Attractor 1
else
return 2 # Attractor 2
end
end## simulate_and_classify (generic function with 1 method)## 100×100 Matrix{Int64}:
## 1 2 1 2 2 2 2 2 2 1 2 1 2 … 2 2 2 2 1 1 2 1 2 2 2 1
## 1 1 2 2 1 2 2 2 2 2 2 2 2 2 2 2 2 1 1 2 2 2 2 1 2
## 1 2 2 2 2 2 2 1 2 2 2 2 1 2 2 2 2 2 2 2 2 2 2 2 2
## 2 2 2 1 2 2 2 1 1 1 2 2 2 2 2 2 2 2 2 1 1 2 2 2 2
## 2 2 1 2 2 1 1 2 1 2 2 2 1 2 2 1 2 2 2 1 1 1 2 2 2
## 1 1 2 2 2 2 2 2 1 2 2 2 2 … 2 1 2 1 2 2 2 2 1 2 2 1
## 1 1 2 2 2 2 2 2 1 2 2 2 2 2 2 2 2 1 2 2 2 2 2 1 2
## 2 2 2 2 2 2 1 1 2 2 2 2 2 2 2 2 2 2 2 1 1 2 1 2 1
## 2 2 1 2 2 2 2 1 1 2 2 1 2 2 2 2 2 2 2 2 1 2 2 2 2
## 2 2 2 1 2 1 2 2 1 2 1 2 2 2 2 2 1 2 2 1 2 2 2 2 1
## ⋮ ⋮ ⋮ ⋱ ⋮ ⋮
## 2 2 2 2 2 1 1 1 1 1 2 2 2 2 2 2 1 1 2 2 2 2 2 1 2
## 2 2 2 2 1 2 2 2 1 2 1 2 2 2 2 2 1 2 2 2 2 1 2 2 2
## 2 1 2 2 1 2 2 2 2 2 2 2 2 2 2 2 2 2 1 2 2 1 2 2 1
## 1 2 2 2 2 2 1 2 2 1 1 2 1 2 2 2 2 2 2 2 1 2 2 2 1
## 2 2 2 2 2 2 1 2 2 1 2 2 1 … 1 2 2 2 2 2 1 2 2 2 2 2
## 2 2 2 2 2 2 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## 2 2 2 1 1 2 2 2 2 2 2 2 1 2 2 2 1 2 2 2 2 2 1 2 2
## 1 2 2 1 2 2 2 1 1 1 2 2 2 2 2 2 2 2 1 2 1 2 2 2 2
## 2 2 1 1 2 2 2 1 1 2 2 2 2 2 2 1 2 1 2 2 2 2 2 2 2
# Plot the basins of attraction
Plots.heatmap(xrange, yrange, basins, xlabel="v (Membrane Potential)", ylabel="w (Recovery Variable)", title="Basins of Attraction for Stochastic FitzHugh-Nagumo Model", color=:viridis)
The mathematical concept of attractors doesn’t translate trivially to stochastic systems, which can intrinsically be highly unstable. However, we can try to improve this protocol to make this system (stochastic FitzHugh-Nagumo model) yield a few fixed point attractors, rather than the above results with almost random noise attractor points in the plane. The goal is to obtain attractors that are more focal, i.e., only mildly perturbed by the noise.
A revised approach may include the following improvements:
- Reduce the Noise Level, as high noise may dominate the system’s dynamics and make it difficult to identify the underlying deterministic attractors.
- Increase Integration Time, as stochastic systems may take longer to converge to their invariant measures (stochastic attractors). By increasing the integration time (\(T\)), the system may to settle to some fixed points (basin attractors).
- Use Robust Solvers, e.g., the Euler-Maruyama (EM) method is simple but may not be accurate for stiff or highly nonlinear systems. Using a more robust solver, e.g., SOSRI or SRIW1, designed for stochastic differential equations may be beneficial.
- Average Over Multiple Realizations, since a single trajectory of a stochastic system may not fully represent the system’s behavior. Averaging over multiple realizations (e.g., using Monte Carlo simulations) may help identify the most probable attractors.
- Use of Clustering to Identify Attractors, instead of classifying attractors based on a simple threshold (e.g., final_state[1] > 0.5), we can use clustering algorithms (e.g., k-means) to identify the most probable attractors from the final states of multiple trajectories
- Improved Results Visualizations.
Running this code natively on a single core will take too long to complete. To significantly speed up the calculations, we will parallelize the code using Julia’s built-in support for multi-threading or distributed computing. Suppose we have \(10\) cores available, we can leverage multi-threading to parallelize the computation of the basins of attraction.
Note that in Julia, the number of cores
(threads) to use is specified outside the code, either by setting
an environment variable before starting Julia, or by using
a command-line argument. Once Julia is started with a
specific number of threads, the Threads.@threads macro
automatically uses all available threads. There are multipl ealteratives
ways to specify the number of cores (threads), e.g.,
- Set the Environment Variable: Set the
JULIA_NUM_THREADSenvironment variable before startingJulia. For example, in a terminal: \(>\ export\ JULIA\_NUM\_THREADS=10; \ julia\), or \(\$env:JULIA\_NUM\_THREADS=10\). - Use a Command-Line Argument: \(julia {\text{ --threads }} 10\).
- Set Threads Programmatically (after starting
Julia), e.g., \(Threads.nthreads() = 10\),
# First time - install Julia package from R
# import Pkg; Pkg.add("Clustering")
# Pkg.add("DynamicalSystems")
# Pkg.status(["DynamicalSystems", "CairoMakie", "OrdinaryDiffEq", "BenchmarkTools"]; mode = Pkg.PKGMODE_MANIFEST)
# import Pkg; Pkg.add("Clustering")
# ################################ Native Singe-Core Protocol ###########################
# using DynamicalSystems
# using Plots
# using OrdinaryDiffEq
# using StochasticDiffEq
# using Clustering
#
# # Define the stochastic FitzHugh-Nagumo model
# function fitzhugh_nagumo_sde(du, u, p, t)
# v, w = u # State variables: v = membrane potential, w = recovery variable
# a, b, c, I = p # Parameters: a, b, c, and external current I
# dv = v - v^3 / 3 - w + I # dv/dt: Membrane potential dynamics
# dw = (v + a - b * w) / c # dw/dt: Recovery variable dynamics
# du[1] = dv
# du[2] = dw
# end
#
# function noise_term(du, u, p, t)
# du[1] = 0.01 # Reduced noise added to v
# du[2] = 0.01 # Reduced noise added to w
# end
#
# # Parameters: a, b, c, external current I
# p = [0.7, 0.8, 3.0, 0.5]
#
# # Initial condition (use a mutable array)
# u0 = [0.0, 0.0] # [v, w]
#
# # Define the SDE problem
# tspan = (0.0, 10000.0) # Longer time span for integration
# prob = SDEProblem(fitzhugh_nagumo_sde, noise_term, u0, tspan, p)
#
# # Define the range for v and w
# xrange = range(-2.0, 2.0, length=100) # Range for v (membrane potential)
# yrange = range(-2.0, 2.0, length=100) # Range for w (recovery variable)
#
# # Create a grid of initial conditions (use mutable arrays)
# grid = [[x, y] for x in xrange, y in yrange]
#
# # Function to simulate the system and determine the attractor
# function simulate_and_classify(u0, prob, p, n_realizations=10)
# attractor_counts = Dict{Int, Int}() # Count how often each attractor is reached
# for _ in 1:n_realizations
# # Define the SDE problem with the new initial condition
# prob_local = remake(prob, u0=u0)
#
# # Solve the SDE
# sol = solve(prob_local, SOSRI(), dt=0.1, saveat=0.1, reltol=1e-8, abstol=1e-8)
#
# # Extract the final state (assumed to be the attractor)
# final_state = sol.u[end]
#
# # Classify the attractor (you can define your own criteria)
# if final_state[1] > 0.5
# attractor = 1 # Attractor 1
# else
# attractor = 2 # Attractor 2
# end
#
# # Update the count for this attractor
# attractor_counts[attractor] = get(attractor_counts, attractor, 0) + 1
# end
#
# # Return the most frequent attractor
# return argmax(attractor_counts)
# end
#
# # Compute basins of attraction
# basins = [simulate_and_classify(u0, prob, p) for u0 in grid]
#
# # Simulate multiple trajectories and collect final states
# final_states = []
# for u0 in grid
# prob_local = remake(prob, u0=u0)
# sol = solve(prob_local, SOSRI(), dt=0.1, saveat=0.1, reltol=1e-8, abstol=1e-8)
# push!(final_states, sol.u[end])
# end
#
# # Identify attractors using clustering
# attractors = identify_attractors(final_states)
# println("Identified attractors: ", attractors)
#
# # Plot the basins of attraction
# heatmap(xrange, yrange, basins, xlabel="v (Membrane Potential)", ylabel="w (Recovery Variable)", title="Basins of Attraction for Stochastic FitzHugh-Nagumo Model", color=:viridis)
#
# # Overlay the identified attractors
# scatter!([a[1] for a in attractors], [a[2] for a in attractors], color=:red, label="Attractors")
################################ Multi-Core Protocol ###########################
using DynamicalSystems
using Plots
using OrdinaryDiffEq
using StochasticDiffEq
using Clustering
using Base.Threads
Threads.nthreads() # Check the current number of threads
Threads.nthreads() = 10 # Set the number of threads (if supported)
# Define the stochastic FitzHugh-Nagumo model
function fitzhugh_nagumo_sde(du, u, p, t)
v, w = u # State variables: v = membrane potential, w = recovery variable
a, b, c, I = p # Parameters: a, b, c, and external current I
dv = v - v^3 / 3 - w + I # dv/dt: Membrane potential dynamics
dw = (v + a - b * w) / c # dw/dt: Recovery variable dynamics
du[1] = dv
du[2] = dw
end
function noise_term(du, u, p, t)
du[1] = 0.01 # Reduced noise added to v
du[2] = 0.01 # Reduced noise added to w
end
# Parameters: a, b, c, external current I
p = [0.7, 0.8, 3.0, 0.5]
# Initial condition (use a mutable array)
u0 = [0.0, 0.0] # [v, w]
# Define the SDE problem
tspan = (0.0, 10000.0) # Longer time span for integration
prob = SDEProblem(fitzhugh_nagumo_sde, noise_term, u0, tspan, p)
# Define the range for v and w
xrange = range(-2.0, 2.0, length=100) # Range for v (membrane potential)
yrange = range(-2.0, 2.0, length=100) # Range for w (recovery variable)
# Create a grid of initial conditions (use mutable arrays)
grid = [[x, y] for x in xrange, y in yrange]
# Function to simulate the system and determine the attractor
function simulate_and_classify(u0, prob, p, n_realizations=10)
attractor_counts = Dict{Int, Int}() # Count how often each attractor is reached
for _ in 1:n_realizations
# Define the SDE problem with the new initial condition
prob_local = remake(prob, u0=u0)
# Solve the SDE
sol = solve(prob_local, SOSRI(), dt=0.1, saveat=0.1, reltol=1e-8, abstol=1e-8)
# Extract the final state (assumed to be the attractor)
final_state = sol.u[end]
# Classify the attractor (you can define your own criteria)
if final_state[1] > 0.5
attractor = 1 # Attractor 1
else
attractor = 2 # Attractor 2
end
# Update the count for this attractor
attractor_counts[attractor] = get(attractor_counts, attractor, 0) + 1
end
# Return the most frequent attractor
return argmax(attractor_counts)
end
# Preallocate the basins array
basins = zeros(Int, length(xrange), length(yrange))
# Parallelize the computation of basins
Threads.@threads for i in eachindex(grid)
u0 = grid[i]
basins[i] = simulate_and_classify(u0, prob, p)
end
# Preallocate the final states array
final_states = Vector{Vector{Float64}}(undef, length(grid))
# Parallelize the collection of final states
Threads.@threads for i in eachindex(grid)
u0 = grid[i]
prob_local = remake(prob, u0=u0)
sol = solve(prob_local, SOSRI(), dt=0.1, saveat=0.1, reltol=1e-8, abstol=1e-8)
final_states[i] = sol.u[end]
end
# Identify attractors using clustering
attractors = identify_attractors(final_states)
println("Identified attractors: ", attractors)
# Plot the basins of attraction
heatmap(xrange, yrange, basins, xlabel="v (Membrane Potential)", ylabel="w (Recovery Variable)", title="Basins of Attraction for Stochastic FitzHugh-Nagumo Model", color=:viridis)
# Overlay the identified attractors
scatter!([a[1] for a in attractors], [a[2] for a in attractors], color=:red, label="Attractors")4 Biocomplexity Examples
Let’s experiment with several concrete hands-on examples of stochastic differential equation (SDE) models of biocomplexity.
4.1 FitzHugh-Nagumo model of the dynamics of neuronal excitability
Let’s extend this model to include stochasticity (noise) to simulate the effects of random fluctuations in the system. The FitzHugh-Nagumo model describes the interaction between the membrane potential (\(v\)) and a recovery variable (\(w\)) in a neuron. The deterministic version of the model is given by \[\frac{dv}{dt} = v - \frac{v^3}{3} - w + I, \] \[ \frac{dw}{dt} = \frac{v + a - b w}{c} ,\] where \(v\) is the membrane potential, \(w\) is the potential recovery variable, \(I\) is the external current input, and the dynamic model parameters are \(a\), \(b\), and \(c\).
We’ll add stochasticity to the model by introducing noise terms to both equations \[dv = \left(v - \frac{v^3}{3} - w + I\right) dt + \sigma_v dW_v ,\] \[dw = \left(\frac{v + a - b w}{c}\right) dt + \sigma_w dW_w ,\] where \(\sigma_v\), \(\sigma_w\) are the noise intensities and \(dW_v\) and \(dW_w\) are Wiener processes (Brownian motion).
We’ll use Julia’s DifferentialEquations.jl package to
solve the SDEs and Plots.jl for visualization. Also, we’ll
define the parameters, initial conditions, and time span for the
simulation before we use the solve() function to
numerically solve the SDE. Then, we’ll display the time series of the
membrane potential (\(v\)) and the
recovery variable (\(w\)), and
visualize the dynamics in the phase space (\(v\) vs. \(w\)). The plot showing the membrane
potential (\(v(t)\)) and recovery
variable (\(w(t)\)) over time will have
natural stochastic fluctuations. The plot showing the trajectory of the
system in the (\(v\), \(w\)) plane will illustrate the effect of
noise on the system’s dynamics. Note that the noise terms (\(\sigma_v\), \(\sigma_w\)) introduce random fluctuations,
making the system’s behavior more realistic and complex. The phase space
plot will reveal how the system’s state evolves over time, with noise
causing deviations from the deterministic trajectory. This model
demonstrates how stochasticity can influence the dynamics of a
biological system, modeling neuronal excitability.
using DifferentialEquations
using Plots
# Define the drift function (deterministic part)
function fitzhugh_nagumo_drift(du, u, p, t)
v, w = u
a, b, c, I = p[1:4]
du[1] = v - v^3 / 3 - w + I # dv/dt
du[2] = (v + a - b * w) / c # dw/dt
end## fitzhugh_nagumo_drift (generic function with 1 method)
# Define the diffusion function (stochastic part)
function fitzhugh_nagumo_diffusion(du, u, p, t)
sigma_v, sigma_w = p[end-1:end] # Extract noise intensities
du[1] = sigma_v # Noise for dv/dt
du[2] = sigma_w # Noise for dw/dt
end## fitzhugh_nagumo_diffusion (generic function with 1 method)
# Parameters: a, b, c, I, sigma_v, sigma_w
p = [0.7, 0.8, 3.0, 0.5, 0.2, 0.2] # a, b, c, I, sigma_v, sigma_w## 6-element Vector{Float64}:
## 0.7
## 0.8
## 3.0
## 0.5
## 0.2
## 0.2## 2-element Vector{Float64}:
## 0.0
## 0.0## (0.0, 100.0)
# Define the SDE problem
sde_prob = SDEProblem(fitzhugh_nagumo_drift, fitzhugh_nagumo_diffusion, u0, tspan, p)## SDEProblem with uType Vector{Float64} and tType Float64. In-place: true
## timespan: (0.0, 100.0)
## u0: 2-element Vector{Float64}:
## 0.0
## 0.0## retcode: Success
## Interpolation: 1st order linear
## t: 17018-element Vector{Float64}:
## 0.0
## 0.006217268199106658
## 0.00746072183892799
## 0.008859607183726987
## 0.01043335319662586
## 0.012203817461137093
## 0.014195589758712228
## 0.016436333593484256
## 0.018957170407602787
## 0.021793111823486137
## ⋮
## 99.86209244523701
## 99.88039063187713
## 99.8984627442468
## 99.9173065353057
## 99.93621242002448
## 99.95548448772126
## 99.9746922972487
## 99.99350587333002
## 100.0
## u: 17018-element Vector{Vector{Float64}}:
## [0.0, 0.0]
## [-0.006260203204282027, 0.146325997573593]
## [-0.005297807738590495, 0.14608115957297643]
## [-0.01052309953157431, 0.15110615406208128]
## [-0.025543172095291174, 0.1542259788128236]
## [-0.03268519060974071, 0.16800683996628119]
## [-0.04805183596400585, 0.16817208884745788]
## [-0.04578915511656342, 0.16308382905278498]
## [-0.040714331985865715, 0.17531147231509975]
## [-0.04872996678350072, 0.1734997356865586]
## ⋮
## [-1.1953833329387102, -0.12039114574985439]
## [-1.1664753740098468, -0.1449194971299323]
## [-1.1850476817405582, -0.1532538867001277]
## [-1.1574319055798818, -0.14798257022641767]
## [-1.1834241872583833, -0.15396871707855364]
## [-1.1976456253873498, -0.14601125127940637]
## [-1.2172444206447839, -0.14269059069909723]
## [-1.1685142685505379, -0.11511445888258762]
## [-1.16161382679322, -0.13505816389556]
# Plot the time series
Plots.plot(sol, vars=(0, 1), label="v(t)", xlabel="Time", ylabel="Membrane Potential (v)", title="FitzHugh-Nagumo Model with Noise")
# Phase space plot
Plots.plot(sol, vars=(1, 2), label="Phase Space", xlabel="v", ylabel="w", title="Phase Space of FitzHugh-Nagumo Model")
4.2 Lotka-Volterra model of the ecosystem dynamics of predator-prey interactions
Let’s explore another example of stochastic differential equation (SDE) modeling, this time using the Lotka-Volterra model, which describes the dynamics of predator-prey interactions in an ecosystem. We’ll add stochasticity to this model to simulate the effects of random fluctuations in the population sizes of predators and prey.
The Lotka-Volterra model describes the interaction between two species: prey (\(x\)) (rabbits) and predators (\(y\)) (wolves). The deterministic version of the model is \[\frac{dx}{dt} = \alpha x - \beta x y ,\] \[\frac{dy}{dt} = \delta x y - \gamma y ,\] where \(x\) is the population size of prey, \(y\) is hte population size of predators, \(\alpha\) is the prey growth rate, \(\beta\) is the predation rate, \(\delta\) is hte predator growth rate due to predation, and \(\gamma\) is the predator death rate.
Adding stochasticity to the model is accomplished by introducing noise terms to both equations \[dx = (\alpha x - \beta x y) dt + \sigma_x x dW_x ,\] \[dy = (\delta x y - \gamma y) dt + \sigma_y y dW_y ,\] where \(\sigma_x\) and \(\sigma_y\) are the noise intensities, and \(dW_x\) and \(dW_y\) are Wiener processes (Brownian motion).
We’ll use the Julia package DifferentialEquations.jl to
solve the SDEs and Plots.jl for to display the
solutions.
Prior to numerically solving the SDE, using solve(),
we’ll need to formulate the drift function (deterministic part
of the SDE) and the diffusion function (stochastic part)
functions of the Lotka-Volterra model, as well as define the parameters,
initial conditions, and time span for the simulation. Finally, we’ll
plot the time series of the prey (\(x\)) and predator (\(y\)) populations and visualize the dynamics
in the phase space (\(x\) vs. \(y\)).
using DifferentialEquations
using Plots
# Define the drift function (deterministic part)
function lotka_volterra_drift(du, u, p, t)
x, y = u
alpha, beta, delta, gamma = p[1:4]
du[1] = alpha * x - beta * x * y # dx/dt
du[2] = delta * x * y - gamma * y # dy/dt
end## lotka_volterra_drift (generic function with 1 method)
# Define the diffusion function (stochastic part)
function lotka_volterra_diffusion(du, u, p, t)
x, y = u
sigma_x, sigma_y = p[end-1:end] # Extract noise intensities
du[1] = sigma_x * x # Noise for dx/dt
du[2] = sigma_y * y # Noise for dy/dt
end## lotka_volterra_diffusion (generic function with 1 method)
# Parameters: alpha, beta, delta, gamma, sigma_x, sigma_y
p = [1.0, 0.1, 0.1, 1.0, 0.1, 0.1] # alpha, beta, delta, gamma, sigma_x, sigma_y## 6-element Vector{Float64}:
## 1.0
## 0.1
## 0.1
## 1.0
## 0.1
## 0.1## 2-element Vector{Float64}:
## 10.0
## 5.0## (0.0, 50.0)
# Define the SDE problem
sde_prob = SDEProblem(lotka_volterra_drift, lotka_volterra_diffusion, u0, tspan, p)## SDEProblem with uType Vector{Float64} and tType Float64. In-place: true
## timespan: (0.0, 50.0)
## u0: 2-element Vector{Float64}:
## 10.0
## 5.0## retcode: Success
## Interpolation: 1st order linear
## t: 17482-element Vector{Float64}:
## 0.0
## 0.007144386769140138
## 0.008573264122968166
## 0.010180751146024698
## 0.011989174046963294
## 0.014023649810519217
## 0.01631243504451963
## 0.018881527101617356
## 0.02170568294826529
## 0.024789259577085315
## ⋮
## 49.99161716469057
## 49.99275919038803
## 49.99389807407973
## 49.99503060966652
## 49.99616193599755
## 49.997290735912685
## 49.998417690784436
## 49.999545572654796
## 50.0
## u: 17482-element Vector{Vector{Float64}}:
## [10.0, 5.0]
## [9.369746664976073, 5.296491632109372]
## [9.364453502590955, 5.301870561355804]
## [9.329116200512367, 5.274024340898211]
## [9.358722825391986, 5.261085829675542]
## [9.420782495371835, 5.242307541518315]
## [9.38646007130484, 5.235650939647902]
## [9.365856138605231, 5.246631435636469]
## [9.495007812888003, 5.269826689964181]
## [9.415420687059804, 5.309723093050032]
## ⋮
## [32.60463720276907, 3.1544318516813745]
## [32.89945992272652, 3.1542066076459903]
## [32.96941753641505, 3.1595968035648316]
## [32.909346053063594, 3.1717444230253]
## [32.975975376696, 3.1589337145288923]
## [32.818848877145186, 3.1689362104558305]
## [32.82730121528162, 3.1814280734560914]
## [32.96380392264878, 3.1980111438061467]
## [33.068568636882, 3.195920160584027]
# Plot the time series
plot(sol, vars=(0, 1), label="Prey (x)", xlabel="Time", ylabel="Population Size", title="Lotka-Volterra Model with Noise")
# Phase space plot
Plots.plot(sol, vars=(1, 2), label="Phase Space", xlabel="Prey (x)", ylabel="Predator (y)", title="Phase Space of Lotka-Volterra Model")
The time series plot shows the population sizes of prey (\(x(t)\)) and predators (\(y(t)\)) over time, with stochastic fluctuations. While the phase space plot illustrates the trajectory of the system in the (\(x\), \(y\)) plane, depicting the effect of noise on the predator-prey dynamics. The noise terms (\(\sigma_x\), \(\sigma_y\)) introduce random fluctuations, making the system’s behavior more realistic and complex. This model demonstrates how stochasticity can influence the dynamics of ecological systems, such as predator-prey interactions.
4.3 SIR model of infection spread in a population
This example of biocomplexity modeles a dynamical system with SDEs using the SIR model, which describes the spread of infectious diseases in a population. We’ll add stochasticity to this model to simulate the effects of random fluctuations in the transmission and recovery processes.
The SIR model divides a population into three compartments:
- Susceptible (\(S\)): Individuals who can contract the disease.
- Infected (\(I\)): Individuals who are currently infected and can spread the disease.
- Recovered (\(R\)): Individuals who have recovered from the disease and are immune.
The deterministic version of the model is \[\frac{dS}{dt} = -\beta S I ,\] \[\frac{dI}{dt} = \beta S I - \gamma I ,\] \[\frac{dR}{dt} = \gamma I ,\] where \(\beta\) is the transmission rate and \(\gamma\) is the recovery rate.
Adding stochasticity to the model requires noise in the transmission and recovery processes. \[dS = -\beta S I \, dt + \sigma_S S I \, dW_S ,\] \[dI = (\beta S I - \gamma I) \, dt + \sigma_I S I \, dW_I - \sigma_R I \, dW_R ,\] \[dR = \gamma I \, dt + \sigma_R I \, dW_R ,\] where \(\sigma_S\), \(\sigma_I\), and \(\sigma_R\) are the noise intensities and \(dW_S\), \(dW_I\), and \(dW_R\) are Wiener processes (Brownian motion).
In Julia, we can use the
DifferentialEquations.jl to solve the SDEs and
Plots.jl for visualization. Also, we’ll define the
drift (deterministic) and diffusion (stochastic)
components of the SIR model, along with the parameters,
initial conditions, and time span for the
simulation.
using DifferentialEquations
using Plots
# Define the drift function (deterministic part)
function sir_drift(du, u, p, t)
S, I, R = u
beta, gamma = p[1:2]
du[1] = -beta * S * I # dS/dt
du[2] = beta * S * I - gamma * I # dI/dt
du[3] = gamma * I # dR/dt
end## sir_drift (generic function with 1 method)
# Define the diffusion function (stochastic part)
function sir_diffusion(du, u, p, t)
S, I, R = u
sigma_S, sigma_I, sigma_R = p[end-2:end] # Extract noise intensities
du[1] = sigma_S * S * I # Noise for dS/dt
du[2] = sigma_I * S * I - sigma_R * I # Noise for dI/dt
du[3] = sigma_R * I # Noise for dR/dt
end## sir_diffusion (generic function with 1 method)
# Parameters: beta, gamma, sigma_S, sigma_I, sigma_R
p = [0.5, 0.1, 0.05, 0.05, 0.05] # beta, gamma, sigma_S, sigma_I, sigma_R## 5-element Vector{Float64}:
## 0.5
## 0.1
## 0.05
## 0.05
## 0.05## 3-element Vector{Float64}:
## 0.99
## 0.01
## 0.0## (0.0, 100.0)## SDEProblem with uType Vector{Float64} and tType Float64. In-place: true
## timespan: (0.0, 100.0)
## u0: 3-element Vector{Float64}:
## 0.99
## 0.01
## 0.0## retcode: Success
## Interpolation: 1st order linear
## t: 581-element Vector{Float64}:
## 0.0
## 0.01
## 0.016334531125769186
## 0.023460878642259517
## 0.03147801959831114
## 0.040497303173869216
## 0.050643997196372054
## 0.06205902797168775
## 0.0749009375939179
## 0.08934808591892683
## ⋮
## 70.16845899451623
## 73.01796630652339
## 76.22366203253145
## 79.83006972429052
## 83.88727837751948
## 88.45163811240205
## 93.58654281414495
## 99.36331060360571
## 100.0
## u: 581-element Vector{Vector{Float64}}:
## [0.99, 0.01, 0.0]
## [0.989896141817806, 0.01004014208168028, 2.0626462469247073e-6]
## [0.9898761018596202, 0.010065087376344858, 4.825307612635316e-5]
## [0.9897386001187528, 0.010093357946073299, 0.00012632086843023744]
## [0.9897056550971711, 0.01012514267136147, 0.00016573655078850542]
## [0.9897350343574617, 0.01016194947922951, 0.00020224341490014295]
## [0.9896654574652254, 0.010202571172162668, 0.0002108074060563141]
## [0.9896963533222479, 0.010249618980172308, 0.00025307044476750737]
## [0.9897220950170164, 0.01030188881355418, 0.00030483909526440016]
## [0.989575123096555, 0.010360763365876942, 0.00024144473818909161]
## ⋮
## [0.0063956261500236495, 0.0034089952453056744, 0.9501997951219978]
## [0.006371649273820325, 0.002156225634143093, 0.9513259455633439]
## [0.006353842518499896, 0.0016672568841753232, 0.9515726007272535]
## [0.006338716286313503, 0.001026243729145374, 0.9520865617605018]
## [0.0063277578700474, 0.0006811132529227019, 0.9524333629607639]
## [0.0063202397292217586, 0.00044528575485144656, 0.9527758874065222]
## [0.006314422954187603, 0.0002726266679622651, 0.9530359413611499]
## [0.006310869221001017, 0.00015467509419389396, 0.953180232810377]
## [0.006310499465243653, 0.00014394942962273725, 0.9531819327221176]
# Plot the time series
Plots.plot(sol, vars=(0, 1), label="Susceptible (S)", xlabel="Time", ylabel="Proportion of Population", title="SIR Model with Noise")
# Phase space plot
Plots.plot(sol, vars=(1, 2), label="Phase Space", xlabel="Susceptible (S)", ylabel="Infected (I)", title="Phase Space of SIR Model")
Again, the time series plot shows the proportions of susceptible (\(S(t)\)), infected (\(I(t)\)), and recovered (\(R(t)\)) individuals over time, with stochastic fluctuations. Whereas the phase space plot depicts the trajectory of the system in the (\(S\), \(I\)) plane, showing the effect of noise on the disease dynamics. The noise terms (\(\sigma_S\), \(\sigma_I\), and \(\sigma_R\)) introduce random fluctuations, making the system’s behavior more realistic and complex. The phase space plot also reveals how the system’s state evolves over time, with noise causing deviations from the deterministic trajectory.
4.4 Stock Market Modeling
This example of stock market modeling also relies on SDEs and uses the Geometric Brownian Motion (GBM) model, which is a classic model for simulating stock prices. The GBM model incorporates both deterministic growth (drift) term and stochastic fluctuations (volatility) to describe the dynamics of stock prices. The Geometric Brownian Motion model describes the evolution of a stock price \(S(t)\) over time. The deterministic version of the model is \[\frac{dS}{S} = \mu \, dt + \sigma \, dW ,\] where \(S\) is the stock price, \(\mu\) is the drift (expected return), \(\sigma\) is the diffusion (stock volatility or standard deviation of returns), and \(dW\) is a Wiener process (Brownian motion).
The SDE form of the model is \(dS = \mu S \, dt + \sigma S \, dW\).
We will use the following parameters, initial conditions, and time span for the simulation.
- Parameters: \(\mu=0.1\) (drift) and \(\sigma=0.2\) (volatility), i.e., \(p = [0.1, 0.2]\).
- Initial condition: \(S(0)=u_o = 100.0\), the initial stock price.
- Time span: \(t\in [0, 10]\), i.e., 10 years.
using DifferentialEquations
using Plots
# Define the drift function (deterministic part)
function gbm_drift(du, u, p, t)
S = u[1]
mu = p[1]
du[1] = mu * S # dS/dt = mu * S
end## gbm_drift (generic function with 1 method)
# Define the diffusion function (stochastic part)
function gbm_diffusion(du, u, p, t)
S = u[1]
sigma = p[2]
du[1] = sigma * S # dS/dt = sigma * S * dW
end## gbm_diffusion (generic function with 1 method)## 2-element Vector{Float64}:
## 0.1
## 0.2## 1-element Vector{Float64}:
## 100.0## (0.0, 10.0)## SDEProblem with uType Vector{Float64} and tType Float64. In-place: true
## timespan: (0.0, 10.0)
## u0: 1-element Vector{Float64}:
## 100.0## retcode: Success
## Interpolation: 1st order linear
## t: 679-element Vector{Float64}:
## 0.0
## 0.01
## 0.012
## 0.01425
## 0.01678125
## 0.01962890625
## 0.02283251953125
## 0.026436584472656253
## 0.030491157531738287
## 0.03505255222320557
## ⋮
## 9.881563073393501
## 9.896501818833192
## 9.911446549730567
## 9.926350554275137
## 9.941459164181662
## 9.956332212976923
## 9.9712690505874
## 9.986171435457463
## 10.0
## u: 679-element Vector{Vector{Float64}}:
## [100.0]
## [99.87102536103765]
## [100.14612833287843]
## [99.38253130455824]
## [99.2461603892151]
## [99.39321257151232]
## [100.0252974569971]
## [99.31123620036281]
## [99.99358031268386]
## [100.1322623308715]
## ⋮
## [161.5457500279727]
## [156.33426540733637]
## [156.25614916420076]
## [161.60352920690087]
## [158.27774659634417]
## [156.77662226493217]
## [156.38067682320113]
## [157.16568095228402]
## [157.25075499492542]
# Plot the stock price over time
Plots.plot(sol, vars=(0, 1), label="Stock Price (S)", xlabel="Time (Years)", ylabel="Stock Price", title="Geometric Brownian Motion Model")
## EnsembleProblem with problem SDEProblem## EnsembleSolution Solution of length 10 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), Nothing, 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(gbm_drift), typeof(gbm_diffusion), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, typeof(gbm_diffusion), Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, Nothing}, SOSRI, StochasticDiffEq.LinearInterpolationData{Vector{Vector{Float64}}, Vector{Float64}}, DiffEqBase.Stats, Nothing}
# Plot the multiple paths - multiple alternative realizations (all possible, with different likelihoods)
Plots.plot(ensemble_sol, vars=(0, 1), xlabel="Time (Years)", ylabel="Stock Price", title="Multiple Paths of GBM Model", legend=false)
The single path plot shows the simulated stock price over time for a single realization of the GBM model. The multiple paths plot shows multiple realizations of the stock price and demonstrates the stochastic nature of the model.
The drift term controls the expected return of the stock. A higher \(\mu\) leads to faster growth in the stock price. The dissusion, volatility term represents the randomness in the stock price. A higher \(\sigma\) leads to larger fluctuations. The Wiener process (\(dW\)) introduces randomness, making the stock price unpredictable in the short term but following a trend in the long term. Simulating multiple paths helps visualize the range of possible outcomes and the impact of randomness on the stock price.
The GBM model is widely used in the Black-Scholes model for pricing financial derivatives. Simulating multiple paths helps assess the risk associated with stock price movements. Finally, the GBM model can be used to simulate the behavior of diverse assets in a portfolio.
4.5 Global carbon dioxide (\(CO_2\)) dynamics in the atmosphere model
This example of global environmental modeling uses stochastic differential equations to develop a simplified model of carbon dioxide (\(CO_2\)) dynamics in the atmosphere. The model incorporates both deterministic processes (e.g., \(CO_2\) emissions and absorption) and stochastic fluctuations (e.g., random variations in emissions or natural processes).
The model describes the \(CO_2\) concentration in the atmosphere, \(C\), over time. The deterministic version of the model is \(\frac{dC}{dt} = E - A C\), where \(E\) is the \(CO_2\) emissions rate (e.g., from human activities), \(A\) is the \(CO_2\) rate of absorption (e.g., by the environment, oceans and forests), and \(C\) is the \(CO_2\) concentration in the atmosphere.
Again, we’ll add stochasticity to the model by introducing noise terms to account for random fluctuations in emissions and absorption \(dC = (E - A C) \, dt + \sigma_E \, dW_E - \sigma_A C \, dW_A ,\) where \(\sigma_E\) is the noise intensity for emissions, \(\sigma_A\) is the noise intensity for absorption, and \(dW_E\) and \(dW_A\) are Wiener processes (Brownian motions).
Define the following parameters, initial conditions, and time span for the simulation.
- Parameters: \(E= 2.0\) (emissions), \(A=0.05\) (absorption), \(\sigma_E=0.5\), \(\sigma_A=0.1\).
- Initial condition: \(C(0)=u_o = 100.0\), the initial \(CO_2\) concentration (in PPM).
- Time span: \(t \in [0, 100]\), i.e., 100 years.
using DifferentialEquations
using Plots
# Define the drift function (deterministic part)
function co2_drift(du, u, p, t)
C = u[1]
E, A = p[1:2]
du[1] = E - A * C # dC/dt = E - A * C
end## co2_drift (generic function with 1 method)
# Define the diffusion function (stochastic part)
function co2_diffusion(du, u, p, t)
C = u[1]
sigma_E, sigma_A = p[end-1:end] # Extract noise intensities
du[1] = sigma_E - sigma_A * C # Noise for dC/dt
end## co2_diffusion (generic function with 1 method)
# Parameters: E (emissions), A (absorption), sigma_E, sigma_A
p = [2.0, 0.05, 0.5, 0.1] # E = 2.0, A = 0.1, sigma_E = 0.5, sigma_A = 0.1## 4-element Vector{Float64}:
## 2.0
## 0.05
## 0.5
## 0.1## 1-element Vector{Float64}:
## 100.0## (0.0, 100.0)## SDEProblem with uType Vector{Float64} and tType Float64. In-place: true
## timespan: (0.0, 100.0)
## u0: 1-element Vector{Float64}:
## 100.0## retcode: Success
## Interpolation: 1st order linear
## t: 1217-element Vector{Float64}:
## 0.0
## 0.08330977077569733
## 0.0999717249308368
## 0.1187164233553687
## 0.1398042090829671
## 0.16352796802651529
## 0.19021040418996943
## 0.21951817830376547
## 0.2512547477463164
## 0.28535365070171215
## ⋮
## 96.53124275464998
## 96.8951827523102
## 97.30461524967795
## 97.76522680921667
## 98.17157642975758
## 98.62871975286609
## 99.08458547185622
## 99.59743440572011
## 100.0
## u: 1217-element Vector{Vector{Float64}}:
## [100.0]
## [100.47160622700267]
## [99.15804703571416]
## [99.35646833586918]
## [99.61201432787996]
## [99.37069054269502]
## [99.48421772748908]
## [99.63084269135705]
## [98.75139441717124]
## [98.60647979595977]
## ⋮
## [38.26166047252316]
## [39.72555936077761]
## [40.22817004044502]
## [39.851763656086305]
## [39.24020007249665]
## [40.27622804648395]
## [39.60644007465196]
## [37.735289310032044]
## [39.56694812331534]
# Plot the CO₂ concentration over time
Plots.plot(sol, vars=(0, 1), label="CO₂ Concentration (C)", xlabel="Time (Years)", ylabel="CO₂ Concentration (ppm)", title="CO₂ Dynamics Model with Noise")
## EnsembleProblem with problem SDEProblem## EnsembleSolution Solution of length 10 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), Nothing, 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(co2_drift), typeof(co2_diffusion), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, typeof(co2_diffusion), Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, Nothing}, SOSRI, StochasticDiffEq.LinearInterpolationData{Vector{Vector{Float64}}, Vector{Float64}}, DiffEqBase.Stats, Nothing}
# Plot the multiple paths
Plots.plot(ensemble_sol, vars=(0, 1), xlabel="Time (Years)", ylabel="CO₂ Concentration (ppm)", title="Multiple Paths of CO₂ Dynamics Model", legend=false)
The single path plot shows the simulated \(CO_2\) concentration over time for a single realization of the model, while the second multiple paths plot demonstrates multiple realizations of the \(CO_2\) concentration, i.e., the stochastic nature of the model.
The emission (\(E\)) term represents the rate at which \(CO_2\) is added to the atmosphere. Higher emissions lead to faster increases in \(CO_2\) concentration. The absorption (\(A\)) term models the rate at which \(CO_2\) is removed from the atmosphere. Higher absorption leads to slower increases in CO₂ concentration. The noise terms (\(\sigma_E\) and \(\sigma_A\)) introduce random fluctuations, making the \(CO_2\) concentration unpredictable in the short term but following a trend in the long term. Simulating multiple paths helps visualize the range of possible outcomes and the impact of randomness on CO₂ concentration. This climate model can be used to simulate the impact of emissions and absorption processes on global \(CO_2\) levels. The model can help assess the effectiveness of policies aimed at reducing emissions or enhancing absorption. Simulating multiple paths helps assess the risk of exceeding critical \(CO_2\) thresholds.
4.5.1 Realistic SDE model of global carbon dioxide (\(CO_2\)) dynamics in the Earth’s atmosphere
Note that the above model predicts constant reduction of the \(CO_2\), since the initial settings of the problem are unrealistic. Let’s refine the model to reflect the Earth’s environmental conditions. To set up a stochastic differential equation (SDE) model for atmospheric \(CO_2\) concentration, \(C(t)\), we need to define realistic initial conditions and parameters that reflect the global carbon cycle for the SDE model \(dC = (E - A C) \, dt + \sigma_E \, dW_E - \sigma_A C \, dW_A,\) where \(C(t)\) is the \(CO_2\) concentration in parts per million (PPM), \(E\) is the global \(CO_2\) emission rate (in PPM/year), \(A\) is the absorption rate coefficient (in 1/year), \(\sigma_E\) and \(\sigma_A\) are the noise intensities for emissions and absorption, respectively, and \(dW_E\) and \(dW_A\) are independent Wiener processes (Brownian motions) representing random fluctuations in emissions and absorption.
- The initial \(CO_2\) concentration (\(C(0) = u_0\)) based on the current atmospheric \(CO_2\) concentration is approximately 420 PPM (as of 2023). For historical context, pre-industrial \(CO_2\) levels were around 280 PPM. Hence, we’ll set \(C(0) = 420 \, \text{PPM}\).
- Global emission rate (\(E\)) of \(CO_2\) are approximately 40 gigatons of \(CO_2\) per year (GtCO₂/year. Converting this to PPM/year uses the conversion factor \(1 Gt\ CO_2 \approx 0.47 PPM\). Thus, we set \(E \approx 40 \times 0.47 = 18.8 \, \text{PPM/year}\).
- The \(CO_2\) absorption rate coefficient (\(A\)) by natural sinks (oceans, forests, etc.) is proportional to the current \(CO_2\) concentration. The global carbon sink removes approximately \(50\%\) of annual emissions, so the absorption rate is roughly half of the emission rate. Thus, \(A \approx \frac{E}{C(0)} = \frac{18.8}{420} \approx 0.045 \, \text{1/year}\).
- The noise terms \(\sigma_E\) and \(\sigma_A\) represent random fluctuations in emissions and absorption, respectively. The emissions noise (\(\sigma_E\)) fluctuate due to economic activity, policy changes, and technological advancements. A reasonable assumption is that emissions fluctuate by about \(10\%\) of the mean emission rate, i.e., \(\sigma_E \approx 0.1 \times E = 0.1 \times 18.8 = 1.88 \, \text{PPM/year}^{1/2}\). The absorption noise (\(\sigma_A\)) also fluctuates due to climate variability (e.g., El Niño, La Niña) and changes in land use. A rational assumption is that absorption fluctuates by about \(5\%\) of the mean absorption rate. Therefore, we set \(\sigma_A \approx 0.05 \times A = 0.05 \times 0.045 = 0.00225 \, \text{1/year}^{1/2}.\)
| Parameter | Symbol | Value | Description |
|---|---|---|---|
| Initial CO₂ | \(C(0)\) | 420 PPM | Initial atmospheric CO₂ concentration. |
| Emission rate | \(E\) | 18.8 PPM/year | Global CO₂ emission rate. |
| Absorption rate | \(A\) | 0.045 1/year | Absorption rate coefficient. |
| Emissions noise | \(\sigma_E\) | 1.88 PPM/year\(^{1/2}\) | Noise intensity for emissions fluctuations. |
| Absorption noise | \(\sigma_A\) | 0.00225 1/year\(^{1/2}\) | Noise intensity for absorption fluctuations. |
Let’s implement a realistic SDE \(C)_2\) prediction model using
Julia’s DifferentialEquations.jl package.
using DifferentialEquations
using Plots
# # Parameters
# CO = 420.0 # Initial CO₂ concentration (PPM)
# E = 18.8 # Emission rate (PPM/year)
# A = 0.045 # Absorption rate coefficient (1/year)
# sigma_E = 1.88 # Emissions noise intensity (PPM/year^0.5)
# sigma_A = 0.00225 # Absorption noise intensity (1/year^0.5)
# Define the drift function (deterministic part)
function co2_drift(du, u, p, t)
C = u[1]
E, A = p[1:2]
du[1] = E - A * C # dC/dt = E - A * C
end## co2_drift (generic function with 1 method)
# Define the diffusion function (stochastic part)
function co2_diffusion(du, u, p, t)
C = u[1]
sigma_E, sigma_A = p[end-1:end] # Extract noise intensities
du[1] = sigma_E - sigma_A * C # Noise for dC/dt
end## co2_diffusion (generic function with 1 method)
# Parameters: E (emissions), A (absorption), sigma_E, sigma_A
p = [18.8, 0.045, 1.88, 0.00225] # E = 18.8, A = 0.045 , sigma_E = 1.88, sigma_A = 0.00225## 4-element Vector{Float64}:
## 18.8
## 0.045
## 1.88
## 0.00225## 1-element Vector{Float64}:
## 420.0## (0.0, 100.0)## SDEProblem with uType Vector{Float64} and tType Float64. In-place: true
## timespan: (0.0, 100.0)
## u0: 1-element Vector{Float64}:
## 420.0## retcode: Success
## Interpolation: 1st order linear
## t: 63-element Vector{Float64}:
## 0.0
## 0.1
## 0.16915405254249202
## 0.24695236165279555
## 0.33447545940188705
## 0.432938944369615
## 0.5437103649583089
## 0.6683282131205895
## 0.8085232923031553
## 0.9662427563835418
## ⋮
## 66.95500123512892
## 71.4032126432156
## 74.71337471645967
## 78.43730704885924
## 82.62673092280876
## 87.15240264021608
## 92.24378332229932
## 97.97158658964297
## 100.0
## u: 63-element Vector{Vector{Float64}}:
## [420.0]
## [419.8869123941257]
## [419.60658016971183]
## [419.8265018686887]
## [419.7325624674787]
## [419.8415087191063]
## [419.22686667926314]
## [418.9194193264947]
## [418.5336006063764]
## [418.3035469902669]
## ⋮
## [421.1462840986261]
## [419.62199204840147]
## [418.7013918557605]
## [418.8691383447791]
## [419.1080395409242]
## [416.6037481923374]
## [418.5218319011145]
## [416.94433903251064]
## [417.2815946967833]
# Plot the CO₂ concentration over time
plot(sol, vars=(0, 1), label="Instance of CO₂ Concentration (C)", xlabel="Time (years)",
ylabel="CO₂ Concentration (PPM)",
title="Realistic Instance of the Earth Atmospheric CO₂ Dynamic SDE\n Model with Stochastic Emissions and Absorption", legend=false)
## EnsembleProblem with problem SDEProblem## EnsembleSolution Solution of length 10 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), Nothing, 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(co2_drift), typeof(co2_diffusion), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, typeof(co2_diffusion), Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, Nothing}, SOSRI, StochasticDiffEq.LinearInterpolationData{Vector{Vector{Float64}}, Vector{Float64}}, DiffEqBase.Stats, Nothing}
# Plot the multiple paths
plot(ensemble_sol, xlabel="Time (years)", ylabel="CO₂ Concentration (PPM)",
title="Multiple Realistic Realizations of the Earth Atmospheric CO₂ Dynamic SDE\n Model with Stochastic Emissions and Absorption", legend=false)
# # # Define the SDE
# # function co2_dynamics(du, u, p, t)
# # E, A, sigma_E, sigma_A = p
# # du[1] = E - A * u[1] # Drift term
# # end
# # function noise_term(du, u, p, t)
# # E, A, sigma_E, sigma_A = p
# # du[1] = sigma_E # Noise for emissions (dW_E term)
# # du[2] = -sigma_A * u[1] # Noise for absorption (dW_A term)
# # end
#
# # Parameters and initial condition
# p = (E, A, sigma_E, sigma_A)
# u0 = [C0] # Initial CO₂ concentration
# tspan = (0.0, 100.0) # Simulate for 100 years
#
# # Define the drift function (deterministic part)
# function co2_drift(du, u, p, t)
# C = u[1]
# E, A = p[1:2]
# du[1] = E - A * C # dC/dt = E - A * C
# end
#
# # Define the diffusion function (stochastic part)
# function co2_diffusion(du, u, p, t)
# C = u[1]
# sigma_E, sigma_A = p[end-1:end] # Extract noise intensities
# du[1] = sigma_E - sigma_A * C # Noise for dC/dt
# end
#
# # Define the SDE problem
# sde_prob = SDEProblem(co2_drift, co2_diffusion, u0, tspan, p)
#
# # Solve the SDE
# sol = solve(sde_prob, SOSRI(), dt=0.1) # Use the SOSRI solver for SDEs
#
# # Plot the solution
# plot(sol, xlabel="Time (years)", ylabel="CO₂ Concentration (PPM)",
# title="Realistic Earth Atmospheric CO₂ Dynamic SDE Model with Stochastic Emissions and Absorption",
# legend=false)The model simulates the evolution of atmospheric \(CO_2\) concentration over time, accounting for both deterministic dynamics (emissions and absorption) and stochastic fluctuations. The noise terms \(\sigma_E\) and \(\sigma_A\) introduce variability in the trajectory, reflecting real-world uncertainties in emissions and absorption processes.