Euler–Maruyama method from scratch

Euler–Maruyama is the simplest explicit method for approximating an SDE on a time grid. It adds one Brownian diffusion contribution to the familiar Euler drift update.

From the SDE to a computable update

\[dX=f(X,t)dt+g(X,t)dW.\]

On the grid \(t_n=n\Delta t\), Euler–Maruyama replaces the infinitesimal expression by

\[X_{n+1}=X_n+f(X_n,t_n)\Delta t+g(X_n,t_n)\Delta W_n,\qquad \Delta W_n=\sqrt{\Delta t}Z_n.\]

Where the update comes from

Over one interval from \(t_n\) to \(t_{n+1}\), the SDE can be written in integral form:

\[ X(t_{n+1})-X(t_n) = \int_{t_n}^{t_{n+1}}f(X(s),s)ds + \int_{t_n}^{t_{n+1}}g(X(s),s)dW(s). \]

Euler–Maruyama holds both coefficients at their known left-endpoint values during this short interval:

\[ \int_{t_n}^{t_{n+1}}f(X(s),s)ds \approx f(X_n,t_n)\Delta t, \qquad \int_{t_n}^{t_{n+1}}g(X(s),s)dW(s) \approx g(X_n,t_n)\Delta W_n. \]

Adding these approximations to \(X_n\) gives the Euler–Maruyama update. This is a numerical approximation, not an algebraic cancellation of \(dt\) or \(dW\).

Four operations in every step

Evaluate drift at the current state.
Evaluate diffusion at the current state.
Generate one new Brownian increment.
Add drift change and diffusion change to the current state.

Both coefficients use \(X_n\), not the already updated value \(X_{n+1}\). This makes the method explicit.

Scenario: fluctuating logistic population

Let \(X(t)\) be a population growing toward carrying capacity \(K\), while environmental variation changes growth continuously:

\[dX=rX\left(1-\frac{X}{K}\right)dt+\eta XdW.\]

The drift is logistic growth. The diffusion scale \(\eta X\) makes environmental fluctuation proportional to current population size.

Calculate the first update by hand

For the values used in the program, \(X_0=40\), \(r=0.8\), \(K=500\), \(\eta=0.12\) and \(\Delta t=0.02\):

\[ f(40)\Delta t = 0.8(40)\left(1-\frac{40}{500}\right)(0.02) = 0.5888. \]

With seed 31, the first generated value is approximately \(Z_0=-0.3953\), so

\[ \Delta W_0=\sqrt{0.02}(-0.3953)\approx-0.0559, \] \[ g(40)\Delta W_0=(0.12)(40)(-0.0559)\approx-0.2683. \]

\[ X_1=40+0.5888-0.2683\approx40.3205. \]

The population increases because the positive drift contribution is larger than the negative diffusion contribution in this particular step.

Interactive Python laboratory

The code implements Euler–Maruyama line by line and compares it with deterministic Euler using the same drift and initial value.

Interactive PythonEuler–Maruyama from scratch

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

def drift(x, r, K):
    return r * x * (1 - x / K)

def diffusion(x, eta):
    return eta * x

# Model and numerical settings
r = 0.8
K = 500
eta = 0.12
X0 = 40
T = 20
dt = 0.02
seed = 31

if dt <= 0:
    raise ValueError("dt must be positive.")

if T <= 0:
    raise ValueError("T must be positive.")

if K <= 0 or eta < 0 or X0 < 0:
    raise ValueError("Use biologically valid inputs.")

steps = int(round(T / dt))
if not np.isclose(steps * dt, T):
    raise ValueError("Choose dt so that T/dt is an integer.")

times = np.linspace(0, T, steps + 1)

rng = np.random.default_rng(seed)
X_sde = np.zeros(steps + 1)
X_deterministic = np.zeros(steps + 1)
X_sde[0] = X0
X_deterministic[0] = X0

first_steps = []

for n in range(steps):
    Z = rng.normal(0, 1)
    dW = np.sqrt(dt) * Z

drift_change = drift(X_sde[n], r, K) * dt
    diffusion_change = diffusion(
        X_sde[n], eta
    ) * dW

X_sde[n + 1] = (
        X_sde[n]
        + drift_change
        + diffusion_change
    )

X_deterministic[n + 1] = (
        X_deterministic[n]
        + drift(X_deterministic[n], r, K) * dt
    )

if n < 8:
        first_steps.append({
            "n": n,
            "time": times[n],
            "X_n": X_sde[n],
            "Z_n": Z,
            "dW_n": dW,
            "drift change": drift_change,
            "diffusion change": diffusion_change,
            "X_next": X_sde[n + 1]
        })

step_table = pd.DataFrame(first_steps)

print("First 8 Euler–Maruyama updates")
print(step_table.round(4).to_string(index=False))
print()
print(f"Final deterministic value: {X_deterministic[-1]:.3f}")
print(f"Final stochastic value: {X_sde[-1]:.3f}")
print(f"Minimum stochastic value: {X_sde.min():.3f}")

plt.figure(figsize=(9.5, 5.2))
plt.plot(times, X_deterministic,
         color="#174c34", linewidth=2.2,
         label="Deterministic Euler")
plt.plot(times, X_sde,
         color="#b42318", linewidth=1.2,
         label="Euler–Maruyama path")
plt.axhline(K, color="black", linestyle="--",
            label="Carrying capacity")
plt.xlabel("Time")
plt.ylabel("Population X")
plt.title("Drift creates the trend; diffusion creates fluctuation")
plt.legend()
plt.grid(alpha=0.25)
plt.tight_layout()

Output

Run the code to see the result.

Read one stored row

Column	Meaning
`X_n`	Current population before the update.
`dW_n`	New Brownian increment for this interval.
`drift change`	Systematic contribution \(f(X_n)\Delta t\).
`diffusion change`	Random contribution \(g(X_n)\Delta W_n\).
`X_next`	Sum of the current state and both contributions.

Euler versus Euler–Maruyama

Euler	Euler–Maruyama
\(x_{n+1}=x_n+f(x_n)\Delta t\)	Adds \(g(x_n)\Delta W_n\) as well.
Same inputs give the same path.	Different Brownian increments give different paths.
Approximates an ODE.	Approximates an SDE.

Numerical cautions

A smaller \(\Delta t\) usually represents the SDE more accurately but needs more steps. Euler–Maruyama may generate biologically impossible negative values for some models. Do not silently replace them without explaining the chosen boundary treatment; a later lesson examines this issue directly.

What this lesson adds

You can now translate an SDE into an explicit Euler–Maruyama loop, separate drift and diffusion changes, store diagnostic rows, and compare a stochastic path with deterministic Euler.