Module: svIntegratorWeakStochasticRungeKutta

template<size_t numberStages> struct SRKCoefficients

#include <svIntegratorWeakStochasticRungeKutta.h>

Stores the coefficients necessary to use the Stochastic Runge-Kutta method.

The Butcher table looks like:

c0 | A0 | B0 | c1 | A1 | B1 | B2 -—|——————— | alpha | beta0 | beta1

See more in the description of svIntegratorWeakStochasticRungeKutta.

Public Types

using StageSizedArray = std::array<double, numberStages>: Array with size = numberStages

using StageSizedMatrix = std::array<StageSizedArray, numberStages>: Square matrix with size = numberStages

Public Members

StageSizedArray alpha = {}: “alpha” coefficient array of the SRK method

StageSizedArray beta0 = {}: “beta0” coefficient array of the SRK method

StageSizedArray beta1 = {}: “beta1” coefficient array of the SRK method

StageSizedMatrix A0 = {}: “A0” coefficient matrix of the SRK method

StageSizedMatrix B0 = {}: “B0” coefficient matrix of the SRK method

StageSizedMatrix A1 = {}: “A1” coefficient matrix of the SRK method

StageSizedMatrix B1 = {}: “B1” coefficient matrix of the SRK method

StageSizedMatrix B2 = {}: “B2” coefficient matrix of the SRK method

StageSizedArray c0 = {}

“c0” coefficient array of the SRK method

c0 is the row sum of the A(0) matrix:

\[ c_0 = \sum_j a^{(0)}_{ij} \]

StageSizedArray c1 = {}

“c1” coefficient array of the SRK method

c1 is the row sum of the A(1) matrix:

\[ c_1 = \sum_j a^{(1)}_{ij} \]

struct SRKRandomVariables

#include <svIntegratorWeakStochasticRungeKutta.h>

Random variables used in the Weak Stochastic Integrator

Public Members

Eigen::VectorXd Ik: See Eq 3.2 in Tang & Xiao

Eigen::MatrixXd Ikl: See Eq 3.3 in Tang & Xiao

Eigen::VectorXd Ikk: Same as Ikl(k,k)

double xi: See Eq 3.2 in Tang & Xiao

class SRKRandomVariableGenerator

#include <svIntegratorWeakStochasticRungeKutta.h>

Random Generator for the integrator described in:

Tang, X., Xiao, A. Efficient weak second-order stochastic Runge–Kutta methods for Itô stochastic differential equations. Bit Numer Math 57, 241–260 (2017). https://doi.org/10.1007/s10543-016-0618-9

Public Functions

inline void setSeed(size_t seed): Sets the seed to the RNG

inline SRKRandomVariables generate(size_t m, double h)

Generates the random values necessary for one step of the integrator.

Parameters:

m – number of noise sources
h – time step, in seconds

Protected Attributes

std::mt19937 rng = {std::random_device{}()}: Random Number Generator

std::uniform_real_distribution<double> uniform_rv: Uniform random variable, from 0 to 1

std::bernoulli_distribution bernoulli_rv: Bernoulli random variable (either 0 or 1, with equal probability)

template<size_t numberStages> class svIntegratorWeakStochasticRungeKutta : public StateVecStochasticIntegrator 

#include <svIntegratorWeakStochasticRungeKutta.h>

The svIntegratorWeakStochasticRungeKutta class implements a state integrator that provides weak solutions to problems with stochastic dynamics (SDEs).

The method is described in:

Tang, X., Xiao, A. Efficient weak second-order stochastic Runge–Kutta methods
for Itô stochastic differential equations. Bit Numer Math 57, 241–260 (2017).
https://doi.org/10.1007/s10543-016-0618-9

The method described in said paper is designed for autonomous systems (where f and g do not depend on time), so we need to modify it for non-autonomous systems. This is simple: we just need to treat the time as a state (with f=1 and g=0).

The resulting pseudocode for the (explicit) method is (note that the order of the sums and some elements has moved around to reflect the code):

