Stochastic solver API

RosslerSRA

Holds the Butcher tableaus for a Rosser SRA method.

RosslerSRI

Holds the Butcher tableaus for a Roessler SRI method.

SOSRA()

SOSRA: Stability-Optimized SRA Method (Nonstiff) - Optimal for Additive Noise

Stability-optimized adaptive Stochastic Runge-Kutta method for additive noise problems. This is the optimal choice for additive noise SDEs.

Method Properties

• Strong Order: 1.5 (for additive noise)
• Weak Order: 2.0
• Time stepping: Adaptive
• Noise types: Additive noise (diagonal, non-diagonal, and scalar)
• SDE interpretation: Both Itô and Stratonovich
• Stability: Optimized for high tolerances and robust to stiffness

When to Use

• Optimal choice for additive noise problems: du = f(u,p,t)dt + σ dW
• When the diffusion term is independent of the solution u
• For problems requiring high accuracy with additive noise
• When using high tolerances (method is stable)
• For both Itô and Stratonovich interpretations

Algorithm Description

SOSRA is a stability-optimized version of the SRA (Stochastic Runge-Kutta for Additive noise) methods. It exploits the special structure of additive noise to achieve better performance and stability.

Additive Noise Structure

Specialized for SDEs of the form:

\[du = f(u,p,t) dt + σ(t) dW\]

where the diffusion σ does not depend on the solution u.

References

• Rößler A., "Runge–Kutta Methods for the Strong Approximation of Solutions of Stochastic Differential Equations", SIAM J. Numer. Anal., 48 (3), pp. 922–952

SOSRA2()

SOSRA2: Stability-Optimized SRA Method Version 2 (Nonstiff)

Alternative stability-optimized adaptive SRA method for additive noise problems.

Method Properties

• Strong Order: 1.5 (for additive noise)
• Weak Order: 2.0
• Time stepping: Adaptive
• Noise types: Additive noise (diagonal, non-diagonal, and scalar)
• SDE interpretation: Both Itô and Stratonovich
• Stability: Optimized for high tolerances and robust to stiffness

When to Use

• Alternative to SOSRA for additive noise problems
• Different stability characteristics may be preferred for specific problems
• When SOSRA performance is unsatisfactory

References

• Rößler A., "Runge–Kutta Methods for the Strong Approximation of Solutions of Stochastic Differential Equations", SIAM J. Numer. Anal., 48 (3), pp. 922–952

SOSRI()

SOSRI: Stability-Optimized SRI Method (Nonstiff) - Recommended

The Stability-Optimized Stochastic Runge-Kutta method. This is the recommended method for general-purpose solving of diagonal/scalar Itô SDEs.

Method Properties

• Strong Order: 1.5 (for diagonal/scalar noise)
• Weak Order: 2.0
• Time stepping: Adaptive
• Noise types: Diagonal and scalar noise only
• SDE interpretation: Itô
• Stability: Optimized for high tolerances and robust to mild stiffness

When to Use

• Recommended as first choice for diagonal/scalar Itô SDEs
• When high accuracy is required (strong order 1.5)
• For problems with mild stiffness
• When using high tolerances (method is stable)
• For most general SDE applications

Algorithm Description

SOSRI is a stability-optimized version of the SRI methods with specially chosen coefficients to improve stability properties. It provides excellent performance for the most common class of SDE problems.

Restrictions

• Only works with diagonal or scalar noise
• For non-diagonal noise, use other methods like RKMilCommute or LambaEM

References

• Rößler A., "Runge–Kutta Methods for the Strong Approximation of Solutions of Stochastic Differential Equations", SIAM J. Numer. Anal., 48 (3), pp. 922–952

SOSRI2()

SOSRI2: Alternative Stability-Optimized SRI Method (Nonstiff)

Alternative stability-optimized adaptive strong order 1.5 method with different stability characteristics than SOSRI.

Method Properties

• Strong Order: 1.5 (for diagonal/scalar noise)
• Weak Order: 2.0
• Time stepping: Adaptive
• Noise types: Diagonal and scalar noise only
• SDE interpretation: Itô
• Stability: Optimized for high tolerances and robust to stiffness

When to Use

• Alternative to SOSRI with different stability properties
• When SOSRI performance is unsatisfactory
• For benchmarking stability-optimized methods
• Problems requiring different stability characteristics

Algorithm Features

• Different stability optimization than SOSRI
• May perform better on certain problem types
• Maintains high tolerance robustness

References

• Stability-optimized SRI methods

SRA(; tableau = constructSRA1())

SRA: Configurable Stochastic Runge-Kutta for Additive Noise (Nonstiff)

Configurable adaptive strong order 1.5 method for additive noise problems with customizable tableaux.

Method Properties

• Strong Order: 1.5 (for additive noise)
• Weak Order: Depends on tableau (typically 2.0)
• Time stepping: Adaptive
• Noise types: Additive noise (diagonal, non-diagonal, and scalar)
• SDE interpretation: Both Itô and Stratonovich

Parameters

• tableau: Tableau specification (default: constructSRA1())

When to Use

• When custom tableaux are needed for additive noise problems
• For research and experimentation with SRA methods
• When default methods don't provide desired characteristics
• For benchmarking different SRA variants

Available Tableaux

• constructSRA1(): Default SRA1 tableau
• Custom tableaux can be constructed for specialized applications

References

• Rößler A., "Runge–Kutta Methods for the Strong Approximation of Solutions of Stochastic Differential Equations", SIAM J. Numer. Anal., 48 (3), pp. 922–952. DOI: 10.1137/09076636X.

SRA1()

SRA1: Stochastic Runge-Kutta A1 Method (Nonstiff)

Adaptive strong order 1.5 method for additive Itô and Stratonovich SDEs with weak order 2.

Method Properties

• Strong Order: 1.5 (for additive noise)
• Weak Order: 2.0
• Time stepping: Adaptive
• Noise types: Additive noise (diagonal, non-diagonal, and scalar)
• SDE interpretation: Both Itô and Stratonovich

When to Use

• Standard choice for additive noise problems
• When proven theoretical properties are important
• Alternative to SOSRA when stability optimization is not needed
• For both Itô and Stratonovich problems with additive noise

Additive Noise Structure

Specialized for SDEs of the form:

\[du = f(u,p,t) dt + σ(p,t) dW\]

where diffusion σ doesn't depend on solution u.

References

• Rößler A., "Runge–Kutta Methods for the Strong Approximation of Solutions of Stochastic Differential Equations", SIAM J. Numer. Anal., 48 (3), pp. 922–952. DOI:10.1137/09076636X

SRA2()

SRA2: Stochastic Runge-Kutta A2 Method (Nonstiff)

Alternative adaptive strong order 1.5 method for additive noise problems with different coefficients.

Method Properties

• Strong Order: 1.5 (for additive noise)
• Weak Order: 2.0
• Time stepping: Adaptive
• Noise types: Additive noise (diagonal, non-diagonal, and scalar)
• SDE interpretation: Both Itô and Stratonovich

When to Use

• Alternative to SRA1 with different stability/accuracy characteristics
• When SRA1 performance is unsatisfactory
• For benchmarking different SRA variants
• Research and comparison studies

