diff --git a/dev/.documenter-siteinfo.json b/dev/.documenter-siteinfo.json index 0d5d5a6d..4955e069 100644 --- a/dev/.documenter-siteinfo.json +++ b/dev/.documenter-siteinfo.json @@ -1 +1 @@ -{"documenter":{"julia_version":"1.10.7","generation_timestamp":"2024-12-31T16:34:38","documenter_version":"1.8.0"}} \ No newline at end of file +{"documenter":{"julia_version":"1.10.7","generation_timestamp":"2024-12-31T16:49:15","documenter_version":"1.8.0"}} \ No newline at end of file diff --git a/dev/examples/gMAM_Maierstein/0a4574e4.png b/dev/examples/gMAM_Maierstein/0a4574e4.png new file mode 100644 index 00000000..4a841ad8 Binary files /dev/null and b/dev/examples/gMAM_Maierstein/0a4574e4.png differ diff --git a/dev/examples/gMAM_Maierstein/5f4d0b05.png b/dev/examples/gMAM_Maierstein/5f4d0b05.png deleted file mode 100644 index 75314020..00000000 Binary files a/dev/examples/gMAM_Maierstein/5f4d0b05.png and /dev/null differ diff --git a/dev/examples/gMAM_Maierstein/66ea99e4.png b/dev/examples/gMAM_Maierstein/66ea99e4.png deleted file mode 100644 index df6c643c..00000000 Binary files a/dev/examples/gMAM_Maierstein/66ea99e4.png and /dev/null differ diff --git a/dev/examples/gMAM_Maierstein/672fe833.png b/dev/examples/gMAM_Maierstein/672fe833.png new file mode 100644 index 00000000..ba5f9842 Binary files /dev/null and b/dev/examples/gMAM_Maierstein/672fe833.png differ diff --git a/dev/examples/gMAM_Maierstein/7eb9dc25.png b/dev/examples/gMAM_Maierstein/7eb9dc25.png deleted file mode 100644 index bcc5e5a5..00000000 Binary files a/dev/examples/gMAM_Maierstein/7eb9dc25.png and /dev/null differ diff --git a/dev/examples/gMAM_Maierstein/cf5c0022.png b/dev/examples/gMAM_Maierstein/cf5c0022.png new file mode 100644 index 00000000..98ea418b Binary files /dev/null and b/dev/examples/gMAM_Maierstein/cf5c0022.png differ diff --git a/dev/examples/gMAM_Maierstein/index.html b/dev/examples/gMAM_Maierstein/index.html index c7770d1d..dc3bc6b7 100644 --- a/dev/examples/gMAM_Maierstein/index.html +++ b/dev/examples/gMAM_Maierstein/index.html @@ -98,14 +98,14 @@ markersize=10) for i in eachindex(fp)] lines!(ax2, reduce(hcat, tr), linewidth=1, color=(:black, 0.2)) -figExample block output

Basins of attraction

Basins of attraction are the regions in the state space that lead to a particular attractor. We can find the basins of attraction using the function basins.

using Attractors
+fig
Example block output

Basins of attraction

Basins of attraction are the regions in the state space that lead to a particular attractor. We can find the basins of attraction using the function basins.

using Attractors
 
 ba = basins(sys, [0.0, 0], [0.0, 1], [1.0, 0], intervals_to_box([-2, -2], [2, 2]), bstep=[0.01, 0.01], ϵ_mapper=0.001, Ttr=100)
 Ur, Vr, atr, M = ba
 heatmap(Ur, Vr, M)
Example block output

The basin boundaries can be quickly extracted using the function basin_boundaries.

bb = basinboundary(ba)
2×401 Matrix{Float64}:
  -2.0  -1.99  -1.98  -1.97  -1.96  …  1.95  1.96  1.97  1.98  1.99  2.0
   0.0   0.0    0.0    0.0    0.0      0.0   0.0   0.0   0.0   0.0   0.0

Transitions