--- (3.1*) stage definitions ---
for i = 1..s:
    // drift‐stage
    H0[i] = y_n
        + h * sum_{j=1..i-1}( a0[i][j] * f( t_n + c0[j]*h,   H0[j] ) )
        + sum_{k=1..m}( Ihat[k] * sum_{j=1..i-1}( b0[i][j] * g[k]( t_n + c1[j]*h,   Hk[j] ) ) )

    // diffusion‐stages, one per noise source k
    for k = 1..m:
        Hk[i] = y_n
            + h  * sum_{j=1..i-1}( a1[i][j] * f(    t_n + c0[j]*h, H0[j] ) )
            + xi * sum_{j=1..i-1}( b1[i][j] * g[k]( t_n + c1[j]*h, Hk[j] ) )
            + sum_{l=1..m, l!=k}( Ihat_kl[k][l] * sum_{j=1..i-1}( b2[i][j] * g[l]( t_n + c1[i]*h, Hl[i] ) ) );

where
  c0[j] = sum_{p=1..s}( a0[j][p] )
  c1[j] = sum_{p=1..s}( a1[j][p] )

--- (3.2) driving‐variable distributions ---
For each k:
    P( Ihat[k] = +sqrt(3*h) ) = 1/6
    P( Ihat[k] = -sqrt(3*h) ) = 1/6
    P( Ihat[k] =  0        ) = 2/3

xi = eta1 * sqrt(h)    // eta1 is {+1,−1} with equal prob.

--- (3.3) mixed‐integral approximations ---
for each pair (k,l):
    if   k < l:
        Ihat_kl[k][l] = 0.5 * ( Ihat[l] - eta2 * Ihat[l] )
    else if k > l:
        Ihat_kl[k][l] = 0.5 * ( Ihat[l] + eta2 * Ihat[l] )
    else // k == l:
        Ihat_kl[k][k] = 0.5 * ( Ihat[k]*Ihat[k] * xi - xi )

--- (3.1*) state evolution ---
y_n+1 = y_n
        + h * sum_{i=1..s}( alpha[i] * f( t_n + c0[i]*h, H0[i] ) )
        + sum_{k=1..m}( Ihat[k]       * sum_{i=1..s}( beta0[i] * g[k]( t_n + c1[i]*h, Hk[i] ) ) )
        + sum_{k=1..m}( Ihat_kl[k][k] * sum_{i=1..s}( beta1[i] * g[k]( t_n + c1[i]*h, Hk[i] ) ) )

Public Functions

svIntegratorWeakStochasticRungeKutta(DynamicObject *dynIn, const SRKCoefficients<numberStages> &coefficients): Creates an explicit RK integrator for the given DynamicObject using the passed coefficients.

inline void setRNGSeed(size_t seed)

Sets the seed for the Random Number Generator used by this integrator.

As a stochastic integrator, random numbers are drawn during each time step. By default, a randomly generated seed is used each time.

If the seed is set, the integrator will always draw the same numbers during time-stepping.

virtual void integrate(double currentTime, double timeStep) override: Performs the integration of the associated dynamic objects up to time currentTime+timeStep

Public Members

SRKRandomVariableGenerator rvGenerator: Random Number Generator for the integrator

Protected Functions

svIntegratorWeakStochasticRungeKutta(DynamicObject *dynIn, std::unique_ptr<SRKCoefficients<numberStages>> &&coefficients): Can be used by subclasses to support passing coefficients that are subclasses of SRKCoefficients

ExtendedStateVector computeDerivatives(double time, double timeStep)

Computes the derivatives of every state given a time and current states.

Internally, this sets the states on the dynamic objects and calls the equationsOfMotion methods.

std::vector<ExtendedStateVector> computeDiffusions(double time, double timeStep, const std::vector<StateIdToIndexMap> &stateIdToNoiseIndexMaps)

Computes the diffusions of every state given a time and current states.

Internally, this sets the states on the dynamic objects and calls the equationsOfMotionDiffusion methods.

ExtendedStateVector computeDiffusions(double time, double timeStep, const StateIdToIndexMap &stateIdToNoiseIndexMap)

Computes the diffusions for the states and noise index in stateIdToNoiseIndexMap given a time and current states.

Internally, this sets the states on the dynamic objects and calls the equationsOfMotionDiffusion methods.

inline ExtendedStateVector scaledSum(const std::array<double, numberStages> &factors, const std::array<ExtendedStateVector, numberStages> &vectors, size_t length)

Utility function, computes:

\[\text{result} = \sum_{i=0}^{\text{length}-1} \text{factors}[i] \cdot \text{vectors}[i] \]

Protected Attributes

const std::unique_ptr<SRKCoefficients<numberStages>> coefficients: Coefficients to be used in the method