Stochastic simulations
Stochastic simulations are supported via the Euler-Maruyama method. Random normal variables \(W_t\sim N(0, dt)\) can be included in the vector field function, and each simulation instance will have an independent stream of random variables.
As an example, consider the Ornstein-Uhlenbeck process with \(\theta = 1\):
\[ dx_t = (\mu - x_t)dt + \sigma W_t \]
To implement this system, we heuristically write it in Langevin form:
\[ \frac{dx}{dt} = \mu - x + \sigma\eta(t)\]
where \(\eta(t)\) represents white noise. When \(\sigma=0\), the system is a first-order ODE with steady state x=\(\mu\). Increasing \(\sigma\) produces solutions with a steady state distribution \(x\sim N(\mu, \frac{\sigma^2}{2})\).
>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> import clode
>>> from typing import List
>>> def wiener_process(
... t: float,
... y: List[float],
... p: List[float],
... dy: List[float],
... aux: List[float],
... wiener: List[float],
... ) -> None:
... x: float = y[0]
... mu: float = p[0]
... sigma: float = p[1]
... weiner_variable: float = wiener[0]
... dx: float = mu - x + sigma * weiner_variable
... dy[0] = dx
...
>>> variables = {"x": 0.0}
>>> parameters = {"mu": 1.0, "sigma": 0.5}
>>> t_span = (0.0, 1000.0)
>>> integrator = clode.Simulator(
... rhs_equation=wiener_process,
... variables=variables,
... parameters=parameters,
... num_noise=1,
... stepper=clode.Stepper.stochastic_euler,
... single_precision=True,
... t_span=t_span,
... dt=0.001,
... )
>>> nPts = 8192
>>> integrator.set_repeat_ensemble(nPts)
>>> integrator.transient()
>>> XF = integrator.get_final_state()
>>> plt.hist(XF, 30)
(array([ 1., 2., 5., 18., 35., 62., 117., 180., 269., 422., 576.,
736., 762., 813., 876., 810., 717., 569., 432., 326., 193., 126.,
69., 31., 25., 6., 6., 5., 1., 2.]), array([-0.3228741 , -0.22877825, -0.1346824 , -0.04058655, 0.0535093 ,
0.14760516, 0.24170101, 0.33579686, 0.42989271, 0.52398856,
0.61808441, 0.71218026, 0.80627611, 0.90037196, 0.99446781,
1.08856367, 1.18265952, 1.27675537, 1.37085122, 1.46494707,
1.55904292, 1.65313877, 1.74723462, 1.84133047, 1.93542632,
2.02952218, 2.12361803, 2.21771388, 2.31180973, 2.40590558,
2.50000143]), <BarContainer object of 30 artists>)
>>> plt.xlabel("x")
Text(0.5, 0, 'x')
>>> plt.show()
>>> print(f"mean xf: {np.mean(XF) :0.5}")
mean xf: 1.0036
>>> print(f"simulation variance: {np.var(XF) :0.5}")
simulation variance: 0.12492
>>> print(f"expected variance: {parameters['sigma'] ** 2 / 2 :0.5}")
expected variance: 0.125
Output:
mean xf: 1.0041
simulation variance: 0.12629
expected variance: 0.125