@@ -19,19 +19,27 @@ Let's take the Lorenz equation and add noise to each of the states.
1919To show the flexibility of ModelingToolkit,
2020we do not use homogeneous noise, with constant variance,
2121but instead use heterogeneous noise,
22- where the magnitude of the noise scales with (0.1 times) the magnitude of each of the states:
22+ where the magnitude of the noise scales with (0.3 times) the magnitude of each of the states:
2323
2424``` math
2525\begin{aligned}
26- dx &= (\sigma (y-x))dt &+ 0.1xdB \\
27- dy &= (x(\rho-z) - y)dt &+ 0.1ydB \\
28- dz &= (xy - \beta z)dt &+ 0.1zdB \\
26+ \frac{dx}{dt} &= (\sigma (y-x)) &+ 0.1x\frac{dB}{dt} \\
27+ \frac{dy}{dt} &= (x(\rho-z) - y) &+ 0.1y\frac{dB}{dt} \\
28+ \frac{dz}{dt} &= (xy - \beta z) &+ 0.1z\frac{dB}{dt} \\
2929\end{aligned}
3030```
3131
3232Where $B$, is standard Brownian motion, also called the
3333[ Wiener process] ( https://en.wikipedia.org/wiki/Wiener_process ) .
34- In ModelingToolkit this can be done by ` @brownian ` variables.
34+ We use notation similar to the
35+ [ Langevin equation] ( https://en.wikipedia.org/wiki/Stochastic_differential_equation#Use_in_physics ) ,
36+ often used in physics.
37+ By "multiplying" the equations by $dt$, the notation used in
38+ [ probability theory] ( https://en.wikipedia.org/wiki/Stochastic_differential_equation#Use_in_probability_and_mathematical_finance )
39+ can be recovered.
40+
41+ We use this Langevin-like notation because it allows us to extend MTK modeling capacity from ODEs to SDEs,
42+ using only a single new concept, ` @brownian ` variables, which represent $\frac{dB}{dt}$ in the above equation.
3543
3644``` @example SDE
3745using ModelingToolkit, StochasticDiffEq
@@ -41,9 +49,9 @@ using Plots
4149@parameters σ=10.0 ρ=2.33 β=26.0
4250@variables x(t)=5.0 y(t)=5.0 z(t)=1.0
4351@brownian B
44- eqs = [D(x) ~ σ * (y - x) + 0.1B * x ,
45- D(y) ~ x * (ρ - z) - y + 0.1B * y ,
46- D(z) ~ x * y - β * z + 0.1B * z ]
52+ eqs = [D(x) ~ σ * (y - x) + 0.3x * B ,
53+ D(y) ~ x * (ρ - z) - y + 0.3y * B ,
54+ D(z) ~ x * y - β * z + 0.3z * B ]
4755
4856@mtkbuild de = System(eqs, t)
4957```
@@ -58,22 +66,29 @@ We continue by solving and plotting the SDE.
5866
5967``` @example SDE
6068prob = SDEProblem(de, [], (0.0, 100.0), [])
61- sol = solve(prob, LambaEulerHeun ())
69+ sol = solve(prob, SRIW1 ())
6270plot(sol, idxs = [(1, 2, 3)])
6371```
6472
6573The noise present in all 3 equations is correlated, as can be seen on the below figure.
66- If you want uncorrelated noise for each equation,
67- multiple ` @brownian ` variables have to be declared.
74+ The figure also shows the multiplicative nature of the noise.
75+ Because states ` x ` and ` y ` generally take on larger values,
76+ the noise also takes on a more pronounced effect on these states compared to the state ` z ` .
6877
6978``` @example SDE
70- @brownian Bx By Bz
79+ plot(sol)
7180```
7281
73- The figure also shows the multiplicative nature of the noise.
74- Because states ` x ` and ` y ` generally take on larger values,
75- the noise also takes on a more pronounced effect on these states compared to the state ` z ` .
82+ If you want uncorrelated noise for each equation,
83+ multiple ` @brownian ` variables have to be declared.
7684
7785``` @example SDE
86+ @brownian Bx By Bz
87+ eqs = [D(x) ~ σ * (y - x) + 0.3x * Bx,
88+ D(y) ~ x * (ρ - z) - y + 0.3y * By,
89+ D(z) ~ x * y - β * z + 0.3z * Bz]
90+ @mtkbuild de = System(eqs, t)
91+ prob = SDEProblem(de, [], (0.0, 100.0), [])
92+ sol = solve(prob, SRIW1())
7893plot(sol)
7994```
0 commit comments