References

• Rößler A., "Runge–Kutta Methods for the Strong Approximation of Solutions of Stochastic Differential Equations", SIAM J. Numer. Anal., 48 (3), pp. 922–952. DOI:10.1137/09076636X

SRA3()

SRA3: Stochastic Runge-Kutta A3 Method (Nonstiff)

Adaptive strong order 1.5 method for additive noise problems with weak order 3.

Method Properties

• Strong Order: 1.5 (for additive noise)
• Weak Order: 3.0
• Time stepping: Adaptive
• Noise types: Additive noise (non-diagonal and scalar)
• SDE interpretation: Both Itô and Stratonovich

When to Use

• When weak order 3.0 convergence is required for additive noise
• For Monte Carlo simulations needing highest weak accuracy
• Problems where weak convergence dominates computational cost
• When computational cost per step is acceptable for higher weak order

Restrictions

• Does not handle diagonal additive noise (use SRA1/SRA2 instead)
• Limited to non-diagonal and scalar additive noise structures

Algorithm Features

• Highest weak order in the SRA family
• More expensive per step than SRA1/SRA2
• Excellent for statistical calculations requiring high weak accuracy

References

• Rößler A., "Runge–Kutta Methods for the Strong Approximation of Solutions of Stochastic Differential Equations", SIAM J. Numer. Anal., 48 (3), pp. 922–952. DOI:10.1137/09076636X

SRI(; tableau = constructSRIW1(), error_terms = 4)

SRI: Configurable Stochastic Runge-Kutta for Itô SDEs (Nonstiff)

Configurable adaptive strong order 1.5 method for diagonal/scalar Itô SDEs with customizable tableaux.

Method Properties

• Strong Order: 1.5 (for diagonal/scalar noise)
• Weak Order: Depends on tableau (typically 2.0)
• Time stepping: Adaptive
• Noise types: Diagonal and scalar noise only
• SDE interpretation: Itô

Parameters

• tableau: Tableau specification (default: constructSRIW1())
• error_terms::Int = 4: Number of error terms for adaptive stepping

When to Use

• When custom tableaux are needed for diagonal/scalar problems
• For research and experimentation with SRI methods
• When default methods don't provide desired characteristics
• For benchmarking different SRI variants

Available Tableaux

• constructSRIW1(): Default SRIW1 tableau
• Custom tableaux can be constructed for specialized applications

References

• Rößler A., "Runge–Kutta Methods for the Strong Approximation of Solutions of Stochastic Differential Equations", SIAM J. Numer. Anal., 48 (3), pp. 922–952. DOI: 10.1137/09076636X.

SRIW1()

SRIW1: Stochastic Runge-Kutta W1 Method (Nonstiff)

Adaptive stochastic Runge-Kutta method with strong order 1.5 and weak order 2.0 for diagonal/scalar Itô SDEs.

Method Properties

• Strong Order: 1.5 (for diagonal/scalar noise)
• Weak Order: 2.0
• Time stepping: Adaptive
• Noise types: Diagonal and scalar noise only
• SDE interpretation: Itô

When to Use

• Standard choice for diagonal/scalar Itô SDEs
• When proven theoretical properties are important
• Alternative to SOSRI when stability optimization is not needed
• For problems requiring exactly weak order 2.0

Algorithm Features

• Embedded error estimation for adaptive stepping
• Well-established theoretical foundation
• Good balance of accuracy and efficiency

References

• Rößler A., "Runge–Kutta Methods for the Strong Approximation of Solutions of Stochastic Differential Equations", SIAM J. Numer. Anal., 48 (3), pp. 922–952. DOI: 10.1137/09076636X.

SRIW2()

SRIW2: Stochastic Runge-Kutta W2 Method (Nonstiff)

Adaptive stochastic Runge-Kutta method with strong order 1.5 and weak order 3.0 for diagonal/scalar Itô SDEs.

Method Properties

• Strong Order: 1.5 (for diagonal/scalar noise)
• Weak Order: 3.0
• Time stepping: Adaptive
• Noise types: Diagonal and scalar noise only
• SDE interpretation: Itô

When to Use

• When weak order 3.0 convergence is required
• For Monte Carlo simulations needing high weak accuracy
• Problems where weak convergence is more important than strong
• When computational cost per step is acceptable for higher weak order

Algorithm Features

• Highest weak order in the SRI family
• More expensive per step than SRIW1
• Excellent for statistical calculations and expectations

References

• Rößler A., "Runge–Kutta Methods for the Strong Approximation of Solutions of Stochastic Differential Equations", SIAM J. Numer. Anal., 48 (3), pp. 922–952. DOI: 10.1137/09076636X.

IIF1M(; nlsolve = NLSOLVEJL_SETUP())

IIF1M: Integrating Factor Method 1 (Semi-Linear)

First-order integrating factor method for semi-linear SDEs with stiff linear parts.

Method Properties

• Strong Order: 1.0
• Weak Order: 1.0
• Time stepping: Fixed or adaptive
• Problem type: Semi-linear SDEs with stiff linear components
• Treatment: Exponential integrator approach

Parameters

• nlsolve: Nonlinear solver configuration

When to Use

• Semi-linear SDEs: du = (L*u + N(u))dt + g(u)dW where L is stiff linear operator
• Problems amenable to integrating factor techniques
• When exponential integrators are appropriate
• Stiff linear parts with nonlinear perturbations

Algorithm Description

Applies integrating factor exp(L*t) to handle stiff linear part exactly while treating nonlinear parts numerically.

References

• Integrating factor methods for stiff SDEs

IIF1Mil(; nlsolve = NLSOLVEJL_SETUP())

IIF1Mil: Integrating Factor Milstein Method (Semi-Linear)

Integrating factor method combined with Milstein correction for semi-linear SDEs.

Method Properties

• Strong Order: 1.0 (with Milstein correction)
• Weak Order: 1.0
• Time stepping: Fixed or adaptive
• Problem type: Semi-linear SDEs with stiff linear components
• Treatment: Exponential integrator with Milstein correction

Parameters

• nlsolve: Nonlinear solver configuration

When to Use

• Semi-linear SDEs requiring Milstein-type accuracy
• When both stiff linear treatment and higher-order stochastic accuracy are needed
• Alternative to IIF1M with enhanced stochastic treatment

References

• Integrating factor methods with Milstein correction

IIF2M(; nlsolve = NLSOLVEJL_SETUP())

IIF2M: Integrating Factor Method 2 (Semi-Linear)

Second-order integrating factor method for semi-linear SDEs.

Method Properties

• Strong Order: 2.0
• Weak Order: 2.0
• Time stepping: Fixed or adaptive
• Problem type: Semi-linear SDEs with stiff linear components
• Treatment: Higher-order exponential integrator

Parameters

• nlsolve: Nonlinear solver configuration

When to Use

• When higher accuracy than IIF1M is needed
• Semi-linear problems requiring second-order convergence
• More expensive but more accurate than IIF1M

References

• Higher-order integrating factor methods for SDEs