We can quickly find a path which computes a transition from one attractor to another using the function `transition.

paths_ends = (fp[stab][1], fp[stab][2])
-path, time, succes = transition(sys, paths_ends...);
(2-dimensional StateSpaceSet{Float64} with 641 points, [517.7171629237718, 517.7376002302709, 517.7494329181949, 517.7599766401947, 517.7703452646747, 517.7792676865303, 517.788877834018, 517.7979336902946, 517.8064944742007, 517.8149114538829  …  523.8363302652367, 523.8493929443385, 523.8625686444242, 523.8765884429681, 523.8920703025415, 523.908016604301, 523.9248026164025, 523.9409609274759, 523.957922755368, 523.957922755368], true)
fig = Figure(size=(600, 400), fontsize=13)
+path, time, succes = transition(sys, paths_ends...);
(2-dimensional StateSpaceSet{Float64} with 900 points, [106.48141489813125, 106.49780026691688, 106.5119462440886, 106.52693095484169, 106.53908156437723, 106.55011244791457, 106.56099909907199, 106.57074227742457, 106.58044624767932, 106.59005455391035  …  115.17352926435885, 115.18420243831719, 115.19534784999497, 115.20647170494696, 115.21768471454715, 115.22925062199255, 115.24101942028322, 115.25378758660449, 115.2681517737159, 115.2681517737159], true)
fig = Figure(size=(600, 400), fontsize=13)
 ax = Axis(fig[1, 1], xlabel="u", ylabel="v", aspect=1.4,
     xgridcolor=:transparent, ygridcolor=:transparent,
     ylabelrotation=0)
@@ -125,11 +125,11 @@
 fig
 
 lines!(ax, path, color=:black)
-fig
Example block output

If we want to compute many: transitions is the function to use.

tt = transitions(sys, paths_ends..., 3; tmax=1e3);
Transition path ensemble of 3 samples
+fig
Example block output

If we want to compute many: transitions is the function to use.

tt = transitions(sys, paths_ends..., 3; tmax=1e3);
Transition path ensemble of 3 samples
            - sampling success rate:      1.0
-           - mean residence time:        470.887
-           - mean transition time:       8.359
-           - normalized transition rate: 56.3
fig = Figure(size=(600, 400), fontsize=13)
+           - mean residence time:        63.288
+           - mean transition time:       10.049
+           - normalized transition rate: 6.3
fig = Figure(size=(600, 400), fontsize=13)
 ax = Axis(fig[1, 1], xlabel="u", ylabel="v", aspect=1.4,
     xgridcolor=:transparent, ygridcolor=:transparent,
     ylabelrotation=0)
@@ -150,7 +150,7 @@
 for i in 1:length(tt.paths)
     lines!(ax, tt.paths[i])
 end
-fig
Example block output

Large deviation theory

In the context of nonlinear dynamics, Large Deviation Theory provides tools to quantify the probability of rare events that deviate significantly from the system's typical behavior. These rare events might be extreme values of a system's output, sudden transitions between different states, or other phenomena that occur with very low probability but can have significant implications for the system's overall behavior.

Large deviation theory applies principles from probability theory and statistical mechanics to develop a rigorous mathematical description of these rare events. It uses the concept of a rate function, which measures the exponential decay rate of the probability of large deviations from the mean or typical behavior. This rate function plays a crucial role in quantifying the likelihood of rare events and understanding their impact on the system.

For example, in a system exhibiting chaotic behavior, LDT can help quantify the probability of sudden large shifts in the system's trajectory. Similarly, in a system with multiple stable states, it can provide insight into the likelihood and pathways of transitions between these states under fluctuations. In the context of the Minimum Action Method (MAM) and the Geometric Minimum Action Method (gMAM), Large Deviation Theory is used to handle the large deviations action functional on the space of curves. This is a key part of how these methods analyze dynamical systems.

The Maier-Stein model is a typical benchmark to test such LDT techniques. Let us try to reproduce the following figure from Tobias Grafke's blog post:

maier_stein

Let us first make an initial path:

xx = range(-1.0, 1.0, length=100)
+fig
Example block output

Large deviation theory

In the context of nonlinear dynamics, Large Deviation Theory provides tools to quantify the probability of rare events that deviate significantly from the system's typical behavior. These rare events might be extreme values of a system's output, sudden transitions between different states, or other phenomena that occur with very low probability but can have significant implications for the system's overall behavior.

Large deviation theory applies principles from probability theory and statistical mechanics to develop a rigorous mathematical description of these rare events. It uses the concept of a rate function, which measures the exponential decay rate of the probability of large deviations from the mean or typical behavior. This rate function plays a crucial role in quantifying the likelihood of rare events and understanding their impact on the system.

For example, in a system exhibiting chaotic behavior, LDT can help quantify the probability of sudden large shifts in the system's trajectory. Similarly, in a system with multiple stable states, it can provide insight into the likelihood and pathways of transitions between these states under fluctuations. In the context of the Minimum Action Method (MAM) and the Geometric Minimum Action Method (gMAM), Large Deviation Theory is used to handle the large deviations action functional on the space of curves. This is a key part of how these methods analyze dynamical systems.

The Maier-Stein model is a typical benchmark to test such LDT techniques. Let us try to reproduce the following figure from Tobias Grafke's blog post:

maier_stein

Let us first make an initial path:

xx = range(-1.0, 1.0, length=100)
 yy = 0.3 .* (-xx .^ 2 .+ 1)
 init = Matrix([xx yy]')
2×100 Matrix{Float64}:
  -1.0  -0.979798   -0.959596   -0.939394   …  0.959596   0.979798   1.0
@@ -195,4 +195,4 @@
 
 lines!(ax, init, linewidth=3, color=:black, linestyle=:dash)
 lines!(ax, MLP, linewidth=3, color=:orange)
-fig
Example block output

Author: Orjan Ameye (orjan.ameye@hotmail.com) Date: 13 Feb 2024

+figExample block output

Author: Orjan Ameye (orjan.ameye@hotmail.com) Date: 13 Feb 2024

diff --git a/dev/examples/sgMAM_KPO/index.html b/dev/examples/sgMAM_KPO/index.html index 0c1118f6..63e10c93 100644 --- a/dev/examples/sgMAM_KPO/index.html +++ b/dev/examples/sgMAM_KPO/index.html @@ -35,7 +35,7 @@ return Matrix([H_pu H_pv]') end -sys = SgmamSystem(H_x, H_p)
SgmamSystem(Main.H_x, Main.H_p)

We saved this function in the SgmamSystem struct. We want to find the optimal path between two attractors in the phase space. We define the initial trajectory as wiggle between the two attractors.

# setup
+sys = SgmamSystem{false, 2}(H_x, H_p)
Doubled 2-dimensional phase space containing out-of-place functions

We saved this function in the SgmamSystem struct. We want to find the optimal path between two attractors in the phase space. We define the initial trajectory as wiggle between the two attractors.

# setup
 Nt = 500  # number of discrete time steps
 s = collect(range(0; stop=1, length=Nt))
 
@@ -78,4 +78,4 @@
     colormap=[:gray, :gray]
 )
 axislegend(ax)
-fig
Example block output +figExample block output diff --git a/dev/examples/tutorial/1fb3990e.svg b/dev/examples/tutorial/b1a44eca.svg similarity index 88% rename from dev/examples/tutorial/1fb3990e.svg rename to dev/examples/tutorial/b1a44eca.svg index 487687cd..b85287b2 100644 --- a/dev/examples/tutorial/1fb3990e.svg +++ b/dev/examples/tutorial/b1a44eca.svg @@ -1,45 +1,45 @@ - + - + - + - + - + - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/examples/tutorial/index.html b/dev/examples/tutorial/index.html index 7e47add9..36c92b93 100644 --- a/dev/examples/tutorial/index.html +++ b/dev/examples/tutorial/index.html @@ -37,4 +37,4 @@ scatter!([fp1[1], fp2[1]], [fp1[2], fp2[2]], color=:red, markersize=4) xlims!(-1.2, 1.2) ylims!(-0.6, 0.6) -pltExample block output

Hopefully, this helped you to get started. For more info, check out the Manual section of these docs.

+pltExample block output

Hopefully, this helped you to get started. For more info, check out the Manual section of these docs.

diff --git a/dev/index.html b/dev/index.html index fccd7e5b..de84c65b 100644 --- a/dev/index.html +++ b/dev/index.html @@ -1,2 +1,2 @@ -Home · CriticalTransitions.jl

CriticalTransitions.jl

A Julia package for the numerical investigation of noise- and rate-induced transitions in dynamical systems.

Building on DynamicalSystems.jl and DifferentialEquations.jl, this package aims to provide a toolbox for dynamical systems under time-dependent forcing, with a focus on tipping phenomena and metastability.

CT.jl infographic

Current features
  • Stochastic simulation made easy: Gaussian noise, uncorrelated and correlated, additive and multiplicative
  • Transition path sampling: Parallelized ensemble rejection sampling
  • Large deviation theory tools: Action functionals and minimization algorithms (MAM, gMAM)
Planned features
  • Rare event simulation: importance sampling, AMS
  • Quasipotentials: Ordered line integral method (OLIM)
  • Rate-induced tipping tools
  • Symbolic differentiation of action functionals
  • ...?

Developers: Reyk Börner, Ryan Deeley, Raphael Römer and Orjan Ameye

Thanks to Jeroen Wouters, Calvin Nesbitt, Tobias Grafke, George Datseris and Oliver Mehling

This work is part of the CriticalEarth project.

+Home · CriticalTransitions.jl

CriticalTransitions.jl

A Julia package for the numerical investigation of noise- and rate-induced transitions in dynamical systems.

Building on DynamicalSystems.jl and DifferentialEquations.jl, this package aims to provide a toolbox for dynamical systems under time-dependent forcing, with a focus on tipping phenomena and metastability.

CT.jl infographic

Current features
  • Stochastic simulation made easy: Gaussian noise, uncorrelated and correlated, additive and multiplicative
  • Transition path sampling: Parallelized ensemble rejection sampling
  • Large deviation theory tools: Action functionals and minimization algorithms (MAM, gMAM)
Planned features
  • Rare event simulation: importance sampling, AMS
  • Quasipotentials: Ordered line integral method (OLIM)
  • Rate-induced tipping tools
  • Symbolic differentiation of action functionals
  • ...?

Developers: Reyk Börner, Ryan Deeley, Raphael Römer and Orjan Ameye

Thanks to Jeroen Wouters, Calvin Nesbitt, Tobias Grafke, George Datseris and Oliver Mehling

This work is part of the CriticalEarth project.

diff --git a/dev/man/CoupledSDEs/index.html b/dev/man/CoupledSDEs/index.html index f025c34b..608126b8 100644 --- a/dev/man/CoupledSDEs/index.html +++ b/dev/man/CoupledSDEs/index.html @@ -120,4 +120,4 @@ p = [1.,3.,1.,1.,1.,0.] sys = CoupledSDEs(fitzhugh_nagumo, zeros(2), p; noise_strength=0.1) -ls = lyapunovspectrum(CoupledODEs(sys), 10000) +ls = lyapunovspectrum(CoupledODEs(sys), 10000) diff --git a/dev/man/largedeviations/index.html b/dev/man/largedeviations/index.html index 8e7748ec..4be57e98 100644 --- a/dev/man/largedeviations/index.html +++ b/dev/man/largedeviations/index.html @@ -7,15 +7,15 @@ λ, generalized_momentum, path_velocity -)

defined at /home/runner/work/CriticalTransitions.jl/CriticalTransitions.jl/src/largedeviations/MinimumActionPath.jl:29.

source

String method

The string method is a technique for finding transition paths between two states in a dynamical system. The method represents the path as a "string" of points that connects the states and evolves it to minimize the drift along the path. The resulating tangent path is parrallel to the drift of the system, i.e., the string method computes the heteroclinic orbit. For non-gradient systems (detailed -balance is broken), the heteroclinic orbit differs from the transition path, it does correctly predict, it correctly captures the deterministic dynamics from the saddle point onward ("downhill" portion of the path).

CriticalTransitions.string_methodFunction
string_method(
+)

defined at /home/runner/work/CriticalTransitions.jl/CriticalTransitions.jl/src/largedeviations/MinimumActionPath.jl:29.

source

String method

The string method is a technique for finding transition paths between two states in a dynamical system. The method represents the path as a "string" of points that connects the states and evolves it to minimize the drift along the path. The resulating tangent path is parrallel to the drift of the system, i.e., the string method computes the heteroclinic orbit. For non-gradient systems (detailed -balance is broken), the heteroclinic orbit differs from the transition path, it does correctly predict, it correctly captures the deterministic dynamics from the saddle point onward ("downhill" portion of the path).

CriticalTransitions.string_methodFunction
string_method(
     sys::Union{Function, SgmamSystem},
     x_initial::Matrix;
     ϵ,
     iterations,
     show_progress
 ) -> Any
-

Compute the string method for a given system using E et al. (2007).

The string method is an iterative algorithm used to find minimum energy path (MEP) between two points in phase space. It works by evolving a discretized path (string) according to the system's drift while maintaining equal arc-length parametrization between points.

This implementation allows for computation between arbitrary points, not just stable fixed points.

Arguments

  • sys::SgmamSystem: The doubled phase space system for which the string method is computed
  • x_initial: Initial path discretized as a matrix where each column represents a point on the path
  • ϵ::Float64: Step size for the evolution step
  • iterations::Int64: Maximum number of iterations for path convergence
  • show_progress::Bool: Whether to display a progress meter during computation

Returns

  • x: The final converged path representing the MEP
source
string_method(sys::CoupledSDEs, init; kwargs...) -> Any
-

Compute the string method for a given system using E et al. (2007).

The string method is an iterative algorithm used to find minimum energy path (MEP) between two points in phase space. It works by evolving a discretized path (string) according to the system's drift while maintaining equal arc-length parametrization between points.

This implementation allows for computation between arbitrary points, not just stable fixed points.

Arguments

  • sys::CoupledSDEs: The system for which the string method is computed
  • x_initial: Initial path discretized as a matrix where each column represents a point on the path
  • ϵ::Float64: Step size for the evolution step
  • iterations::Int64: Maximum number of iterations for path convergence
  • show_progress::Bool: Whether to display a progress meter during computation

Returns

  • x: The final converged path representing the MEP
source

Minimum action methods

The minimum action method (MAM) is a technique for calculating the most probable transition path between two (meta)stable states in a stochastic dynamical system. In the limit of small noise, this path corresponds to the minimizer of an action functional. The action functional typically takes into account both the deterministic drift and the noise intensity of the system. By discretizing this path and using optimization techniques, MAM finds the trajectory that requires the least "effort" to transition between states in phase space.

Minimum action method (MAM)

Minimization of the action functioal using the optimization algorithm of Optimization.jl.

CriticalTransitions.min_action_methodFunction
min_action_method(
+

Compute the string method for a given system using E et al. (2007).

The string method is an iterative algorithm used to find minimum energy path (MEP) between two points in phase space. It works by evolving a discretized path (string) according to the system's drift while maintaining equal arc-length parametrization between points.

This implementation allows for computation between arbitrary points, not just stable fixed points.

Arguments

  • sys::SgmamSystem: The doubled phase space system for which the string method is computed
  • x_initial: Initial path discretized as a matrix where each column represents a point on the path
  • ϵ::Float64: Step size for the evolution step
  • iterations::Int64: Maximum number of iterations for path convergence
  • show_progress::Bool: Whether to display a progress meter during computation

Returns

  • x: The final converged path representing the MEP
source
string_method(sys::CoupledSDEs, init; kwargs...) -> Any
+

Compute the string method for a given system using E et al. (2007).

The string method is an iterative algorithm used to find minimum energy path (MEP) between two points in phase space. It works by evolving a discretized path (string) according to the system's drift while maintaining equal arc-length parametrization between points.

This implementation allows for computation between arbitrary points, not just stable fixed points.

Arguments

  • sys::CoupledSDEs: The system for which the string method is computed
  • x_initial: Initial path discretized as a matrix where each column represents a point on the path
  • ϵ::Float64: Step size for the evolution step
  • iterations::Int64: Maximum number of iterations for path convergence
  • show_progress::Bool: Whether to display a progress meter during computation

Returns

  • x: The final converged path representing the MEP
source

Minimum action methods

The minimum action method (MAM) is a technique for calculating the most probable transition path between two (meta)stable states in a stochastic dynamical system. In the limit of small noise, this path corresponds to the minimizer of an action functional. The action functional typically takes into account both the deterministic drift and the noise intensity of the system. By discretizing this path and using optimization techniques, MAM finds the trajectory that requires the least "effort" to transition between states in phase space.

Minimum action method (MAM)

Minimization of the action functioal using the optimization algorithm of Optimization.jl.

CriticalTransitions.min_action_methodFunction
min_action_method(
     sys::CoupledSDEs,
     x_i,
     x_f,
@@ -23,7 +23,7 @@
     T::Real;
     kwargs...
 ) -> MinimumActionPath{_A, _B, _C, Nothing, Nothing, Nothing, Nothing, Nothing} where {_A, _B<:Real, _C}
-

Runs the Minimum Action Method (MAM) to find the minimum action path (instanton) between an initial state x_i and final state x_f in phase space.

This algorithm uses the minimizers of the Optimization.jl package to minimize the action functional (see fw_action) of a path for the given CoupledSDEs sys. The path is initialized as a straight line between x_i and x_f, parameterized in time via N equidistant points and total time T. Thus, the time step between discretized path points is $\Delta t = T/N$. To set an initial path different from a straight line, see the multiple dispatch method

  • min_action_method(sys::CoupledSDEs, init::Matrix, T::Real; kwargs...).

The minimization can be performed in blocks to save intermediate results.

Keyword arguments

  • functional = "FW": type of action functional to minimize. Defaults to fw_action, alternative: om_action.
  • maxiter = 100: maximum number of iterations before the algorithm stops.
  • abstol=1e-8: absolute tolerance of action gradient to determine convergence
  • reltol=1e-8: relative tolerance of action gradient to determine convergence
  • method = Adam(): minimization algorithm (see Optimization.jl)
  • verbose = true: whether to print Optimization information during the run
  • show_progress = false: whether to print a progress bar

Output

If save_iterations, returns Optim.OptimizationResults. Else, returns only the optimizer (path). If blocks > 1, a vector of results/optimizers is returned.

source
min_action_method(
+

Runs the Minimum Action Method (MAM) to find the minimum action path (instanton) between an initial state x_i and final state x_f in phase space.

This algorithm uses the minimizers of the Optimization.jl package to minimize the action functional (see fw_action) of a path for the given CoupledSDEs sys. The path is initialized as a straight line between x_i and x_f, parameterized in time via N equidistant points and total time T. Thus, the time step between discretized path points is $\Delta t = T/N$. To set an initial path different from a straight line, see the multiple dispatch method

  • min_action_method(sys::CoupledSDEs, init::Matrix, T::Real; kwargs...).

The minimization can be performed in blocks to save intermediate results.

Keyword arguments

  • functional = "FW": type of action functional to minimize. Defaults to fw_action, alternative: om_action.
  • maxiter = 100: maximum number of iterations before the algorithm stops.
  • abstol=1e-8: absolute tolerance of action gradient to determine convergence
  • reltol=1e-8: relative tolerance of action gradient to determine convergence
  • method = Adam(): minimization algorithm (see Optimization.jl)
  • verbose = true: whether to print Optimization information during the run
  • show_progress = false: whether to print a progress bar

Output

If save_iterations, returns Optim.OptimizationResults. Else, returns only the optimizer (path). If blocks > 1, a vector of results/optimizers is returned.

source
min_action_method(
     sys::CoupledSDEs,
     init::Matrix,
     T::Real;
@@ -36,14 +36,14 @@
     verbose,
     show_progress
 ) -> MinimumActionPath{_A, _B, _C, Nothing, Nothing, Nothing, Nothing, Nothing} where {_A, _B<:Real, _C}
-

Runs the Minimum Action Method (MAM) to find the minimum action path (instanton) from an initial condition init, given a system sys and total path time T.

The initial path init must be a matrix of size (D, N), where D is the dimension of the system and N is the number of path points. The physical time of the path is specified by T, such that the time step between consecutive path points is $\Delta t = T/N$.

For more information see the main method, min_action_method(sys::CoupledSDEs, x_i, x_f, N::Int, T::Real; kwargs...).

source

Geometric minimum action method (gMAM)

Minimization of the geometric action following Heymann and Vanden-Eijnden, PRL (2008). The gMAM reformulates MAM to avoid double optimisation of both the action and the transition time. It achieves this by using a geometric action functional that is independent of the time parametrization of the path. This reparametrization invariance makes the method more robust and computationally efficient, particularly for systems with metastable states separated by large barriers.

CriticalTransitions.geometric_min_action_methodFunction
geometric_min_action_method(
+

Runs the Minimum Action Method (MAM) to find the minimum action path (instanton) from an initial condition init, given a system sys and total path time T.

The initial path init must be a matrix of size (D, N), where D is the dimension of the system and N is the number of path points. The physical time of the path is specified by T, such that the time step between consecutive path points is $\Delta t = T/N$.

For more information see the main method, min_action_method(sys::CoupledSDEs, x_i, x_f, N::Int, T::Real; kwargs...).

source

Geometric minimum action method (gMAM)

Minimization of the geometric action following Heymann and Vanden-Eijnden, PRL (2008). The gMAM reformulates MAM to avoid double optimisation of both the action and the transition time. It achieves this by using a geometric action functional that is independent of the time parametrization of the path. This reparametrization invariance makes the method more robust and computationally efficient, particularly for systems with metastable states separated by large barriers.

CriticalTransitions.geometric_min_action_methodFunction
geometric_min_action_method(
     sys::CoupledSDEs,
     x_i,
     x_f;
     N,
     kwargs...
 ) -> MinimumActionPath{_A, _B, _C, Nothing, Nothing, Nothing, Nothing, Nothing} where {_A, _B<:Real, _C}
-

Computes the minimizer of the geometric Freidlin-Wentzell action based on the geometric minimum action method (gMAM), using optimizers of OPtimization.jl or the original formulation by Heymann and Vanden-Eijnden[1]. Only Freidlin-Wentzell action has a geometric formulation.

To set an initial path different from a straight line, see the multiple dispatch method

  • geometric_min_action_method(sys::CoupledSDEs, init::Matrix, arclength::Float64; kwargs...).

Keyword arguments

  • maxiter::Int=100: maximum number of optimization iterations before the alogrithm stops
  • abstol=1e-8: absolute tolerance of action gradient to determine convergence
  • reltol=1e-8: relative tolerance of action gradient to determine convergence
  • method = Adam(): minimization algorithm (see Optimization.jl)
  • =0.1: step size parameter in gradient descent HeymannVandenEijnden implementation.
  • verbose=false: if true, print additional output
  • show_progress=true: if true, display a progress bar

Optimization algorithms

The method keyword argument takes solver methods of the Optimization.jl package; alternatively, the option solver = "HeymannVandenEijnden" uses the original gMAM algorithm[1].

source
geometric_min_action_method(
+

Computes the minimizer of the geometric Freidlin-Wentzell action based on the geometric minimum action method (gMAM), using optimizers of OPtimization.jl or the original formulation by Heymann and Vanden-Eijnden[1]. Only Freidlin-Wentzell action has a geometric formulation.

To set an initial path different from a straight line, see the multiple dispatch method

  • geometric_min_action_method(sys::CoupledSDEs, init::Matrix, arclength::Float64; kwargs...).

Keyword arguments

  • maxiter::Int=100: maximum number of optimization iterations before the alogrithm stops
  • abstol=1e-8: absolute tolerance of action gradient to determine convergence
  • reltol=1e-8: relative tolerance of action gradient to determine convergence
  • method = Adam(): minimization algorithm (see Optimization.jl)
  • =0.1: step size parameter in gradient descent HeymannVandenEijnden implementation.
  • verbose=false: if true, print additional output
  • show_progress=true: if true, display a progress bar

Optimization algorithms

The method keyword argument takes solver methods of the Optimization.jl package; alternatively, the option solver = "HeymannVandenEijnden" uses the original gMAM algorithm[1].

source
geometric_min_action_method(
     sys::CoupledSDEs,
     init::Matrix;
     maxiter,
@@ -55,7 +55,7 @@
     verbose,
     show_progress
 ) -> MinimumActionPath{_A, _B, _C, Nothing, Nothing, Nothing, Nothing, Nothing} where {_A, _B<:Real, _C}
-

Runs the geometric Minimum Action Method (gMAM) to find the minimum action path (instanton) from an initial condition init, given a system sys and total arc length arclength.

The initial path init must be a matrix of size (D, N), where D is the dimension of the system and N is the number of path points.

For more information see the main method, geometric_min_action_method(sys::CoupledSDEs, x_i, x_f, arclength::Float64; kwargs...).

source

Simple Geometric minimum action method (sgMAM)

Simplified minimization of the geometric action following Grafke et al. (2017). The simple gMAM reduces the complexity of the original gMAM by requiring only first-order derivatives of the underlying Hamiltonian optimisation formulation. This simplifies the numerical treatment and the computational complexity.

The implementation below perform a constrained gradient descent where it assumes an autonomous system with additive Gaussian noise.

CriticalTransitions.sgmamFunction
sgmam(
+

Runs the geometric Minimum Action Method (gMAM) to find the minimum action path (instanton) from an initial condition init, given a system sys and total arc length arclength.

The initial path init must be a matrix of size (D, N), where D is the dimension of the system and N is the number of path points.

For more information see the main method, geometric_min_action_method(sys::CoupledSDEs, x_i, x_f, arclength::Float64; kwargs...).

source

Simple Geometric minimum action method (sgMAM)

Simplified minimization of the geometric action following Grafke et al. (2017). The simple gMAM reduces the complexity of the original gMAM by requiring only first-order derivatives of the underlying Hamiltonian optimisation formulation. This simplifies the numerical treatment and the computational complexity.

The implementation below perform a constrained gradient descent where it assumes an autonomous system with additive Gaussian noise.

CriticalTransitions.sgmamFunction
sgmam(
     sys::SgmamSystem,
     x_initial::Matrix{<:Real};
     ϵ,
@@ -63,14 +63,14 @@
     show_progress,
     reltol
 ) -> MinimumActionPath{_A, _B, V, Nothing, Nothing, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}} where {_A, _B<:Real, V<:(AbstractVector)}
-

Performs the simplified geometric Minimal Action Method (sgMAM) on the given system sys. Our implementation is only valid for additive noise.

This method computes the optimal path in the phase space of a Hamiltonian system that minimizes the Freidlin–Wentzell action. The Hamiltonian functions H_x and H_p define the system's dynamics in a doubled phase. The initial state x_initial is evolved iteratively using constrained gradient descent with step size parameter ϵ over a specified number of iterations. The method can display a progress meter and will stop early if the relative tolerance reltol is achieved.

The function returns a tuple containing the final state, the action value, the Lagrange multipliers, the momentum, and the state derivatives. The implementation is based on the work of Grafke et al. (2019).

source
CriticalTransitions.SgmamSystemType

A structure representing a system with Hamiltonian functions Hx and Hp.

This system operates in an extended phase space where the Hamiltonian is assumed to be quadratic in the extended momentum. The phase space coordinates x are doubled to form a 2n-dimensional extended phase space.

source

Action functionals

Freidlin-Wentzell action

CriticalTransitions.fw_actionFunction
fw_action(sys::CoupledSDEs, path, time) -> Any
-

Calculates the Freidlin-Wentzell action of a given path with time points time in a drift field specified by the deterministic dynamics f = dynamic_rule(sys) and (normalized) noise covariance matrix covariance_matrix(sys).

The path must be a (D x N) matrix, where D is the dimensionality of the system sys and N is the number of path points. The time array must have length N.

Returns a single number, which is the value of the action functional

$S_T[\phi_t] = \frac{1}{2} \int_0^T || \dot \phi_t - f(\phi_t) ||^2_Q \text{d}t$

where $\phi_t$ denotes the path in state space, $b$ the drift field, and $T$ the total time of the path. The subscript $Q$ refers to the generalized norm $||a||_Q^2 := \langle a, Q^{-1} b \rangle$ (see anorm). Here $Q$ is the noise covariance matrix normalized by $D/L_1(Q)$, with $L_1$ being the L1 matrix norm.

source

Geometric Freidlin-Wentzell action

CriticalTransitions.geometric_actionFunction
geometric_action(sys::CoupledSDEs, path) -> Any
+

Performs the simplified geometric Minimal Action Method (sgMAM) on the given system sys. Our implementation is only valid for additive noise.

This method computes the optimal path in the phase space of a Hamiltonian system that minimizes the Freidlin–Wentzell action. The Hamiltonian functions H_x and H_p define the system's dynamics in a doubled phase. The initial state x_initial is evolved iteratively using constrained gradient descent with step size parameter ϵ over a specified number of iterations. The method can display a progress meter and will stop early if the relative tolerance reltol is achieved.

The function returns a tuple containing the final state, the action value, the Lagrange multipliers, the momentum, and the state derivatives. The implementation is based on the work of Grafke et al. (2019).

source
CriticalTransitions.SgmamSystemType

A structure representing a system with Hamiltonian functions Hx and Hp.

This system operates in an extended phase space where the Hamiltonian is assumed to be quadratic in the extended momentum. The phase space coordinates x are doubled to form a 2n-dimensional extended phase space.

source

Action functionals

Freidlin-Wentzell action

CriticalTransitions.fw_actionFunction
fw_action(sys::CoupledSDEs, path, time) -> Any
+

Calculates the Freidlin-Wentzell action of a given path with time points time in a drift field specified by the deterministic dynamics f = dynamic_rule(sys) and (normalized) noise covariance matrix covariance_matrix(sys).

The path must be a (D x N) matrix, where D is the dimensionality of the system sys and N is the number of path points. The time array must have length N.

Returns a single number, which is the value of the action functional

$S_T[\phi_t] = \frac{1}{2} \int_0^T || \dot \phi_t - f(\phi_t) ||^2_Q \text{d}t$

where $\phi_t$ denotes the path in state space, $b$ the drift field, and $T$ the total time of the path. The subscript $Q$ refers to the generalized norm $||a||_Q^2 := \langle a, Q^{-1} b \rangle$ (see anorm). Here $Q$ is the noise covariance matrix normalized by $D/L_1(Q)$, with $L_1$ being the L1 matrix norm.

source

Geometric Freidlin-Wentzell action

CriticalTransitions.geometric_actionFunction
geometric_action(sys::CoupledSDEs, path) -> Any
 geometric_action(sys::CoupledSDEs, path, arclength) -> Any
-

Calculates the geometric action of a given path with specified arclength for the drift field specified by the deterministic dynamics f = dynamic_rule(sys) and (normalized) noise covariance matrix covariance_matrix(sys).

For a given path $\varphi$, the geometric action $\bar S$ corresponds to the minimum of the Freidlin-Wentzell action $S_T(\varphi)$ over all travel times $T>0$, where $\varphi$ denotes the path's parameterization in physical time (see fw_action). It is given by the integral

$\bar S[\varphi] = \int_0^L \left( ||\varphi'||_Q \, ||f(\varphi)||_Q - \langle \varphi', \, f(\varphi) \rangle_Q \right) \, \text{d}s$

where $s$ is the arclength coordinate, $L$ the arclength, $f$ the drift field, and the subscript $Q$ refers to the generalized dot product $\langle a, b \rangle_Q := a^{\top} \cdot Q^{-1} b$ (see anorm). Here $Q$ is the noise covariance matrix normalized by $D/L_1(Q)$, with $L_1$ being the L1 matrix norm.

Returns the value of the geometric action $\bar S$.

source

Onsager-Machlup action

CriticalTransitions.om_actionFunction
om_action(sys::CoupledSDEs, path, time, noise_strength)

Calculates the Onsager-Machlup action of a given path with time points time for the drift field f = dynamic_rule(sys) at given noise_strength.

The path must be a (D x N) matrix, where D is the dimensionality of the system sys and N is the number of path points. The time array must have length N.

Returns a single number, which is the value of the action functional

$S^{\sigma}_T[\phi_t] = \frac{1}{2} \int_0^T \left( || \dot \phi_t - f(\phi_t) ||^2_Q + \sigma^2 \nabla \cdot f \right) \, \text{d} t$

where $\phi_t$ denotes the path in state space, $b$ the drift field, $T$ the total time of the path, and $\sigma$ the noise strength. The subscript $Q$ refers to the generalized norm $||a||_Q^2 := \langle a, Q^{-1} b \rangle$ (see anorm). Here $Q$ is the noise covariance matrix normalized by $D/L_1(Q)$, with $L_1$ being the L1 matrix norm.

source

For convenience, a general action function is available where the type of functional is set as an argument:

CriticalTransitions.actionFunction
action(
+

Calculates the geometric action of a given path with specified arclength for the drift field specified by the deterministic dynamics f = dynamic_rule(sys) and (normalized) noise covariance matrix covariance_matrix(sys).

For a given path $\varphi$, the geometric action $\bar S$ corresponds to the minimum of the Freidlin-Wentzell action $S_T(\varphi)$ over all travel times $T>0$, where $\varphi$ denotes the path's parameterization in physical time (see fw_action). It is given by the integral

$\bar S[\varphi] = \int_0^L \left( ||\varphi'||_Q \, ||f(\varphi)||_Q - \langle \varphi', \, f(\varphi) \rangle_Q \right) \, \text{d}s$

where $s$ is the arclength coordinate, $L$ the arclength, $f$ the drift field, and the subscript $Q$ refers to the generalized dot product $\langle a, b \rangle_Q := a^{\top} \cdot Q^{-1} b$ (see anorm). Here $Q$ is the noise covariance matrix normalized by $D/L_1(Q)$, with $L_1$ being the L1 matrix norm.

Returns the value of the geometric action $\bar S$.

source

Onsager-Machlup action

CriticalTransitions.om_actionFunction
om_action(sys::CoupledSDEs, path, time, noise_strength)

Calculates the Onsager-Machlup action of a given path with time points time for the drift field f = dynamic_rule(sys) at given noise_strength.

The path must be a (D x N) matrix, where D is the dimensionality of the system sys and N is the number of path points. The time array must have length N.

Returns a single number, which is the value of the action functional

$S^{\sigma}_T[\phi_t] = \frac{1}{2} \int_0^T \left( || \dot \phi_t - f(\phi_t) ||^2_Q + \sigma^2 \nabla \cdot f \right) \, \text{d} t$

where $\phi_t$ denotes the path in state space, $b$ the drift field, $T$ the total time of the path, and $\sigma$ the noise strength. The subscript $Q$ refers to the generalized norm $||a||_Q^2 := \langle a, Q^{-1} b \rangle$ (see anorm). Here $Q$ is the noise covariance matrix normalized by $D/L_1(Q)$, with $L_1$ being the L1 matrix norm.

source

For convenience, a general action function is available where the type of functional is set as an argument:

CriticalTransitions.actionFunction
action(
     sys::CoupledSDEs,
     path::Matrix,
     time,
     functional;
     noise_strength
 ) -> Any
-

Computes the action functional specified by functional for a given CoupledSDEs sys and path parameterized by time.

  • functional = "FW": Returns the Freidlin-Wentzell action (fw_action)
  • functional = "OM": Returns the Onsager-Machlup action (om_action)
source
+

Computes the action functional specified by functional for a given CoupledSDEs sys and path parameterized by time.

source
diff --git a/dev/man/sampling/index.html b/dev/man/sampling/index.html index cdc6733a..d7fc62ec 100644 --- a/dev/man/sampling/index.html +++ b/dev/man/sampling/index.html @@ -6,7 +6,7 @@ residence_time, transition_time, sciml_ensemble -)

defined at /home/runner/work/CriticalTransitions.jl/CriticalTransitions.jl/src/trajectories/transition.jl:14.

source

... by direct simulation

These functions generate noise-induced transitions between an initial and final state.

CriticalTransitions.transitionFunction
transition(
+)

defined at /home/runner/work/CriticalTransitions.jl/CriticalTransitions.jl/src/trajectories/transition.jl:14.

source

... by direct simulation

These functions generate noise-induced transitions between an initial and final state.

CriticalTransitions.transitionFunction
transition(
     sys::CoupledSDEs,
     x_i,
     x_f;
@@ -16,4 +16,4 @@
     cut_start,
     kwargs...
 ) -> Tuple{Any, Any, Bool}
-

Generates a sample transition from point x_i to point x_f.

This function simulates sys in time, starting from initial condition x_i, until entering a ball of given radius around x_f.

Keyword arguments

  • radii=(0.1, 0.1): radius of the ball around x_i and x_f, respectively
  • tmax=1e3: maximum time until the simulation stops even x_f has not been reached
  • radius_directions=1:length(sys.u): the directions in phase space to consider when calculating the radii rad_i and rad_f. Defaults to all directions. To consider only a subspace of state space, insert a vector of indices of the dimensions to be included.
  • cut_start=true: if false, returns the whole trajectory up to the transition
  • kwargs...: keyword arguments passed to CommonSolve.solve

Output

[path, times, success]

  • path (Matrix): transition path (size [dim × N], where N is the number of time points)
  • times (Vector): time values (since start of simulation) of the path points (size N)
  • success (bool): if true, a transition occured (i.e. the ball around x_f has been reached), else false

See also transitions, trajectory.

source
CriticalTransitions.transitionsFunction
function transitions(sys::CoupledSDEs, x_i, x_f, N=1; kwargs...)

Generates an ensemble of N transition samples of sys from point x_i to point x_f.

This function repeatedly calls the transition function to efficiently generate an ensemble of transitions. Multi-threading is enabled.

Keyword arguments

  • radii=(0.1, 0.1): radius of the ball around x_i and x_f, respectively
  • Nmax: number of attempts before the algorithm stops even if less than N transitions occurred.
  • tmax=1e3: maximum time when the simulation stops even x_f has not been reached
  • radius_directions=1:length(sys.u): the directions in phase space to consider when calculating the radii rad_i and rad_f. Defaults to all directions. To consider only a subspace of state space, insert a vector of indices of the dimensions to be included.
  • cut_start=true: if false, returns the whole trajectory up to the transition
  • show_progress=true: shows a progress bar with respect to Nmax
  • kwargs...: keyword arguments passed to CommonSolve.solve

See also transition.

Returns a TransitionEnsemble object.

source

... in pathspace

Coming soon...

+

Generates a sample transition from point x_i to point x_f.

This function simulates sys in time, starting from initial condition x_i, until entering a ball of given radius around x_f.

Keyword arguments

Output

[path, times, success]

See also transitions, trajectory.

source
CriticalTransitions.transitionsFunction
function transitions(sys::CoupledSDEs, x_i, x_f, N=1; kwargs...)

Generates an ensemble of N transition samples of sys from point x_i to point x_f.

This function repeatedly calls the transition function to efficiently generate an ensemble of transitions. Multi-threading is enabled.

Keyword arguments

  • radii=(0.1, 0.1): radius of the ball around x_i and x_f, respectively
  • Nmax: number of attempts before the algorithm stops even if less than N transitions occurred.
  • tmax=1e3: maximum time when the simulation stops even x_f has not been reached
  • radius_directions=1:length(sys.u): the directions in phase space to consider when calculating the radii rad_i and rad_f. Defaults to all directions. To consider only a subspace of state space, insert a vector of indices of the dimensions to be included.
  • cut_start=true: if false, returns the whole trajectory up to the transition
  • show_progress=true: shows a progress bar with respect to Nmax
  • kwargs...: keyword arguments passed to CommonSolve.solve

See also transition.

Returns a TransitionEnsemble object.

source

... in pathspace

Coming soon...

diff --git a/dev/man/simulation/index.html b/dev/man/simulation/index.html index be95c701..2de02519 100644 --- a/dev/man/simulation/index.html +++ b/dev/man/simulation/index.html @@ -11,4 +11,4 @@ diffeq, kwargs... ) -> Tuple{Any, Any} -

Simulates the deterministic (noise-free) dynamics of CoupledSDEs sys in time for a duration T, starting at initial condition init.

This function is equivalent to calling trajectory on the deterministic part of the CoupledSDEs (with noise_strength=0). It works with the ODE solvers used for CoupledODEs.

Keyword arguments

For more info, see ODEProblem. For stochastic integration, see trajectory.

source +

Simulates the deterministic (noise-free) dynamics of CoupledSDEs sys in time for a duration T, starting at initial condition init.

This function is equivalent to calling trajectory on the deterministic part of the CoupledSDEs (with noise_strength=0). It works with the ODE solvers used for CoupledODEs.

Keyword arguments

For more info, see ODEProblem. For stochastic integration, see trajectory.

source diff --git a/dev/man/systemanalysis/index.html b/dev/man/systemanalysis/index.html index 7cfb27f8..6a83e3e4 100644 --- a/dev/man/systemanalysis/index.html +++ b/dev/man/systemanalysis/index.html @@ -4,7 +4,7 @@ bmin::Vector, bmax::Vector ) -> Tuple{StateSpaceSet, Vector{Vector{ComplexF64}}, Vector{Bool}} -

Returns fixed points, their eigenvalues and stability of the system sys within the state space volume defined by bmin and bmax.

This is a wrapper around the fixedpoints function of DynamicalSystems.jl.

Input

Example: fixedpoints(sys, [-2,-1,0], [2,1,1]) finds the fixed points of the 3D system sys in a cube defined by the intervals [-2,2] × [-1,1] × [0,1].

Output

[fp, eigs, stable]

Additional methods

source

Edge tracking

The edge tracking algorithm is a simple numerical method to find the edge state or (possibly chaotic) saddle on the boundary between two basins of attraction. It is first introduced by Battelino et al. (1988) and further described by Skufca et al. (2006).

Attractors.edgetrackingFunction
edgetracking(
+

Returns fixed points, their eigenvalues and stability of the system sys within the state space volume defined by bmin and bmax.

This is a wrapper around the fixedpoints function of DynamicalSystems.jl.

Input

  • bmin (Vector): lower limits of the state space box to be considered, as a vector of coordinates
  • bmax (Vector): upper limits
  • alternatively box (IntervalBox) can replace bmin and bmax

Example: fixedpoints(sys, [-2,-1,0], [2,1,1]) finds the fixed points of the 3D system sys in a cube defined by the intervals [-2,2] × [-1,1] × [0,1].

Output

[fp, eigs, stable]

  • fp: StateSpaceSet of fixed points
  • eigs: vector of Jacobian eigenvalues of each fixed point
  • stable: vector of booleans indicating the stability of each fixed point (true=stable, false=unstable)

Additional methods

  • fixedpoints(sys::CoupedSDEs, box)
source

Edge tracking

The edge tracking algorithm is a simple numerical method to find the edge state or (possibly chaotic) saddle on the boundary between two basins of attraction. It is first introduced by Battelino et al. (1988) and further described by Skufca et al. (2006).

Attractors.edgetrackingFunction
edgetracking(
     sys::CoupledSDEs,
     u1,
     u2,
@@ -22,7 +22,7 @@
     output_all,
     kwargs...
 )
-

Runs the edge tracking algorithm.

Input

  • sys: dynamical system of type CoupledSDEs
  • u1, u2: initial states; must belong to different basins of attraction
  • attractors: vector of state vectors corresponding to the stable fixed points of sys

Keyword arguments

  • eps1 = 1e-9: tolerance for bisection distance
  • eps2 = 1e-8: tolerance for divergence of trajectories before re-bisecting
  • converge = 1e-5: convergence criterion for M state accuracy (Euclidean distance)
  • dt = 0.01: integration time step
  • tmax = Inf: maximum integration time of parallel trajectories until re-bisection
  • ϵ_mapper = 0.1: distance threshold for AttractorMapper
  • dt_mapper = 0.01: time step for AttractorMapper (keyword argument Δt)
  • solver = Vern9(): ODE solver from DifferentialEquations.jl
  • maxit = 100: maximum number of iterations before algorithm stops
  • verbose = true: print status updates during run
  • output_all = false: if false, returns M state, else returns all points of the track
  • kwargs...: additional keyword arguments of AttractorsViaProximity may be passed

Returns

If output_all, a single state vector corresponding to the found edge state is returned. Else, a triple edge, track1, track2 is returned, where track1 and track2 are the tracks along the edge within the basin of attraction of u1 and u2, respectively; edge is the track along the edge derived from (track1 + track2)/2.

Warning

May behave erroneously when run with solver = SimpleATsit5(), which is the default solver for AttractorsViaProximity. The recommended solver is Vern9().

source
Attractors.bisect_to_edgeFunction
bisect_to_edge(
+

Runs the edge tracking algorithm.

Input

  • sys: dynamical system of type CoupledSDEs
  • u1, u2: initial states; must belong to different basins of attraction
  • attractors: vector of state vectors corresponding to the stable fixed points of sys

Keyword arguments

  • eps1 = 1e-9: tolerance for bisection distance
  • eps2 = 1e-8: tolerance for divergence of trajectories before re-bisecting
  • converge = 1e-5: convergence criterion for M state accuracy (Euclidean distance)
  • dt = 0.01: integration time step
  • tmax = Inf: maximum integration time of parallel trajectories until re-bisection
  • ϵ_mapper = 0.1: distance threshold for AttractorMapper
  • dt_mapper = 0.01: time step for AttractorMapper (keyword argument Δt)
  • solver = Vern9(): ODE solver from DifferentialEquations.jl
  • maxit = 100: maximum number of iterations before algorithm stops
  • verbose = true: print status updates during run
  • output_all = false: if false, returns M state, else returns all points of the track
  • kwargs...: additional keyword arguments of AttractorsViaProximity may be passed

Returns

If output_all, a single state vector corresponding to the found edge state is returned. Else, a triple edge, track1, track2 is returned, where track1 and track2 are the tracks along the edge within the basin of attraction of u1 and u2, respectively; edge is the track along the edge derived from (track1 + track2)/2.

Warning

May behave erroneously when run with solver = SimpleATsit5(), which is the default solver for AttractorsViaProximity. The recommended solver is Vern9().

source
Attractors.bisect_to_edgeFunction
bisect_to_edge(
     sys::CoupledSDEs,
     u1,
     u2,
@@ -36,9 +36,9 @@
     relto,
     kwargs...
 )
-

Bisects to the basin boundary between two initial points u1 and u2. Returns the two final points, one on each side of the basin boundary, that are less than eps1 apart from each other.

Keyword arguments

  • eps1=1e-9: tolerance for final distance between the two states
  • ϵ_mapper=0.1: ϵ parameter for AttractorMapper
  • dt_mapper = 0.01: time step for AttractorMapper (keyword argument Δt)
  • solver=Vern9(): solver for AttractorMapper
  • kwargs...: additional kwargs that can be passed to AttractorMapper
source

Baisin of attraction

CriticalTransitions.basinboundaryFunction
basinboundary(X, Y, h; coords, A, B, C) -> Any
-

Computes the basin boundary for given output X, Y, h of the basins function.

To be further documented.

source
basinboundary(boa) -> Any
-

Computes the basin boundary for given output boa of the basins function.

To be further documented.

source
CriticalTransitions.basboundaryFunction
basboundary(
+

Bisects to the basin boundary between two initial points u1 and u2. Returns the two final points, one on each side of the basin boundary, that are less than eps1 apart from each other.

Keyword arguments

  • eps1=1e-9: tolerance for final distance between the two states
  • ϵ_mapper=0.1: ϵ parameter for AttractorMapper
  • dt_mapper = 0.01: time step for AttractorMapper (keyword argument Δt)
  • solver=Vern9(): solver for AttractorMapper
  • kwargs...: additional kwargs that can be passed to AttractorMapper
source

Baisin of attraction

CriticalTransitions.basinboundaryFunction
basinboundary(X, Y, h; coords, A, B, C) -> Any
+

Computes the basin boundary for given output X, Y, h of the basins function.

To be further documented.

source
basinboundary(boa) -> Any
+

Computes the basin boundary for given output boa of the basins function.

To be further documented.

source
CriticalTransitions.basboundaryFunction
basboundary(
     sys::CoupledSDEs,
     xrange::Vector,
     yrange::Vector,
@@ -50,7 +50,7 @@
     solver,
     maxit
 ) -> Tuple{Any, Any}
-

This function computes the basin boundary.

source
CriticalTransitions.basinsFunction
basins(
+

This function computes the basin boundary.

source
CriticalTransitions.basinsFunction
basins(
     sys::CoupledSDEs,
     A,
     B,
@@ -62,4 +62,4 @@
     show_progress,
     kwargs...
 ) -> Any
-

Computes the basins of attraction of CoupledSDEs sys on a plane spanned by the distinct points A, B, C and limited by the box H. Uses the AttractorsViaProximity function from DynamicalSystems.jl to compute the basins of attraction.

A, B, C are elements of $\mathbb{R}^d$ (where $d$ is the dimension of the sys) and H is a hyperrectangle in $\mathbb{R}^d$.

The plane is given by $P_{U,V} := \{A+u(B-A)+v(C-A)\in\mathbb{R}^d: u \in U,\, v\in V\}$ for some closed and bounded real intervals $U$ and $V$ which are selected such that both

  • $P_{U,\,V} \subseteq H$, and
  • $U\times V\subseteq\mathbb{R}^2$ has maximal area,

i.e. $P_{U,\,V}$ is the "largest" possible plane contained within H. This plane is determined behind the scenes.

This function returns a four-dimensional vector. The first two entries are discretised versions of the interval $U$ and $V$ (as defined above, of lengths $\ell_U,\,\ell_V$ respectively); the third entry is a dictionary of the attractors (stable equilibria) of the system within H, and the final entry is an $\ell_V\times\ell_U$ matrix of integers that group the initial conditions (written in terms of $A+u(B-A)+v(C-A)$ where $u\in U$ and $v\in V$) by which attractor they will in time converge to.

Keyword arguments

  • bstep = [0.01, 0.01]: a vector of length two whose elements respectively specify the length of the incremental steps taken across each dimension in the discretisation of your plane
  • pstep = [0.1, 0.1]: a vector of length two whose elements give the increments of the mesh that the maximisation process of finding a plane from a box is taken over (for more information see the source code of the function plane in the src/systemanalysis/planeofbox.jl file)
  • ϵ_mapper = 0.01: ϵ parameter of AttractorsViaProximity
  • kwargs...: keyword arguments passed to the AttractorsViaProximity function (namely, Ttr, Δt, horizon_limit, mx_chk_lost)
  • show_progress = false: if true, displays a progress bar
source
+

Computes the basins of attraction of CoupledSDEs sys on a plane spanned by the distinct points A, B, C and limited by the box H. Uses the AttractorsViaProximity function from DynamicalSystems.jl to compute the basins of attraction.

A, B, C are elements of $\mathbb{R}^d$ (where $d$ is the dimension of the sys) and H is a hyperrectangle in $\mathbb{R}^d$.

The plane is given by $P_{U,V} := \{A+u(B-A)+v(C-A)\in\mathbb{R}^d: u \in U,\, v\in V\}$ for some closed and bounded real intervals $U$ and $V$ which are selected such that both

i.e. $P_{U,\,V}$ is the "largest" possible plane contained within H. This plane is determined behind the scenes.

This function returns a four-dimensional vector. The first two entries are discretised versions of the interval $U$ and $V$ (as defined above, of lengths $\ell_U,\,\ell_V$ respectively); the third entry is a dictionary of the attractors (stable equilibria) of the system within H, and the final entry is an $\ell_V\times\ell_U$ matrix of integers that group the initial conditions (written in terms of $A+u(B-A)+v(C-A)$ where $u\in U$ and $v\in V$) by which attractor they will in time converge to.

Keyword arguments

source diff --git a/dev/man/utils/index.html b/dev/man/utils/index.html index 41108b34..e4b2d5ca 100644 --- a/dev/man/utils/index.html +++ b/dev/man/utils/index.html @@ -1,8 +1,8 @@ Utilities · CriticalTransitions.jl

Utility functions

CoupledSDEs utility functions

StochasticSystemsBase.diffusion_matrixFunction
diffusion_matrix(ds::CoupledSDEs)

Returns the diffusion matrix of the stochastic term of the CoupledSDEs ds, provided that the diffusion function g can be expressed as a constant invertible matrix. If this is not the case, returns nothing.

Note: The diffusion matrix $Σ$ is the square root of the noise covariance matrix $Q$ (see covariance_matrix), defined via the Cholesky decomposition $Q = Σ Σ^\top$.

source
CriticalTransitions.noise_processFunction
noise_process(sys::CoupledSDEs) -> Any
-

Fetches the stochastic process $\mathcal{N}$ specified in the intergrator of sys. Returns the type DiffEqNoiseProcess.NoiseProcess.

source
CriticalTransitions.driftFunction
drift(sys::CoupledSDEs{IIP}, x; t) -> Any
-

Returns the deterministic drift $f(x)$ of the CoupledSDEs sys at state x. For time-dependent systems, the time can be specified as a keyword argument t (by default t=0).

source
StochasticSystemsBase.diffusion_matrixFunction
diffusion_matrix(ds::CoupledSDEs)

Returns the diffusion matrix of the stochastic term of the CoupledSDEs ds, provided that the diffusion function g can be expressed as a constant invertible matrix. If this is not the case, returns nothing.

Note: The diffusion matrix $Σ$ is the square root of the noise covariance matrix $Q$ (see covariance_matrix), defined via the Cholesky decomposition $Q = Σ Σ^\top$.

source
CriticalTransitions.noise_processFunction
noise_process(sys::CoupledSDEs) -> Any
+

Fetches the stochastic process $\mathcal{N}$ specified in the intergrator of sys. Returns the type DiffEqNoiseProcess.NoiseProcess.

source
CriticalTransitions.driftFunction
drift(sys::CoupledSDEs{IIP}, x; t) -> Any
+

Returns the deterministic drift $f(x)$ of the CoupledSDEs sys at state x. For time-dependent systems, the time can be specified as a keyword argument t (by default t=0).

source
CriticalTransitions.div_driftFunction
div_drift(sys::CoupledSDEs, x) -> Any
 div_drift(sys::CoupledSDEs, x, t) -> Any
-

Computes the divergence of the drift field $f(x)$ at state x. For time- dependent systems, the time can be specified as a keyword argument t (by default t=0).

source

More

CriticalTransitions.intervals_to_boxFunction
intervals_to_box(bmin::Vector, bmax::Vector) -> Any
-

Generates a box from specifying the interval limits in each dimension.

  • bmin (Vector): lower limit of the box in each dimension
  • bmax (Vector): upper limit

Example

intervals_to_box([-2,-1,0], [2,1,1]) returns a 3D box of dimensions [-2,2] × [-1,1] × [0,1].

source
+

Computes the divergence of the drift field $f(x)$ at state x. For time- dependent systems, the time can be specified as a keyword argument t (by default t=0).

source

More

CriticalTransitions.intervals_to_boxFunction
intervals_to_box(bmin::Vector, bmax::Vector) -> Any
+

Generates a box from specifying the interval limits in each dimension.

  • bmin (Vector): lower limit of the box in each dimension
  • bmax (Vector): upper limit

Example

intervals_to_box([-2,-1,0], [2,1,1]) returns a 3D box of dimensions [-2,2] × [-1,1] × [0,1].

source
diff --git a/dev/quickstart/index.html b/dev/quickstart/index.html index b4ccead1..cd13e90f 100644 --- a/dev/quickstart/index.html +++ b/dev/quickstart/index.html @@ -1,2 +1,2 @@ -Quickstart · CriticalTransitions.jl

Quickstart

Installation

As this module is not published yet, there are two ways to access it:

  • Option 1 (recommended): Install from GitHub
    1. Enter the Julia package manager by typing ] in the REPL: julia> ]
    2. type add https://github.com/juliadynamics/CriticalTransitions.jl.git
  • Option 2: Load module locally
    1. Clone the repo: git clone https://github.com/juliadynamics/CriticalTransitions.jl.git
    2. In Julia, include the module file: include("PATH/src/CriticalTransitions.jl"), where PATH is the relative path to the repo you just cloned
    3. Load the module: using .CriticalTransitions

Basic usage

The general workflow of CriticalTransitions essentially follows two steps:

  1. Define your system (see Define a CoupledSDEs system)
  2. Investigate the system by calling methods (see Methods)
New system type: RateSystem

We are planning to introduce the the struct RateSystem along CoupledSDEs. In a RateSystem, the time dependence of parameters can conveniently be specified, laying the foundation for a toolbox to study rate-induced tipping, or R-tipping.

Methods

Currently the following functions are implemented to analyze a CoupledSDEs and corresponding sample transition paths.

+Quickstart · CriticalTransitions.jl

Quickstart

Installation

As this module is not published yet, there are two ways to access it:

  • Option 1 (recommended): Install from GitHub
    1. Enter the Julia package manager by typing ] in the REPL: julia> ]
    2. type add https://github.com/juliadynamics/CriticalTransitions.jl.git
  • Option 2: Load module locally
    1. Clone the repo: git clone https://github.com/juliadynamics/CriticalTransitions.jl.git
    2. In Julia, include the module file: include("PATH/src/CriticalTransitions.jl"), where PATH is the relative path to the repo you just cloned
    3. Load the module: using .CriticalTransitions

Basic usage

The general workflow of CriticalTransitions essentially follows two steps:

  1. Define your system (see Define a CoupledSDEs system)
  2. Investigate the system by calling methods (see Methods)
New system type: RateSystem

We are planning to introduce the the struct RateSystem along CoupledSDEs. In a RateSystem, the time dependence of parameters can conveniently be specified, laying the foundation for a toolbox to study rate-induced tipping, or R-tipping.

Methods

Currently the following functions are implemented to analyze a CoupledSDEs and corresponding sample transition paths.

diff --git a/dev/search_index.js b/dev/search_index.js index 05b1e60d..d8718a2f 100644 --- a/dev/search_index.js +++ b/dev/search_index.js @@ -1,3 +1,3 @@ var documenterSearchIndex = {"docs": -[{"location":"examples/tutorial/#Tutorial","page":"Tutorial","title":"Tutorial","text":"","category":"section"},{"location":"examples/tutorial/","page":"Tutorial","title":"Tutorial","text":"To give you an idea of how our package works, this tutorial provides some example code with explanations.","category":"page"},{"location":"examples/tutorial/#Example:-FitzHugh-Nagumo-model","page":"Tutorial","title":"Example: FitzHugh-Nagumo model","text":"","category":"section"},{"location":"examples/tutorial/","page":"Tutorial","title":"Tutorial","text":"Let's consider a simple 2-dimensional dynamical system - the FitzHugh-Nagumo model:","category":"page"},{"location":"examples/tutorial/","page":"Tutorial","title":"Tutorial","text":"beginaligned\nfracdudt = frac1epsilon left( -alpha u^3 + gamma u - kappa v + I right) \nfracdvdt = -beta v + u \nendaligned","category":"page"},{"location":"examples/tutorial/","page":"Tutorial","title":"Tutorial","text":"where epsilon is the parameter of time scale separation between the state variables u and v. The parameters alpha 0, beta 1, gamma0, and kappa0 are real constants, and I denotes a driving term.","category":"page"},{"location":"examples/tutorial/","page":"Tutorial","title":"Tutorial","text":"Let's investigate this system under stochastic forcing.","category":"page"},{"location":"examples/tutorial/#System-definition","page":"Tutorial","title":"System definition","text":"","category":"section"},{"location":"examples/tutorial/","page":"Tutorial","title":"Tutorial","text":"First, we need to translate the system equations above into Julia code.","category":"page"},{"location":"examples/tutorial/","page":"Tutorial","title":"Tutorial","text":"This works exactly as in DynamicalSystems.jl by defining a function f(u,p,t) which takes as input a vector u of state variables (u,v), a vector p of parameters, and time t. The function must return an array of flow increments (textdu, textdv). For performance reasons, it is advisable to return a StaticArray SA[du, dv] rather than just a Vector [du, dv].","category":"page"},{"location":"examples/tutorial/","page":"Tutorial","title":"Tutorial","text":"using CriticalTransitions\nimport Random # hide\nRandom.seed!(1) # hide\n\nfunction fitzhugh_nagumo(u,p,t)\n u, v = u\n ϵ, β, α, γ, κ, I = p\n\n du = (-α*u^3 + γ*u - κ*v + I)/ϵ\n dv = -β*v + u\n\n SA[du, dv]\nend","category":"page"},{"location":"examples/tutorial/","page":"Tutorial","title":"Tutorial","text":"tip: In-place vs. out-of-place\nThe function fitzhugh_nagumo(u,p,t) is defined out-of-place. It is also possible to define the system in-place as fitzhugh_nagumo!(du,u,p,t). For more info, see here.","category":"page"},{"location":"examples/tutorial/#CoupledSDE","page":"Tutorial","title":"CoupledSDE","text":"","category":"section"},{"location":"examples/tutorial/","page":"Tutorial","title":"Tutorial","text":"Next, we construct a stochastic system with the fitzhugh_nagumo equation as the deterministic part. Suppose we would like to force both state variables u and v with additive, uncorrelated Gaussian noise of intensity sigma. This is the default case. We simply write","category":"page"},{"location":"examples/tutorial/","page":"Tutorial","title":"Tutorial","text":"p = [1., 3., 1., 1., 1., 0.] # Parameters (ϵ, β, α, γ, κ, I)\nσ = 0.2 # noise strength\n\n# CoupledSDE\nsys = CoupledSDEs(fitzhugh_nagumo, zeros(2), p; noise_strength=σ)","category":"page"},{"location":"examples/tutorial/","page":"Tutorial","title":"Tutorial","text":"Here the first field fitzhugh_nagumo specifies the deterministic dynamics f (see Define a CoupledSDEs system). We have chosen zeros(2) as the initial state of the system, which is the second field. The length of this vector must match the system's dimensionality. In the (optional) third field, we specify the parameter vector p, which includes the parameters of f followed by the parameters of g (in this case, there are no parameters for g). Lastly, noise_strength sets the noise strength. Since we have not specified a noise process, the default case of an uncorrelated Wiener process is used.","category":"page"},{"location":"examples/tutorial/","page":"Tutorial","title":"Tutorial","text":"note: Multiplicative and/or correlated noise\nOf course, it is also possible to define more complicated noise processes than simple additive white noise. This is done by specifying a custom noise function and covariance matrix in the CoupledSDEs definition. For more info, see Define a CoupledSDEs system.","category":"page"},{"location":"examples/tutorial/","page":"Tutorial","title":"Tutorial","text":"That's it! Now we can apply the toolbox of CriticalTransitions to our stochastic FitzHugh-Nagumo system sys.","category":"page"},{"location":"examples/tutorial/#Find-stable-equilibria","page":"Tutorial","title":"Find stable equilibria","text":"","category":"section"},{"location":"examples/tutorial/","page":"Tutorial","title":"Tutorial","text":"For the parameters chosen above, the FitzHugh-Nagumo system is bistable. Let's compute the fixed points using the fixedpoints function. This function is borrowed from ChaosTools.jl and is loaded as an extension when we write using ChaosTools.","category":"page"},{"location":"examples/tutorial/","page":"Tutorial","title":"Tutorial","text":"using ChaosTools\n# Calculate fixed points and store the stable ones\neqs, eigs, stab = fixedpoints(sys, [-2,-2], [2,2])\nfp1, fp2 = eqs[stab]","category":"page"},{"location":"examples/tutorial/#Stochastic-simulation","page":"Tutorial","title":"Stochastic simulation","text":"","category":"section"},{"location":"examples/tutorial/","page":"Tutorial","title":"Tutorial","text":"Using the trajectory function, we now run a simulation of our system for 100 time units starting out from the fixed point fp1:","category":"page"},{"location":"examples/tutorial/","page":"Tutorial","title":"Tutorial","text":"sim = trajectory(sys, 100, fp1)","category":"page"},{"location":"examples/tutorial/","page":"Tutorial","title":"Tutorial","text":"In the keyword arguments, we have specified at which interval the solution is saved. Further keyword arguments can be used to change the solver (the default is SOSRA() for stochastic integration) and other settings.","category":"page"},{"location":"examples/tutorial/","page":"Tutorial","title":"Tutorial","text":"The simulated trajectory is stored in sim in the usual output format of the solve method of DifferentialEquations.jl, including the solution sim.u and the vector of time points sim.t. The solution can also be accessed as a matrix sim[i, t], where i is the i-th component of u and t the time index.","category":"page"},{"location":"examples/tutorial/","page":"Tutorial","title":"Tutorial","text":"Let's plot the result. Did the trajectory transition to the other attractor?","category":"page"},{"location":"examples/tutorial/","page":"Tutorial","title":"Tutorial","text":"using Plots\nplt = plot(sim[1][:,1], sim[1][:,2]; xlabel=\"u\", ylabel=\"v\", legend=false)\nscatter!([fp1[1], fp2[1]], [fp1[2], fp2[2]], color=:red, markersize=4)\nxlims!(-1.2, 1.2)\nylims!(-0.6, 0.6)\nplt","category":"page"},{"location":"examples/tutorial/","page":"Tutorial","title":"Tutorial","text":"Hopefully, this helped you to get started. For more info, check out the Manual section of these docs.","category":"page"},{"location":"man/simulation/#Simulating-the-system","page":"Simulating the system","title":"Simulating the system","text":"","category":"section"},{"location":"man/simulation/","page":"Simulating the system","title":"Simulating the system","text":"We provide two main functions to simulate a CoupledSDEs forward in time:","category":"page"},{"location":"man/simulation/","page":"Simulating the system","title":"Simulating the system","text":"DynamicalSystemsBase.trajectory, which integrates the stochastic CoupledSDEs system forward in time\ndeterministic_orbit, which integrates only the deterministic part of the CoupledSDEs system ","category":"page"},{"location":"man/simulation/","page":"Simulating the system","title":"Simulating the system","text":"deterministic_orbit","category":"page"},{"location":"man/simulation/#CriticalTransitions.deterministic_orbit","page":"Simulating the system","title":"CriticalTransitions.deterministic_orbit","text":"deterministic_orbit(\n sys::CoupledSDEs,\n T;\n ...\n) -> Tuple{Any, Any}\ndeterministic_orbit(\n sys::CoupledSDEs,\n T,\n init;\n diffeq,\n kwargs...\n) -> Tuple{Any, Any}\n\n\nSimulates the deterministic (noise-free) dynamics of CoupledSDEs sys in time for a duration T, starting at initial condition init.\n\nThis function is equivalent to calling trajectory on the deterministic part of the CoupledSDEs (with noise_strength=0). It works with the ODE solvers used for CoupledODEs.\n\nKeyword arguments\n\ndiffeq=(alg=Tsit5(), abstol = 1e-6, reltol = 1e-6): ODE solver settings (see CoupledODEs)\nkwargs...: keyword arguments passed to trajectory\n\nFor more info, see ODEProblem. For stochastic integration, see trajectory.\n\n\n\n\n\n","category":"function"},{"location":"man/utils/#Utility-functions","page":"Utilities","title":"Utility functions","text":"","category":"section"},{"location":"man/utils/#CoupledSDEs-utility-functions","page":"Utilities","title":"CoupledSDEs utility functions","text":"","category":"section"},{"location":"man/utils/","page":"Utilities","title":"Utilities","text":"solver\nStochasticSystemsBase.covariance_matrix\nStochasticSystemsBase.diffusion_matrix\nnoise_process\ndrift\ndiv_drift","category":"page"},{"location":"man/utils/#CriticalTransitions.solver","page":"Utilities","title":"CriticalTransitions.solver","text":"solver(ds::CoupledSDEs) -> Any\n\n\nReturns the SDE solver specified in the diffeq settings of the CoupledSDEs.\n\n\n\n\n\n","category":"function"},{"location":"man/utils/#StochasticSystemsBase.covariance_matrix","page":"Utilities","title":"StochasticSystemsBase.covariance_matrix","text":"covariance_matrix(ds::CoupledSDEs)\n\nReturns the covariance matrix of the stochastic term of the CoupledSDEs ds, provided that the diffusion function g can be expressed as a constant invertible matrix. If this is not the case, returns nothing.\n\nSee also diffusion_matrix.\n\n\n\n\n\n","category":"function"},{"location":"man/utils/#StochasticSystemsBase.diffusion_matrix","page":"Utilities","title":"StochasticSystemsBase.diffusion_matrix","text":"diffusion_matrix(ds::CoupledSDEs)\n\nReturns the diffusion matrix of the stochastic term of the CoupledSDEs ds, provided that the diffusion function g can be expressed as a constant invertible matrix. If this is not the case, returns nothing.\n\nNote: The diffusion matrix Σ is the square root of the noise covariance matrix Q (see covariance_matrix), defined via the Cholesky decomposition Q = Σ Σ^top.\n\n\n\n\n\n","category":"function"},{"location":"man/utils/#CriticalTransitions.noise_process","page":"Utilities","title":"CriticalTransitions.noise_process","text":"noise_process(sys::CoupledSDEs) -> Any\n\n\nFetches the stochastic process mathcalN specified in the intergrator of sys. Returns the type DiffEqNoiseProcess.NoiseProcess.\n\n\n\n\n\n","category":"function"},{"location":"man/utils/#CriticalTransitions.drift","page":"Utilities","title":"CriticalTransitions.drift","text":"drift(sys::CoupledSDEs{IIP}, x; t) -> Any\n\n\nReturns the deterministic drift f(x) of the CoupledSDEs sys at state x. For time-dependent systems, the time can be specified as a keyword argument t (by default t=0).\n\n\n\n\n\n","category":"function"},{"location":"man/utils/#CriticalTransitions.div_drift","page":"Utilities","title":"CriticalTransitions.div_drift","text":"div_drift(sys::CoupledSDEs, x) -> Any\ndiv_drift(sys::CoupledSDEs, x, t) -> Any\n\n\nComputes the divergence of the drift field f(x) at state x. For time- dependent systems, the time can be specified as a keyword argument t (by default t=0).\n\n\n\n\n\n","category":"function"},{"location":"man/utils/#More","page":"Utilities","title":"More","text":"","category":"section"},{"location":"man/utils/","page":"Utilities","title":"Utilities","text":"intervals_to_box","category":"page"},{"location":"man/utils/#CriticalTransitions.intervals_to_box","page":"Utilities","title":"CriticalTransitions.intervals_to_box","text":"intervals_to_box(bmin::Vector, bmax::Vector) -> Any\n\n\nGenerates a box from specifying the interval limits in each dimension.\n\nbmin (Vector): lower limit of the box in each dimension\nbmax (Vector): upper limit\n\nExample\n\nintervals_to_box([-2,-1,0], [2,1,1]) returns a 3D box of dimensions [-2,2] × [-1,1] × [0,1].\n\n\n\n\n\n","category":"function"},{"location":"examples/gMAM_Maierstein/#The-Maier-Stein-model.","page":"Anlyses for the Maier-Stein system","title":"The Maier-Stein model.","text":"","category":"section"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"using CriticalTransitions\n\nusing CairoMakie\nusing CairoMakie.Makie.MathTeXEngine: get_font\nfont = (;\n regular=get_font(:regular), bold=get_font(:bold),\n italic=get_font(:italic), bold_italic=get_font(:bolditalic)\n);","category":"page"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"Let us explore the features of CriticalTransitions.jl with Maier-Stein model.","category":"page"},{"location":"examples/gMAM_Maierstein/#Maier-stein-model","page":"Anlyses for the Maier-Stein system","title":"Maier-stein model","text":"","category":"section"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"The Maier-Stein model (J. Stat. Phys 83, 3–4 (1996)) is commonly used in the field of nonlinear dynamics for benchmarking Large Deviation Theory (LDT) techniques, e.g., stoachastic transitions between different stable states. It is a simple model that describes the dynamics of a system with two degrees of freedom u and v, and is given by the following set of ordinary differential equations:","category":"page"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"beginaligned\n dotu = u-u^3 - beta*u*v^2\n dotv = -alpha (1+u^2)*v\nendaligned","category":"page"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"The parameter alpha0 controls the strength of the drift field and beta0 represents the softening of that drift field.","category":"page"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"function meier_stein!(du, u, p, t) # in-place\n x, y = u\n du[1] = x - x^3 - 10 * x * y^2\n du[2] = -(1 + x^2) * y\nend\nfunction meier_stein(u, p, t) # out-of-place\n x, y = u\n dx = x - x^3 - 10 * x * y^2\n dy = -(1 + x^2) * y\n SA[dx, dy]\nend\nσ = 0.25\nsys = CoupledSDEs(meier_stein, zeros(2), (); noise_strength=σ)","category":"page"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"A good reference to read about the large deviations methods is this or this blog post by Tobias Grafke.","category":"page"},{"location":"examples/gMAM_Maierstein/#Attractors","page":"Anlyses for the Maier-Stein system","title":"Attractors","text":"","category":"section"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"We start by investigating the deterministic dynamics of the Maier-Stein model.","category":"page"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"The function fixed points return the attractors, their eigenvalues and stability within the state space volume defined by bmin and bmax.","category":"page"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"using ChaosTools\n\nu_min = -1.1;\nu_max = 1.1;\nv_min = -0.4;\nv_max = 0.4;\nbmin = [u_min, v_min];\nbmax = [u_max, v_max];\nfp, eig, stab = fixedpoints(sys, bmin, bmax)\nstable_fp = fp[stab]","category":"page"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"using LinearAlgebra: norm\nres = 100\nu_range = range(u_min, u_max, length=res)\nv_range = range(v_min, v_max, length=res)\n\ndu(u, v) = u - u^3 - 10 * u * v^2\ndv(u, v) = -(1 + u^2) * v\nodeSol(u, v) = Point2f(du(u, v), dv(u, v))\n\nz = [norm([du(x, y), dv(x, y)]) for x in u_range, y in v_range]\nzmin, zmax = minimum(z), maximum(z)\n\nfig = Figure(size=(600, 400), fontsize=13)\nax = Axis(fig[1, 1], xlabel=\"u\", ylabel=\"v\", aspect=1.4,\n xgridcolor=:transparent, ygridcolor=:transparent,\n ylabelrotation=0)\n\nhm = heatmap!(ax, u_range, v_range, z, colormap=:Blues, colorrange=(zmin, zmax))\nColorbar(fig[1, 2], hm; label=\"\", width=15, ticksize=15, tickalign=1)\nstreamplot!(ax, odeSol, (u_min, u_max), (v_min, v_max);\n gridsize=(20, 20), arrow_size=10, stepsize=0.01,\n colormap=[:black, :black]\n)\ncolgap!(fig.layout, 7)\nlimits!(u_min, u_max, v_min, v_max)\nfig\n\n[scatter!(ax, Point(fp[i]), color=stab[i] > 0 ? :red : :dodgerblue,\n markersize=10) for i in eachindex(fp)]\nfig","category":"page"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"We can simulate a stochastic trajectory using the function trajectory.","category":"page"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"tr, ts = trajectory(sys, 1000)\n\nfig = Figure(size=(1000, 400), fontsize=13)\nax1 = Axis(fig[1, 1], xlabel=\"t\", ylabel=\"u\", aspect=1.2,\n xgridcolor=:transparent, ygridcolor=:transparent,\n ylabelrotation=0)\nax2 = Axis(fig[1, 2], xlabel=\"u\", ylabel=\"v\", aspect=1.2,\n xgridcolor=:transparent, ygridcolor=:transparent,\n ylabelrotation=0)\n\nlines!(ax1, ts, first.(tr); linewidth=2, color=:black)\n\nhm = heatmap!(ax2, u_range, v_range, z, colormap=:Blues, colorrange=(zmin, zmax))\nColorbar(fig[1, 3], hm; label=\"\", width=15, ticksize=15, tickalign=1)\nstreamplot!(ax2, odeSol, (u_min, u_max), (v_min, v_max);\n gridsize=(20, 20), arrow_size=10, stepsize=0.01,\n colormap=[:white, :white]\n)\ncolgap!(fig.layout, 7)\nlimits!(u_min, u_max, v_min, v_max)\nfig\n\n[scatter!(ax2, Point(fp[i]), color=stab[i] > 0 ? :red : :dodgerblue,\n markersize=10) for i in eachindex(fp)]\n\nlines!(ax2, reduce(hcat, tr), linewidth=1, color=(:black, 0.2))\nfig","category":"page"},{"location":"examples/gMAM_Maierstein/#Basins-of-attraction","page":"Anlyses for the Maier-Stein system","title":"Basins of attraction","text":"","category":"section"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"Basins of attraction are the regions in the state space that lead to a particular attractor. We can find the basins of attraction using the function basins.","category":"page"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"using Attractors\n\nba = basins(sys, [0.0, 0], [0.0, 1], [1.0, 0], intervals_to_box([-2, -2], [2, 2]), bstep=[0.01, 0.01], ϵ_mapper=0.001, Ttr=100)\nUr, Vr, atr, M = ba\nheatmap(Ur, Vr, M)","category":"page"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"The basin boundaries can be quickly extracted using the function basin_boundaries.","category":"page"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"bb = basinboundary(ba)","category":"page"},{"location":"examples/gMAM_Maierstein/#Transitions","page":"Anlyses for the Maier-Stein system","title":"Transitions","text":"","category":"section"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"We can quickly find a path which computes a transition from one attractor to another using the function `transition.","category":"page"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"paths_ends = (fp[stab][1], fp[stab][2])\npath, time, succes = transition(sys, paths_ends...);","category":"page"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"fig = Figure(size=(600, 400), fontsize=13)\nax = Axis(fig[1, 1], xlabel=\"u\", ylabel=\"v\", aspect=1.4,\n xgridcolor=:transparent, ygridcolor=:transparent,\n ylabelrotation=0)\n\nhm = heatmap!(ax, u_range, v_range, z, colormap=:Blues, colorrange=(zmin, zmax))\nColorbar(fig[1, 2], hm; label=\"\", width=15, ticksize=15, tickalign=1)\nstreamplot!(ax, odeSol, (u_min, u_max), (v_min, v_max);\n gridsize=(20, 20), arrow_size=10, stepsize=0.01,\n colormap=[:white, :white]\n)\ncolgap!(fig.layout, 7)\nlimits!(u_min, u_max, v_min, v_max)\nfig\n\n[scatter!(ax, Point(fp[i]), color=stab[i] > 0 ? :red : :dodgerblue,\n markersize=10) for i in eachindex(fp)]\nfig\n\nlines!(ax, path, color=:black)\nfig","category":"page"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"If we want to compute many: transitions is the function to use.","category":"page"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"tt = transitions(sys, paths_ends..., 3; tmax=1e3);","category":"page"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"fig = Figure(size=(600, 400), fontsize=13)\nax = Axis(fig[1, 1], xlabel=\"u\", ylabel=\"v\", aspect=1.4,\n xgridcolor=:transparent, ygridcolor=:transparent,\n ylabelrotation=0)\n\nhm = heatmap!(ax, u_range, v_range, z, colormap=:Blues, colorrange=(zmin, zmax))\nColorbar(fig[1, 2], hm; label=\"\", width=15, ticksize=15, tickalign=1)\nstreamplot!(ax, odeSol, (u_min, u_max), (v_min, v_max);\n gridsize=(20, 20), arrow_size=10, stepsize=0.01,\n colormap=[:black, :black]\n)\ncolgap!(fig.layout, 7)\nlimits!(u_min, u_max, v_min, v_max)\nfig\n\n[scatter!(ax, Point(fp[i]), color=stab[i] > 0 ? :red : :dodgerblue,\n markersize=10) for i in eachindex(fp)]\n\nfor i in 1:length(tt.paths)\n lines!(ax, tt.paths[i])\nend\nfig","category":"page"},{"location":"examples/gMAM_Maierstein/#Large-deviation-theory","page":"Anlyses for the Maier-Stein system","title":"Large deviation theory","text":"","category":"section"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"In the context of nonlinear dynamics, Large Deviation Theory provides tools to quantify the probability of rare events that deviate significantly from the system's typical behavior. These rare events might be extreme values of a system's output, sudden transitions between different states, or other phenomena that occur with very low probability but can have significant implications for the system's overall behavior.","category":"page"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"Large deviation theory applies principles from probability theory and statistical mechanics to develop a rigorous mathematical description of these rare events. It uses the concept of a rate function, which measures the exponential decay rate of the probability of large deviations from the mean or typical behavior. This rate function plays a crucial role in quantifying the likelihood of rare events and understanding their impact on the system.","category":"page"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"For example, in a system exhibiting chaotic behavior, LDT can help quantify the probability of sudden large shifts in the system's trajectory. Similarly, in a system with multiple stable states, it can provide insight into the likelihood and pathways of transitions between these states under fluctuations. In the context of the Minimum Action Method (MAM) and the Geometric Minimum Action Method (gMAM), Large Deviation Theory is used to handle the large deviations action functional on the space of curves. This is a key part of how these methods analyze dynamical systems.","category":"page"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"The Maier-Stein model is a typical benchmark to test such LDT techniques. Let us try to reproduce the following figure from Tobias Grafke's blog post:","category":"page"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"(Image: maier_stein)","category":"page"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"Let us first make an initial path:","category":"page"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"xx = range(-1.0, 1.0, length=100)\nyy = 0.3 .* (-xx .^ 2 .+ 1)\ninit = Matrix([xx yy]')","category":"page"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"geometric_min_action_method computes the minimizer of the Freidlin-Wentzell action using the geometric minimum action method (gMAM), to find the minimum action path (instanton) between an initial state xi and final state xf. The Minimum Action Method (MAM) is a more traditional approach, while the Geometric Minimum Action Method (gMAM) is a blend of the original MAM and the string method.","category":"page"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"gm = geometric_min_action_method(sys, init, maxiter=500; show_progress=false)\nMLP = gm.path","category":"page"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"fig = Figure(size=(600, 400), fontsize=13)\nax = Axis(fig[1, 1], xlabel=\"u\", ylabel=\"v\", aspect=1.4,\n xgridcolor=:transparent, ygridcolor=:transparent,\n ylabelrotation=0)\n\nhm = heatmap!(ax, u_range, v_range, z, colormap=:Blues, colorrange=(zmin, zmax))\nColorbar(fig[1, 2], hm; label=\"\", width=15, ticksize=15, tickalign=1)\nstreamplot!(ax, odeSol, (u_min, u_max), (v_min, v_max);\n gridsize=(20, 20), arrow_size=10, stepsize=0.01,\n colormap=[:black, :black]\n)\ncolgap!(fig.layout, 7)\nlimits!(u_min, u_max, v_min, v_max)\nfig\n\n[scatter!(ax, Point(fp[i]), color=stab[i] > 0 ? :red : :dodgerblue,\n markersize=10) for i in eachindex(fp)]\n\nlines!(ax, init, linewidth=3, color=:black, linestyle=:dash)\nlines!(ax, MLP, linewidth=3, color=:orange)\nfig","category":"page"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"","category":"page"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"Author: Orjan Ameye (orjan.ameye@hotmail.com) Date: 13 Feb 2024","category":"page"},{"location":"man/largedeviations/#Large-deviation-theory","page":"Large deviation theory","title":"Large deviation theory","text":"","category":"section"},{"location":"man/largedeviations/#Mimumum-action-path","page":"Large deviation theory","title":"Mimumum action path","text":"","category":"section"},{"location":"man/largedeviations/","page":"Large deviation theory","title":"Large deviation theory","text":"MinimumActionPath","category":"page"},{"location":"man/largedeviations/#CriticalTransitions.MinimumActionPath","page":"Large deviation theory","title":"CriticalTransitions.MinimumActionPath","text":"struct MinimumActionPath{D, T<:Real, V, Phis, Ahis, Lambda, PV, GPV}\n\nThe minimum action path between two points in a D-dimensional phase space.\n\nFields\n\npath::StateSpaceSet{D, T, V} where {D, T<:Real, V}: The path matrix.\naction::Real: The action value associated to the path.\npath_history::Any: The history of action of the paths in the optimisation algorithm (optional).\naction_history::Any: The history of action of the paths in the optimisation algorithm (optional).\nλ::Any: The Lagrange multiplier parameter for the minimum action path (optional).\ngeneralized_momentum::Any: The generalized momentum of the phase space variables (optional).\npath_velocity::Any: The path velocity (optional).\n\nConstructors\n\nMinimumActionPath(\n path,\n action;\n path_history,\n action_history,\n λ,\n generalized_momentum,\n path_velocity\n)\n\ndefined at /home/runner/work/CriticalTransitions.jl/CriticalTransitions.jl/src/largedeviations/MinimumActionPath.jl:29.\n\n\n\n\n\n","category":"type"},{"location":"man/largedeviations/#String-method","page":"Large deviation theory","title":"String method","text":"","category":"section"},{"location":"man/largedeviations/","page":"Large deviation theory","title":"Large deviation theory","text":"The string method is a technique for finding transition paths between two states in a dynamical system. The method represents the path as a \"string\" of points that connects the states and evolves it to minimize the drift along the path. The resulating tangent path is parrallel to the drift of the system, i.e., the string method computes the heteroclinic orbit. For non-gradient systems (detailed -balance is broken), the heteroclinic orbit differs from the transition path, it does correctly predict, it correctly captures the deterministic dynamics from the saddle point onward (\"downhill\" portion of the path).","category":"page"},{"location":"man/largedeviations/","page":"Large deviation theory","title":"Large deviation theory","text":"string_method","category":"page"},{"location":"man/largedeviations/#CriticalTransitions.string_method","page":"Large deviation theory","title":"CriticalTransitions.string_method","text":"string_method(\n sys::Union{Function, SgmamSystem},\n x_initial::Matrix;\n ϵ,\n iterations,\n show_progress\n) -> Any\n\n\nCompute the string method for a given system using E et al. (2007).\n\nThe string method is an iterative algorithm used to find minimum energy path (MEP) between two points in phase space. It works by evolving a discretized path (string) according to the system's drift while maintaining equal arc-length parametrization between points.\n\nThis implementation allows for computation between arbitrary points, not just stable fixed points.\n\nArguments\n\nsys::SgmamSystem: The doubled phase space system for which the string method is computed\nx_initial: Initial path discretized as a matrix where each column represents a point on the path\nϵ::Float64: Step size for the evolution step\niterations::Int64: Maximum number of iterations for path convergence\nshow_progress::Bool: Whether to display a progress meter during computation\n\nReturns\n\nx: The final converged path representing the MEP\n\n\n\n\n\nstring_method(sys::CoupledSDEs, init; kwargs...) -> Any\n\n\nCompute the string method for a given system using E et al. (2007).\n\nThe string method is an iterative algorithm used to find minimum energy path (MEP) between two points in phase space. It works by evolving a discretized path (string) according to the system's drift while maintaining equal arc-length parametrization between points.\n\nThis implementation allows for computation between arbitrary points, not just stable fixed points.\n\nArguments\n\nsys::CoupledSDEs: The system for which the string method is computed\nx_initial: Initial path discretized as a matrix where each column represents a point on the path\nϵ::Float64: Step size for the evolution step\niterations::Int64: Maximum number of iterations for path convergence\nshow_progress::Bool: Whether to display a progress meter during computation\n\nReturns\n\nx: The final converged path representing the MEP\n\n\n\n\n\n","category":"function"},{"location":"man/largedeviations/#Minimum-action-methods","page":"Large deviation theory","title":"Minimum action methods","text":"","category":"section"},{"location":"man/largedeviations/","page":"Large deviation theory","title":"Large deviation theory","text":"The minimum action method (MAM) is a technique for calculating the most probable transition path between two (meta)stable states in a stochastic dynamical system. In the limit of small noise, this path corresponds to the minimizer of an action functional. The action functional typically takes into account both the deterministic drift and the noise intensity of the system. By discretizing this path and using optimization techniques, MAM finds the trajectory that requires the least \"effort\" to transition between states in phase space.","category":"page"},{"location":"man/largedeviations/#Minimum-action-method-(MAM)","page":"Large deviation theory","title":"Minimum action method (MAM)","text":"","category":"section"},{"location":"man/largedeviations/","page":"Large deviation theory","title":"Large deviation theory","text":"Minimization of the action functioal using the optimization algorithm of Optimization.jl.","category":"page"},{"location":"man/largedeviations/","page":"Large deviation theory","title":"Large deviation theory","text":"min_action_method","category":"page"},{"location":"man/largedeviations/#CriticalTransitions.min_action_method","page":"Large deviation theory","title":"CriticalTransitions.min_action_method","text":"min_action_method(\n sys::CoupledSDEs,\n x_i,\n x_f,\n N::Int64,\n T::Real;\n kwargs...\n) -> MinimumActionPath{_A, _B, _C, Nothing, Nothing, Nothing, Nothing, Nothing} where {_A, _B<:Real, _C}\n\n\nRuns the Minimum Action Method (MAM) to find the minimum action path (instanton) between an initial state x_i and final state x_f in phase space.\n\nThis algorithm uses the minimizers of the Optimization.jl package to minimize the action functional (see fw_action) of a path for the given CoupledSDEs sys. The path is initialized as a straight line between x_i and x_f, parameterized in time via N equidistant points and total time T. Thus, the time step between discretized path points is Delta t = TN. To set an initial path different from a straight line, see the multiple dispatch method\n\nmin_action_method(sys::CoupledSDEs, init::Matrix, T::Real; kwargs...).\n\nThe minimization can be performed in blocks to save intermediate results.\n\nKeyword arguments\n\nfunctional = \"FW\": type of action functional to minimize. Defaults to fw_action, alternative: om_action.\nmaxiter = 100: maximum number of iterations before the algorithm stops.\nabstol=1e-8: absolute tolerance of action gradient to determine convergence\nreltol=1e-8: relative tolerance of action gradient to determine convergence\nmethod = Adam(): minimization algorithm (see Optimization.jl)\nverbose = true: whether to print Optimization information during the run\nshow_progress = false: whether to print a progress bar\n\nOutput\n\nIf save_iterations, returns Optim.OptimizationResults. Else, returns only the optimizer (path). If blocks > 1, a vector of results/optimizers is returned.\n\n\n\n\n\nmin_action_method(\n sys::CoupledSDEs,\n init::Matrix,\n T::Real;\n functional,\n maxiter,\n abstol,\n reltol,\n method,\n AD,\n verbose,\n show_progress\n) -> MinimumActionPath{_A, _B, _C, Nothing, Nothing, Nothing, Nothing, Nothing} where {_A, _B<:Real, _C}\n\n\nRuns the Minimum Action Method (MAM) to find the minimum action path (instanton) from an initial condition init, given a system sys and total path time T.\n\nThe initial path init must be a matrix of size (D, N), where D is the dimension of the system and N is the number of path points. The physical time of the path is specified by T, such that the time step between consecutive path points is Delta t = TN.\n\nFor more information see the main method, min_action_method(sys::CoupledSDEs, x_i, x_f, N::Int, T::Real; kwargs...).\n\n\n\n\n\n","category":"function"},{"location":"man/largedeviations/#Geometric-minimum-action-method-(gMAM)","page":"Large deviation theory","title":"Geometric minimum action method (gMAM)","text":"","category":"section"},{"location":"man/largedeviations/","page":"Large deviation theory","title":"Large deviation theory","text":"Minimization of the geometric action following Heymann and Vanden-Eijnden, PRL (2008). The gMAM reformulates MAM to avoid double optimisation of both the action and the transition time. It achieves this by using a geometric action functional that is independent of the time parametrization of the path. This reparametrization invariance makes the method more robust and computationally efficient, particularly for systems with metastable states separated by large barriers.","category":"page"},{"location":"man/largedeviations/","page":"Large deviation theory","title":"Large deviation theory","text":"geometric_min_action_method","category":"page"},{"location":"man/largedeviations/#CriticalTransitions.geometric_min_action_method","page":"Large deviation theory","title":"CriticalTransitions.geometric_min_action_method","text":"geometric_min_action_method(\n sys::CoupledSDEs,\n x_i,\n x_f;\n N,\n kwargs...\n) -> MinimumActionPath{_A, _B, _C, Nothing, Nothing, Nothing, Nothing, Nothing} where {_A, _B<:Real, _C}\n\n\nComputes the minimizer of the geometric Freidlin-Wentzell action based on the geometric minimum action method (gMAM), using optimizers of OPtimization.jl or the original formulation by Heymann and Vanden-Eijnden[1]. Only Freidlin-Wentzell action has a geometric formulation.\n\nTo set an initial path different from a straight line, see the multiple dispatch method\n\ngeometric_min_action_method(sys::CoupledSDEs, init::Matrix, arclength::Float64; kwargs...).\n\nKeyword arguments\n\nmaxiter::Int=100: maximum number of optimization iterations before the alogrithm stops\nabstol=1e-8: absolute tolerance of action gradient to determine convergence\nreltol=1e-8: relative tolerance of action gradient to determine convergence\nmethod = Adam(): minimization algorithm (see Optimization.jl)\n=0.1: step size parameter in gradient descent HeymannVandenEijnden implementation.\nverbose=false: if true, print additional output\nshow_progress=true: if true, display a progress bar\n\nOptimization algorithms\n\nThe method keyword argument takes solver methods of the Optimization.jl package; alternatively, the option solver = \"HeymannVandenEijnden\" uses the original gMAM algorithm[1].\n\n[1]: Heymann and Vanden-Eijnden, PRL (2008)\n\n\n\n\n\ngeometric_min_action_method(\n sys::CoupledSDEs,\n init::Matrix;\n maxiter,\n abstol,\n reltol,\n method,\n AD,\n ϵ,\n verbose,\n show_progress\n) -> MinimumActionPath{_A, _B, _C, Nothing, Nothing, Nothing, Nothing, Nothing} where {_A, _B<:Real, _C}\n\n\nRuns the geometric Minimum Action Method (gMAM) to find the minimum action path (instanton) from an initial condition init, given a system sys and total arc length arclength.\n\nThe initial path init must be a matrix of size (D, N), where D is the dimension of the system and N is the number of path points.\n\nFor more information see the main method, geometric_min_action_method(sys::CoupledSDEs, x_i, x_f, arclength::Float64; kwargs...).\n\n\n\n\n\n","category":"function"},{"location":"man/largedeviations/#Simple-Geometric-minimum-action-method-(sgMAM)","page":"Large deviation theory","title":"Simple Geometric minimum action method (sgMAM)","text":"","category":"section"},{"location":"man/largedeviations/","page":"Large deviation theory","title":"Large deviation theory","text":"Simplified minimization of the geometric action following Grafke et al. (2017). The simple gMAM reduces the complexity of the original gMAM by requiring only first-order derivatives of the underlying Hamiltonian optimisation formulation. This simplifies the numerical treatment and the computational complexity.","category":"page"},{"location":"man/largedeviations/","page":"Large deviation theory","title":"Large deviation theory","text":"The implementation below perform a constrained gradient descent where it assumes an autonomous system with additive Gaussian noise.","category":"page"},{"location":"man/largedeviations/","page":"Large deviation theory","title":"Large deviation theory","text":"sgmam\nSgmamSystem","category":"page"},{"location":"man/largedeviations/#CriticalTransitions.sgmam","page":"Large deviation theory","title":"CriticalTransitions.sgmam","text":"sgmam(\n sys::SgmamSystem,\n x_initial::Matrix{<:Real};\n ϵ,\n iterations,\n show_progress,\n reltol\n) -> MinimumActionPath{_A, _B, V, Nothing, Nothing, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}} where {_A, _B<:Real, V<:(AbstractVector)}\n\n\nPerforms the simplified geometric Minimal Action Method (sgMAM) on the given system sys. Our implementation is only valid for additive noise.\n\nThis method computes the optimal path in the phase space of a Hamiltonian system that minimizes the Freidlin–Wentzell action. The Hamiltonian functions H_x and H_p define the system's dynamics in a doubled phase. The initial state x_initial is evolved iteratively using constrained gradient descent with step size parameter ϵ over a specified number of iterations. The method can display a progress meter and will stop early if the relative tolerance reltol is achieved.\n\nThe function returns a tuple containing the final state, the action value, the Lagrange multipliers, the momentum, and the state derivatives. The implementation is based on the work of Grafke et al. (2019).\n\n\n\n\n\n","category":"function"},{"location":"man/largedeviations/#CriticalTransitions.SgmamSystem","page":"Large deviation theory","title":"CriticalTransitions.SgmamSystem","text":"A structure representing a system with Hamiltonian functions Hx and Hp.\n\nThis system operates in an extended phase space where the Hamiltonian is assumed to be quadratic in the extended momentum. The phase space coordinates x are doubled to form a 2n-dimensional extended phase space.\n\n\n\n\n\n","category":"type"},{"location":"man/largedeviations/#Action-functionals","page":"Large deviation theory","title":"Action functionals","text":"","category":"section"},{"location":"man/largedeviations/#Freidlin-Wentzell-action","page":"Large deviation theory","title":"Freidlin-Wentzell action","text":"","category":"section"},{"location":"man/largedeviations/","page":"Large deviation theory","title":"Large deviation theory","text":"fw_action","category":"page"},{"location":"man/largedeviations/#CriticalTransitions.fw_action","page":"Large deviation theory","title":"CriticalTransitions.fw_action","text":"fw_action(sys::CoupledSDEs, path, time) -> Any\n\n\nCalculates the Freidlin-Wentzell action of a given path with time points time in a drift field specified by the deterministic dynamics f = dynamic_rule(sys) and (normalized) noise covariance matrix covariance_matrix(sys).\n\nThe path must be a (D x N) matrix, where D is the dimensionality of the system sys and N is the number of path points. The time array must have length N.\n\nReturns a single number, which is the value of the action functional\n\nS_Tphi_t = frac12 int_0^T dot phi_t - f(phi_t) ^2_Q textdt\n\nwhere phi_t denotes the path in state space, b the drift field, and T the total time of the path. The subscript Q refers to the generalized norm a_Q^2 = langle a Q^-1 b rangle (see anorm). Here Q is the noise covariance matrix normalized by DL_1(Q), with L_1 being the L1 matrix norm.\n\n\n\n\n\n","category":"function"},{"location":"man/largedeviations/#Geometric-Freidlin-Wentzell-action","page":"Large deviation theory","title":"Geometric Freidlin-Wentzell action","text":"","category":"section"},{"location":"man/largedeviations/","page":"Large deviation theory","title":"Large deviation theory","text":"geometric_action","category":"page"},{"location":"man/largedeviations/#CriticalTransitions.geometric_action","page":"Large deviation theory","title":"CriticalTransitions.geometric_action","text":"geometric_action(sys::CoupledSDEs, path) -> Any\ngeometric_action(sys::CoupledSDEs, path, arclength) -> Any\n\n\nCalculates the geometric action of a given path with specified arclength for the drift field specified by the deterministic dynamics f = dynamic_rule(sys) and (normalized) noise covariance matrix covariance_matrix(sys).\n\nFor a given path varphi, the geometric action bar S corresponds to the minimum of the Freidlin-Wentzell action S_T(varphi) over all travel times T0, where varphi denotes the path's parameterization in physical time (see fw_action). It is given by the integral\n\nbar Svarphi = int_0^L left( varphi_Q f(varphi)_Q - langle varphi f(varphi) rangle_Q right) textds\n\nwhere s is the arclength coordinate, L the arclength, f the drift field, and the subscript Q refers to the generalized dot product langle a b rangle_Q = a^top cdot Q^-1 b (see anorm). Here Q is the noise covariance matrix normalized by DL_1(Q), with L_1 being the L1 matrix norm.\n\nReturns the value of the geometric action bar S.\n\n\n\n\n\n","category":"function"},{"location":"man/largedeviations/#Onsager-Machlup-action","page":"Large deviation theory","title":"Onsager-Machlup action","text":"","category":"section"},{"location":"man/largedeviations/","page":"Large deviation theory","title":"Large deviation theory","text":"om_action","category":"page"},{"location":"man/largedeviations/#CriticalTransitions.om_action","page":"Large deviation theory","title":"CriticalTransitions.om_action","text":"om_action(sys::CoupledSDEs, path, time, noise_strength)\n\nCalculates the Onsager-Machlup action of a given path with time points time for the drift field f = dynamic_rule(sys) at given noise_strength.\n\nThe path must be a (D x N) matrix, where D is the dimensionality of the system sys and N is the number of path points. The time array must have length N.\n\nReturns a single number, which is the value of the action functional\n\nS^sigma_Tphi_t = frac12 int_0^T left( dot phi_t - f(phi_t) ^2_Q + sigma^2 nabla cdot f right) textd t\n\nwhere phi_t denotes the path in state space, b the drift field, T the total time of the path, and sigma the noise strength. The subscript Q refers to the generalized norm a_Q^2 = langle a Q^-1 b rangle (see anorm). Here Q is the noise covariance matrix normalized by DL_1(Q), with L_1 being the L1 matrix norm.\n\n\n\n\n\n","category":"function"},{"location":"man/largedeviations/","page":"Large deviation theory","title":"Large deviation theory","text":"For convenience, a general action function is available where the type of functional is set as an argument:","category":"page"},{"location":"man/largedeviations/","page":"Large deviation theory","title":"Large deviation theory","text":"action","category":"page"},{"location":"man/largedeviations/#CriticalTransitions.action","page":"Large deviation theory","title":"CriticalTransitions.action","text":"action(\n sys::CoupledSDEs,\n path::Matrix,\n time,\n functional;\n noise_strength\n) -> Any\n\n\nComputes the action functional specified by functional for a given CoupledSDEs sys and path parameterized by time.\n\nfunctional = \"FW\": Returns the Freidlin-Wentzell action (fw_action)\nfunctional = \"OM\": Returns the Onsager-Machlup action (om_action)\n\n\n\n\n\n","category":"function"},{"location":"man/systemanalysis/#Analyzing-a-system's-stability-properties","page":"Stability analysis","title":"Analyzing a system's stability properties","text":"","category":"section"},{"location":"man/systemanalysis/","page":"Stability analysis","title":"Stability analysis","text":"To use the following functionalities, you need to load ChoasTools.jl and Attractors.jl.","category":"page"},{"location":"man/systemanalysis/#Fixed-points","page":"Stability analysis","title":"Fixed points","text":"","category":"section"},{"location":"man/systemanalysis/","page":"Stability analysis","title":"Stability analysis","text":"fixedpoints","category":"page"},{"location":"man/systemanalysis/#ChaosTools.fixedpoints","page":"Stability analysis","title":"ChaosTools.fixedpoints","text":"fixedpoints(\n sys::CoupledSDEs,\n bmin::Vector,\n bmax::Vector\n) -> Tuple{StateSpaceSet, Vector{Vector{ComplexF64}}, Vector{Bool}}\n\n\nReturns fixed points, their eigenvalues and stability of the system sys within the state space volume defined by bmin and bmax.\n\nThis is a wrapper around the fixedpoints function of DynamicalSystems.jl.\n\nInput\n\nbmin (Vector): lower limits of the state space box to be considered, as a vector of coordinates\nbmax (Vector): upper limits\nalternatively box (IntervalBox) can replace bmin and bmax\n\nExample: fixedpoints(sys, [-2,-1,0], [2,1,1]) finds the fixed points of the 3D system sys in a cube defined by the intervals [-2,2] × [-1,1] × [0,1].\n\nOutput\n\n[fp, eigs, stable]\n\nfp: StateSpaceSet of fixed points\neigs: vector of Jacobian eigenvalues of each fixed point\nstable: vector of booleans indicating the stability of each fixed point (true=stable, false=unstable)\n\nAdditional methods\n\nfixedpoints(sys::CoupedSDEs, box)\n\n\n\n\n\n","category":"function"},{"location":"man/systemanalysis/#Edge-tracking","page":"Stability analysis","title":"Edge tracking","text":"","category":"section"},{"location":"man/systemanalysis/","page":"Stability analysis","title":"Stability analysis","text":"The edge tracking algorithm is a simple numerical method to find the edge state or (possibly chaotic) saddle on the boundary between two basins of attraction. It is first introduced by Battelino et al. (1988) and further described by Skufca et al. (2006).","category":"page"},{"location":"man/systemanalysis/","page":"Stability analysis","title":"Stability analysis","text":"edgetracking\nbisect_to_edge","category":"page"},{"location":"man/systemanalysis/#Attractors.edgetracking","page":"Stability analysis","title":"Attractors.edgetracking","text":"edgetracking(\n sys::CoupledSDEs,\n u1,\n u2,\n attractors::Vector;\n eps1,\n eps2,\n converge,\n dt,\n tmax,\n ϵ_mapper,\n dt_mapper,\n solver,\n maxit,\n verbose,\n output_all,\n kwargs...\n)\n\n\nRuns the edge tracking algorithm.\n\nInput\n\nsys: dynamical system of type CoupledSDEs\nu1, u2: initial states; must belong to different basins of attraction\nattractors: vector of state vectors corresponding to the stable fixed points of sys\n\nKeyword arguments\n\neps1 = 1e-9: tolerance for bisection distance\neps2 = 1e-8: tolerance for divergence of trajectories before re-bisecting\nconverge = 1e-5: convergence criterion for M state accuracy (Euclidean distance)\ndt = 0.01: integration time step\ntmax = Inf: maximum integration time of parallel trajectories until re-bisection\nϵ_mapper = 0.1: distance threshold for AttractorMapper\ndt_mapper = 0.01: time step for AttractorMapper (keyword argument Δt)\nsolver = Vern9(): ODE solver from DifferentialEquations.jl\nmaxit = 100: maximum number of iterations before algorithm stops\nverbose = true: print status updates during run\noutput_all = false: if false, returns M state, else returns all points of the track\nkwargs...: additional keyword arguments of AttractorsViaProximity may be passed\n\nReturns\n\nIf output_all, a single state vector corresponding to the found edge state is returned. Else, a triple edge, track1, track2 is returned, where track1 and track2 are the tracks along the edge within the basin of attraction of u1 and u2, respectively; edge is the track along the edge derived from (track1 + track2)/2.\n\nwarning: Warning\nMay behave erroneously when run with solver = SimpleATsit5(), which is the default solver for AttractorsViaProximity. The recommended solver is Vern9().\n\n\n\n\n\n","category":"function"},{"location":"man/systemanalysis/#Attractors.bisect_to_edge","page":"Stability analysis","title":"Attractors.bisect_to_edge","text":"bisect_to_edge(\n sys::CoupledSDEs,\n u1,\n u2,\n attractors::Vector;\n eps1,\n ϵ_mapper,\n dt_mapper,\n solver,\n maxit,\n absto,\n relto,\n kwargs...\n)\n\n\nBisects to the basin boundary between two initial points u1 and u2. Returns the two final points, one on each side of the basin boundary, that are less than eps1 apart from each other.\n\nKeyword arguments\n\neps1=1e-9: tolerance for final distance between the two states\nϵ_mapper=0.1: ϵ parameter for AttractorMapper\ndt_mapper = 0.01: time step for AttractorMapper (keyword argument Δt)\nsolver=Vern9(): solver for AttractorMapper\nkwargs...: additional kwargs that can be passed to AttractorMapper\n\n\n\n\n\n","category":"function"},{"location":"man/systemanalysis/#Baisin-of-attraction","page":"Stability analysis","title":"Baisin of attraction","text":"","category":"section"},{"location":"man/systemanalysis/","page":"Stability analysis","title":"Stability analysis","text":"basinboundary\nbasboundary\nbasins","category":"page"},{"location":"man/systemanalysis/#CriticalTransitions.basinboundary","page":"Stability analysis","title":"CriticalTransitions.basinboundary","text":"basinboundary(X, Y, h; coords, A, B, C) -> Any\n\n\nComputes the basin boundary for given output X, Y, h of the basins function.\n\nTo be further documented.\n\n\n\n\n\nbasinboundary(boa) -> Any\n\n\nComputes the basin boundary for given output boa of the basins function.\n\nTo be further documented.\n\n\n\n\n\n","category":"function"},{"location":"man/systemanalysis/#CriticalTransitions.basboundary","page":"Stability analysis","title":"CriticalTransitions.basboundary","text":"basboundary(\n sys::CoupledSDEs,\n xrange::Vector,\n yrange::Vector,\n xspacing::Float64,\n attractors::Vector;\n eps1,\n ϵ_mapper,\n dt_mapper,\n solver,\n maxit\n) -> Tuple{Any, Any}\n\n\nThis function computes the basin boundary.\n\n\n\n\n\n","category":"function"},{"location":"man/systemanalysis/#CriticalTransitions.basins","page":"Stability analysis","title":"CriticalTransitions.basins","text":"basins(\n sys::CoupledSDEs,\n A,\n B,\n C,\n H;\n bstep,\n pstep,\n ϵ_mapper,\n show_progress,\n kwargs...\n) -> Any\n\n\nComputes the basins of attraction of CoupledSDEs sys on a plane spanned by the distinct points A, B, C and limited by the box H. Uses the AttractorsViaProximity function from DynamicalSystems.jl to compute the basins of attraction.\n\nA, B, C are elements of mathbbR^d (where d is the dimension of the sys) and H is a hyperrectangle in mathbbR^d.\n\nThe plane is given by P_UV = A+u(B-A)+v(C-A)inmathbbR^d u in U vin V for some closed and bounded real intervals U and V which are selected such that both\n\nP_UV subseteq H, and\nUtimes VsubseteqmathbbR^2 has maximal area,\n\ni.e. P_UV is the \"largest\" possible plane contained within H. This plane is determined behind the scenes.\n\nThis function returns a four-dimensional vector. The first two entries are discretised versions of the interval U and V (as defined above, of lengths ell_Uell_V respectively); the third entry is a dictionary of the attractors (stable equilibria) of the system within H, and the final entry is an ell_Vtimesell_U matrix of integers that group the initial conditions (written in terms of A+u(B-A)+v(C-A) where uin U and vin V) by which attractor they will in time converge to.\n\nKeyword arguments\n\nbstep = [0.01, 0.01]: a vector of length two whose elements respectively specify the length of the incremental steps taken across each dimension in the discretisation of your plane\npstep = [0.1, 0.1]: a vector of length two whose elements give the increments of the mesh that the maximisation process of finding a plane from a box is taken over (for more information see the source code of the function plane in the src/systemanalysis/planeofbox.jl file)\nϵ_mapper = 0.01: ϵ parameter of AttractorsViaProximity\nkwargs...: keyword arguments passed to the AttractorsViaProximity function (namely, Ttr, Δt, horizon_limit, mx_chk_lost)\nshow_progress = false: if true, displays a progress bar\n\n\n\n\n\n","category":"function"},{"location":"examples/sgMAM_KPO/","page":"sgMAM for the Kerr Parametric Oscillator","title":"sgMAM for the Kerr Parametric Oscillator","text":"We demonstrate the simple geometric minimum action method (sgMAM) on the Kerr parametric oscillator (KPO). The method computes the optimal path between two attractors in the phase space that minimizes the action of the system. It is a simplification of the geometric minimum action method (gMAM) by avoiding the computation of the second order derivatives of the extended Hamiltonian of the optimisation problem.","category":"page"},{"location":"examples/sgMAM_KPO/","page":"sgMAM for the Kerr Parametric Oscillator","title":"sgMAM for the Kerr Parametric Oscillator","text":"using CriticalTransitions, CairoMakie","category":"page"},{"location":"examples/sgMAM_KPO/","page":"sgMAM for the Kerr Parametric Oscillator","title":"sgMAM for the Kerr Parametric Oscillator","text":"The KPO equation is a nonlinear ordinary differential equation that describes the response of the nonlinear parametrically driven resonator at its dominant resonant condition. The equation of motion are of the form:","category":"page"},{"location":"examples/sgMAM_KPO/","page":"sgMAM for the Kerr Parametric Oscillator","title":"sgMAM for the Kerr Parametric Oscillator","text":"dotmathbfx = mathbff(mathbfx) + sigmamathbfξ(t)","category":"page"},{"location":"examples/sgMAM_KPO/","page":"sgMAM for the Kerr Parametric Oscillator","title":"sgMAM for the Kerr Parametric Oscillator","text":"where f is an autonemous drift function and and we have brownian noise ξ with intensity σ.","category":"page"},{"location":"examples/sgMAM_KPO/","page":"sgMAM for the Kerr Parametric Oscillator","title":"sgMAM for the Kerr Parametric Oscillator","text":"Here we define the define the drift of each seperable variable u and v. In addition, we hard-code the Jacobian of the drift function. ","category":"page"},{"location":"examples/sgMAM_KPO/","page":"sgMAM for the Kerr Parametric Oscillator","title":"sgMAM for the Kerr Parametric Oscillator","text":"const λ = 3 / 1.21 * 2 / 295\nconst ω0 = 1.000\nconst ω = 1.000\nconst γ = 1 / 295\nconst η = 0\nconst α = -1\n\nfunction fu(u, v)\n return (-4 * γ * ω * u - 2 * λ * v - 4 * (ω0 - ω^2) * v - 3 * α * v * (u^2 + v^2)) /\n (8 * ω)\nend\nfunction fv(u, v)\n return (-4 * γ * ω * v - 2 * λ * u + 4 * (ω0 - ω^2) * u + 3 * α * u * (u^2 + v^2)) /\n (8 * ω)\nend\nstream(u,v) = Point2f(fu(u, v), fv(u, v))\ndfvdv(u, v) = (-4 * γ * ω + 6 * α * u * v) / (8 * ω)\ndfudu(u, v) = (-4 * γ * ω - 6 * α * u * v) / (8 * ω)\ndfvdu(u, v) = (-2 * λ + 4 * (ω0 - ω^2) + 9 * α * u^2 + 3 * α * v^2) / (8 * ω)\ndfudv(u, v) = (-2 * λ - 4 * (ω0 - ω^2) - 3 * α * u^2 - 9 * α * v^2) / (8 * ω)","category":"page"},{"location":"examples/sgMAM_KPO/","page":"sgMAM for the Kerr Parametric Oscillator","title":"sgMAM for the Kerr Parametric Oscillator","text":"The optimisation is performed in a doubled phase space, i.e., every variable of the SDE system is considered as a generelised coordinate mathbfx and gets a corresponding generalised momentum mathbfp. The makes it that also systems with dissipative flow can be solved. As such, we extend the phase space by defining the hamiltionian","category":"page"},{"location":"examples/sgMAM_KPO/","page":"sgMAM for the Kerr Parametric Oscillator","title":"sgMAM for the Kerr Parametric Oscillator","text":"H = sum_i fracp_i^22 + f_i(mathbfx)p_i","category":"page"},{"location":"examples/sgMAM_KPO/","page":"sgMAM for the Kerr Parametric Oscillator","title":"sgMAM for the Kerr Parametric Oscillator","text":"Hence, to use the sgMAM method, we need to define the derivatives of the Hamiltonian with respect to the phase space coordinates and the generalised momentum: ","category":"page"},{"location":"examples/sgMAM_KPO/","page":"sgMAM for the Kerr Parametric Oscillator","title":"sgMAM for the Kerr Parametric Oscillator","text":"function H_x(x, p) # ℜ² → ℜ²\n u, v = eachrow(x)\n pu, pv = eachrow(p)\n\n H_u = @. pu * dfudu(u, v) + pv * dfvdu(u, v)\n H_v = @. pu * dfudv(u, v) + pv * dfvdv(u, v)\n return Matrix([H_u H_v]')\nend\nfunction H_p(x, p) # ℜ² → ℜ²\n u, v = eachrow(x)\n pu, pv = eachrow(p)\n\n H_pu = @. pu + fu(u, v)\n H_pv = @. pv + fv(u, v)\n return Matrix([H_pu H_pv]')\nend\n\nsys = SgmamSystem(H_x, H_p)","category":"page"},{"location":"examples/sgMAM_KPO/","page":"sgMAM for the Kerr Parametric Oscillator","title":"sgMAM for the Kerr Parametric Oscillator","text":"We saved this function in the SgmamSystem struct. We want to find the optimal path between two attractors in the phase space. We define the initial trajectory as wiggle between the two attractors.","category":"page"},{"location":"examples/sgMAM_KPO/","page":"sgMAM for the Kerr Parametric Oscillator","title":"sgMAM for the Kerr Parametric Oscillator","text":"# setup\nNt = 500 # number of discrete time steps\ns = collect(range(0; stop=1, length=Nt))\n\nxa = [-0.0208, 0.0991]\nxb = -xa\nxsaddle = [0.0, 0.0]\n\n# Initial trajectory\nxx = @. (xb[1] - xa[1]) * s + xa[1] + 4 * s * (1 - s) * xsaddle[1]\nyy = @. (xb[2] - xa[2]) * s + xa[2] + 4 * s * (1 - s) * xsaddle[2] + 0.01 * sin(2π * s)\nx_initial = Matrix([xx yy]')","category":"page"},{"location":"examples/sgMAM_KPO/","page":"sgMAM for the Kerr Parametric Oscillator","title":"sgMAM for the Kerr Parametric Oscillator","text":"The optimisation is the performed by the sgmam function:","category":"page"},{"location":"examples/sgMAM_KPO/","page":"sgMAM for the Kerr Parametric Oscillator","title":"sgMAM for the Kerr Parametric Oscillator","text":"MLP = sgmam(\n sys, x_initial; iterations=1_000, ϵ=10e2, show_progress=false\n)\nx_min = MLP.path;","category":"page"},{"location":"examples/sgMAM_KPO/","page":"sgMAM for the Kerr Parametric Oscillator","title":"sgMAM for the Kerr Parametric Oscillator","text":"The function returns the optimal path x_min, the minmal action S_min, the Lagrange multipliers lambda associated with the optimal path, the optimal generalised momentum p, and the time derivative of the optimal path xdot. We can plot the initial trajectory and the optimal path:","category":"page"},{"location":"examples/sgMAM_KPO/","page":"sgMAM for the Kerr Parametric Oscillator","title":"sgMAM for the Kerr Parametric Oscillator","text":"fig, ax, _ = lines(x_initial[1, :], x_initial[2, :]; label=\"init\", linewidth=3, color=:black)\nlines!(x_min[1, :], x_min[2, :]; label=\"MLP\", linewidth=3, color=:red)\nstreamplot!(ax, stream, (-0.08, 0.08), (-0.15, 0.15);\n gridsize=(20, 20), arrow_size=10, stepsize=0.001,\n colormap=[:gray, :gray]\n)\naxislegend(ax)\nfig","category":"page"},{"location":"man/sampling/#Sampling-transitions","page":"Sampling transitions","title":"Sampling transitions","text":"","category":"section"},{"location":"man/sampling/","page":"Sampling transitions","title":"Sampling transitions","text":"TransitionEnsemble","category":"page"},{"location":"man/sampling/#CriticalTransitions.TransitionEnsemble","page":"Sampling transitions","title":"CriticalTransitions.TransitionEnsemble","text":"struct TransitionEnsemble{SSS, T, Tstat, ES}\n\nEnsemble of transition paths between two points in a state space.\n\nFields\n\npaths::Vector\ntimes::Vector\nsuccess_rate::Any\nresidence_time::Any\ntransition_time::Any\nsciml_ensemble::Any\n\nConstructors\n\nTransitionEnsemble(\n paths,\n times,\n success_rate,\n residence_time,\n transition_time,\n sciml_ensemble\n)\n\ndefined at /home/runner/work/CriticalTransitions.jl/CriticalTransitions.jl/src/trajectories/transition.jl:14.\n\n\n\n\n\n","category":"type"},{"location":"man/sampling/#...-by-direct-simulation","page":"Sampling transitions","title":"... by direct simulation","text":"","category":"section"},{"location":"man/sampling/","page":"Sampling transitions","title":"Sampling transitions","text":"These functions generate noise-induced transitions between an initial and final state.","category":"page"},{"location":"man/sampling/","page":"Sampling transitions","title":"Sampling transitions","text":"transition\ntransitions","category":"page"},{"location":"man/sampling/#CriticalTransitions.transition","page":"Sampling transitions","title":"CriticalTransitions.transition","text":"transition(\n sys::CoupledSDEs,\n x_i,\n x_f;\n radii,\n tmax,\n radius_directions,\n cut_start,\n kwargs...\n) -> Tuple{Any, Any, Bool}\n\n\nGenerates a sample transition from point x_i to point x_f.\n\nThis function simulates sys in time, starting from initial condition x_i, until entering a ball of given radius around x_f.\n\nKeyword arguments\n\nradii=(0.1, 0.1): radius of the ball around x_i and x_f, respectively\ntmax=1e3: maximum time until the simulation stops even x_f has not been reached\nradius_directions=1:length(sys.u): the directions in phase space to consider when calculating the radii rad_i and rad_f. Defaults to all directions. To consider only a subspace of state space, insert a vector of indices of the dimensions to be included.\ncut_start=true: if false, returns the whole trajectory up to the transition\nkwargs...: keyword arguments passed to CommonSolve.solve\n\nOutput\n\n[path, times, success]\n\npath (Matrix): transition path (size [dim × N], where N is the number of time points)\ntimes (Vector): time values (since start of simulation) of the path points (size N)\nsuccess (bool): if true, a transition occured (i.e. the ball around x_f has been reached), else false\n\nSee also transitions, trajectory.\n\n\n\n\n\n","category":"function"},{"location":"man/sampling/#CriticalTransitions.transitions","page":"Sampling transitions","title":"CriticalTransitions.transitions","text":"function transitions(sys::CoupledSDEs, x_i, x_f, N=1; kwargs...)\n\nGenerates an ensemble of N transition samples of sys from point x_i to point x_f.\n\nThis function repeatedly calls the transition function to efficiently generate an ensemble of transitions. Multi-threading is enabled.\n\nKeyword arguments\n\nradii=(0.1, 0.1): radius of the ball around x_i and x_f, respectively\nNmax: number of attempts before the algorithm stops even if less than N transitions occurred.\ntmax=1e3: maximum time when the simulation stops even x_f has not been reached\nradius_directions=1:length(sys.u): the directions in phase space to consider when calculating the radii rad_i and rad_f. Defaults to all directions. To consider only a subspace of state space, insert a vector of indices of the dimensions to be included.\ncut_start=true: if false, returns the whole trajectory up to the transition\nshow_progress=true: shows a progress bar with respect to Nmax\nkwargs...: keyword arguments passed to CommonSolve.solve\n\nSee also transition.\n\nReturns a TransitionEnsemble object.\n\n\n\n\n\n","category":"function"},{"location":"man/sampling/#...-in-pathspace","page":"Sampling transitions","title":"... in pathspace","text":"","category":"section"},{"location":"man/sampling/","page":"Sampling transitions","title":"Sampling transitions","text":"Coming soon...","category":"page"},{"location":"quickstart/#Quickstart","page":"Quickstart","title":"Quickstart","text":"","category":"section"},{"location":"quickstart/#Installation","page":"Quickstart","title":"Installation","text":"","category":"section"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"As this module is not published yet, there are two ways to access it:","category":"page"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"Option 1 (recommended): Install from GitHub\nEnter the Julia package manager by typing ] in the REPL: julia> ]\ntype add https://github.com/juliadynamics/CriticalTransitions.jl.git\nOption 2: Load module locally\nClone the repo: git clone https://github.com/juliadynamics/CriticalTransitions.jl.git\nIn Julia, include the module file: include(\"PATH/src/CriticalTransitions.jl\"), where PATH is the relative path to the repo you just cloned\nLoad the module: using .CriticalTransitions","category":"page"},{"location":"quickstart/#Basic-usage","page":"Quickstart","title":"Basic usage","text":"","category":"section"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"The general workflow of CriticalTransitions essentially follows two steps:","category":"page"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"Define your system (see Define a CoupledSDEs system)\nInvestigate the system by calling methods (see Methods)","category":"page"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"info: New system type: RateSystem\nWe are planning to introduce the the struct RateSystem along CoupledSDEs. In a RateSystem, the time dependence of parameters can conveniently be specified, laying the foundation for a toolbox to study rate-induced tipping, or R-tipping.","category":"page"},{"location":"quickstart/#Methods","page":"Quickstart","title":"Methods","text":"","category":"section"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"Currently the following functions are implemented to analyze a CoupledSDEs and corresponding sample transition paths.","category":"page"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"Pages = [\"man/simulation.md\", \"man/sampling.md\", \"man/largedeviations.md\"]","category":"page"},{"location":"#CriticalTransitions.jl","page":"Home","title":"CriticalTransitions.jl","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"A Julia package for the numerical investigation of noise- and rate-induced transitions in dynamical systems.","category":"page"},{"location":"","page":"Home","title":"Home","text":"Building on DynamicalSystems.jl and DifferentialEquations.jl, this package aims to provide a toolbox for dynamical systems under time-dependent forcing, with a focus on tipping phenomena and metastability.","category":"page"},{"location":"","page":"Home","title":"Home","text":"(Image: CT.jl infographic)","category":"page"},{"location":"","page":"Home","title":"Home","text":"info: Current features\nStochastic simulation made easy: Gaussian noise, uncorrelated and correlated, additive and multiplicative\nTransition path sampling: Parallelized ensemble rejection sampling\nLarge deviation theory tools: Action functionals and minimization algorithms (MAM, gMAM)","category":"page"},{"location":"","page":"Home","title":"Home","text":"ukw: Planned features\nRare event simulation: importance sampling, AMS\nQuasipotentials: Ordered line integral method (OLIM)\nRate-induced tipping tools\nSymbolic differentiation of action functionals\n...?","category":"page"},{"location":"","page":"Home","title":"Home","text":"Developers: Reyk Börner, Ryan Deeley, Raphael Römer and Orjan Ameye","category":"page"},{"location":"","page":"Home","title":"Home","text":"Thanks to Jeroen Wouters, Calvin Nesbitt, Tobias Grafke, George Datseris and Oliver Mehling","category":"page"},{"location":"","page":"Home","title":"Home","text":"This work is part of the CriticalEarth project.","category":"page"},{"location":"man/CoupledSDEs/#Define-a-CoupledSDEs-system","page":"Define a CoupledSDEs system","title":"Define a CoupledSDEs system","text":"","category":"section"},{"location":"man/CoupledSDEs/","page":"Define a CoupledSDEs system","title":"Define a CoupledSDEs system","text":"CoupledSDEs","category":"page"},{"location":"man/CoupledSDEs/#DynamicalSystemsBase.CoupledSDEs","page":"Define a CoupledSDEs system","title":"DynamicalSystemsBase.CoupledSDEs","text":"CoupledSDEs(ds::CoupledODEs, p; kwargs...)\n\nConverts a CoupledODEs system into a CoupledSDEs.\n\n\n\n\n\n","category":"type"},{"location":"man/CoupledSDEs/","page":"Define a CoupledSDEs system","title":"Define a CoupledSDEs system","text":"info: Info\nNote that nonlinear mixings of the Noise Process mathcalW fall into the class of random ordinary differential equations (RODEs) which have a separate set of solvers. See this example of DifferentialEquations.jl.","category":"page"},{"location":"man/CoupledSDEs/#defining-stochastic-dynamics","page":"Define a CoupledSDEs system","title":"Examples: Defining stochastic dynamics","text":"","category":"section"},{"location":"man/CoupledSDEs/","page":"Define a CoupledSDEs system","title":"Define a CoupledSDEs system","text":"Let's look at some examples of the different types of stochastic systems that can be defined.","category":"page"},{"location":"man/CoupledSDEs/","page":"Define a CoupledSDEs system","title":"Define a CoupledSDEs system","text":"For simplicity, we choose a slow exponential growth in 2 dimensions as the deterministic dynamics f:","category":"page"},{"location":"man/CoupledSDEs/","page":"Define a CoupledSDEs system","title":"Define a CoupledSDEs system","text":"using DynamicalSystemsBase, StochasticDiffEq, DiffEqNoiseProcess\nusing CairoMakie\nimport Random # hide\nRandom.seed!(10) # hide\nf!(du, u, p, t) = du .= 1.01u # deterministic part\n\nfunction plot_trajectory(Y, t)\n fig = Figure()\n ax = Axis(fig[1,1]; xlabel = \"time\", ylabel = \"variable\")\n for var in columns(Y)\n lines!(ax, t, var)\n end\n fig\nend;","category":"page"},{"location":"man/CoupledSDEs/#Additive-noise","page":"Define a CoupledSDEs system","title":"Additive noise","text":"","category":"section"},{"location":"man/CoupledSDEs/","page":"Define a CoupledSDEs system","title":"Define a CoupledSDEs system","text":"When g(u p t) is independent of the state u, the noise is called additive; otherwise, it is multiplicative. We can define a simple additive noise system as follows:","category":"page"},{"location":"man/CoupledSDEs/","page":"Define a CoupledSDEs system","title":"Define a CoupledSDEs system","text":"sde = CoupledSDEs(f!, zeros(2));","category":"page"},{"location":"man/CoupledSDEs/","page":"Define a CoupledSDEs system","title":"Define a CoupledSDEs system","text":"which is equivalent to","category":"page"},{"location":"man/CoupledSDEs/","page":"Define a CoupledSDEs system","title":"Define a CoupledSDEs system","text":"t0 = 0.0; W0 = zeros(2);\nW = WienerProcess(t0, W0, 0.0)\nsde = CoupledSDEs(f!, zeros(2);\n noise_process=W, covariance=[1 0; 0 1], noise_strength=1.0\n );","category":"page"},{"location":"man/CoupledSDEs/","page":"Define a CoupledSDEs system","title":"Define a CoupledSDEs system","text":"We defined a Wiener process W, whose increments are vectors of normally distributed random numbers of length matching the output of g. The noise is applied element-wise, i.e., g.*dW. Since the noise processes are uncorrelated, meaning the covariance matrix is diagonal, this type of noise is referred to as diagonal.","category":"page"},{"location":"man/CoupledSDEs/","page":"Define a CoupledSDEs system","title":"Define a CoupledSDEs system","text":"We can sample a trajectory from this system using the trajectory function also used for the deterministic systems:","category":"page"},{"location":"man/CoupledSDEs/","page":"Define a CoupledSDEs system","title":"Define a CoupledSDEs system","text":"tr = trajectory(sde, 1.0)\nplot_trajectory(tr...)","category":"page"},{"location":"man/CoupledSDEs/#Correlated-noise","page":"Define a CoupledSDEs system","title":"Correlated noise","text":"","category":"section"},{"location":"man/CoupledSDEs/","page":"Define a CoupledSDEs system","title":"Define a CoupledSDEs system","text":"In the case of correlated noise, the random numbers in a vector increment dW are correlated. This can be achieved by specifying the covariance matrix Sigma via the covariance keyword:","category":"page"},{"location":"man/CoupledSDEs/","page":"Define a CoupledSDEs system","title":"Define a CoupledSDEs system","text":"ρ = 0.3\nΣ = [1 ρ; ρ 1]\ndiffeq = (alg = LambaEM(), dt=0.1)\nsde = CoupledSDEs(f!, zeros(2); covariance=Σ, diffeq=diffeq)","category":"page"},{"location":"man/CoupledSDEs/","page":"Define a CoupledSDEs system","title":"Define a CoupledSDEs system","text":"Alternatively, we can parametrise the covariance matrix by defining the diffusion function g ourselves:","category":"page"},{"location":"man/CoupledSDEs/","page":"Define a CoupledSDEs system","title":"Define a CoupledSDEs system","text":"g!(du, u, p, t) = (du .= [1 p[1]; p[1] 1]; return nothing) \nsde = CoupledSDEs(f!, zeros(2), (ρ); g=g!, noise_prototype=zeros(2, 2))","category":"page"},{"location":"man/CoupledSDEs/","page":"Define a CoupledSDEs system","title":"Define a CoupledSDEs system","text":"Here, we had to provide noise_prototype to indicate that the diffusion function g will output a 2x2 matrix.","category":"page"},{"location":"man/CoupledSDEs/#Scalar-noise","page":"Define a CoupledSDEs system","title":"Scalar noise","text":"","category":"section"},{"location":"man/CoupledSDEs/","page":"Define a CoupledSDEs system","title":"Define a CoupledSDEs system","text":"If all state variables are forced by the same single random variable, we have scalar noise. To define scalar noise, one has to give an one-dimensional noise process to the noise_process keyword of the CoupledSDEs constructor. ","category":"page"},{"location":"man/CoupledSDEs/","page":"Define a CoupledSDEs system","title":"Define a CoupledSDEs system","text":"t0 = 0.0; W0 = 0.0;\nnoise = WienerProcess(t0, W0, 0.0)\nsde = CoupledSDEs(f!, rand(2)/10; noise_process=noise)\n\ntr = trajectory(sde, 1.0)\nplot_trajectory(tr...)","category":"page"},{"location":"man/CoupledSDEs/","page":"Define a CoupledSDEs system","title":"Define a CoupledSDEs system","text":"We can see that noise applied to each variable is the same.","category":"page"},{"location":"man/CoupledSDEs/#Multiplicative-and-time-dependent-noise","page":"Define a CoupledSDEs system","title":"Multiplicative and time-dependent noise","text":"","category":"section"},{"location":"man/CoupledSDEs/","page":"Define a CoupledSDEs system","title":"Define a CoupledSDEs system","text":"In the SciML ecosystem, multiplicative noise is defined through the condition g_i(t u)=a_i u. However, in the literature the name is more broadly used for any situation where the noise is non-additive and depends on the state u, possibly also in a non-linear way. When defining a CoupledSDEs, we can make the noise term time- and state-dependent by specifying an explicit time- or state-dependence in the noise function g, just like we would define f. For example, we can define a system with temporally decreasing multiplicative noise as follows:","category":"page"},{"location":"man/CoupledSDEs/","page":"Define a CoupledSDEs system","title":"Define a CoupledSDEs system","text":"function g!(du, u, p, t)\n du .= u ./ (1+t)\n return nothing\nend\nsde = CoupledSDEs(f!, rand(2)./10; g=g!)","category":"page"},{"location":"man/CoupledSDEs/#Non-diagonal-noise","page":"Define a CoupledSDEs system","title":"Non-diagonal noise","text":"","category":"section"},{"location":"man/CoupledSDEs/","page":"Define a CoupledSDEs system","title":"Define a CoupledSDEs system","text":"Non-diagonal noise allows for the terms to be linearly mixed (correlated) via g being a matrix. Suppose we have two Wiener processes and two state variables such that the output of g is a 2x2 matrix. Therefore, we have","category":"page"},{"location":"man/CoupledSDEs/","page":"Define a CoupledSDEs system","title":"Define a CoupledSDEs system","text":"du_1 = f_1(upt)dt + g_11(upt)dW_1 + g_12(upt)dW_2 \ndu_2 = f_2(upt)dt + g_21(upt)dW_1 + g_22(upt)dW_2","category":"page"},{"location":"man/CoupledSDEs/","page":"Define a CoupledSDEs system","title":"Define a CoupledSDEs system","text":"To indicate the structure that g should have, we must use the noise_prototype keyword. Let us define a special type of non-diagonal noise called commutative noise. For this we can utilize the RKMilCommute algorithm which is designed to utilize the structure of commutative noise.","category":"page"},{"location":"man/CoupledSDEs/","page":"Define a CoupledSDEs system","title":"Define a CoupledSDEs system","text":"σ = 0.25 # noise strength\nfunction g!(du, u, p, t)\n du[1,1] = σ*u[1]\n du[2,1] = σ*u[2]\n du[1,2] = σ*u[1]\n du[2,2] = σ*u[2]\n return nothing\nend\ndiffeq = (alg = RKMilCommute(), reltol = 1e-3, abstol = 1e-3, dt=0.1)\nsde = CoupledSDEs(f!, rand(2)./10; g=g!, noise_prototype = zeros(2, 2), diffeq = diffeq)","category":"page"},{"location":"man/CoupledSDEs/","page":"Define a CoupledSDEs system","title":"Define a CoupledSDEs system","text":"warning: Warning\nNon-diagonal problems need specific solvers. See the SciML recommendations.","category":"page"},{"location":"man/CoupledSDEs/#Interface-to-DynamicalSystems.jl","page":"Define a CoupledSDEs system","title":"Interface to DynamicalSystems.jl","text":"","category":"section"},{"location":"man/CoupledSDEs/#Converting-between-CoupledSDEs-and-CoupledODEs","page":"Define a CoupledSDEs system","title":"Converting between CoupledSDEs and CoupledODEs","text":"","category":"section"},{"location":"man/CoupledSDEs/","page":"Define a CoupledSDEs system","title":"Define a CoupledSDEs system","text":"tip: Analyzing deterministic dynamics with DynamicalSystems.jl\nThe deterministic part of a CoupledSDEs system can easily be extracted as a CoupledODEs, a common subtype of a ContinuousTimeDynamicalSystem in DynamicalSystems.jl.","category":"page"},{"location":"man/CoupledSDEs/","page":"Define a CoupledSDEs system","title":"Define a CoupledSDEs system","text":"CoupledODEs(sde::CoupledSDEs) extracts the deterministic part of sde as a CoupledODEs.\nCoupledSDEs(ode::CoupledODEs; kwargs) turns ode into a CoupledSDEs.","category":"page"},{"location":"man/CoupledSDEs/","page":"Define a CoupledSDEs system","title":"Define a CoupledSDEs system","text":"CoupledODEs","category":"page"},{"location":"man/CoupledSDEs/#DynamicalSystemsBase.CoupledODEs","page":"Define a CoupledSDEs system","title":"DynamicalSystemsBase.CoupledODEs","text":"CoupledODEs(ds::CoupledSDEs; kwargs...)\n\nConverts a CoupledSDEs into a CoupledODEs by extracting the deterministic part of ds.\n\n\n\n\n\n","category":"type"},{"location":"man/CoupledSDEs/","page":"Define a CoupledSDEs system","title":"Define a CoupledSDEs system","text":"For example, the Lyapunov spectrum of a CoupledSDEs in the absence of noise, here exemplified by the FitzHugh-Nagumo model, can be computed by typing:","category":"page"},{"location":"man/CoupledSDEs/","page":"Define a CoupledSDEs system","title":"Define a CoupledSDEs system","text":"using CriticalTransitions\nusing DynamicalSystems: lyapunovspectrum\n\nfunction fitzhugh_nagumo(u, p, t)\n x, y = u\n ϵ, β, α, γ, κ, I = p\n\n dx = (-α * x^3 + γ * x - κ * y + I) / ϵ\n dy = -β * y + x\n\n return SA[dx, dy]\nend\n\np = [1.,3.,1.,1.,1.,0.]\n\nsys = CoupledSDEs(fitzhugh_nagumo, zeros(2), p; noise_strength=0.1)\nls = lyapunovspectrum(CoupledODEs(sys), 10000)","category":"page"}] +[{"location":"examples/tutorial/#Tutorial","page":"Tutorial","title":"Tutorial","text":"","category":"section"},{"location":"examples/tutorial/","page":"Tutorial","title":"Tutorial","text":"To give you an idea of how our package works, this tutorial provides some example code with explanations.","category":"page"},{"location":"examples/tutorial/#Example:-FitzHugh-Nagumo-model","page":"Tutorial","title":"Example: FitzHugh-Nagumo model","text":"","category":"section"},{"location":"examples/tutorial/","page":"Tutorial","title":"Tutorial","text":"Let's consider a simple 2-dimensional dynamical system - the FitzHugh-Nagumo model:","category":"page"},{"location":"examples/tutorial/","page":"Tutorial","title":"Tutorial","text":"beginaligned\nfracdudt = frac1epsilon left( -alpha u^3 + gamma u - kappa v + I right) \nfracdvdt = -beta v + u \nendaligned","category":"page"},{"location":"examples/tutorial/","page":"Tutorial","title":"Tutorial","text":"where epsilon is the parameter of time scale separation between the state variables u and v. The parameters alpha 0, beta 1, gamma0, and kappa0 are real constants, and I denotes a driving term.","category":"page"},{"location":"examples/tutorial/","page":"Tutorial","title":"Tutorial","text":"Let's investigate this system under stochastic forcing.","category":"page"},{"location":"examples/tutorial/#System-definition","page":"Tutorial","title":"System definition","text":"","category":"section"},{"location":"examples/tutorial/","page":"Tutorial","title":"Tutorial","text":"First, we need to translate the system equations above into Julia code.","category":"page"},{"location":"examples/tutorial/","page":"Tutorial","title":"Tutorial","text":"This works exactly as in DynamicalSystems.jl by defining a function f(u,p,t) which takes as input a vector u of state variables (u,v), a vector p of parameters, and time t. The function must return an array of flow increments (textdu, textdv). For performance reasons, it is advisable to return a StaticArray SA[du, dv] rather than just a Vector [du, dv].","category":"page"},{"location":"examples/tutorial/","page":"Tutorial","title":"Tutorial","text":"using CriticalTransitions\nimport Random # hide\nRandom.seed!(1) # hide\n\nfunction fitzhugh_nagumo(u,p,t)\n u, v = u\n ϵ, β, α, γ, κ, I = p\n\n du = (-α*u^3 + γ*u - κ*v + I)/ϵ\n dv = -β*v + u\n\n SA[du, dv]\nend","category":"page"},{"location":"examples/tutorial/","page":"Tutorial","title":"Tutorial","text":"tip: In-place vs. out-of-place\nThe function fitzhugh_nagumo(u,p,t) is defined out-of-place. It is also possible to define the system in-place as fitzhugh_nagumo!(du,u,p,t). For more info, see here.","category":"page"},{"location":"examples/tutorial/#CoupledSDE","page":"Tutorial","title":"CoupledSDE","text":"","category":"section"},{"location":"examples/tutorial/","page":"Tutorial","title":"Tutorial","text":"Next, we construct a stochastic system with the fitzhugh_nagumo equation as the deterministic part. Suppose we would like to force both state variables u and v with additive, uncorrelated Gaussian noise of intensity sigma. This is the default case. We simply write","category":"page"},{"location":"examples/tutorial/","page":"Tutorial","title":"Tutorial","text":"p = [1., 3., 1., 1., 1., 0.] # Parameters (ϵ, β, α, γ, κ, I)\nσ = 0.2 # noise strength\n\n# CoupledSDE\nsys = CoupledSDEs(fitzhugh_nagumo, zeros(2), p; noise_strength=σ)","category":"page"},{"location":"examples/tutorial/","page":"Tutorial","title":"Tutorial","text":"Here the first field fitzhugh_nagumo specifies the deterministic dynamics f (see Define a CoupledSDEs system). We have chosen zeros(2) as the initial state of the system, which is the second field. The length of this vector must match the system's dimensionality. In the (optional) third field, we specify the parameter vector p, which includes the parameters of f followed by the parameters of g (in this case, there are no parameters for g). Lastly, noise_strength sets the noise strength. Since we have not specified a noise process, the default case of an uncorrelated Wiener process is used.","category":"page"},{"location":"examples/tutorial/","page":"Tutorial","title":"Tutorial","text":"note: Multiplicative and/or correlated noise\nOf course, it is also possible to define more complicated noise processes than simple additive white noise. This is done by specifying a custom noise function and covariance matrix in the CoupledSDEs definition. For more info, see Define a CoupledSDEs system.","category":"page"},{"location":"examples/tutorial/","page":"Tutorial","title":"Tutorial","text":"That's it! Now we can apply the toolbox of CriticalTransitions to our stochastic FitzHugh-Nagumo system sys.","category":"page"},{"location":"examples/tutorial/#Find-stable-equilibria","page":"Tutorial","title":"Find stable equilibria","text":"","category":"section"},{"location":"examples/tutorial/","page":"Tutorial","title":"Tutorial","text":"For the parameters chosen above, the FitzHugh-Nagumo system is bistable. Let's compute the fixed points using the fixedpoints function. This function is borrowed from ChaosTools.jl and is loaded as an extension when we write using ChaosTools.","category":"page"},{"location":"examples/tutorial/","page":"Tutorial","title":"Tutorial","text":"using ChaosTools\n# Calculate fixed points and store the stable ones\neqs, eigs, stab = fixedpoints(sys, [-2,-2], [2,2])\nfp1, fp2 = eqs[stab]","category":"page"},{"location":"examples/tutorial/#Stochastic-simulation","page":"Tutorial","title":"Stochastic simulation","text":"","category":"section"},{"location":"examples/tutorial/","page":"Tutorial","title":"Tutorial","text":"Using the trajectory function, we now run a simulation of our system for 100 time units starting out from the fixed point fp1:","category":"page"},{"location":"examples/tutorial/","page":"Tutorial","title":"Tutorial","text":"sim = trajectory(sys, 100, fp1)","category":"page"},{"location":"examples/tutorial/","page":"Tutorial","title":"Tutorial","text":"In the keyword arguments, we have specified at which interval the solution is saved. Further keyword arguments can be used to change the solver (the default is SOSRA() for stochastic integration) and other settings.","category":"page"},{"location":"examples/tutorial/","page":"Tutorial","title":"Tutorial","text":"The simulated trajectory is stored in sim in the usual output format of the solve method of DifferentialEquations.jl, including the solution sim.u and the vector of time points sim.t. The solution can also be accessed as a matrix sim[i, t], where i is the i-th component of u and t the time index.","category":"page"},{"location":"examples/tutorial/","page":"Tutorial","title":"Tutorial","text":"Let's plot the result. Did the trajectory transition to the other attractor?","category":"page"},{"location":"examples/tutorial/","page":"Tutorial","title":"Tutorial","text":"using Plots\nplt = plot(sim[1][:,1], sim[1][:,2]; xlabel=\"u\", ylabel=\"v\", legend=false)\nscatter!([fp1[1], fp2[1]], [fp1[2], fp2[2]], color=:red, markersize=4)\nxlims!(-1.2, 1.2)\nylims!(-0.6, 0.6)\nplt","category":"page"},{"location":"examples/tutorial/","page":"Tutorial","title":"Tutorial","text":"Hopefully, this helped you to get started. For more info, check out the Manual section of these docs.","category":"page"},{"location":"man/simulation/#Simulating-the-system","page":"Simulating the system","title":"Simulating the system","text":"","category":"section"},{"location":"man/simulation/","page":"Simulating the system","title":"Simulating the system","text":"We provide two main functions to simulate a CoupledSDEs forward in time:","category":"page"},{"location":"man/simulation/","page":"Simulating the system","title":"Simulating the system","text":"DynamicalSystemsBase.trajectory, which integrates the stochastic CoupledSDEs system forward in time\ndeterministic_orbit, which integrates only the deterministic part of the CoupledSDEs system ","category":"page"},{"location":"man/simulation/","page":"Simulating the system","title":"Simulating the system","text":"deterministic_orbit","category":"page"},{"location":"man/simulation/#CriticalTransitions.deterministic_orbit","page":"Simulating the system","title":"CriticalTransitions.deterministic_orbit","text":"deterministic_orbit(\n sys::CoupledSDEs,\n T;\n ...\n) -> Tuple{Any, Any}\ndeterministic_orbit(\n sys::CoupledSDEs,\n T,\n init;\n diffeq,\n kwargs...\n) -> Tuple{Any, Any}\n\n\nSimulates the deterministic (noise-free) dynamics of CoupledSDEs sys in time for a duration T, starting at initial condition init.\n\nThis function is equivalent to calling trajectory on the deterministic part of the CoupledSDEs (with noise_strength=0). It works with the ODE solvers used for CoupledODEs.\n\nKeyword arguments\n\ndiffeq=(alg=Tsit5(), abstol = 1e-6, reltol = 1e-6): ODE solver settings (see CoupledODEs)\nkwargs...: keyword arguments passed to trajectory\n\nFor more info, see ODEProblem. For stochastic integration, see trajectory.\n\n\n\n\n\n","category":"function"},{"location":"man/utils/#Utility-functions","page":"Utilities","title":"Utility functions","text":"","category":"section"},{"location":"man/utils/#CoupledSDEs-utility-functions","page":"Utilities","title":"CoupledSDEs utility functions","text":"","category":"section"},{"location":"man/utils/","page":"Utilities","title":"Utilities","text":"solver\nStochasticSystemsBase.covariance_matrix\nStochasticSystemsBase.diffusion_matrix\nnoise_process\ndrift\ndiv_drift","category":"page"},{"location":"man/utils/#CriticalTransitions.solver","page":"Utilities","title":"CriticalTransitions.solver","text":"solver(ds::CoupledSDEs) -> Any\n\n\nReturns the SDE solver specified in the diffeq settings of the CoupledSDEs.\n\n\n\n\n\n","category":"function"},{"location":"man/utils/#StochasticSystemsBase.covariance_matrix","page":"Utilities","title":"StochasticSystemsBase.covariance_matrix","text":"covariance_matrix(ds::CoupledSDEs)\n\nReturns the covariance matrix of the stochastic term of the CoupledSDEs ds, provided that the diffusion function g can be expressed as a constant invertible matrix. If this is not the case, returns nothing.\n\nSee also diffusion_matrix.\n\n\n\n\n\n","category":"function"},{"location":"man/utils/#StochasticSystemsBase.diffusion_matrix","page":"Utilities","title":"StochasticSystemsBase.diffusion_matrix","text":"diffusion_matrix(ds::CoupledSDEs)\n\nReturns the diffusion matrix of the stochastic term of the CoupledSDEs ds, provided that the diffusion function g can be expressed as a constant invertible matrix. If this is not the case, returns nothing.\n\nNote: The diffusion matrix Σ is the square root of the noise covariance matrix Q (see covariance_matrix), defined via the Cholesky decomposition Q = Σ Σ^top.\n\n\n\n\n\n","category":"function"},{"location":"man/utils/#CriticalTransitions.noise_process","page":"Utilities","title":"CriticalTransitions.noise_process","text":"noise_process(sys::CoupledSDEs) -> Any\n\n\nFetches the stochastic process mathcalN specified in the intergrator of sys. Returns the type DiffEqNoiseProcess.NoiseProcess.\n\n\n\n\n\n","category":"function"},{"location":"man/utils/#CriticalTransitions.drift","page":"Utilities","title":"CriticalTransitions.drift","text":"drift(sys::CoupledSDEs{IIP}, x; t) -> Any\n\n\nReturns the deterministic drift f(x) of the CoupledSDEs sys at state x. For time-dependent systems, the time can be specified as a keyword argument t (by default t=0).\n\n\n\n\n\n","category":"function"},{"location":"man/utils/#CriticalTransitions.div_drift","page":"Utilities","title":"CriticalTransitions.div_drift","text":"div_drift(sys::CoupledSDEs, x) -> Any\ndiv_drift(sys::CoupledSDEs, x, t) -> Any\n\n\nComputes the divergence of the drift field f(x) at state x. For time- dependent systems, the time can be specified as a keyword argument t (by default t=0).\n\n\n\n\n\n","category":"function"},{"location":"man/utils/#More","page":"Utilities","title":"More","text":"","category":"section"},{"location":"man/utils/","page":"Utilities","title":"Utilities","text":"intervals_to_box","category":"page"},{"location":"man/utils/#CriticalTransitions.intervals_to_box","page":"Utilities","title":"CriticalTransitions.intervals_to_box","text":"intervals_to_box(bmin::Vector, bmax::Vector) -> Any\n\n\nGenerates a box from specifying the interval limits in each dimension.\n\nbmin (Vector): lower limit of the box in each dimension\nbmax (Vector): upper limit\n\nExample\n\nintervals_to_box([-2,-1,0], [2,1,1]) returns a 3D box of dimensions [-2,2] × [-1,1] × [0,1].\n\n\n\n\n\n","category":"function"},{"location":"examples/gMAM_Maierstein/#The-Maier-Stein-model.","page":"Anlyses for the Maier-Stein system","title":"The Maier-Stein model.","text":"","category":"section"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"using CriticalTransitions\n\nusing CairoMakie\nusing CairoMakie.Makie.MathTeXEngine: get_font\nfont = (;\n regular=get_font(:regular), bold=get_font(:bold),\n italic=get_font(:italic), bold_italic=get_font(:bolditalic)\n);","category":"page"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"Let us explore the features of CriticalTransitions.jl with Maier-Stein model.","category":"page"},{"location":"examples/gMAM_Maierstein/#Maier-stein-model","page":"Anlyses for the Maier-Stein system","title":"Maier-stein model","text":"","category":"section"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"The Maier-Stein model (J. Stat. Phys 83, 3–4 (1996)) is commonly used in the field of nonlinear dynamics for benchmarking Large Deviation Theory (LDT) techniques, e.g., stoachastic transitions between different stable states. It is a simple model that describes the dynamics of a system with two degrees of freedom u and v, and is given by the following set of ordinary differential equations:","category":"page"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"beginaligned\n dotu = u-u^3 - beta*u*v^2\n dotv = -alpha (1+u^2)*v\nendaligned","category":"page"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"The parameter alpha0 controls the strength of the drift field and beta0 represents the softening of that drift field.","category":"page"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"function meier_stein!(du, u, p, t) # in-place\n x, y = u\n du[1] = x - x^3 - 10 * x * y^2\n du[2] = -(1 + x^2) * y\nend\nfunction meier_stein(u, p, t) # out-of-place\n x, y = u\n dx = x - x^3 - 10 * x * y^2\n dy = -(1 + x^2) * y\n SA[dx, dy]\nend\nσ = 0.25\nsys = CoupledSDEs(meier_stein, zeros(2), (); noise_strength=σ)","category":"page"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"A good reference to read about the large deviations methods is this or this blog post by Tobias Grafke.","category":"page"},{"location":"examples/gMAM_Maierstein/#Attractors","page":"Anlyses for the Maier-Stein system","title":"Attractors","text":"","category":"section"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"We start by investigating the deterministic dynamics of the Maier-Stein model.","category":"page"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"The function fixed points return the attractors, their eigenvalues and stability within the state space volume defined by bmin and bmax.","category":"page"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"using ChaosTools\n\nu_min = -1.1;\nu_max = 1.1;\nv_min = -0.4;\nv_max = 0.4;\nbmin = [u_min, v_min];\nbmax = [u_max, v_max];\nfp, eig, stab = fixedpoints(sys, bmin, bmax)\nstable_fp = fp[stab]","category":"page"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"using LinearAlgebra: norm\nres = 100\nu_range = range(u_min, u_max, length=res)\nv_range = range(v_min, v_max, length=res)\n\ndu(u, v) = u - u^3 - 10 * u * v^2\ndv(u, v) = -(1 + u^2) * v\nodeSol(u, v) = Point2f(du(u, v), dv(u, v))\n\nz = [norm([du(x, y), dv(x, y)]) for x in u_range, y in v_range]\nzmin, zmax = minimum(z), maximum(z)\n\nfig = Figure(size=(600, 400), fontsize=13)\nax = Axis(fig[1, 1], xlabel=\"u\", ylabel=\"v\", aspect=1.4,\n xgridcolor=:transparent, ygridcolor=:transparent,\n ylabelrotation=0)\n\nhm = heatmap!(ax, u_range, v_range, z, colormap=:Blues, colorrange=(zmin, zmax))\nColorbar(fig[1, 2], hm; label=\"\", width=15, ticksize=15, tickalign=1)\nstreamplot!(ax, odeSol, (u_min, u_max), (v_min, v_max);\n gridsize=(20, 20), arrow_size=10, stepsize=0.01,\n colormap=[:black, :black]\n)\ncolgap!(fig.layout, 7)\nlimits!(u_min, u_max, v_min, v_max)\nfig\n\n[scatter!(ax, Point(fp[i]), color=stab[i] > 0 ? :red : :dodgerblue,\n markersize=10) for i in eachindex(fp)]\nfig","category":"page"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"We can simulate a stochastic trajectory using the function trajectory.","category":"page"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"tr, ts = trajectory(sys, 1000)\n\nfig = Figure(size=(1000, 400), fontsize=13)\nax1 = Axis(fig[1, 1], xlabel=\"t\", ylabel=\"u\", aspect=1.2,\n xgridcolor=:transparent, ygridcolor=:transparent,\n ylabelrotation=0)\nax2 = Axis(fig[1, 2], xlabel=\"u\", ylabel=\"v\", aspect=1.2,\n xgridcolor=:transparent, ygridcolor=:transparent,\n ylabelrotation=0)\n\nlines!(ax1, ts, first.(tr); linewidth=2, color=:black)\n\nhm = heatmap!(ax2, u_range, v_range, z, colormap=:Blues, colorrange=(zmin, zmax))\nColorbar(fig[1, 3], hm; label=\"\", width=15, ticksize=15, tickalign=1)\nstreamplot!(ax2, odeSol, (u_min, u_max), (v_min, v_max);\n gridsize=(20, 20), arrow_size=10, stepsize=0.01,\n colormap=[:white, :white]\n)\ncolgap!(fig.layout, 7)\nlimits!(u_min, u_max, v_min, v_max)\nfig\n\n[scatter!(ax2, Point(fp[i]), color=stab[i] > 0 ? :red : :dodgerblue,\n markersize=10) for i in eachindex(fp)]\n\nlines!(ax2, reduce(hcat, tr), linewidth=1, color=(:black, 0.2))\nfig","category":"page"},{"location":"examples/gMAM_Maierstein/#Basins-of-attraction","page":"Anlyses for the Maier-Stein system","title":"Basins of attraction","text":"","category":"section"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"Basins of attraction are the regions in the state space that lead to a particular attractor. We can find the basins of attraction using the function basins.","category":"page"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"using Attractors\n\nba = basins(sys, [0.0, 0], [0.0, 1], [1.0, 0], intervals_to_box([-2, -2], [2, 2]), bstep=[0.01, 0.01], ϵ_mapper=0.001, Ttr=100)\nUr, Vr, atr, M = ba\nheatmap(Ur, Vr, M)","category":"page"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"The basin boundaries can be quickly extracted using the function basin_boundaries.","category":"page"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"bb = basinboundary(ba)","category":"page"},{"location":"examples/gMAM_Maierstein/#Transitions","page":"Anlyses for the Maier-Stein system","title":"Transitions","text":"","category":"section"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"We can quickly find a path which computes a transition from one attractor to another using the function `transition.","category":"page"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"paths_ends = (fp[stab][1], fp[stab][2])\npath, time, succes = transition(sys, paths_ends...);","category":"page"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"fig = Figure(size=(600, 400), fontsize=13)\nax = Axis(fig[1, 1], xlabel=\"u\", ylabel=\"v\", aspect=1.4,\n xgridcolor=:transparent, ygridcolor=:transparent,\n ylabelrotation=0)\n\nhm = heatmap!(ax, u_range, v_range, z, colormap=:Blues, colorrange=(zmin, zmax))\nColorbar(fig[1, 2], hm; label=\"\", width=15, ticksize=15, tickalign=1)\nstreamplot!(ax, odeSol, (u_min, u_max), (v_min, v_max);\n gridsize=(20, 20), arrow_size=10, stepsize=0.01,\n colormap=[:white, :white]\n)\ncolgap!(fig.layout, 7)\nlimits!(u_min, u_max, v_min, v_max)\nfig\n\n[scatter!(ax, Point(fp[i]), color=stab[i] > 0 ? :red : :dodgerblue,\n markersize=10) for i in eachindex(fp)]\nfig\n\nlines!(ax, path, color=:black)\nfig","category":"page"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"If we want to compute many: transitions is the function to use.","category":"page"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"tt = transitions(sys, paths_ends..., 3; tmax=1e3);","category":"page"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"fig = Figure(size=(600, 400), fontsize=13)\nax = Axis(fig[1, 1], xlabel=\"u\", ylabel=\"v\", aspect=1.4,\n xgridcolor=:transparent, ygridcolor=:transparent,\n ylabelrotation=0)\n\nhm = heatmap!(ax, u_range, v_range, z, colormap=:Blues, colorrange=(zmin, zmax))\nColorbar(fig[1, 2], hm; label=\"\", width=15, ticksize=15, tickalign=1)\nstreamplot!(ax, odeSol, (u_min, u_max), (v_min, v_max);\n gridsize=(20, 20), arrow_size=10, stepsize=0.01,\n colormap=[:black, :black]\n)\ncolgap!(fig.layout, 7)\nlimits!(u_min, u_max, v_min, v_max)\nfig\n\n[scatter!(ax, Point(fp[i]), color=stab[i] > 0 ? :red : :dodgerblue,\n markersize=10) for i in eachindex(fp)]\n\nfor i in 1:length(tt.paths)\n lines!(ax, tt.paths[i])\nend\nfig","category":"page"},{"location":"examples/gMAM_Maierstein/#Large-deviation-theory","page":"Anlyses for the Maier-Stein system","title":"Large deviation theory","text":"","category":"section"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"In the context of nonlinear dynamics, Large Deviation Theory provides tools to quantify the probability of rare events that deviate significantly from the system's typical behavior. These rare events might be extreme values of a system's output, sudden transitions between different states, or other phenomena that occur with very low probability but can have significant implications for the system's overall behavior.","category":"page"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"Large deviation theory applies principles from probability theory and statistical mechanics to develop a rigorous mathematical description of these rare events. It uses the concept of a rate function, which measures the exponential decay rate of the probability of large deviations from the mean or typical behavior. This rate function plays a crucial role in quantifying the likelihood of rare events and understanding their impact on the system.","category":"page"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"For example, in a system exhibiting chaotic behavior, LDT can help quantify the probability of sudden large shifts in the system's trajectory. Similarly, in a system with multiple stable states, it can provide insight into the likelihood and pathways of transitions between these states under fluctuations. In the context of the Minimum Action Method (MAM) and the Geometric Minimum Action Method (gMAM), Large Deviation Theory is used to handle the large deviations action functional on the space of curves. This is a key part of how these methods analyze dynamical systems.","category":"page"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"The Maier-Stein model is a typical benchmark to test such LDT techniques. Let us try to reproduce the following figure from Tobias Grafke's blog post:","category":"page"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"(Image: maier_stein)","category":"page"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"Let us first make an initial path:","category":"page"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"xx = range(-1.0, 1.0, length=100)\nyy = 0.3 .* (-xx .^ 2 .+ 1)\ninit = Matrix([xx yy]')","category":"page"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"geometric_min_action_method computes the minimizer of the Freidlin-Wentzell action using the geometric minimum action method (gMAM), to find the minimum action path (instanton) between an initial state xi and final state xf. The Minimum Action Method (MAM) is a more traditional approach, while the Geometric Minimum Action Method (gMAM) is a blend of the original MAM and the string method.","category":"page"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"gm = geometric_min_action_method(sys, init, maxiter=500; show_progress=false)\nMLP = gm.path","category":"page"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"fig = Figure(size=(600, 400), fontsize=13)\nax = Axis(fig[1, 1], xlabel=\"u\", ylabel=\"v\", aspect=1.4,\n xgridcolor=:transparent, ygridcolor=:transparent,\n ylabelrotation=0)\n\nhm = heatmap!(ax, u_range, v_range, z, colormap=:Blues, colorrange=(zmin, zmax))\nColorbar(fig[1, 2], hm; label=\"\", width=15, ticksize=15, tickalign=1)\nstreamplot!(ax, odeSol, (u_min, u_max), (v_min, v_max);\n gridsize=(20, 20), arrow_size=10, stepsize=0.01,\n colormap=[:black, :black]\n)\ncolgap!(fig.layout, 7)\nlimits!(u_min, u_max, v_min, v_max)\nfig\n\n[scatter!(ax, Point(fp[i]), color=stab[i] > 0 ? :red : :dodgerblue,\n markersize=10) for i in eachindex(fp)]\n\nlines!(ax, init, linewidth=3, color=:black, linestyle=:dash)\nlines!(ax, MLP, linewidth=3, color=:orange)\nfig","category":"page"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"","category":"page"},{"location":"examples/gMAM_Maierstein/","page":"Anlyses for the Maier-Stein system","title":"Anlyses for the Maier-Stein system","text":"Author: Orjan Ameye (orjan.ameye@hotmail.com) Date: 13 Feb 2024","category":"page"},{"location":"man/largedeviations/#Large-deviation-theory","page":"Large deviation theory","title":"Large deviation theory","text":"","category":"section"},{"location":"man/largedeviations/#Mimumum-action-path","page":"Large deviation theory","title":"Mimumum action path","text":"","category":"section"},{"location":"man/largedeviations/","page":"Large deviation theory","title":"Large deviation theory","text":"MinimumActionPath","category":"page"},{"location":"man/largedeviations/#CriticalTransitions.MinimumActionPath","page":"Large deviation theory","title":"CriticalTransitions.MinimumActionPath","text":"struct MinimumActionPath{D, T<:Real, V, Phis, Ahis, Lambda, PV, GPV}\n\nThe minimum action path between two points in a D-dimensional phase space.\n\nFields\n\npath::StateSpaceSet{D, T, V} where {D, T<:Real, V}: The path matrix.\naction::Real: The action value associated to the path.\npath_history::Any: The history of action of the paths in the optimisation algorithm (optional).\naction_history::Any: The history of action of the paths in the optimisation algorithm (optional).\nλ::Any: The Lagrange multiplier parameter for the minimum action path (optional).\ngeneralized_momentum::Any: The generalized momentum of the phase space variables (optional).\npath_velocity::Any: The path velocity (optional).\n\nConstructors\n\nMinimumActionPath(\n path,\n action;\n path_history,\n action_history,\n λ,\n generalized_momentum,\n path_velocity\n)\n\ndefined at /home/runner/work/CriticalTransitions.jl/CriticalTransitions.jl/src/largedeviations/MinimumActionPath.jl:29.\n\n\n\n\n\n","category":"type"},{"location":"man/largedeviations/#String-method","page":"Large deviation theory","title":"String method","text":"","category":"section"},{"location":"man/largedeviations/","page":"Large deviation theory","title":"Large deviation theory","text":"The string method is a technique for finding transition paths between two states in a dynamical system. The method represents the path as a \"string\" of points that connects the states and evolves it to minimize the drift along the path. The resulating tangent path is parrallel to the drift of the system, i.e., the string method computes the heteroclinic orbit. For non-gradient systems (detailed -balance is broken), the heteroclinic orbit differs from the transition path, it does correctly predict, it correctly captures the deterministic dynamics from the saddle point onward (\"downhill\" portion of the path).","category":"page"},{"location":"man/largedeviations/","page":"Large deviation theory","title":"Large deviation theory","text":"string_method","category":"page"},{"location":"man/largedeviations/#CriticalTransitions.string_method","page":"Large deviation theory","title":"CriticalTransitions.string_method","text":"string_method(\n sys::Union{Function, SgmamSystem},\n x_initial::Matrix;\n ϵ,\n iterations,\n show_progress\n) -> Any\n\n\nCompute the string method for a given system using E et al. (2007).\n\nThe string method is an iterative algorithm used to find minimum energy path (MEP) between two points in phase space. It works by evolving a discretized path (string) according to the system's drift while maintaining equal arc-length parametrization between points.\n\nThis implementation allows for computation between arbitrary points, not just stable fixed points.\n\nArguments\n\nsys::SgmamSystem: The doubled phase space system for which the string method is computed\nx_initial: Initial path discretized as a matrix where each column represents a point on the path\nϵ::Float64: Step size for the evolution step\niterations::Int64: Maximum number of iterations for path convergence\nshow_progress::Bool: Whether to display a progress meter during computation\n\nReturns\n\nx: The final converged path representing the MEP\n\n\n\n\n\nstring_method(sys::CoupledSDEs, init; kwargs...) -> Any\n\n\nCompute the string method for a given system using E et al. (2007).\n\nThe string method is an iterative algorithm used to find minimum energy path (MEP) between two points in phase space. It works by evolving a discretized path (string) according to the system's drift while maintaining equal arc-length parametrization between points.\n\nThis implementation allows for computation between arbitrary points, not just stable fixed points.\n\nArguments\n\nsys::CoupledSDEs: The system for which the string method is computed\nx_initial: Initial path discretized as a matrix where each column represents a point on the path\nϵ::Float64: Step size for the evolution step\niterations::Int64: Maximum number of iterations for path convergence\nshow_progress::Bool: Whether to display a progress meter during computation\n\nReturns\n\nx: The final converged path representing the MEP\n\n\n\n\n\n","category":"function"},{"location":"man/largedeviations/#Minimum-action-methods","page":"Large deviation theory","title":"Minimum action methods","text":"","category":"section"},{"location":"man/largedeviations/","page":"Large deviation theory","title":"Large deviation theory","text":"The minimum action method (MAM) is a technique for calculating the most probable transition path between two (meta)stable states in a stochastic dynamical system. In the limit of small noise, this path corresponds to the minimizer of an action functional. The action functional typically takes into account both the deterministic drift and the noise intensity of the system. By discretizing this path and using optimization techniques, MAM finds the trajectory that requires the least \"effort\" to transition between states in phase space.","category":"page"},{"location":"man/largedeviations/#Minimum-action-method-(MAM)","page":"Large deviation theory","title":"Minimum action method (MAM)","text":"","category":"section"},{"location":"man/largedeviations/","page":"Large deviation theory","title":"Large deviation theory","text":"Minimization of the action functioal using the optimization algorithm of Optimization.jl.","category":"page"},{"location":"man/largedeviations/","page":"Large deviation theory","title":"Large deviation theory","text":"min_action_method","category":"page"},{"location":"man/largedeviations/#CriticalTransitions.min_action_method","page":"Large deviation theory","title":"CriticalTransitions.min_action_method","text":"min_action_method(\n sys::CoupledSDEs,\n x_i,\n x_f,\n N::Int64,\n T::Real;\n kwargs...\n) -> MinimumActionPath{_A, _B, _C, Nothing, Nothing, Nothing, Nothing, Nothing} where {_A, _B<:Real, _C}\n\n\nRuns the Minimum Action Method (MAM) to find the minimum action path (instanton) between an initial state x_i and final state x_f in phase space.\n\nThis algorithm uses the minimizers of the Optimization.jl package to minimize the action functional (see fw_action) of a path for the given CoupledSDEs sys. The path is initialized as a straight line between x_i and x_f, parameterized in time via N equidistant points and total time T. Thus, the time step between discretized path points is Delta t = TN. To set an initial path different from a straight line, see the multiple dispatch method\n\nmin_action_method(sys::CoupledSDEs, init::Matrix, T::Real; kwargs...).\n\nThe minimization can be performed in blocks to save intermediate results.\n\nKeyword arguments\n\nfunctional = \"FW\": type of action functional to minimize. Defaults to fw_action, alternative: om_action.\nmaxiter = 100: maximum number of iterations before the algorithm stops.\nabstol=1e-8: absolute tolerance of action gradient to determine convergence\nreltol=1e-8: relative tolerance of action gradient to determine convergence\nmethod = Adam(): minimization algorithm (see Optimization.jl)\nverbose = true: whether to print Optimization information during the run\nshow_progress = false: whether to print a progress bar\n\nOutput\n\nIf save_iterations, returns Optim.OptimizationResults. Else, returns only the optimizer (path). If blocks > 1, a vector of results/optimizers is returned.\n\n\n\n\n\nmin_action_method(\n sys::CoupledSDEs,\n init::Matrix,\n T::Real;\n functional,\n maxiter,\n abstol,\n reltol,\n method,\n AD,\n verbose,\n show_progress\n) -> MinimumActionPath{_A, _B, _C, Nothing, Nothing, Nothing, Nothing, Nothing} where {_A, _B<:Real, _C}\n\n\nRuns the Minimum Action Method (MAM) to find the minimum action path (instanton) from an initial condition init, given a system sys and total path time T.\n\nThe initial path init must be a matrix of size (D, N), where D is the dimension of the system and N is the number of path points. The physical time of the path is specified by T, such that the time step between consecutive path points is Delta t = TN.\n\nFor more information see the main method, min_action_method(sys::CoupledSDEs, x_i, x_f, N::Int, T::Real; kwargs...).\n\n\n\n\n\n","category":"function"},{"location":"man/largedeviations/#Geometric-minimum-action-method-(gMAM)","page":"Large deviation theory","title":"Geometric minimum action method (gMAM)","text":"","category":"section"},{"location":"man/largedeviations/","page":"Large deviation theory","title":"Large deviation theory","text":"Minimization of the geometric action following Heymann and Vanden-Eijnden, PRL (2008). The gMAM reformulates MAM to avoid double optimisation of both the action and the transition time. It achieves this by using a geometric action functional that is independent of the time parametrization of the path. This reparametrization invariance makes the method more robust and computationally efficient, particularly for systems with metastable states separated by large barriers.","category":"page"},{"location":"man/largedeviations/","page":"Large deviation theory","title":"Large deviation theory","text":"geometric_min_action_method","category":"page"},{"location":"man/largedeviations/#CriticalTransitions.geometric_min_action_method","page":"Large deviation theory","title":"CriticalTransitions.geometric_min_action_method","text":"geometric_min_action_method(\n sys::CoupledSDEs,\n x_i,\n x_f;\n N,\n kwargs...\n) -> MinimumActionPath{_A, _B, _C, Nothing, Nothing, Nothing, Nothing, Nothing} where {_A, _B<:Real, _C}\n\n\nComputes the minimizer of the geometric Freidlin-Wentzell action based on the geometric minimum action method (gMAM), using optimizers of OPtimization.jl or the original formulation by Heymann and Vanden-Eijnden[1]. Only Freidlin-Wentzell action has a geometric formulation.\n\nTo set an initial path different from a straight line, see the multiple dispatch method\n\ngeometric_min_action_method(sys::CoupledSDEs, init::Matrix, arclength::Float64; kwargs...).\n\nKeyword arguments\n\nmaxiter::Int=100: maximum number of optimization iterations before the alogrithm stops\nabstol=1e-8: absolute tolerance of action gradient to determine convergence\nreltol=1e-8: relative tolerance of action gradient to determine convergence\nmethod = Adam(): minimization algorithm (see Optimization.jl)\n=0.1: step size parameter in gradient descent HeymannVandenEijnden implementation.\nverbose=false: if true, print additional output\nshow_progress=true: if true, display a progress bar\n\nOptimization algorithms\n\nThe method keyword argument takes solver methods of the Optimization.jl package; alternatively, the option solver = \"HeymannVandenEijnden\" uses the original gMAM algorithm[1].\n\n[1]: Heymann and Vanden-Eijnden, PRL (2008)\n\n\n\n\n\ngeometric_min_action_method(\n sys::CoupledSDEs,\n init::Matrix;\n maxiter,\n abstol,\n reltol,\n method,\n AD,\n ϵ,\n verbose,\n show_progress\n) -> MinimumActionPath{_A, _B, _C, Nothing, Nothing, Nothing, Nothing, Nothing} where {_A, _B<:Real, _C}\n\n\nRuns the geometric Minimum Action Method (gMAM) to find the minimum action path (instanton) from an initial condition init, given a system sys and total arc length arclength.\n\nThe initial path init must be a matrix of size (D, N), where D is the dimension of the system and N is the number of path points.\n\nFor more information see the main method, geometric_min_action_method(sys::CoupledSDEs, x_i, x_f, arclength::Float64; kwargs...).\n\n\n\n\n\n","category":"function"},{"location":"man/largedeviations/#Simple-Geometric-minimum-action-method-(sgMAM)","page":"Large deviation theory","title":"Simple Geometric minimum action method (sgMAM)","text":"","category":"section"},{"location":"man/largedeviations/","page":"Large deviation theory","title":"Large deviation theory","text":"Simplified minimization of the geometric action following Grafke et al. (2017). The simple gMAM reduces the complexity of the original gMAM by requiring only first-order derivatives of the underlying Hamiltonian optimisation formulation. This simplifies the numerical treatment and the computational complexity.","category":"page"},{"location":"man/largedeviations/","page":"Large deviation theory","title":"Large deviation theory","text":"The implementation below perform a constrained gradient descent where it assumes an autonomous system with additive Gaussian noise.","category":"page"},{"location":"man/largedeviations/","page":"Large deviation theory","title":"Large deviation theory","text":"sgmam\nSgmamSystem","category":"page"},{"location":"man/largedeviations/#CriticalTransitions.sgmam","page":"Large deviation theory","title":"CriticalTransitions.sgmam","text":"sgmam(\n sys::SgmamSystem,\n x_initial::Matrix{<:Real};\n ϵ,\n iterations,\n show_progress,\n reltol\n) -> MinimumActionPath{_A, _B, V, Nothing, Nothing, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}} where {_A, _B<:Real, V<:(AbstractVector)}\n\n\nPerforms the simplified geometric Minimal Action Method (sgMAM) on the given system sys. Our implementation is only valid for additive noise.\n\nThis method computes the optimal path in the phase space of a Hamiltonian system that minimizes the Freidlin–Wentzell action. The Hamiltonian functions H_x and H_p define the system's dynamics in a doubled phase. The initial state x_initial is evolved iteratively using constrained gradient descent with step size parameter ϵ over a specified number of iterations. The method can display a progress meter and will stop early if the relative tolerance reltol is achieved.\n\nThe function returns a tuple containing the final state, the action value, the Lagrange multipliers, the momentum, and the state derivatives. The implementation is based on the work of Grafke et al. (2019).\n\n\n\n\n\n","category":"function"},{"location":"man/largedeviations/#CriticalTransitions.SgmamSystem","page":"Large deviation theory","title":"CriticalTransitions.SgmamSystem","text":"A structure representing a system with Hamiltonian functions Hx and Hp.\n\nThis system operates in an extended phase space where the Hamiltonian is assumed to be quadratic in the extended momentum. The phase space coordinates x are doubled to form a 2n-dimensional extended phase space.\n\n\n\n\n\n","category":"type"},{"location":"man/largedeviations/#Action-functionals","page":"Large deviation theory","title":"Action functionals","text":"","category":"section"},{"location":"man/largedeviations/#Freidlin-Wentzell-action","page":"Large deviation theory","title":"Freidlin-Wentzell action","text":"","category":"section"},{"location":"man/largedeviations/","page":"Large deviation theory","title":"Large deviation theory","text":"fw_action","category":"page"},{"location":"man/largedeviations/#CriticalTransitions.fw_action","page":"Large deviation theory","title":"CriticalTransitions.fw_action","text":"fw_action(sys::CoupledSDEs, path, time) -> Any\n\n\nCalculates the Freidlin-Wentzell action of a given path with time points time in a drift field specified by the deterministic dynamics f = dynamic_rule(sys) and (normalized) noise covariance matrix covariance_matrix(sys).\n\nThe path must be a (D x N) matrix, where D is the dimensionality of the system sys and N is the number of path points. The time array must have length N.\n\nReturns a single number, which is the value of the action functional\n\nS_Tphi_t = frac12 int_0^T dot phi_t - f(phi_t) ^2_Q textdt\n\nwhere phi_t denotes the path in state space, b the drift field, and T the total time of the path. The subscript Q refers to the generalized norm a_Q^2 = langle a Q^-1 b rangle (see anorm). Here Q is the noise covariance matrix normalized by DL_1(Q), with L_1 being the L1 matrix norm.\n\n\n\n\n\n","category":"function"},{"location":"man/largedeviations/#Geometric-Freidlin-Wentzell-action","page":"Large deviation theory","title":"Geometric Freidlin-Wentzell action","text":"","category":"section"},{"location":"man/largedeviations/","page":"Large deviation theory","title":"Large deviation theory","text":"geometric_action","category":"page"},{"location":"man/largedeviations/#CriticalTransitions.geometric_action","page":"Large deviation theory","title":"CriticalTransitions.geometric_action","text":"geometric_action(sys::CoupledSDEs, path) -> Any\ngeometric_action(sys::CoupledSDEs, path, arclength) -> Any\n\n\nCalculates the geometric action of a given path with specified arclength for the drift field specified by the deterministic dynamics f = dynamic_rule(sys) and (normalized) noise covariance matrix covariance_matrix(sys).\n\nFor a given path varphi, the geometric action bar S corresponds to the minimum of the Freidlin-Wentzell action S_T(varphi) over all travel times T0, where varphi denotes the path's parameterization in physical time (see fw_action). It is given by the integral\n\nbar Svarphi = int_0^L left( varphi_Q f(varphi)_Q - langle varphi f(varphi) rangle_Q right) textds\n\nwhere s is the arclength coordinate, L the arclength, f the drift field, and the subscript Q refers to the generalized dot product langle a b rangle_Q = a^top cdot Q^-1 b (see anorm). Here Q is the noise covariance matrix normalized by DL_1(Q), with L_1 being the L1 matrix norm.\n\nReturns the value of the geometric action bar S.\n\n\n\n\n\n","category":"function"},{"location":"man/largedeviations/#Onsager-Machlup-action","page":"Large deviation theory","title":"Onsager-Machlup action","text":"","category":"section"},{"location":"man/largedeviations/","page":"Large deviation theory","title":"Large deviation theory","text":"om_action","category":"page"},{"location":"man/largedeviations/#CriticalTransitions.om_action","page":"Large deviation theory","title":"CriticalTransitions.om_action","text":"om_action(sys::CoupledSDEs, path, time, noise_strength)\n\nCalculates the Onsager-Machlup action of a given path with time points time for the drift field f = dynamic_rule(sys) at given noise_strength.\n\nThe path must be a (D x N) matrix, where D is the dimensionality of the system sys and N is the number of path points. The time array must have length N.\n\nReturns a single number, which is the value of the action functional\n\nS^sigma_Tphi_t = frac12 int_0^T left( dot phi_t - f(phi_t) ^2_Q + sigma^2 nabla cdot f right) textd t\n\nwhere phi_t denotes the path in state space, b the drift field, T the total time of the path, and sigma the noise strength. The subscript Q refers to the generalized norm a_Q^2 = langle a Q^-1 b rangle (see anorm). Here Q is the noise covariance matrix normalized by DL_1(Q), with L_1 being the L1 matrix norm.\n\n\n\n\n\n","category":"function"},{"location":"man/largedeviations/","page":"Large deviation theory","title":"Large deviation theory","text":"For convenience, a general action function is available where the type of functional is set as an argument:","category":"page"},{"location":"man/largedeviations/","page":"Large deviation theory","title":"Large deviation theory","text":"action","category":"page"},{"location":"man/largedeviations/#CriticalTransitions.action","page":"Large deviation theory","title":"CriticalTransitions.action","text":"action(\n sys::CoupledSDEs,\n path::Matrix,\n time,\n functional;\n noise_strength\n) -> Any\n\n\nComputes the action functional specified by functional for a given CoupledSDEs sys and path parameterized by time.\n\nfunctional = \"FW\": Returns the Freidlin-Wentzell action (fw_action)\nfunctional = \"OM\": Returns the Onsager-Machlup action (om_action)\n\n\n\n\n\n","category":"function"},{"location":"man/systemanalysis/#Analyzing-a-system's-stability-properties","page":"Stability analysis","title":"Analyzing a system's stability properties","text":"","category":"section"},{"location":"man/systemanalysis/","page":"Stability analysis","title":"Stability analysis","text":"To use the following functionalities, you need to load ChoasTools.jl and Attractors.jl.","category":"page"},{"location":"man/systemanalysis/#Fixed-points","page":"Stability analysis","title":"Fixed points","text":"","category":"section"},{"location":"man/systemanalysis/","page":"Stability analysis","title":"Stability analysis","text":"fixedpoints","category":"page"},{"location":"man/systemanalysis/#ChaosTools.fixedpoints","page":"Stability analysis","title":"ChaosTools.fixedpoints","text":"fixedpoints(\n sys::CoupledSDEs,\n bmin::Vector,\n bmax::Vector\n) -> Tuple{StateSpaceSet, Vector{Vector{ComplexF64}}, Vector{Bool}}\n\n\nReturns fixed points, their eigenvalues and stability of the system sys within the state space volume defined by bmin and bmax.\n\nThis is a wrapper around the fixedpoints function of DynamicalSystems.jl.\n\nInput\n\nbmin (Vector): lower limits of the state space box to be considered, as a vector of coordinates\nbmax (Vector): upper limits\nalternatively box (IntervalBox) can replace bmin and bmax\n\nExample: fixedpoints(sys, [-2,-1,0], [2,1,1]) finds the fixed points of the 3D system sys in a cube defined by the intervals [-2,2] × [-1,1] × [0,1].\n\nOutput\n\n[fp, eigs, stable]\n\nfp: StateSpaceSet of fixed points\neigs: vector of Jacobian eigenvalues of each fixed point\nstable: vector of booleans indicating the stability of each fixed point (true=stable, false=unstable)\n\nAdditional methods\n\nfixedpoints(sys::CoupedSDEs, box)\n\n\n\n\n\n","category":"function"},{"location":"man/systemanalysis/#Edge-tracking","page":"Stability analysis","title":"Edge tracking","text":"","category":"section"},{"location":"man/systemanalysis/","page":"Stability analysis","title":"Stability analysis","text":"The edge tracking algorithm is a simple numerical method to find the edge state or (possibly chaotic) saddle on the boundary between two basins of attraction. It is first introduced by Battelino et al. (1988) and further described by Skufca et al. (2006).","category":"page"},{"location":"man/systemanalysis/","page":"Stability analysis","title":"Stability analysis","text":"edgetracking\nbisect_to_edge","category":"page"},{"location":"man/systemanalysis/#Attractors.edgetracking","page":"Stability analysis","title":"Attractors.edgetracking","text":"edgetracking(\n sys::CoupledSDEs,\n u1,\n u2,\n attractors::Vector;\n eps1,\n eps2,\n converge,\n dt,\n tmax,\n ϵ_mapper,\n dt_mapper,\n solver,\n maxit,\n verbose,\n output_all,\n kwargs...\n)\n\n\nRuns the edge tracking algorithm.\n\nInput\n\nsys: dynamical system of type CoupledSDEs\nu1, u2: initial states; must belong to different basins of attraction\nattractors: vector of state vectors corresponding to the stable fixed points of sys\n\nKeyword arguments\n\neps1 = 1e-9: tolerance for bisection distance\neps2 = 1e-8: tolerance for divergence of trajectories before re-bisecting\nconverge = 1e-5: convergence criterion for M state accuracy (Euclidean distance)\ndt = 0.01: integration time step\ntmax = Inf: maximum integration time of parallel trajectories until re-bisection\nϵ_mapper = 0.1: distance threshold for AttractorMapper\ndt_mapper = 0.01: time step for AttractorMapper (keyword argument Δt)\nsolver = Vern9(): ODE solver from DifferentialEquations.jl\nmaxit = 100: maximum number of iterations before algorithm stops\nverbose = true: print status updates during run\noutput_all = false: if false, returns M state, else returns all points of the track\nkwargs...: additional keyword arguments of AttractorsViaProximity may be passed\n\nReturns\n\nIf output_all, a single state vector corresponding to the found edge state is returned. Else, a triple edge, track1, track2 is returned, where track1 and track2 are the tracks along the edge within the basin of attraction of u1 and u2, respectively; edge is the track along the edge derived from (track1 + track2)/2.\n\nwarning: Warning\nMay behave erroneously when run with solver = SimpleATsit5(), which is the default solver for AttractorsViaProximity. The recommended solver is Vern9().\n\n\n\n\n\n","category":"function"},{"location":"man/systemanalysis/#Attractors.bisect_to_edge","page":"Stability analysis","title":"Attractors.bisect_to_edge","text":"bisect_to_edge(\n sys::CoupledSDEs,\n u1,\n u2,\n attractors::Vector;\n eps1,\n ϵ_mapper,\n dt_mapper,\n solver,\n maxit,\n absto,\n relto,\n kwargs...\n)\n\n\nBisects to the basin boundary between two initial points u1 and u2. Returns the two final points, one on each side of the basin boundary, that are less than eps1 apart from each other.\n\nKeyword arguments\n\neps1=1e-9: tolerance for final distance between the two states\nϵ_mapper=0.1: ϵ parameter for AttractorMapper\ndt_mapper = 0.01: time step for AttractorMapper (keyword argument Δt)\nsolver=Vern9(): solver for AttractorMapper\nkwargs...: additional kwargs that can be passed to AttractorMapper\n\n\n\n\n\n","category":"function"},{"location":"man/systemanalysis/#Baisin-of-attraction","page":"Stability analysis","title":"Baisin of attraction","text":"","category":"section"},{"location":"man/systemanalysis/","page":"Stability analysis","title":"Stability analysis","text":"basinboundary\nbasboundary\nbasins","category":"page"},{"location":"man/systemanalysis/#CriticalTransitions.basinboundary","page":"Stability analysis","title":"CriticalTransitions.basinboundary","text":"basinboundary(X, Y, h; coords, A, B, C) -> Any\n\n\nComputes the basin boundary for given output X, Y, h of the basins function.\n\nTo be further documented.\n\n\n\n\n\nbasinboundary(boa) -> Any\n\n\nComputes the basin boundary for given output boa of the basins function.\n\nTo be further documented.\n\n\n\n\n\n","category":"function"},{"location":"man/systemanalysis/#CriticalTransitions.basboundary","page":"Stability analysis","title":"CriticalTransitions.basboundary","text":"basboundary(\n sys::CoupledSDEs,\n xrange::Vector,\n yrange::Vector,\n xspacing::Float64,\n attractors::Vector;\n eps1,\n ϵ_mapper,\n dt_mapper,\n solver,\n maxit\n) -> Tuple{Any, Any}\n\n\nThis function computes the basin boundary.\n\n\n\n\n\n","category":"function"},{"location":"man/systemanalysis/#CriticalTransitions.basins","page":"Stability analysis","title":"CriticalTransitions.basins","text":"basins(\n sys::CoupledSDEs,\n A,\n B,\n C,\n H;\n bstep,\n pstep,\n ϵ_mapper,\n show_progress,\n kwargs...\n) -> Any\n\n\nComputes the basins of attraction of CoupledSDEs sys on a plane spanned by the distinct points A, B, C and limited by the box H. Uses the AttractorsViaProximity function from DynamicalSystems.jl to compute the basins of attraction.\n\nA, B, C are elements of mathbbR^d (where d is the dimension of the sys) and H is a hyperrectangle in mathbbR^d.\n\nThe plane is given by P_UV = A+u(B-A)+v(C-A)inmathbbR^d u in U vin V for some closed and bounded real intervals U and V which are selected such that both\n\nP_UV subseteq H, and\nUtimes VsubseteqmathbbR^2 has maximal area,\n\ni.e. P_UV is the \"largest\" possible plane contained within H. This plane is determined behind the scenes.\n\nThis function returns a four-dimensional vector. The first two entries are discretised versions of the interval U and V (as defined above, of lengths ell_Uell_V respectively); the third entry is a dictionary of the attractors (stable equilibria) of the system within H, and the final entry is an ell_Vtimesell_U matrix of integers that group the initial conditions (written in terms of A+u(B-A)+v(C-A) where uin U and vin V) by which attractor they will in time converge to.\n\nKeyword arguments\n\nbstep = [0.01, 0.01]: a vector of length two whose elements respectively specify the length of the incremental steps taken across each dimension in the discretisation of your plane\npstep = [0.1, 0.1]: a vector of length two whose elements give the increments of the mesh that the maximisation process of finding a plane from a box is taken over (for more information see the source code of the function plane in the src/systemanalysis/planeofbox.jl file)\nϵ_mapper = 0.01: ϵ parameter of AttractorsViaProximity\nkwargs...: keyword arguments passed to the AttractorsViaProximity function (namely, Ttr, Δt, horizon_limit, mx_chk_lost)\nshow_progress = false: if true, displays a progress bar\n\n\n\n\n\n","category":"function"},{"location":"examples/sgMAM_KPO/","page":"sgMAM for the Kerr Parametric Oscillator","title":"sgMAM for the Kerr Parametric Oscillator","text":"We demonstrate the simple geometric minimum action method (sgMAM) on the Kerr parametric oscillator (KPO). The method computes the optimal path between two attractors in the phase space that minimizes the action of the system. It is a simplification of the geometric minimum action method (gMAM) by avoiding the computation of the second order derivatives of the extended Hamiltonian of the optimisation problem.","category":"page"},{"location":"examples/sgMAM_KPO/","page":"sgMAM for the Kerr Parametric Oscillator","title":"sgMAM for the Kerr Parametric Oscillator","text":"using CriticalTransitions, CairoMakie","category":"page"},{"location":"examples/sgMAM_KPO/","page":"sgMAM for the Kerr Parametric Oscillator","title":"sgMAM for the Kerr Parametric Oscillator","text":"The KPO equation is a nonlinear ordinary differential equation that describes the response of the nonlinear parametrically driven resonator at its dominant resonant condition. The equation of motion are of the form:","category":"page"},{"location":"examples/sgMAM_KPO/","page":"sgMAM for the Kerr Parametric Oscillator","title":"sgMAM for the Kerr Parametric Oscillator","text":"dotmathbfx = mathbff(mathbfx) + sigmamathbfξ(t)","category":"page"},{"location":"examples/sgMAM_KPO/","page":"sgMAM for the Kerr Parametric Oscillator","title":"sgMAM for the Kerr Parametric Oscillator","text":"where f is an autonemous drift function and and we have brownian noise ξ with intensity σ.","category":"page"},{"location":"examples/sgMAM_KPO/","page":"sgMAM for the Kerr Parametric Oscillator","title":"sgMAM for the Kerr Parametric Oscillator","text":"Here we define the define the drift of each seperable variable u and v. In addition, we hard-code the Jacobian of the drift function. ","category":"page"},{"location":"examples/sgMAM_KPO/","page":"sgMAM for the Kerr Parametric Oscillator","title":"sgMAM for the Kerr Parametric Oscillator","text":"const λ = 3 / 1.21 * 2 / 295\nconst ω0 = 1.000\nconst ω = 1.000\nconst γ = 1 / 295\nconst η = 0\nconst α = -1\n\nfunction fu(u, v)\n return (-4 * γ * ω * u - 2 * λ * v - 4 * (ω0 - ω^2) * v - 3 * α * v * (u^2 + v^2)) /\n (8 * ω)\nend\nfunction fv(u, v)\n return (-4 * γ * ω * v - 2 * λ * u + 4 * (ω0 - ω^2) * u + 3 * α * u * (u^2 + v^2)) /\n (8 * ω)\nend\nstream(u,v) = Point2f(fu(u, v), fv(u, v))\ndfvdv(u, v) = (-4 * γ * ω + 6 * α * u * v) / (8 * ω)\ndfudu(u, v) = (-4 * γ * ω - 6 * α * u * v) / (8 * ω)\ndfvdu(u, v) = (-2 * λ + 4 * (ω0 - ω^2) + 9 * α * u^2 + 3 * α * v^2) / (8 * ω)\ndfudv(u, v) = (-2 * λ - 4 * (ω0 - ω^2) - 3 * α * u^2 - 9 * α * v^2) / (8 * ω)","category":"page"},{"location":"examples/sgMAM_KPO/","page":"sgMAM for the Kerr Parametric Oscillator","title":"sgMAM for the Kerr Parametric Oscillator","text":"The optimisation is performed in a doubled phase space, i.e., every variable of the SDE system is considered as a generelised coordinate mathbfx and gets a corresponding generalised momentum mathbfp. The makes it that also systems with dissipative flow can be solved. As such, we extend the phase space by defining the hamiltionian","category":"page"},{"location":"examples/sgMAM_KPO/","page":"sgMAM for the Kerr Parametric Oscillator","title":"sgMAM for the Kerr Parametric Oscillator","text":"H = sum_i fracp_i^22 + f_i(mathbfx)p_i","category":"page"},{"location":"examples/sgMAM_KPO/","page":"sgMAM for the Kerr Parametric Oscillator","title":"sgMAM for the Kerr Parametric Oscillator","text":"Hence, to use the sgMAM method, we need to define the derivatives of the Hamiltonian with respect to the phase space coordinates and the generalised momentum: ","category":"page"},{"location":"examples/sgMAM_KPO/","page":"sgMAM for the Kerr Parametric Oscillator","title":"sgMAM for the Kerr Parametric Oscillator","text":"function H_x(x, p) # ℜ² → ℜ²\n u, v = eachrow(x)\n pu, pv = eachrow(p)\n\n H_u = @. pu * dfudu(u, v) + pv * dfvdu(u, v)\n H_v = @. pu * dfudv(u, v) + pv * dfvdv(u, v)\n return Matrix([H_u H_v]')\nend\nfunction H_p(x, p) # ℜ² → ℜ²\n u, v = eachrow(x)\n pu, pv = eachrow(p)\n\n H_pu = @. pu + fu(u, v)\n H_pv = @. pv + fv(u, v)\n return Matrix([H_pu H_pv]')\nend\n\nsys = SgmamSystem{false, 2}(H_x, H_p)","category":"page"},{"location":"examples/sgMAM_KPO/","page":"sgMAM for the Kerr Parametric Oscillator","title":"sgMAM for the Kerr Parametric Oscillator","text":"We saved this function in the SgmamSystem struct. We want to find the optimal path between two attractors in the phase space. We define the initial trajectory as wiggle between the two attractors.","category":"page"},{"location":"examples/sgMAM_KPO/","page":"sgMAM for the Kerr Parametric Oscillator","title":"sgMAM for the Kerr Parametric Oscillator","text":"# setup\nNt = 500 # number of discrete time steps\ns = collect(range(0; stop=1, length=Nt))\n\nxa = [-0.0208, 0.0991]\nxb = -xa\nxsaddle = [0.0, 0.0]\n\n# Initial trajectory\nxx = @. (xb[1] - xa[1]) * s + xa[1] + 4 * s * (1 - s) * xsaddle[1]\nyy = @. (xb[2] - xa[2]) * s + xa[2] + 4 * s * (1 - s) * xsaddle[2] + 0.01 * sin(2π * s)\nx_initial = Matrix([xx yy]')","category":"page"},{"location":"examples/sgMAM_KPO/","page":"sgMAM for the Kerr Parametric Oscillator","title":"sgMAM for the Kerr Parametric Oscillator","text":"The optimisation is the performed by the sgmam function:","category":"page"},{"location":"examples/sgMAM_KPO/","page":"sgMAM for the Kerr Parametric Oscillator","title":"sgMAM for the Kerr Parametric Oscillator","text":"MLP = sgmam(\n sys, x_initial; iterations=1_000, ϵ=10e2, show_progress=false\n)\nx_min = MLP.path;","category":"page"},{"location":"examples/sgMAM_KPO/","page":"sgMAM for the Kerr Parametric Oscillator","title":"sgMAM for the Kerr Parametric Oscillator","text":"The function returns the optimal path x_min, the minmal action S_min, the Lagrange multipliers lambda associated with the optimal path, the optimal generalised momentum p, and the time derivative of the optimal path xdot. We can plot the initial trajectory and the optimal path:","category":"page"},{"location":"examples/sgMAM_KPO/","page":"sgMAM for the Kerr Parametric Oscillator","title":"sgMAM for the Kerr Parametric Oscillator","text":"fig, ax, _ = lines(x_initial[1, :], x_initial[2, :]; label=\"init\", linewidth=3, color=:black)\nlines!(x_min[1, :], x_min[2, :]; label=\"MLP\", linewidth=3, color=:red)\nstreamplot!(ax, stream, (-0.08, 0.08), (-0.15, 0.15);\n gridsize=(20, 20), arrow_size=10, stepsize=0.001,\n colormap=[:gray, :gray]\n)\naxislegend(ax)\nfig","category":"page"},{"location":"man/sampling/#Sampling-transitions","page":"Sampling transitions","title":"Sampling transitions","text":"","category":"section"},{"location":"man/sampling/","page":"Sampling transitions","title":"Sampling transitions","text":"TransitionEnsemble","category":"page"},{"location":"man/sampling/#CriticalTransitions.TransitionEnsemble","page":"Sampling transitions","title":"CriticalTransitions.TransitionEnsemble","text":"struct TransitionEnsemble{SSS, T, Tstat, ES}\n\nEnsemble of transition paths between two points in a state space.\n\nFields\n\npaths::Vector\ntimes::Vector\nsuccess_rate::Any\nresidence_time::Any\ntransition_time::Any\nsciml_ensemble::Any\n\nConstructors\n\nTransitionEnsemble(\n paths,\n times,\n success_rate,\n residence_time,\n transition_time,\n sciml_ensemble\n)\n\ndefined at /home/runner/work/CriticalTransitions.jl/CriticalTransitions.jl/src/trajectories/transition.jl:14.\n\n\n\n\n\n","category":"type"},{"location":"man/sampling/#...-by-direct-simulation","page":"Sampling transitions","title":"... by direct simulation","text":"","category":"section"},{"location":"man/sampling/","page":"Sampling transitions","title":"Sampling transitions","text":"These functions generate noise-induced transitions between an initial and final state.","category":"page"},{"location":"man/sampling/","page":"Sampling transitions","title":"Sampling transitions","text":"transition\ntransitions","category":"page"},{"location":"man/sampling/#CriticalTransitions.transition","page":"Sampling transitions","title":"CriticalTransitions.transition","text":"transition(\n sys::CoupledSDEs,\n x_i,\n x_f;\n radii,\n tmax,\n radius_directions,\n cut_start,\n kwargs...\n) -> Tuple{Any, Any, Bool}\n\n\nGenerates a sample transition from point x_i to point x_f.\n\nThis function simulates sys in time, starting from initial condition x_i, until entering a ball of given radius around x_f.\n\nKeyword arguments\n\nradii=(0.1, 0.1): radius of the ball around x_i and x_f, respectively\ntmax=1e3: maximum time until the simulation stops even x_f has not been reached\nradius_directions=1:length(sys.u): the directions in phase space to consider when calculating the radii rad_i and rad_f. Defaults to all directions. To consider only a subspace of state space, insert a vector of indices of the dimensions to be included.\ncut_start=true: if false, returns the whole trajectory up to the transition\nkwargs...: keyword arguments passed to CommonSolve.solve\n\nOutput\n\n[path, times, success]\n\npath (Matrix): transition path (size [dim × N], where N is the number of time points)\ntimes (Vector): time values (since start of simulation) of the path points (size N)\nsuccess (bool): if true, a transition occured (i.e. the ball around x_f has been reached), else false\n\nSee also transitions, trajectory.\n\n\n\n\n\n","category":"function"},{"location":"man/sampling/#CriticalTransitions.transitions","page":"Sampling transitions","title":"CriticalTransitions.transitions","text":"function transitions(sys::CoupledSDEs, x_i, x_f, N=1; kwargs...)\n\nGenerates an ensemble of N transition samples of sys from point x_i to point x_f.\n\nThis function repeatedly calls the transition function to efficiently generate an ensemble of transitions. Multi-threading is enabled.\n\nKeyword arguments\n\nradii=(0.1, 0.1): radius of the ball around x_i and x_f, respectively\nNmax: number of attempts before the algorithm stops even if less than N transitions occurred.\ntmax=1e3: maximum time when the simulation stops even x_f has not been reached\nradius_directions=1:length(sys.u): the directions in phase space to consider when calculating the radii rad_i and rad_f. Defaults to all directions. To consider only a subspace of state space, insert a vector of indices of the dimensions to be included.\ncut_start=true: if false, returns the whole trajectory up to the transition\nshow_progress=true: shows a progress bar with respect to Nmax\nkwargs...: keyword arguments passed to CommonSolve.solve\n\nSee also transition.\n\nReturns a TransitionEnsemble object.\n\n\n\n\n\n","category":"function"},{"location":"man/sampling/#...-in-pathspace","page":"Sampling transitions","title":"... in pathspace","text":"","category":"section"},{"location":"man/sampling/","page":"Sampling transitions","title":"Sampling transitions","text":"Coming soon...","category":"page"},{"location":"quickstart/#Quickstart","page":"Quickstart","title":"Quickstart","text":"","category":"section"},{"location":"quickstart/#Installation","page":"Quickstart","title":"Installation","text":"","category":"section"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"As this module is not published yet, there are two ways to access it:","category":"page"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"Option 1 (recommended): Install from GitHub\nEnter the Julia package manager by typing ] in the REPL: julia> ]\ntype add https://github.com/juliadynamics/CriticalTransitions.jl.git\nOption 2: Load module locally\nClone the repo: git clone https://github.com/juliadynamics/CriticalTransitions.jl.git\nIn Julia, include the module file: include(\"PATH/src/CriticalTransitions.jl\"), where PATH is the relative path to the repo you just cloned\nLoad the module: using .CriticalTransitions","category":"page"},{"location":"quickstart/#Basic-usage","page":"Quickstart","title":"Basic usage","text":"","category":"section"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"The general workflow of CriticalTransitions essentially follows two steps:","category":"page"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"Define your system (see Define a CoupledSDEs system)\nInvestigate the system by calling methods (see Methods)","category":"page"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"info: New system type: RateSystem\nWe are planning to introduce the the struct RateSystem along CoupledSDEs. In a RateSystem, the time dependence of parameters can conveniently be specified, laying the foundation for a toolbox to study rate-induced tipping, or R-tipping.","category":"page"},{"location":"quickstart/#Methods","page":"Quickstart","title":"Methods","text":"","category":"section"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"Currently the following functions are implemented to analyze a CoupledSDEs and corresponding sample transition paths.","category":"page"},{"location":"quickstart/","page":"Quickstart","title":"Quickstart","text":"Pages = [\"man/simulation.md\", \"man/sampling.md\", \"man/largedeviations.md\"]","category":"page"},{"location":"#CriticalTransitions.jl","page":"Home","title":"CriticalTransitions.jl","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"A Julia package for the numerical investigation of noise- and rate-induced transitions in dynamical systems.","category":"page"},{"location":"","page":"Home","title":"Home","text":"Building on DynamicalSystems.jl and DifferentialEquations.jl, this package aims to provide a toolbox for dynamical systems under time-dependent forcing, with a focus on tipping phenomena and metastability.","category":"page"},{"location":"","page":"Home","title":"Home","text":"(Image: CT.jl infographic)","category":"page"},{"location":"","page":"Home","title":"Home","text":"info: Current features\nStochastic simulation made easy: Gaussian noise, uncorrelated and correlated, additive and multiplicative\nTransition path sampling: Parallelized ensemble rejection sampling\nLarge deviation theory tools: Action functionals and minimization algorithms (MAM, gMAM)","category":"page"},{"location":"","page":"Home","title":"Home","text":"ukw: Planned features\nRare event simulation: importance sampling, AMS\nQuasipotentials: Ordered line integral method (OLIM)\nRate-induced tipping tools\nSymbolic differentiation of action functionals\n...?","category":"page"},{"location":"","page":"Home","title":"Home","text":"Developers: Reyk Börner, Ryan Deeley, Raphael Römer and Orjan Ameye","category":"page"},{"location":"","page":"Home","title":"Home","text":"Thanks to Jeroen Wouters, Calvin Nesbitt, Tobias Grafke, George Datseris and Oliver Mehling","category":"page"},{"location":"","page":"Home","title":"Home","text":"This work is part of the CriticalEarth project.","category":"page"},{"location":"man/CoupledSDEs/#Define-a-CoupledSDEs-system","page":"Define a CoupledSDEs system","title":"Define a CoupledSDEs system","text":"","category":"section"},{"location":"man/CoupledSDEs/","page":"Define a CoupledSDEs system","title":"Define a CoupledSDEs system","text":"CoupledSDEs","category":"page"},{"location":"man/CoupledSDEs/#DynamicalSystemsBase.CoupledSDEs","page":"Define a CoupledSDEs system","title":"DynamicalSystemsBase.CoupledSDEs","text":"CoupledSDEs(ds::CoupledODEs, p; kwargs...)\n\nConverts a CoupledODEs system into a CoupledSDEs.\n\n\n\n\n\n","category":"type"},{"location":"man/CoupledSDEs/","page":"Define a CoupledSDEs system","title":"Define a CoupledSDEs system","text":"info: Info\nNote that nonlinear mixings of the Noise Process mathcalW fall into the class of random ordinary differential equations (RODEs) which have a separate set of solvers. See this example of DifferentialEquations.jl.","category":"page"},{"location":"man/CoupledSDEs/#defining-stochastic-dynamics","page":"Define a CoupledSDEs system","title":"Examples: Defining stochastic dynamics","text":"","category":"section"},{"location":"man/CoupledSDEs/","page":"Define a CoupledSDEs system","title":"Define a CoupledSDEs system","text":"Let's look at some examples of the different types of stochastic systems that can be defined.","category":"page"},{"location":"man/CoupledSDEs/","page":"Define a CoupledSDEs system","title":"Define a CoupledSDEs system","text":"For simplicity, we choose a slow exponential growth in 2 dimensions as the deterministic dynamics f:","category":"page"},{"location":"man/CoupledSDEs/","page":"Define a CoupledSDEs system","title":"Define a CoupledSDEs system","text":"using DynamicalSystemsBase, StochasticDiffEq, DiffEqNoiseProcess\nusing CairoMakie\nimport Random # hide\nRandom.seed!(10) # hide\nf!(du, u, p, t) = du .= 1.01u # deterministic part\n\nfunction plot_trajectory(Y, t)\n fig = Figure()\n ax = Axis(fig[1,1]; xlabel = \"time\", ylabel = \"variable\")\n for var in columns(Y)\n lines!(ax, t, var)\n end\n fig\nend;","category":"page"},{"location":"man/CoupledSDEs/#Additive-noise","page":"Define a CoupledSDEs system","title":"Additive noise","text":"","category":"section"},{"location":"man/CoupledSDEs/","page":"Define a CoupledSDEs system","title":"Define a CoupledSDEs system","text":"When g(u p t) is independent of the state u, the noise is called additive; otherwise, it is multiplicative. We can define a simple additive noise system as follows:","category":"page"},{"location":"man/CoupledSDEs/","page":"Define a CoupledSDEs system","title":"Define a CoupledSDEs system","text":"sde = CoupledSDEs(f!, zeros(2));","category":"page"},{"location":"man/CoupledSDEs/","page":"Define a CoupledSDEs system","title":"Define a CoupledSDEs system","text":"which is equivalent to","category":"page"},{"location":"man/CoupledSDEs/","page":"Define a CoupledSDEs system","title":"Define a CoupledSDEs system","text":"t0 = 0.0; W0 = zeros(2);\nW = WienerProcess(t0, W0, 0.0)\nsde = CoupledSDEs(f!, zeros(2);\n noise_process=W, covariance=[1 0; 0 1], noise_strength=1.0\n );","category":"page"},{"location":"man/CoupledSDEs/","page":"Define a CoupledSDEs system","title":"Define a CoupledSDEs system","text":"We defined a Wiener process W, whose increments are vectors of normally distributed random numbers of length matching the output of g. The noise is applied element-wise, i.e., g.*dW. Since the noise processes are uncorrelated, meaning the covariance matrix is diagonal, this type of noise is referred to as diagonal.","category":"page"},{"location":"man/CoupledSDEs/","page":"Define a CoupledSDEs system","title":"Define a CoupledSDEs system","text":"We can sample a trajectory from this system using the trajectory function also used for the deterministic systems:","category":"page"},{"location":"man/CoupledSDEs/","page":"Define a CoupledSDEs system","title":"Define a CoupledSDEs system","text":"tr = trajectory(sde, 1.0)\nplot_trajectory(tr...)","category":"page"},{"location":"man/CoupledSDEs/#Correlated-noise","page":"Define a CoupledSDEs system","title":"Correlated noise","text":"","category":"section"},{"location":"man/CoupledSDEs/","page":"Define a CoupledSDEs system","title":"Define a CoupledSDEs system","text":"In the case of correlated noise, the random numbers in a vector increment dW are correlated. This can be achieved by specifying the covariance matrix Sigma via the covariance keyword:","category":"page"},{"location":"man/CoupledSDEs/","page":"Define a CoupledSDEs system","title":"Define a CoupledSDEs system","text":"ρ = 0.3\nΣ = [1 ρ; ρ 1]\ndiffeq = (alg = LambaEM(), dt=0.1)\nsde = CoupledSDEs(f!, zeros(2); covariance=Σ, diffeq=diffeq)","category":"page"},{"location":"man/CoupledSDEs/","page":"Define a CoupledSDEs system","title":"Define a CoupledSDEs system","text":"Alternatively, we can parametrise the covariance matrix by defining the diffusion function g ourselves:","category":"page"},{"location":"man/CoupledSDEs/","page":"Define a CoupledSDEs system","title":"Define a CoupledSDEs system","text":"g!(du, u, p, t) = (du .= [1 p[1]; p[1] 1]; return nothing) \nsde = CoupledSDEs(f!, zeros(2), (ρ); g=g!, noise_prototype=zeros(2, 2))","category":"page"},{"location":"man/CoupledSDEs/","page":"Define a CoupledSDEs system","title":"Define a CoupledSDEs system","text":"Here, we had to provide noise_prototype to indicate that the diffusion function g will output a 2x2 matrix.","category":"page"},{"location":"man/CoupledSDEs/#Scalar-noise","page":"Define a CoupledSDEs system","title":"Scalar noise","text":"","category":"section"},{"location":"man/CoupledSDEs/","page":"Define a CoupledSDEs system","title":"Define a CoupledSDEs system","text":"If all state variables are forced by the same single random variable, we have scalar noise. To define scalar noise, one has to give an one-dimensional noise process to the noise_process keyword of the CoupledSDEs constructor. ","category":"page"},{"location":"man/CoupledSDEs/","page":"Define a CoupledSDEs system","title":"Define a CoupledSDEs system","text":"t0 = 0.0; W0 = 0.0;\nnoise = WienerProcess(t0, W0, 0.0)\nsde = CoupledSDEs(f!, rand(2)/10; noise_process=noise)\n\ntr = trajectory(sde, 1.0)\nplot_trajectory(tr...)","category":"page"},{"location":"man/CoupledSDEs/","page":"Define a CoupledSDEs system","title":"Define a CoupledSDEs system","text":"We can see that noise applied to each variable is the same.","category":"page"},{"location":"man/CoupledSDEs/#Multiplicative-and-time-dependent-noise","page":"Define a CoupledSDEs system","title":"Multiplicative and time-dependent noise","text":"","category":"section"},{"location":"man/CoupledSDEs/","page":"Define a CoupledSDEs system","title":"Define a CoupledSDEs system","text":"In the SciML ecosystem, multiplicative noise is defined through the condition g_i(t u)=a_i u. However, in the literature the name is more broadly used for any situation where the noise is non-additive and depends on the state u, possibly also in a non-linear way. When defining a CoupledSDEs, we can make the noise term time- and state-dependent by specifying an explicit time- or state-dependence in the noise function g, just like we would define f. For example, we can define a system with temporally decreasing multiplicative noise as follows:","category":"page"},{"location":"man/CoupledSDEs/","page":"Define a CoupledSDEs system","title":"Define a CoupledSDEs system","text":"function g!(du, u, p, t)\n du .= u ./ (1+t)\n return nothing\nend\nsde = CoupledSDEs(f!, rand(2)./10; g=g!)","category":"page"},{"location":"man/CoupledSDEs/#Non-diagonal-noise","page":"Define a CoupledSDEs system","title":"Non-diagonal noise","text":"","category":"section"},{"location":"man/CoupledSDEs/","page":"Define a CoupledSDEs system","title":"Define a CoupledSDEs system","text":"Non-diagonal noise allows for the terms to be linearly mixed (correlated) via g being a matrix. Suppose we have two Wiener processes and two state variables such that the output of g is a 2x2 matrix. Therefore, we have","category":"page"},{"location":"man/CoupledSDEs/","page":"Define a CoupledSDEs system","title":"Define a CoupledSDEs system","text":"du_1 = f_1(upt)dt + g_11(upt)dW_1 + g_12(upt)dW_2 \ndu_2 = f_2(upt)dt + g_21(upt)dW_1 + g_22(upt)dW_2","category":"page"},{"location":"man/CoupledSDEs/","page":"Define a CoupledSDEs system","title":"Define a CoupledSDEs system","text":"To indicate the structure that g should have, we must use the noise_prototype keyword. Let us define a special type of non-diagonal noise called commutative noise. For this we can utilize the RKMilCommute algorithm which is designed to utilize the structure of commutative noise.","category":"page"},{"location":"man/CoupledSDEs/","page":"Define a CoupledSDEs system","title":"Define a CoupledSDEs system","text":"σ = 0.25 # noise strength\nfunction g!(du, u, p, t)\n du[1,1] = σ*u[1]\n du[2,1] = σ*u[2]\n du[1,2] = σ*u[1]\n du[2,2] = σ*u[2]\n return nothing\nend\ndiffeq = (alg = RKMilCommute(), reltol = 1e-3, abstol = 1e-3, dt=0.1)\nsde = CoupledSDEs(f!, rand(2)./10; g=g!, noise_prototype = zeros(2, 2), diffeq = diffeq)","category":"page"},{"location":"man/CoupledSDEs/","page":"Define a CoupledSDEs system","title":"Define a CoupledSDEs system","text":"warning: Warning\nNon-diagonal problems need specific solvers. See the SciML recommendations.","category":"page"},{"location":"man/CoupledSDEs/#Interface-to-DynamicalSystems.jl","page":"Define a CoupledSDEs system","title":"Interface to DynamicalSystems.jl","text":"","category":"section"},{"location":"man/CoupledSDEs/#Converting-between-CoupledSDEs-and-CoupledODEs","page":"Define a CoupledSDEs system","title":"Converting between CoupledSDEs and CoupledODEs","text":"","category":"section"},{"location":"man/CoupledSDEs/","page":"Define a CoupledSDEs system","title":"Define a CoupledSDEs system","text":"tip: Analyzing deterministic dynamics with DynamicalSystems.jl\nThe deterministic part of a CoupledSDEs system can easily be extracted as a CoupledODEs, a common subtype of a ContinuousTimeDynamicalSystem in DynamicalSystems.jl.","category":"page"},{"location":"man/CoupledSDEs/","page":"Define a CoupledSDEs system","title":"Define a CoupledSDEs system","text":"CoupledODEs(sde::CoupledSDEs) extracts the deterministic part of sde as a CoupledODEs.\nCoupledSDEs(ode::CoupledODEs; kwargs) turns ode into a CoupledSDEs.","category":"page"},{"location":"man/CoupledSDEs/","page":"Define a CoupledSDEs system","title":"Define a CoupledSDEs system","text":"CoupledODEs","category":"page"},{"location":"man/CoupledSDEs/#DynamicalSystemsBase.CoupledODEs","page":"Define a CoupledSDEs system","title":"DynamicalSystemsBase.CoupledODEs","text":"CoupledODEs(ds::CoupledSDEs; kwargs...)\n\nConverts a CoupledSDEs into a CoupledODEs by extracting the deterministic part of ds.\n\n\n\n\n\n","category":"type"},{"location":"man/CoupledSDEs/","page":"Define a CoupledSDEs system","title":"Define a CoupledSDEs system","text":"For example, the Lyapunov spectrum of a CoupledSDEs in the absence of noise, here exemplified by the FitzHugh-Nagumo model, can be computed by typing:","category":"page"},{"location":"man/CoupledSDEs/","page":"Define a CoupledSDEs system","title":"Define a CoupledSDEs system","text":"using CriticalTransitions\nusing DynamicalSystems: lyapunovspectrum\n\nfunction fitzhugh_nagumo(u, p, t)\n x, y = u\n ϵ, β, α, γ, κ, I = p\n\n dx = (-α * x^3 + γ * x - κ * y + I) / ϵ\n dy = -β * y + x\n\n return SA[dx, dy]\nend\n\np = [1.,3.,1.,1.,1.,0.]\n\nsys = CoupledSDEs(fitzhugh_nagumo, zeros(2), p; noise_strength=0.1)\nls = lyapunovspectrum(CoupledODEs(sys), 10000)","category":"page"}] }