In this post I want to explore how randomness affects deterministic systems — visually and mathematically. The setting is the 3-dimensional Lorenz system with the classic parameters $(\sigma,\rho,\beta)=(10,28,8/3)$, which produces the well-known butterfly attractor.

The question is: what happens when we include stochastic forcing on top of that already-chaotic dynamics? We’ll work through additive noise, multiplicative noise, and Ornstein–Uhlenbeck (colored) noise, visualize how the attractor changes under each. If the stochastic calculus feels unfamiliar, my earlier post on Itô’s Lemma covers the foundations — quadratic variation, stochastic integrals, and the SDE variants we’ll be using here.

The deterministic system

The state $s=(x, y, z)^\top \in \mathbb{R}^3$ evolves as:

\[\begin{aligned} \frac{dx}{dt} &= \sigma(y - x) \\ \frac{dy}{dt} &= x(\rho - z) - y \\ \frac{dz}{dt} &= xy - \beta z \end{aligned} \tag{1}\]

With the parameters above, the system has three unstable fixed points and trajectories wind aperiodically between two lobes. All simulations in this post start from $(x_0, y_0, z_0) = (1, 1, 1)$ — the exact initial condition doesn’t matter much since the trajectory quickly converges to the attractor.

Here’s an interactive 3D view (Euler-integrated in your browser, $dt = 0.005$, $9000$ steps). The attractor lies around the two unstable fixed points $C^\pm = (\pm\sqrt{\beta(\rho-1)},\, \pm\sqrt{\beta(\rho-1)},\, \rho-1) \approx (\pm 8.49,\, \pm 8.49,\, 27)$, with one lobe spiraling around each. Rotate and zoom around:

Adding noise

There are basically three standard ways to inject stochasticity into an ODE like this, and they behave quite differently. Let’s go through each of them and check how the system and visualisation changes.

Additive noise

The simplest version — just rely on i.i.d. Brownian increments:

\[\begin{aligned} dx &= \sigma(y - x)\,dt + \sqrt{2D}\,dW_t^{(x)} \\ dy &= [x(\rho - z) - y]\,dt + \sqrt{2D}\,dW_t^{(y)} \\ dz &= (xy - \beta z)\,dt + \sqrt{2D}\,dW_t^{(z)} \end{aligned} \tag{2}\]

$dW_t^{(x)}, dW_t^{(y)}, dW_t^{(z)}$ are independent standard Brownian increments, $D$ sets the noise intensity. The noise amplitude is state-independent — it doesn’t care if $x$ is 0.1 or 100, the fluctuations are the same size. This is appropriate when the perturbation is genuinely external:

  • Thermal noise — at any nonzero temperature, molecules move randomly. This shows up as tiny random forces on anything they bump into — a pollen grain in water (Brownian motion).
  • Wind — turbulent air currents pushing on a physical system. The forcing doesn’t know or depend on the system’s state.
  • Measurement error — sensor imprecision, digitization errors, timing uncertainty. The true trajectory is clean, but what you record has noise on top.

For the interactive simulation below we use $dt = 0.005$, $9000$ steps, and a default noise intensity of $D = 3.0$ (adjustable via the slider). The per-step noise coefficient is $\sqrt{2D\,dt}$. We turned the deterministic system into a stochastic differential equation (SDE) with additive noise, i.e. state-independent diffusion coeffcient $\sqrt{2D}$.

Drag the slider to change $D$ (each change draws a fresh realization):

3.0

Notice that each coordinate of the state has its own randomness injected via $(dW_t^{x}, dW_t^{y}, dW_t^{z})^\top$.

Here is the same system but with shared Brownian motion — one $dW_t$ applied to all three coordinates simultaneously. Compare the shape of the perturbations to the independent case above:

3.0

The trajectory still wiggles, but the perturbations look smoother in 3D: because all three coordinates receive the same random kick at each step, the noise can only push the state along one direction rather than scattering it freely in three dimensions. The overall shape of the attractor is better preserved — the trajectory stays closer to the deterministic butterfly.

Multiplicative noise

More realistic in many settings — the noise scales with the state:

\[\begin{aligned} dx &= \sigma(y - x)\,dt + x\cdot\sqrt{2D}\,dW_t^{(x)} \\ dy &= [x(\rho - z) - y]\,dt + y\cdot\sqrt{2D}\,dW_t^{(y)} \\ dz &= (xy - \beta z)\,dt + z\cdot\sqrt{2D}\,dW_t^{(z)} \end{aligned} \tag{3}\]