ISSEM(;
    theta = 1, symplectic = false, linsolve = nothing,
    nlsolve = NLNewton(), extrapolant = :constant,
    new_jac_conv_bound = 1.0e-3,
    autodiff = AutoForwardDiff(), concrete_jac = nothing
)

ISSEM: Implicit Split-Step Euler-Maruyama (Fully Stiff)

Implicit split-step Euler-Maruyama method for Itô SDEs. Unlike ImplicitEM, which is implicit only in the drift, the split-step formulation applies the implicit solve to an intermediate state and then adds the diffusion, which gives stability for problems where the diffusion term is stiff as well.

Method Properties

• Strong Order: 0.5
• Weak Order: 1.0
• Time stepping: Adaptive
• Noise types: All (diagonal, non-diagonal, scalar)
• SDE interpretation: Itô

Mathematical Formulation

The step is split into a purely deterministic implicit solve followed by the diffusion update, so the noise never enters the nonlinear system:

\[u^* = u_n + dt (θ f(u^*, t_{n+1}) + (1 - θ) f(u_n, t_n))\]

\[u_{n+1} = u^* + g(u_n, t_n) ΔW_n\]

When to Use

• SDEs that are stiff in both the drift and the diffusion
• Problems where ImplicitEM still requires very small time steps
• Mean-square stability sensitive problems (e.g. large multiplicative noise)

For Stratonovich problems, use ISSEulerHeun.

Keyword Arguments

• theta: Implicitness parameter of the drift solve (default: 1).
• symplectic: If true, uses the symplectic midpoint discretization (forces theta = 1/2).
• linsolve: LinearSolve.jl algorithm used for the Newton linear systems (default: nothing, i.e. the LinearSolve.jl default).
• nlsolve: Nonlinear solver used for the implicit stage (default: NLNewton()).
• extrapolant: Initial guess method for the nonlinear solve (default: :constant).
• new_jac_conv_bound: Convergence bound below which the Jacobian is reused instead of recomputed (default: 1e-3).
• autodiff: ADTypes.jl backend used for the Jacobian (default: AutoForwardDiff()).
• concrete_jac: Whether to build a concrete Jacobian matrix even when a matrix-free linear solver is used (default: nothing, i.e. decided automatically).

References

• Higham, D.J., Mao, X., Stuart, A.M., "Strong Convergence of Euler-Type Methods for Nonlinear Stochastic Differential Equations", SIAM J. Numer. Anal., 40 (3), pp. 1041–1063 (2002). DOI: 10.1137/S0036142901389530.

ISSEulerHeun(;
    theta = 1, symplectic = false, linsolve = nothing,
    nlsolve = NLNewton(), extrapolant = :constant,
    new_jac_conv_bound = 1.0e-3,
    autodiff = AutoForwardDiff(), concrete_jac = nothing
)

ISSEulerHeun: Implicit Split-Step Euler-Heun (Fully Stiff, Stratonovich)

Implicit split-step Euler-Heun method, the Stratonovich counterpart of ISSEM. The drift is advanced with an implicit solve that contains no noise terms, and the Stratonovich diffusion is then applied with the Heun predictor-corrector, giving stability for problems that are stiff in both the drift and the diffusion.

Method Properties

• Strong Order: 0.5
• Weak Order: 1.0
• Time stepping: Adaptive
• Noise types: All (diagonal, non-diagonal, scalar)
• SDE interpretation: Stratonovich

When to Use

• Stratonovich SDEs that are stiff in both drift and diffusion
• Problems where ImplicitEulerHeun still requires very small time steps
• Physical models with large multiplicative noise in the Stratonovich sense

Keyword Arguments

• theta: Implicitness parameter of the drift solve (default: 1).
• symplectic: If true, uses the symplectic midpoint discretization (forces theta = 1/2).
• linsolve: LinearSolve.jl algorithm used for the Newton linear systems (default: nothing, i.e. the LinearSolve.jl default).
• nlsolve: Nonlinear solver used for the implicit stage (default: NLNewton()).
• extrapolant: Initial guess method for the nonlinear solve (default: :constant).
• new_jac_conv_bound: Convergence bound below which the Jacobian is reused instead of recomputed (default: 1e-3).
• autodiff: ADTypes.jl backend used for the Jacobian (default: AutoForwardDiff()).
• concrete_jac: Whether to build a concrete Jacobian matrix even when a matrix-free linear solver is used (default: nothing, i.e. decided automatically).

References

• Higham, D.J., Mao, X., Stuart, A.M., "Strong Convergence of Euler-Type Methods for Nonlinear Stochastic Differential Equations", SIAM J. Numer. Anal., 40 (3), pp. 1041–1063 (2002). DOI: 10.1137/S0036142901389530.

ImplicitEM(;
    theta = 1, symplectic = false, linsolve = nothing,
    nlsolve = NLNewton(), extrapolant = :constant,
    new_jac_conv_bound = 1.0e-3,
    autodiff = AutoForwardDiff(), concrete_jac = nothing
)

ImplicitEM: Drift-Implicit Euler-Maruyama (Stiff)

Drift-implicit theta-method version of Euler-Maruyama for Itô SDEs. The drift is integrated with the theta method while the diffusion is kept explicit, which gives stability for problems with a stiff drift term.

Method Properties

• Strong Order: 0.5
• Weak Order: 1.0
• Time stepping: Adaptive
• Noise types: All (diagonal, non-diagonal, scalar)
• SDE interpretation: Itô

Mathematical Formulation

\[u_{n+1} = u_n + dt (θ f(u_{n+1}, t_{n+1}) + (1 - θ) f(u_n, t_n)) + g(u_n, t_n) ΔW_n\]

When to Use

• SDEs with a stiff drift term and non-stiff diffusion
• As a robust fallback when explicit methods require impractically small steps
• Problems where only a low-order method is needed but stability is essential

If the diffusion term is also stiff, prefer the split-step method ISSEM. For Stratonovich problems, use ImplicitEulerHeun.

Keyword Arguments

• theta: Implicitness parameter of the drift discretization (default: 1, i.e. drift-implicit Euler). theta = 1/2 gives the trapezoidal rule (STrapezoid).
• symplectic: If true, uses the symplectic midpoint discretization (forces theta = 1/2). See SImplicitMidpoint.
• linsolve: LinearSolve.jl algorithm used for the Newton linear systems (default: nothing, i.e. the LinearSolve.jl default).
• nlsolve: Nonlinear solver used for the implicit stage (default: NLNewton()).
• extrapolant: Initial guess method for the nonlinear solve (default: :constant).
• new_jac_conv_bound: Convergence bound below which the Jacobian is reused instead of recomputed (default: 1e-3).
• autodiff: ADTypes.jl backend used for the Jacobian (default: AutoForwardDiff()).
• concrete_jac: Whether to build a concrete Jacobian matrix even when a matrix-free linear solver is used (default: nothing, i.e. decided automatically).

References

• Kloeden, P.E., Platen, E., "Numerical Solution of Stochastic Differential Equations", Springer (1992).
• Milstein, G.N., "Numerical Integration of Stochastic Differential Equations", Kluwer (1995).

