utils.solver.milstein_method

utils.solver.milstein_method(
    f,
    t_span,
    y0,
    dt,
    noise_func=None,
    rng=None,
    args=(),
)

Fixed-step Milstein integrator for stochastic differential equations.

Achieves strong order 1.0 for both additive and multiplicative noise (diagonal diffusion), making it the right choice when sde_noise depends on the state. For additive noise it reduces to Euler-Maruyama plus a zero correction; prefer :func:heun_maruyama_method in that case as it also improves the drift approximation.

Solves:

dy = f(t, y) dt + g(t, y) dW

using the Milstein correction::

y' = y + f dt + g dW + ½ g (∂g/∂y)(dW² - dt)

The derivative ∂g/∂y is approximated element-wise by a forward finite difference with step h = √ε · max(1, |y|), so no analytical Jacobian is required from the model author.

Parameters

f : callable

Drift function with signature f(t, y, *args).

t_span : tuple of float

(t0, tf) integration bounds.

y0 : array - like

Initial state vector.

dt : float

Fixed timestep.

noise_func : callable or None = None

Diffusion function with signature noise_func(t, y), returning a vector of per-state diffusion scales. If None, the stochastic term is zero and forward Euler is recovered.

rng : numpy.random.Generator or None = None

Random generator for Wiener increments. A fresh generator is created if None.

args : tuple = ()

Extra positional arguments forwarded to f.

Returns

solution : Solution

Object with attributes t (time axis) and y (state trajectory).

Raises

: ValueError

If noise_func returns a vector whose shape does not match the state vector.

Notes

The finite-difference step uses h = √(machine_eps) · max(1, |yᵢ|) per component, which balances truncation and floating-point cancellation errors. For SDEs where an exact ∂g/∂y is available, accuracy can be improved by subclassing and overriding the correction directly, but the finite-difference approximation is adequate for all current ClimateCritters use cases.