Skip to content

Commit

Permalink
Merge pull request #3084 from ArnoStrouwen/docs2
Browse files Browse the repository at this point in the history
Rewrite of SDE tutorial
  • Loading branch information
ChrisRackauckas authored Oct 12, 2024
2 parents 84ae858 + 87cbf0e commit aaf6893
Showing 1 changed file with 80 additions and 24 deletions.
104 changes: 80 additions & 24 deletions docs/src/tutorials/stochastic_diffeq.md
Original file line number Diff line number Diff line change
@@ -1,38 +1,94 @@
# Modeling with Stochasticity

All models with `ODESystem` are deterministic. `SDESystem` adds another element
to the model: randomness. This is a
All previous differential equations tutorials deal with deterministic `ODESystem`s.
In this tutorial, we add randomness.
In particular, we show how to represent a
[stochastic differential equation](https://en.wikipedia.org/wiki/Stochastic_differential_equation)
which has a deterministic (drift) component and a stochastic (diffusion)
component. Let's take the Lorenz equation from the first tutorial and extend
it to have multiplicative noise by creating `@brownian` variables in the equations.
as a `SDESystem`.

!!! note

The high level `@mtkmodel` macro used in the
[getting started tutorial](@ref getting_started)
is not yet compatible with `SDESystem`.
We thus have to use a lower level interface to define stochastic differential equations.
For an introduction to this interface, read the
[programmatically generating ODESystems tutorial](@ref programmatically).

Let's take the Lorenz equation and add noise to each of the states.
To show the flexibility of ModelingToolkit,
we do not use homogeneous noise, with constant variance,
but instead use heterogeneous noise,
where the magnitude of the noise scales with (0.3 times) the magnitude of each of the states:

```math
\begin{aligned}
\frac{dx}{dt} &= (\sigma (y-x)) &+ 0.1x\frac{dB}{dt} \\
\frac{dy}{dt} &= (x(\rho-z) - y) &+ 0.1y\frac{dB}{dt} \\
\frac{dz}{dt} &= (xy - \beta z) &+ 0.1z\frac{dB}{dt} \\
\end{aligned}
```

Where $B$, is standard Brownian motion, also called the
[Wiener process](https://en.wikipedia.org/wiki/Wiener_process).
We use notation similar to the
[Langevin equation](https://en.wikipedia.org/wiki/Stochastic_differential_equation#Use_in_physics),
often used in physics.
By "multiplying" the equations by $dt$, the notation used in
[probability theory](https://en.wikipedia.org/wiki/Stochastic_differential_equation#Use_in_probability_and_mathematical_finance)
can be recovered.

We use this Langevin-like notation because it allows us to extend MTK modeling capacity from ODEs to SDEs,
using only a single new concept, `@brownian` variables, which represent $\frac{dB}{dt}$ in the above equation.

```@example SDE
using ModelingToolkit, StochasticDiffEq
using ModelingToolkit: t_nounits as t, D_nounits as D
using Plots
# Define some variables
@parameters σ ρ β
@variables x(t) y(t) z(t)
@brownian a
eqs = [D(x) ~ σ * (y - x) + 0.1a * x,
D(y) ~ x * (ρ - z) - y + 0.1a * y,
D(z) ~ x * y - β * z + 0.1a * z]
@parameters σ=10.0 ρ=2.33 β=26.0
@variables x(t)=5.0 y(t)=5.0 z(t)=1.0
@brownian B
eqs = [D(x) ~ σ * (y - x) + 0.3x * B,
D(y) ~ x * (ρ - z) - y + 0.3y * B,
D(z) ~ x * y - β * z + 0.3z * B]
@mtkbuild de = System(eqs, t)
```

Even though we did not explicitly use `SDESystem`, ModelingToolkit can still infer this from the equations.

```@example SDE
typeof(de)
```

We continue by solving and plotting the SDE.

```@example SDE
prob = SDEProblem(de, [], (0.0, 100.0), [])
sol = solve(prob, SRIW1())
plot(sol, idxs = [(1, 2, 3)])
```

u0map = [
x => 1.0,
y => 0.0,
z => 0.0
]
The noise present in all 3 equations is correlated, as can be seen on the below figure.
The figure also shows the multiplicative nature of the noise.
Because states `x` and `y` generally take on larger values,
the noise also takes on a more pronounced effect on these states compared to the state `z`.

parammap = [
σ => 10.0,
β => 26.0,
ρ => 2.33
]
```@example SDE
plot(sol)
```

If you want uncorrelated noise for each equation,
multiple `@brownian` variables have to be declared.

prob = SDEProblem(de, u0map, (0.0, 100.0), parammap)
sol = solve(prob, LambaEulerHeun())
```@example SDE
@brownian Bx By Bz
eqs = [D(x) ~ σ * (y - x) + 0.3x * Bx,
D(y) ~ x * (ρ - z) - y + 0.3y * By,
D(z) ~ x * y - β * z + 0.3z * Bz]
@mtkbuild de = System(eqs, t)
prob = SDEProblem(de, [], (0.0, 100.0), [])
sol = solve(prob, SRIW1())
plot(sol)
```

0 comments on commit aaf6893

Please sign in to comment.