ImplicitEulerHeun(;
    theta = 1, symplectic = false, linsolve = nothing,
    nlsolve = NLNewton(), extrapolant = :constant,
    new_jac_conv_bound = 1.0e-3,
    autodiff = AutoForwardDiff(), concrete_jac = nothing
)

ImplicitEulerHeun: Drift-Implicit Euler-Heun (Stiff, Stratonovich)

Drift-implicit theta-method version of the Euler-Heun scheme, the Stratonovich counterpart of ImplicitEM. The drift is treated implicitly while the Stratonovich diffusion is handled with the explicit Heun predictor-corrector.

Method Properties

• Strong Order: 0.5
• Weak Order: 1.0
• Time stepping: Adaptive
• Noise types: All (diagonal, non-diagonal, scalar)
• SDE interpretation: Stratonovich

When to Use

• Stratonovich SDEs with a stiff drift term
• Physical models where the noise is a smooth approximation limit (Wong-Zakai), which is naturally Stratonovich
• As the Stratonovich analogue of ImplicitEM()

If the diffusion is also stiff, prefer the split-step method ISSEulerHeun.

Keyword Arguments

• theta: Implicitness parameter of the drift discretization (default: 1, i.e. drift-implicit Euler). theta = 1/2 gives the trapezoidal rule.
• symplectic: If true, uses the symplectic midpoint discretization (forces theta = 1/2).
• linsolve: LinearSolve.jl algorithm used for the Newton linear systems (default: nothing, i.e. the LinearSolve.jl default).
• nlsolve: Nonlinear solver used for the implicit stage (default: NLNewton()).
• extrapolant: Initial guess method for the nonlinear solve (default: :constant).
• new_jac_conv_bound: Convergence bound below which the Jacobian is reused instead of recomputed (default: 1e-3).
• autodiff: ADTypes.jl backend used for the Jacobian (default: AutoForwardDiff()).
• concrete_jac: Whether to build a concrete Jacobian matrix even when a matrix-free linear solver is used (default: nothing, i.e. decided automatically).

References

• Kloeden, P.E., Platen, E., "Numerical Solution of Stochastic Differential Equations", Springer (1992).

ImplicitRKMil(;
    theta = 1, symplectic = false, linsolve = nothing,
    nlsolve = NLNewton(), extrapolant = :constant,
    new_jac_conv_bound = 1.0e-3,
    interpretation = AlgorithmInterpretation.Ito,
    autodiff = AutoForwardDiff(), concrete_jac = nothing
)

ImplicitRKMil: Drift-Implicit Runge-Kutta Milstein (Stiff)

Drift-implicit theta-method version of the derivative-free Runge-Kutta Milstein scheme. The Milstein correction gives strong order 1.0, twice the strong order of ImplicitEM, at the cost of being restricted to diagonal and scalar noise.

Method Properties

• Strong Order: 1.0
• Weak Order: 1.0
• Time stepping: Adaptive
• Noise types: Diagonal and scalar noise only
• SDE interpretation: Itô or Stratonovich (selected by interpretation)

When to Use

• Stiff SDEs with diagonal or scalar noise where strong order 1.0 is needed
• When ImplicitEM converges too slowly in the strong sense
• Pathwise-accurate simulations (e.g. multilevel Monte Carlo) of stiff problems

Note that ImplicitRKMil is only compatible with diagonal-noise problems; for non-diagonal noise use ImplicitEM or ISSEM.

Keyword Arguments

• theta: Implicitness parameter of the drift discretization (default: 1, i.e. drift-implicit Euler). theta = 1/2 gives the trapezoidal rule.
• symplectic: If true, uses the symplectic midpoint discretization (forces theta = 1/2).
• interpretation: AlgorithmInterpretation.Ito (default) or AlgorithmInterpretation.Stratonovich, selecting which SDE interpretation the Milstein correction is built for.
• linsolve: LinearSolve.jl algorithm used for the Newton linear systems (default: nothing, i.e. the LinearSolve.jl default).
• nlsolve: Nonlinear solver used for the implicit stage (default: NLNewton()).
• extrapolant: Initial guess method for the nonlinear solve (default: :constant).
• new_jac_conv_bound: Convergence bound below which the Jacobian is reused instead of recomputed (default: 1e-3).
• autodiff: ADTypes.jl backend used for the Jacobian (default: AutoForwardDiff()).
• concrete_jac: Whether to build a concrete Jacobian matrix even when a matrix-free linear solver is used (default: nothing, i.e. decided automatically).

References

• Milstein, G.N., "Numerical Integration of Stochastic Differential Equations", Kluwer (1995).
• Kloeden, P.E., Platen, E., "Numerical Solution of Stochastic Differential Equations", Springer (1992).

SImplicitMidpoint(; kwargs...)

SImplicitMidpoint: Stochastic Implicit Midpoint Rule (Stiff, Symplectic)

Convenience constructor for the symplectic midpoint member of the ImplicitEM family, i.e. ImplicitEM(theta = 1/2, symplectic = true). The midpoint discretization of the drift makes the deterministic part of the integrator symplectic, so it preserves quadratic invariants and has good long-time energy behavior.

Method Properties

• Strong Order: 0.5
• Weak Order: 1.0
• Time stepping: Adaptive
• Noise types: All (diagonal, non-diagonal, scalar)
• SDE interpretation: Itô
• Symplectic: Yes (deterministic part)

When to Use

• Stiff Hamiltonian or otherwise structure-preserving SDEs
• Long-time integrations where drift in conserved quantities must be avoided
• When the non-damping behavior of the midpoint rule is preferable to the dissipativity of backward Euler

Keyword Arguments

All keyword arguments of ImplicitEM are accepted and forwarded, except that theta is fixed to 1/2 and symplectic is fixed to true.

References

• Milstein, G.N., Repin, Yu.M., Tretyakov, M.V., "Numerical Methods for Stochastic Systems Preserving Symplectic Structure", SIAM J. Numer. Anal., 40 (4), pp. 1583–1604 (2002). DOI: 10.1137/S0036142901395588.

SKenCarp(;
    linsolve = nothing, nlsolve = NLNewton(),
    smooth_est = true, extrapolant = :min_correct,
    new_jac_conv_bound = 1.0e-3, ode_error_est = true,
    autodiff = AutoForwardDiff(), concrete_jac = nothing
)

SKenCarp: Stochastic KenCarp for Additive Noise (Stiff, Recommended)

Stochastic extension of the KenCarp4 additive Runge-Kutta (ESDIRK) method for SDEs with additive noise. The deterministic part is integrated with the L-stable, stiffly-accurate KenCarp4 tableau, while the additive noise is handled with the higher-order stochastic quadrature of Rößler's SRA schemes, so the method retains strong order 1.5 on additive-noise problems.

Method Properties

• Strong Order: 1.5 (for additive noise)
• Deterministic Order: 4 (KenCarp4)
• Time stepping: Adaptive
• Noise types: Additive noise (diagonal, non-diagonal, and scalar)
• SDE interpretation: Itô (Itô and Stratonovich coincide for additive noise)
• Stability: L-stable, stiffly accurate

When to Use

