utils.solver.heun_maruyama_method
utils.solver.heun_maruyama_method(
f,
t_span,
y0,
dt,
noise_func=None,
rng=None,
args=(),
)Fixed-step Heun-Maruyama integrator for stochastic differential equations.
A predictor-corrector scheme that achieves strong order 1.0 for SDEs with additive noise (diffusion independent of state) and weak order 2.0. This is a meaningful improvement over :func:euler_maruyama_method (strong order 0.5) when transition timing is the quantity of interest, as in bistable climate models.
Solves:
dy = f(t, y) dt + g(t, y) dW
using the two-stage Heun scheme::
ỹ = y + f(t, y) dt + g(t, y) dW (Euler predictor)
y' = y + ½[f(t, y) + f(t+dt, ỹ)] dt
+ ½[g(t, y) + g(t+dt, ỹ)] dW (Heun corrector)
A single Wiener increment dW is shared between predictor and corrector steps, which is the standard approach for strong convergence.
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. IfNone, the stochastic term is zero and the Heun deterministic ODE solver 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) andy(state trajectory).
Raises
: ValueError-
If
noise_funcreturns a vector whose shape does not match the state vector.
Notes
For purely additive noise (g constant or state-independent) this scheme achieves strong order 1.0. For multiplicative noise (g depends on y) strong order drops back toward 0.5 but weak order 2.0 is retained, still outperforming Euler-Maruyama in distribution-level statistics. See Rößler (2010) for a full convergence analysis.