Skip to content

SDE Examples: minor rewording of explanation of "diagonal vs non-diagonal noise" #779

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 4 commits into from
May 5, 2025
Merged

SDE Examples: minor rewording of explanation of "diagonal vs non-diagonal noise" #779

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
0bf17df
Minor rewording of explanation of "diagonal vs non-diagonal noise" in…
hurak May 3, 2025
72675bc
Update docs/src/tutorials/sde_example.md
ChrisRackauckas May 4, 2025
a291a26
Update docs/src/tutorials/sde_example.md
ChrisRackauckas May 4, 2025
0de11d4
Update make.jl
ChrisRackauckas May 5, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 3 additions & 1 deletion docs/make.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -39,7 +39,9 @@ makedocs(
"https://www.wolframalpha.com/input/?i=u%27%3D-sqrt%28u%29",
"https://www.mathworks.com/help/simulink/gui/absolutetolerance.html",
"https://www.mathworks.com/help/matlab/math/choose-an-ode-solver.html",
"https://journals.ametsoc.org/view/journals/atsc/20/2/1520-0469_1963_020_0130_dnf_2_0_co_2.xml"
"https://journals.ametsoc.org/view/journals/atsc/20/2/1520-0469_1963_020_0130_dnf_2_0_co_2.xml",
"https://mathematica.stackexchange.com/questions/40122/help-to-plot-poincar%C3%A9-section-for-double-pendulum",
"https://github.com/SciML/DiffEqProblemLibrary.jl/blob/master/lib/SDEProblemLibrary/src/SDEProblemLibrary.jl",
],
doctest = false, clean = true,
warnonly = [:missing_docs],
Expand Down
64 changes: 30 additions & 34 deletions docs/src/tutorials/sde_example.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -142,49 +142,45 @@ In general, a system of SDEs
du = f(u,p,t)dt + g(u,p,t)dW,
```

where `g` is now a matrix of values, is numerically integrated in the
same way as ODEs. A common scenario is when we have diagonal noise, which
is the default for DifferentialEquations.jl. Physically this means that
every variable in the system gets a different random kick. Consequently, `g` is a
diagonal matrix and we can handle this in a simple manner by defining
the deterministic part `f(du,u,p,t)` and the stochastic part
`g(du2,u,p,t)` as in-place functions.

Consider for example a stochastic variant of the Lorenz equations, where we introduce a
simple additive noise to each of `x,y,z`, which is simply `3*N(0,dt)`. Here, `N` is the normal
distribution and `dt` is the time step. This is done as follows:
where `u` is now a vector of variables, `f` is a vector, and `g` is a matrix, is numerically integrated in the same way as ODEs. A common scenario, which is the default for DifferentialEquations.jl, is that every variable in the system gets a different random kick. This is the case when `g` is a diagonal matrix. Correspondingly, we say that we have a diagonal noise.

We handle this in a simple manner by defining the deterministic part `f!(du,u,p,t)` and the stochastic part
`g!(du2,u,p,t)` as in-place functions, but note that our convention is that the function `g!` only defines and modifies the diagonal entries of `g` matrix.

As an example, we consider a stochastic variant of the Lorenz equations, where we add to each of `u₁, u₂, u₃` their own simple noise `3*N(0,dt)`. Here, `N` is the normal distribution and `dt` is the time step. This is done as follows:

```@example sde2
using DifferentialEquations
using Plots
function lorenz!(du, u, p, t)

function f!(du, u, p, t)
du[1] = 10.0 * (u[2] - u[1])
du[2] = u[1] * (28.0 - u[3]) - u[2]
du[3] = u[1] * u[2] - (8 / 3) * u[3]
end

function σ_lorenz!(du, u, p, t)
function g!(du, u, p, t) # It actually represents a diagonal matrix [3.0 0 0; 0 3.0 0; 0 0 3.0]
du[1] = 3.0
du[2] = 3.0
du[3] = 3.0
end

prob_sde_lorenz = SDEProblem(lorenz!, σ_lorenz!, [1.0, 0.0, 0.0], (0.0, 10.0))
prob_sde_lorenz = SDEProblem(f!, g!, [1.0, 0.0, 0.0], (0.0, 10.0))
sol = solve(prob_sde_lorenz)
plot(sol, idxs = (1, 2, 3))
```

Note that it's okay for the noise function to mix terms. For example

```@example sde2
function σ_lorenz!(du, u, p, t)
function g!(du, u, p, t)
du[1] = sin(u[3]) * 3.0
du[2] = u[2] * u[1] * 3.0
du[3] = 3.0
end
```

is a valid noise function, which will once again give diagonal noise by `du2.*W`.
is a valid noise function.

## Example 3: Systems of SDEs with Scalar Noise

Expand All @@ -200,12 +196,12 @@ let's solve a linear SDE with scalar noise using a high order algorithm:
```@example sde3
using DifferentialEquations
using Plots
f(du, u, p, t) = (du .= u)
g(du, u, p, t) = (du .= u)
f!(du, u, p, t) = (du .= u)
g!(du, u, p, t) = (du .= u)
u0 = rand(4, 2)