• The recommended default for stiff SDEs with additive noise
• Problems where the deterministic dynamics dominate and need high accuracy
• Chemical, mechanical, or PDE-discretization models driven by additive forcing

Additive noise means the diffusion function does not depend on the state (g(u,p,t) is independent of u); for state-dependent noise use ImplicitEM, ImplicitRKMil, or ISSEM.

Keyword Arguments

• smooth_est: Whether to use the smoothed error estimate, which is more robust for stiff problems (default: true).
• ode_error_est: Whether the deterministic (ODE) embedded error estimate is included in the adaptivity controller (default: true).
• linsolve: LinearSolve.jl algorithm used for the Newton linear systems (default: nothing, i.e. the LinearSolve.jl default).
• nlsolve: Nonlinear solver used for the implicit stages (default: NLNewton()).
• extrapolant: Initial guess method for the nonlinear solves (default: :min_correct).
• new_jac_conv_bound: Convergence bound below which the Jacobian is reused instead of recomputed (default: 1e-3).
• autodiff: ADTypes.jl backend used for the Jacobian (default: AutoForwardDiff()).
• concrete_jac: Whether to build a concrete Jacobian matrix even when a matrix-free linear solver is used (default: nothing, i.e. decided automatically).

References

• Kennedy, C.A., Carpenter, M.H., "Additive Runge-Kutta schemes for convection-diffusion-reaction equations", Applied Numerical Mathematics, 44 (1-2), pp. 139–181 (2003). DOI: 10.1016/S0168-9274(02)00138-1.
• Rackauckas, C., Nie, Q., "Stability-Optimized High Order Methods and Stiffness Detection for Pathwise Stiff Stochastic Differential Equations", 2020 IEEE High Performance Extreme Computing Conference (HPEC). DOI: 10.1109/HPEC43674.2020.9286178.
• Rößler A., "Runge-Kutta Methods for the Strong Approximation of Solutions of Stochastic Differential Equations", SIAM J. Numer. Anal., 48 (3), pp. 922–952 (2010). DOI: 10.1137/09076636X.

STrapezoid(; kwargs...)

STrapezoid: Stochastic Trapezoidal Rule (Stiff)

Convenience constructor for the trapezoidal-rule member of the ImplicitEM family, i.e. ImplicitEM(theta = 1/2). The drift is discretized with the trapezoidal rule and the diffusion is kept explicit.

Method Properties

• Strong Order: 0.5
• Weak Order: 1.0
• Time stepping: Adaptive
• Noise types: All (diagonal, non-diagonal, scalar)
• SDE interpretation: Itô

When to Use

• Stiff Itô SDEs where the second-order accurate (in the deterministic part) trapezoidal drift discretization is preferred over backward Euler
• Problems where backward Euler's numerical damping is undesirable

Note that the trapezoidal rule is only A-stable, not L-stable, so for very stiff drifts the default ImplicitEM() (theta = 1) may be more robust.

Keyword Arguments

All keyword arguments of ImplicitEM are accepted and forwarded, except that theta is fixed to 1/2.

CaoTauLeaping()

CaoTauLeaping: Cao's Adaptive Tau-Leaping Method (Jump-Diffusion)

Advanced tau-leaping method with adaptive tau selection and improved error control.

Method Properties

• Problem type: Jump-diffusion processes
• Approach: Adaptive tau selection with error control
• Time stepping: Adaptive tau based on error estimates
• Accuracy: Superior to basic tau-leaping

When to Use

• Production jump-diffusion simulations requiring reliability
• When adaptive tau selection is needed
• Problems where basic TauLeaping gives poor accuracy
• Chemical reaction networks requiring precise control

Algorithm Features

• Adaptive tau selection based on error estimates
• Better stability and accuracy than basic tau-leaping
• Automatic step size control
• More sophisticated error estimation

Tau Selection

Automatically adjusts tau based on:

• Local error estimates
• Jump rate variations
• Solution stability requirements

References

• Cao, Y., Gillespie, D.T., Petzold, L.R., "Efficient step size selection for the tau-leaping method"

ImplicitTauLeaping(; nlsolve = NLFunctional())

ImplicitTauLeaping: First Order Implicit Tau-Leaping Method (Jump-Diffusion)

An implicit (backward Euler) tau-leaping method for stiff chemical kinetic systems. Uses backward Euler discretization to provide improved stability for systems with fast reversible reactions or stiff rate constants.

Method Properties

Mathematical Formulation

The method solves the implicit equation:

\[X_{n+1} = X_n + ν ⋅ Poisson(dt ⋅ a(X_{n+1}))\]

which is approximated by:

\[X_{n+1} = X_n + ν ⋅ k + dt ⋅ (drift(X_{n+1}) - drift(X_n))\]

where k ~ Poisson(dt * a(X_n)) and drift(u) = ν * a(u).

This corresponds to ThetaTrapezoidalTauLeaping with θ = 1 (fully implicit).

Keyword Arguments

• nlsolve: Nonlinear solver algorithm (default: NLFunctional()). Options include NLFunctional(), NLAnderson(), and NLNewton().

Example

using StochasticDiffEq, JumpProcesses

# Define rate function and stoichiometry
rate(out, u, p, t) = (out[1] = 0.1 * u[1]; out[2] = 0.05 * u[2])
c(du, u, p, t, counts, mark) = (du[1] = -counts[1] + counts[2]; du[2] = counts[1] - counts[2])

rj = RegularJump(rate, c, 2)
prob = DiscreteProblem([100.0, 0.0], (0.0, 10.0))
jprob = JumpProblem(prob, Direct(), rj)

sol = solve(jprob, ImplicitTauLeaping(); dt = 0.1)

References

• Rathinam, M., Petzold, L.R., Cao, Y., Gillespie, D.T., "Stiffness in stochastic chemically reacting systems: The implicit tau-leaping method", J. Chem. Phys. 119, 12784 (2003)

TauLeaping()

TauLeaping: Basic Tau-Leaping Method (Jump-Diffusion)

Basic tau-leaping method for approximating jump-diffusion processes by "leaping" over multiple potential jump events.

Method Properties

• Problem type: Jump-diffusion processes
• Approach: Approximate multiple jumps per time step
• Time stepping: Fixed tau approach
• Accuracy: Depends on tau selection

When to Use

• Jump-diffusion systems with many small jumps
• When exact jump simulation is computationally prohibitive
• Chemical reaction networks with fast reactions
• Population models with high birth-death rates
• Initial exploration of jump-diffusion problems

Algorithm Description

Approximates Poisson processes by assuming constant propensities over time interval tau, then sampling number of jumps from Poisson distribution.

Tau Selection

Critical parameter: tau should be small enough that jump rates don't change significantly over [t, t+tau].

References

• Gillespie, D.T., "Approximate accelerated stochastic simulation of chemically reacting systems"

ThetaTrapezoidalTauLeaping(; theta = 0.5, max_iters = 10, abstol = 1.0e-8, reltol = 1.0e-6)

ThetaTrapezoidalTauLeaping: Implicit Weak Second Order Tau-Leaping Method (Jump-Diffusion)

An implicit tau-leaping method achieving weak second order accuracy in the large volume scaling. Uses fixed-point iteration to solve the implicit equation. Based on the work of Hu, Li, and Min (2011) and Anderson and Mattingly (2011).

