Time-step choice, negative values and model conditions

A trustworthy SDE result requires more than code that runs: the time grid, state boundaries, coefficient conditions and sensitivity to numerical choices must be examined.

A valid equation does not guarantee a valid numerical path

Euler–Maruyama uses normal increments with unbounded support. Even when the true process should remain in \([0,N]\), a finite-step proposal can be negative or exceed \(N\). The risk grows near boundaries, with large diffusion, or with a large time step.

Time-step requirements

dt <= 0 is invalid exactly because a step length cannot be zero or negative.
\(T/\Delta t\) should be handled consistently so the grid ends at the intended horizon.
Results should be repeated with smaller steps to check numerical stability.
A smaller step reduces discretisation error but increases computation.

This is different from checking whether a floating-point calculation is approximately zero; exact input conditions and computed numerical residuals require different tests.

Four boundary approaches

Approach	Consequence
Stop and report	Exposes invalidity but produces incomplete paths.
Projection or clipping	Simple and bounded, but changes the Euler–Maruyama scheme.
Absorbing rule at zero	Represents extinction once zero is reached, but crossing behaviour remains numerical.
Alternative transformed or positivity-preserving method	Can respect structure better but requires a different numerical method.

No correction should be applied silently.

Interactive sensitivity experiment

Run 500 SIS paths at three time steps. The code records every raw out-of-range proposal before applying the same declared projection rule.

Interactive PythonTime-step and boundary diagnostics

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

N = 200
I0 = 3
beta = 0.25
gamma = 0.15
T = 40
M = 500
time_steps = [0.50, 0.10, 0.02]
base_seed = 150

rows = []
mean_paths = {}

for experiment, dt in enumerate(time_steps):
    steps = int(round(T / dt))
    times = np.linspace(0, T, steps + 1)
    rng = np.random.default_rng(base_seed + experiment)

I = np.zeros((M, steps + 1))
    I[:, 0] = I0
    invalid_proposals = 0
    paths_ever_corrected = np.zeros(M, dtype=bool)

for n in range(steps):
        current = I[:, n]
        infection_rate = beta * (N - current) * current / N
        recovery_rate = gamma * current
        drift = infection_rate - recovery_rate
        diffusion = np.sqrt(np.maximum(
            infection_rate + recovery_rate, 0
        ))

dW = np.sqrt(dt) * rng.normal(0, 1, size=M)
        proposals = current + drift * dt + diffusion * dW

invalid = (proposals < 0) | (proposals > N)
        invalid_proposals += invalid.sum()
        paths_ever_corrected |= invalid

I[:, n + 1] = np.clip(proposals, 0, N)

mean_paths[dt] = (times, I.mean(axis=0))

rows.append({
        "dt": dt,
        "steps": steps,
        "invalid proposals": invalid_proposals,
        "fraction paths corrected": paths_ever_corrected.mean(),
        "mean final I": I[:, -1].mean(),
        "fraction final I = 0": (I[:, -1] == 0).mean()
    })

results = pd.DataFrame(rows)
print(results.round(4).to_string(index=False))

plt.figure(figsize=(9.5, 5.2))
for dt, (times, mean_I) in mean_paths.items():
    plt.plot(times, mean_I, label=f"dt = {dt}")

plt.xlabel("Time (days)")
plt.ylabel("Mean infectious population")
plt.title("Mean paths under three numerical time steps")
plt.legend()
plt.grid(alpha=0.25)
plt.tight_layout()

plt.figure(figsize=(7.8, 4.8))
plt.plot(
    results["dt"],
    results["fraction paths corrected"],
    "o-", color="#b42318"
)
plt.xscale("log")
plt.xlabel("Time step dt (log scale)")
plt.ylabel("Fraction of paths needing correction")
plt.title("Boundary diagnostics must be reported")
plt.grid(alpha=0.25)
plt.tight_layout()

Output

Run the code to see the result.

How to choose a practical time step

Choose a plausible initial \(\Delta t\) from the biological time scale.
Run the same outcome experiment with \(\Delta t/2\) and \(\Delta t/4\).
Compare means, quantiles, extinction measures and correction frequency.
Reduce the step until conclusions are acceptably stable for the question.
Report the chosen step and diagnostic results.

Mathematical coefficient conditions

Standard Euler–Maruyama convergence theory commonly assumes sufficient regularity such as Lipschitz continuity and controlled growth of drift and diffusion. Square-root epidemic diffusion coefficients are not globally Lipschitz at zero, so boundary behaviour deserves particular care. A numerical plot alone is not a proof of convergence.

Validation checklist

Initial compartments are non-negative and sum to \(N\).
Parameters and event rates are non-negative.
The time step is positive and the horizon is reached as intended.
Population conservation errors are checked for multivariable models.
Boundary corrections are counted and reported.
Important outcomes remain stable under step refinement.

What this lesson adds

You can now distinguish model validity from numerical validity, diagnose raw invalid proposals, evaluate a declared boundary rule, perform step-refinement experiments and explain why epidemic square-root diffusions require care near extinction.