W = WienerProcess(0.0, 0.0, 0.0)
prob = SDEProblem(f, g, u0, (0.0, 1.0), noise = W)
prob = SDEProblem(f!, g!, u0, (0.0, 1.0), noise = W)
sol = solve(prob, SRIW1())
plot(sol)
```
Expand All @@ -216,7 +212,7 @@ In the previous examples we had diagonal noise, that is a vector of random numbe
`dW` whose size matches the output of `g` where the noise is applied element-wise,
and scalar noise where a single random variable is applied to all dependent variables.
However, a more general type of noise allows for the terms to linearly mixed via `g`
being a matrix.
being a general nondiagonal matrix.

Note that nonlinear mixings are not SDEs but fall under the more general class of
random ordinary differential equations (RODEs) which have a
Expand All @@ -228,8 +224,8 @@ is `g(u,p,t)*dW`, the matrix multiplication. For example, we can do the followin

```@example sde4
using DifferentialEquations
f(du, u, p, t) = du .= 1.01u
function g(du, u, p, t)
f!(du, u, p, t) = du .= 1.01u
function g!(du, u, p, t)
du[1, 1] = 0.3u[1]
du[1, 2] = 0.6u[1]
du[1, 3] = 0.9u[1]
Expand All @@ -239,10 +235,10 @@ function g(du, u, p, t)
du[2, 3] = 0.3u[2]
du[2, 4] = 1.8u[2]
end
prob = SDEProblem(f, g, ones(2), (0.0, 1.0), noise_rate_prototype = zeros(2, 4))
prob = SDEProblem(f!, g!, ones(2), (0.0, 1.0), noise_rate_prototype = zeros(2, 4))
```

In our `g` we define the functions for computing the values of the matrix.
In our `g!` we define the functions for computing the values of the matrix.
We can now think of the SDE that this solves as the system of equations

```math
Expand All @@ -259,7 +255,7 @@ the same random number in the first and second SDEs.
noise. This is discussed [in the SDE solvers page](@ref sde_solve).

The matrix itself is determined by the keyword argument `noise_rate_prototype` in the `SDEProblem`
constructor. This is a prototype for the type that `du` will be in `g`. This can
constructor. This is a prototype for the type that `du` will be in `g!`. This can
be any `AbstractMatrix` type. Thus, we can define the problem as

```@example sde4
Expand All @@ -271,18 +267,18 @@ A[1, 4] = 1
A[2, 4] = 1
A = sparse(A)

# Make `g` write the sparse matrix values
function g(du, u, p, t)
# Make `g!` write the sparse matrix values
function g!(du, u, p, t)
du[1, 1] = 0.3u[1]
du[1, 4] = 0.12u[2]
du[2, 4] = 1.8u[2]
end

# Make `g` use the sparse matrix
prob = SDEProblem(f, g, ones(2), (0.0, 1.0), noise_rate_prototype = A)
# Make `g!` use the sparse matrix
prob = SDEProblem(f!, g!, ones(2), (0.0, 1.0), noise_rate_prototype = A)
```

and now `g(u,p,t)` writes into a sparse matrix, and `g(u,p,t)*dW` is sparse matrix
and now `g!(u,p,t)` writes into a sparse matrix, and `g!(u,p,t)*dW` is sparse matrix
multiplication.

## Example 4: Colored Noise
Expand All @@ -292,7 +288,7 @@ In that portion of the docs, it is shown how to define your own noise process
`my_noise`, which can be passed to the SDEProblem

```julia
SDEProblem(f, g, u0, tspan, noise = my_noise)
SDEProblem(f!, g!, u0, tspan, noise = my_noise)
```

Note that general colored noise problems are only compatible with the `EM` and `EulerHeun` methods.
Expand All @@ -311,11 +307,11 @@ dW_1 dW_2 = ρ dt
In this problem, we have a diagonal noise problem given by:

```@example sde4
function f(du, u, p, t)
function f!(du, u, p, t)
du[1] = μ * u[1]
du[2] = κ * (Θ - u[2])
end
function g(du, u, p, t)
function g!(du, u, p, t)
du[1] = √u[2] * u[1]
du[2] = σ * √u[2]
end
Expand All @@ -339,7 +335,7 @@ heston_noise = CorrelatedWienerProcess!(Γ, tspan[1], zeros(2), zeros(2))
This is then used to build the SDE:

```@example sde4
SDEProblem(f, g, ones(2), tspan, noise = heston_noise)
SDEProblem(f!, g!, ones(2), tspan, noise = heston_noise)
```

Of course, to fully define this problem, we need to define our constants. Constructors
Expand Down
Loading