Method Properties

• Problem type: Jump-diffusion processes
• Order: Weak order 2 (in the large volume scaling)
• Time stepping: Fixed or adaptive tau
• Accuracy: Superior to both Euler tau-leaping and midpoint tau-leaping
• Implicit treatment: Uses nonlinear solver for implicit rate equations

Parameters

• theta::Float64: Implicitness parameter (default: 0.5)
  • Must be in range (0, 1)
  • theta = 0.5 gives trapezoidal method (recommended for balanced accuracy/stability)
  • theta = 1.0 gives backward Euler (maximum stability)
• max_iters::Int: Maximum iterations for nonlinear solver (default: 10)
• abstol::Float64: Absolute tolerance for convergence (default: 1e-8)
• reltol::Float64: Relative tolerance for convergence (default: 1e-6)

When to Use

• When higher accuracy is needed compared to standard tau-leaping methods
• Chemical reaction networks requiring weak second order accuracy
• Systems where accurate mean and covariance estimates are important
• When both Euler and midpoint tau-leaping provide insufficient accuracy
• Stiff chemical systems where implicit treatment provides stability

Algorithm Description

The method solves the implicit equation:

\[X_{n+1} = X_n + ν⋅k + θ⋅dt⋅ν⋅(a(X_{n+1}) - a(X_n))\]

where k ~ Poisson(dt⋅a(X_n)) are the jump counts.

This is solved using fixed-point iteration:

1. Generate Poisson jumps with rate dt·a(X_n)
2. Initialize z = 0
3. Iterate: z_{new} = θ·dt·(drift(X_n + ν·k + z) - drift(X_n))
4. Final state: X_{n+1} = X_n + ν·k + z

Convergence Properties

The local truncation error for covariance is O(τ³V⁻¹) when τ = V^(-β) for 0 < β < 1 and system size V → ∞, which is higher order than both Euler and midpoint methods.

References

• Hu, Y., Li, T., Min, B., "A weak second order tau-leaping method for chemical kinetic systems", J. Chem. Phys. 135, 024113 (2011)
• Anderson, D.F., Mattingly, J.C., "A weak trapezoidal method for a class of stochastic differential equations", Comm. Math. Sci. 9, 301 (2011)

EM(; split = true)

Select the fixed-step Euler-Maruyama method for an Itô SDE.

Keywords

• split = true: Whether to use split-step handling.

Examples

sol = solve(prob, EM(); dt = 0.01)

Use split = false when the problem should be advanced without split-step handling.

EulerHeun()

EulerHeun: Stochastic Euler-Heun Method

A two-stage predictor-corrector (Heun-type) method for Stratonovich SDEs, the Stratonovich analogue of Euler-Maruyama. An Euler step forms the predictor ũ, then the drift and diffusion are re-evaluated there and trapezoidally averaged to form the update:

• predictor: ũ = uₙ + f(uₙ)·Δt + g(uₙ)·ΔW
• corrector: uₙ₊₁ = uₙ + ½(f(uₙ) + f(ũ))·Δt + ½(g(uₙ) + g(ũ))·ΔW

This corresponds to the improved-Euler scheme of Roberts (2012). Note this is a genuine two-stage scheme (two drift and two diffusion evaluations per step), not a single-stage predictor-corrector.

Method Properties

• Strong Order: 1.0 for Stratonovich SDEs with commutative noise (e.g. scalar, diagonal, or additive); 1/2 for general non-commutative noise, which is the value reported by alg_order (the same convention as EM)
• Weak Order: 1.0
• Time stepping: Fixed step size
• Noise types: General (scalar, diagonal, non-diagonal)
• SDE interpretation: Stratonovich

When to Use

• For Stratonovich SDEs where a simple, low-cost method is sufficient
• As the Stratonovich counterpart to EM for Itô problems
• When adaptive time stepping is not required; see LambaEulerHeun for an adaptive variant

References

• Roberts, A.J., "Modify the improved Euler scheme to integrate stochastic differential equations", arXiv:1210.0933 (2012)
• Kloeden, P.E., Platen, E., Numerical Solution of Stochastic Differential Equations, Springer, Berlin Heidelberg, p. 373 (1992)

LambaEM(; split = true)

Select the adaptive Lamba Euler-Maruyama method for an Itô SDE.

Keywords

• split = true: Whether to use split-step handling.

Examples

sol = solve(prob, LambaEM(); reltol = 1.0e-3, abstol = 1.0e-3)

LambaEulerHeun()

Select the adaptive Lamba Euler-Heun method for a Stratonovich SDE.

Examples

sol = solve(prob, LambaEulerHeun(); reltol = 1.0e-3, abstol = 1.0e-3)

PCEuler(ggprime; theta = 1 / 2, eta = 1 / 2)

Select a fixed-step predictor-corrector Euler method for an Itô SDE, supplying the derivative of the diffusion coefficient as ggprime.

Arguments

• ggprime: Derivative of the diffusion coefficient.

Keywords

• theta = 1 / 2: Drift evaluation weight.
• eta = 1 / 2: Diffusion evaluation weight.

Examples

sol = solve(prob, PCEuler(ggprime); dt = 0.01)

RKMil(; interpretation = SciMLBase.AlgorithmInterpretation.Ito)

Select the adaptive Runge-Kutta Milstein method for scalar or diagonal-noise SDEs.

Keywords

• interpretation = SciMLBase.AlgorithmInterpretation.Ito: Stochastic interpretation used by the method.

Examples

sol = solve(prob, RKMil(); reltol = 1.0e-3, abstol = 1.0e-3)

RKMilCommute(;
    interpretation = SciMLBase.AlgorithmInterpretation.Ito,
    ii_approx = IICommutative()
)

Select the adaptive Runge-Kutta Milstein method for commutative-noise SDEs.

Fields

• interpretation: Stochastic interpretation used by the method.
• ii_approx: Iterated-integral approximation for commutative noise.

Keywords

• interpretation = SciMLBase.AlgorithmInterpretation.Ito: Stochastic interpretation used by the method.
• ii_approx = IICommutative(): Iterated-integral approximation for commutative noise.

Examples

sol = solve(prob, RKMilCommute(); reltol = 1.0e-3, abstol = 1.0e-3)

SimplifiedEM()

Select the simplified fixed-step Euler-Maruyama method for weak SDE approximations.

Examples

sol = solve(prob, SimplifiedEM(); dt = 0.01)

SplitEM()

Select the fixed-step split Euler-Maruyama method for an Itô SDE with a splittable drift.

Examples

sol = solve(prob, SplitEM(); dt = 0.01)

RKMilGeneral(; interpretation = Ito, ii_approx = IILevyArea(), c = 1, p = nothing, dt = nothing)

RKMilGeneral: Generalized Runge-Kutta Milstein Method

Adaptive step-size Milstein method supporting general (non-diagonal, non-commutative) noise via iterated stochastic integrals (Levy area approximation).

Method Properties