Now $x$, $y$, $z$ multiply the noise terms directly. Near the origin the noise vanishes; far away it blows up. This is the same structure as geometric Brownian motion $dX_t = \mu X_t\,dt + \sigma X_t\,dW_t$, which was solved using Itô’s Lemma in my stochastic calculus post. The key insight there: applying $\log$ and Itô’s Lemma produces the $-\sigma^2/2$ drift correction that doesn’t exist in deterministic calculus. This is the natural model whenever fluctuations are proportional — stock returns or population dynamics.

For the simulation below: $dt = 0.005$, $9000$ steps, default $D = 0.02$ (adjustable). Because the state values on the Lorenz system can become large, even a small $D$ produces visible perturbations once multiplied by $x$, $y$, or $z$. Drag the slider:

0.02

With increasing $D$, the global Lorenz attractor shape is lost and mostly influenced by noise.

Colored noise (Ornstein–Uhlenbeck)

In the first additive case above, the noise at each time step is completely uncorrelated with the noise at any other step — it has no memory. But real perturbations usually don’t work that way: a current of wind persists for a while, an instrument drifts slowly, a temperature fluctuation decays over seconds or minutes. To capture that, we replace the raw Brownian motion fluctuations with an Ornstein–Uhlenbeck process that has a tunable correlation time:

\[d\eta_i = -\frac{1}{\tau}\eta_i\,dt + \sqrt{\frac{2D}{\tau}}\,dW_t\]

$\tau$ is the correlation time. $\tau \to 0$ recovers white noise; large $\tau$ gives slowly-wandering forcing. The autocorrelation is exponential:

\[\mathbb{E}[\eta_i(t)\,\eta_i(t')] \propto e^{-|t-t'|/\tau}\]

This follows from the explicit OU solution and Itô isometry — see the derivation below.

Coupling this to the Lorenz drift gives a 6-dimensional system:

\[\begin{aligned} dx &= \sigma(y - x)\,dt + \eta_x\,dt \\ dy &= [x(\rho - z) - y]\,dt + \eta_y\,dt \\ dz &= (xy - \beta z)\,dt + \eta_z\,dt \\ d\eta_x &= -\frac{\eta_x}{\tau}\,dt + \sqrt{\frac{2D}{\tau}}\,dW_t^{(x)} \\ d\eta_y &= -\frac{\eta_y}{\tau}\,dt + \sqrt{\frac{2D}{\tau}}\,dW_t^{(y)} \\ d\eta_z &= -\frac{\eta_z}{\tau}\,dt + \sqrt{\frac{2D}{\tau}}\,dW_t^{(z)} \end{aligned} \tag{4}\]

This matters when the perturbation source itself has inertia — climate variability, slow instrument drift, biological signaling with finite response times. The noise has to “relax” toward new values rather than teleporting.

The OU process has a nice closed-form stationary distribution. In my Itô’s Lemma post I derived the explicit solution of the general OU process $d\eta = -\theta\eta\,dt + \sigma\,dW$ using an integrating-factor trick and Itô’s Lemma, and showed that the stationary variance is $\sigma^2/(2\theta)$. Plugging in $\theta = 1/\tau$ and $\sigma^2 = 2D/\tau$:

\[\eta_i \;\xrightarrow{t\to\infty}\; \mathcal{N}\!\left(0,\; \frac{2D/\tau}{2/\tau}\right) = \mathcal{N}(0,\, D)\]

So the two parameters control different things: $D$ determines how large the noise is (the stationary variance), while $\tau$ determines how fast the noise changes. With a large $\tau$, the noise drifts slowly — it picks a value and sticks with it for a while before moving on. With a small $\tau$, the noise jumps around rapidly from step to step, behaving more and more like the uncorrelated additive noise from earlier.

Interactive OU-driven trajectory (note the smoother deviations compared to the additive case above):

5.0
5.0

What does the noisy attractor look like?

Once we add stochastic fluctuations, the trajectory no longer stays on that exact curve. It wanders around it, tracing out a thicker cloud. The butterfly shape is still recognizable, but now it has width: some regions get visited more densely than others, reflecting the underlying deterministic dynamics.

Below are 2D projections from an OU-noise simulation that show this broadening:

XY-projection of the stochastic Lorenz system with Ornstein-Uhlenbeck noise.
XZ-projection of the stochastic Lorenz system with Ornstein-Uhlenbeck noise.
YZ-projection of the stochastic Lorenz system with Ornstein-Uhlenbeck noise.

Color = time. You can see the lobe-switching dynamics are preserved even with non-trivial OU forcing.

Numerics and Simulation

The visualisations shown earlier implement as baseline method an Euler–Maruyama SDE discretization:

\[X_{n+1} = X_n + f(X_n)\,\Delta t + g(X_n)\,\sqrt{\Delta t}\,\xi_n~,\]

where $f$ is drift, $g$ is diffusion coefficient, $\xi_n \sim \mathcal{N}(0,1)$. The $\sqrt{\Delta t}$ factor comes from the quadratic variation of Brownian motion. (As explained in my previous Itô’s Lemma post.)

For additive noise the code would look like this:

def euler_maruyama_additive(dt, n_steps, x0, y0, z0, D):
    x, y, z = np.zeros(n_steps), np.zeros(n_steps), np.zeros(n_steps)
    x[0], y[0], z[0] = x0, y0, z0
    
    noise_coeff = np.sqrt(2.0 * D * dt)
    
    for i in range(1, n_steps):
        # Deterministic drift
        dx_det = sigma * (y[i-1] - x[i-1])
        dy_det = x[i-1] * (rho - z[i-1]) - y[i-1]
        dz_det = x[i-1] * y[i-1] - beta * z[i-1]
        
        # Stochastic increment
        dW = np.random.randn(3)
        
        x[i] = x[i-1] + dx_det * dt + noise_coeff * dW[0]
        y[i] = y[i-1] + dy_det * dt + noise_coeff * dW[1]
        z[i] = z[i-1] + dz_det * dt + noise_coeff * dW[2]
    
    return x, y, z

(For multiplicative noise, replace noise_coeff * dW[0] with noise_coeff * x[i-1] * dW[0], etc. For OU you need to carry the $\eta$ state too.

Code & notebooks

I put together a small Python module for playing with all of this: deterministic, additive, multiplicative, and OU-noise on the Lorenz system, plus some plotting utilities.

  • Jupyter notebook — interactive exploration (might be too large to display on GitHub. Recommended to download and open locally)
  • Python module — importable
from lorenz_system import LorenzSystem, LorenzPlotter, NoiseType

# Create system with default parameters
lorenz = LorenzSystem()

# Solve with Ornstein-Uhlenbeck noise
result = lorenz.solve(
    noise_type=NoiseType.ORNSTEIN_UHLENBECK,
    D=5.0,
    correlation_time=5.0,
    t_span=(0, 50),
    dt=0.001
)

# Visualize
plotter = LorenzPlotter(gradient_name="pastel")
plotter.plot_3d(result, title="Stochastic Lorenz Attractor")
plotter.plot_2d_projections(result)

Appendix: OU autocorrelation derivation

Start from the explicit solution of the OU SDE, obtained via the integrating-factor trick (derived in my Itô’s Lemma post):

\[\eta(t) = \eta(0)\,e^{-t/\tau} + \sqrt{\frac{2D}{\tau}}\int_0^t e^{-(t-s)/\tau}\,dW_s\]

To compute $\mathbb{E}[\eta(t)\,\eta(t’)]$ in the stationary regime, we drop the transient $\eta(0)\,e^{-t/\tau}$ term as it vanishes as $t \to \infty$. Next, we apply Itô isometry:

\[\mathbb{E}\!\left[\int f\,dW \cdot \int g\,dW\right] = \int fg\,ds\]

to get (assuming $t’ \ge t$):

\[\mathbb{E}[\eta(t)\,\eta(t')] = \frac{2D}{\tau}\int_0^t e^{-(t-s)/\tau}\,e^{-(t'-s)/\tau}\,ds = \frac{2D}{\tau}\,e^{-(t+t')/\tau}\int_0^t e^{2s/\tau}\,ds = D\,e^{-|t'-t|/\tau}\]

This confirms the exponential autocorrelation with correlation time $\tau$ and stationary variance $D$.

Further reading

  • My post on Itô’s Lemma for more details on stochastic calculus.
  • Øksendal, Stochastic Differential Equations — the standard complete reference

All code is on GitHub.