• Strong Order: 1.0
• Weak Order: 1.0
• Time stepping: Adaptive step size
• Noise types: General (scalar, diagonal, non-diagonal, commutative)
• SDE interpretation: Configurable (Ito or Stratonovich)

Keyword Arguments

• interpretation: SciMLBase.AlgorithmInterpretation.Ito (default) or SciMLBase.AlgorithmInterpretation.Stratonovich
• ii_approx: Iterated integral approximation method (default: IILevyArea())
• c: Constant in the Levy area accuracy target ε = c dt^(3/2) (default: 1)
• p: Number of Fourier terms retained in the Levy area series. nothing (the default) picks it from the noise dimension and the solver's step size so that the truncation error meets ε; pass an Integer to fix it instead.
• dt: Ignored, accepted for backwards compatibility.

When to Use

• For SDEs with general (non-commutative) noise requiring strong order 1.0
• When adaptive step size is needed with Milstein-type accuracy
• When Levy area computation is acceptable for accuracy improvement
• For non-diagonal noise SDEs where RKMil/RKMilCommute are not applicable

Performance Notes

• Levy area computation scales with noise dimension
• Adaptive truncation balances accuracy and efficiency

References

• Kloeden, P.E., Platen, E., Numerical Solution of Stochastic Differential Equations. Springer. Berlin Heidelberg (2011)
• Kastner, F. and Roessler, A., "LevyArea.jl: A Julia package for Levy area computation", arXiv:2201.08424
• LevyArea.jl: https://github.com/stochastics-uni-luebeck/LevyArea.jl

WangLi3SMil_A()

WangLi3SMil_A: 3-Stage Milstein Method A (Nonstiff)

Fixed step-size explicit 3-stage Milstein method with strong and weak order 1.0 for Ito SDEs.

Method Properties

• Strong Order: 1.0
• Weak Order: 1.0
• Time stepping: Fixed step size
• Noise types: Depends on tableau (typically diagonal/scalar)
• SDE interpretation: Ito

When to Use

• When fixed step size is preferred
• For Ito SDEs requiring order 1.0 accuracy
• Part of WangLi family - compare performance with other variants
• When computational cost per step is less important than simplicity

References

• Wang and Li, "Three-stage stochastic Runge-Kutta methods for stochastic differential equations"

WangLi3SMil_B()

WangLi3SMil_B: 3-Stage Milstein Method B (Nonstiff)

Alternative 3-stage Milstein method with different stability and accuracy characteristics.

Method Properties

• Strong Order: 1.0
• Weak Order: 1.0
• Time stepping: Fixed step size
• Noise types: Depends on tableau (typically diagonal/scalar)
• SDE interpretation: Ito

References

• Wang and Li, "Three-stage stochastic Runge-Kutta methods for stochastic differential equations"

WangLi3SMil_C()

WangLi3SMil_C: 3-Stage Milstein Method C (Nonstiff)

Third variant in the WangLi 3-stage Milstein family.

Method Properties

• Strong Order: 1.0
• Weak Order: 1.0
• Time stepping: Fixed step size
• Noise types: Depends on tableau (typically diagonal/scalar)
• SDE interpretation: Ito

References

• Wang and Li, "Three-stage stochastic Runge-Kutta methods for stochastic differential equations"

WangLi3SMil_D()

WangLi3SMil_D: 3-Stage Milstein Method D (Nonstiff)

Fourth variant in the WangLi 3-stage Milstein family.

Method Properties

• Strong Order: 1.0
• Weak Order: 1.0
• Time stepping: Fixed step size
• Noise types: Depends on tableau (typically diagonal/scalar)
• SDE interpretation: Ito

References

• Wang and Li, "Three-stage stochastic Runge-Kutta methods for stochastic differential equations"

WangLi3SMil_E()

WangLi3SMil_E: 3-Stage Milstein Method E (Nonstiff)

Fifth variant in the WangLi 3-stage Milstein family.

Method Properties

• Strong Order: 1.0
• Weak Order: 1.0
• Time stepping: Fixed step size
• Noise types: Depends on tableau (typically diagonal/scalar)
• SDE interpretation: Ito

References

• Wang and Li, "Three-stage stochastic Runge-Kutta methods for stochastic differential equations"

WangLi3SMil_F()

WangLi3SMil_F: 3-Stage Milstein Method F (Nonstiff)

Sixth and final variant in the WangLi 3-stage Milstein family.

Method Properties

• Strong Order: 1.0
• Weak Order: 1.0
• Time stepping: Fixed step size
• Noise types: Depends on tableau (typically diagonal/scalar)
• SDE interpretation: Ito

When to Use (WangLi Family)

• Compare all variants (A-F) to find best performance for your problem
• Fixed step applications where step size is predetermined
• Benchmarking against adaptive methods
• When Milstein accuracy is needed with explicit fixed steps

References

• Wang and Li, "Three-stage stochastic Runge-Kutta methods for stochastic differential equations"

KomBurSROCK2(; eigen_est = nothing)

KomBurSROCK2: Komori-Burrage Second-Order SROCK Method

Alternative second-order stabilized method with different coefficients and stability properties.

Method Properties

• Strong Order: 1.0
• Weak Order: 2.0
• Time stepping: Fixed step size with extended stability
• Stability: Extended along negative real axis

When to Use

• Alternative to SROCK2 with different stability characteristics
• When SROCK2 performance is unsatisfactory
• Benchmarking against other second-order ROCK methods

References

• Komori and Burrage stabilized methods

SKSROCK(; post_processing = false, eigen_est = nothing)

SKSROCK: SK-SROCK Stabilized Method

Fixed step stabilized explicit method for stiff Ito problems with enhanced stability domain and optional post-processing.

Method Properties

• Strong Order: 0.5 (up to 2.0 with post-processing)
• Weak Order: 1.0 (up to 2.0 with post-processing)
• Time stepping: Fixed step size with enhanced stability
• Noise types: 1-dimensional, diagonal, and multi-dimensional noise
• SDE interpretation: Ito only
• Stability: Better stability domain than SROCK1

Parameters

• post_processing::Bool = false: Enable post-processing for higher accuracy (experimental)
• eigen_est: Eigenvalue estimation for stability (automatic if nothing)

When to Use

• Stiff Ito problems requiring better stability than SROCK1
• Ergodic dynamical systems (with post-processing)
• Problems where enhanced stability domain is crucial
• When experimenting with post-processing techniques

Post-Processing (Experimental)

• Can achieve order 2 accuracy for ergodic systems
• Particularly useful for Brownian dynamics
• Currently under development - use with caution

Algorithm Features

• Enhanced stability compared to SROCK1
• Handles various noise structures
• Optional post-processing for specialized applications

References

• SK-SROCK methods for stochastic problems
• Post-processing techniques for ergodic systems

SROCK1(;
    interpretation = SciMLBase.AlgorithmInterpretation.Ito,
    eigen_est = nothing
)

Select the first-order stabilized Runge-Kutta Chebyshev method for a mildly stiff SDE.

Keywords

• interpretation = SciMLBase.AlgorithmInterpretation.Ito: Stochastic interpretation used by the method.
• eigen_est = nothing: Optional estimate of the stiff eigenvalues.

Examples

sol = solve(prob, SROCK1(); dt = 0.01)

Set eigen_est when a problem-specific estimate of the stiff eigenvalues is available.

SROCK2(; eigen_est = nothing)

SROCK2: Second-Order Stabilized Runge-Kutta Chebyshev Method

Second-order stabilized explicit method with weak order 2.0 for mildly stiff SDE problems.

Method Properties

• Strong Order: 1.0
• Weak Order: 2.0
• Time stepping: Fixed step size with extended stability
• Stability: Extended along negative real axis

When to Use

• When higher accuracy than SROCK1 is needed
• Parabolic PDEs requiring better weak convergence
• Problems where second-order accuracy justifies increased cost

References

• Second-order ROCK methods for stochastic problems

SROCKC2(; eigen_est = nothing)

SROCKC2: Conservative Second-Order SROCK Method

Conservative second-order stabilized method designed for robust performance.

Method Properties

• Strong Order: 1.0
• Weak Order: 2.0
• Time stepping: Fixed step size with extended stability
• Stability: Conservative stability region, more robust

When to Use

• When robustness is more important than efficiency
• For difficult problems where other ROCK methods fail
• As a fallback option for problematic cases

References

• Conservative ROCK methods for stochastic problems

SROCKEM(; strong_order_1 = true, eigen_est = nothing)

SROCKEM: ROCK-Stabilized Euler-Maruyama Method

Fixed step Euler-Maruyama method with first-order ROCK stabilization for handling stiff problems.

Method Properties

• Strong Order: 1.0 (default) or 0.5 (if strong_order_1=false)
• Weak Order: 1.0 (default) or 0.5 (if strong_order_1=false)
• Time stepping: Fixed step size with ROCK stabilization
• Noise types: 1-dimensional, diagonal, and multi-dimensional noise
• SDE interpretation: Ito only
• Stability: ROCK stabilization for moderate stiffness

Parameters

• strong_order_1::Bool = true: Use strong/weak order 1.0 (true) or 0.5 (false)
• eigen_est: Eigenvalue estimation for stability (automatic if nothing)

When to Use

• Stiff problems where standard EM fails
• When ROCK stabilization is preferred over full implicit treatment
• Problems requiring Euler-Maruyama structure with enhanced stability
• Multi-dimensional stiff SDEs

Algorithm Description

Combines Euler-Maruyama discretization with ROCK stabilization techniques to extend the stability region without requiring linear solves.

References

• ROCK stabilization techniques applied to SDEs
• Stabilized Euler methods for stiff problems

TangXiaoSROCK2(; version_num = 5, eigen_est = nothing)

TangXiaoSROCK2: Tang-Xiao Second-Order SROCK Method

Fixed step size stabilized explicit method with multiple variants offering different stability domains.

Method Properties

• Strong Order: 1.0
• Weak Order: 2.0
• Time stepping: Fixed step size with extended stability
• Noise types: Various (depends on version)
• SDE interpretation: Ito only
• Stability: Version-dependent stability domains

Parameters

• version_num::Int = 5: Choose version 1-5 with different stability characteristics
• eigen_est: Eigenvalue estimation for stability (automatic if nothing)

When to Use

• When experimenting with different stability domains
• Problems requiring weak order 2.0 with fixed steps
• Benchmarking different ROCK variants
• Note: Currently under development

Versions

• Versions 1-5 offer different stability domains
• Version 5 (default) typically provides good general performance
• Choose version based on problem-specific stability requirements

Development Status

• Method is under active development
• Use with caution in production code
• Consider more established ROCK methods for critical applications

References

• Tang and Xiao, "Second-order SROCK methods for stochastic problems"

BAOAB(; gamma = 1.0, scale_noise = true)

BAOAB: Langevin Dynamics Integrator (Specialized)

Specialized integrator for Langevin dynamics in molecular dynamics simulations, particularly effective for configurational sampling.

Method Properties

• Problem type: Langevin dynamics (second-order SDEs)
• Structure: Position-velocity formulation
• Sampling: Designed for equilibrium sampling
• Time stepping: Fixed step size
• Conservation: Preserves equilibrium distributions

Parameters

• gamma::Real = 1.0: Friction coefficient
• scale_noise::Bool = true: Whether to scale noise appropriately

System Structure

Designed for Langevin systems:

\[\begin{align*} du &= v \, dt \\ dv &= f(v,u) \, dt - γv \, dt + g(u) \sqrt{2γ} \, dW \end{align*}\]

where:

• $u$: position coordinates
• $v$: velocity coordinates
• $γ$: friction coefficient
• $f(v,u)$: force function
• $g(u)$: noise scaling function

When to Use

• Molecular dynamics simulations with Langevin thermostat
• Configurational sampling of molecular systems
• Equilibrium sampling from canonical ensemble
• Second-order SDEs with damping and noise

Algorithm Features

• BAOAB splitting: B(kick) - A(drift) - O(Ornstein-Uhlenbeck) - A(drift) - B(kick)
• Preserves correct equilibrium distribution
• Robust and efficient for molecular sampling
• Well-suited for long-time integration

References

• Leimkuhler B., Matthews C., "Robust and efficient configurational molecular sampling via Langevin dynamics", J. Chem. Phys. 138, 174102 (2013)

RandomEM()

RandomEM: Random Euler Method (RODE)

Euler method for Random Ordinary Differential Equations (RODEs) with random parameters.

Method Properties

• Problem type: Random ODEs (RODEs)
• Strong Order: 1.0 (for deterministic part)
• Randomness: Handles random parameters, not Brownian motion
• Time stepping: Fixed step size

When to Use

• Random ODEs with random parameters but no Brownian motion
• Uncertainty quantification with parameter randomness
• Problems with random coefficients or initial conditions
• Monte Carlo simulation of deterministic systems with random inputs

RODE vs SDE

• RODE: Random parameters, deterministic evolution
• SDE: Fixed parameters, stochastic (Brownian) evolution

References

• Random ordinary differential equation methods

RandomHeun()

RandomHeun: Random Heun Method (RODE)

Heun method for Random Ordinary Differential Equations with improved accuracy.

Method Properties

• Problem type: Random ODEs (RODEs)
• Strong Order: 2.0 (for deterministic part)
• Randomness: Handles random parameters
• Time stepping: Fixed step size

When to Use

• RODEs requiring higher accuracy than RandomEM
• When computational cost per step is acceptable
• Random parameter problems needing second-order accuracy

References

• Higher-order methods for random ODEs

RandomTamedEM()

RandomTamedEM: Tamed Random Euler Method (RODE)

Tamed Euler method for RODEs with potentially explosive behavior.

Method Properties

• Problem type: Random ODEs with potential blow-up
• Approach: Taming to prevent numerical explosion
• Stability: Enhanced stability for unstable random systems
• Time stepping: Fixed step size with taming

When to Use

• RODEs that may exhibit explosive growth
• When RandomEM gives unstable or explosive solutions
• Random systems with strong nonlinearities
• Problems requiring enhanced numerical stability

Taming Mechanism

Applies taming technique to prevent numerical blow-up while maintaining accuracy for well-behaved solutions.

References

• Tamed methods for random differential equations