Theory

Stochastic Hamiltonian systems

The systems this package integrates are stochastic differential equations in the Stratonovich sense, written with $\circ$:

\[dq(t) = v(t, q) \, dt + \sum_{r=1}^{m} B^{\cdot r}(t, q) \circ dW^r(t) .\]

The Stratonovich convention is what makes the geometry available. Its differential obeys the ordinary chain rule, so a symplectic form, a coordinate change or a conserved quantity behaves as it does in the deterministic case; in the Itô calculus each picks up a second-order correction and the structure these integrators preserve would not be there to preserve.

The case of real interest is a stochastic Hamiltonian system, with a Hamiltonian $H$ and one stochastic Hamiltonian $h_r$ per noise dimension:

\[\begin{aligned} dq &= \frac{\partial H}{\partial p} \, dt + \sum_{r} \frac{\partial h_r}{\partial p} \circ dW^r , \\ dp &= -\frac{\partial H}{\partial q} \, dt - \sum_{r} \frac{\partial h_r}{\partial q} \circ dW^r . \end{aligned}\]

Its flow is symplectic for every realisation of the noise, pathwise, exactly as in the deterministic case. A forced system adds a deterministic force $F$ and stochastic forces $f_r$:

\[\begin{aligned} dq &= \frac{\partial H}{\partial p} \, dt + \sum_{r} \frac{\partial h_r}{\partial p} \circ dW^r , \\ dp &= \left[ -\frac{\partial H}{\partial q} + F \right] dt + \sum_{r} \left[ -\frac{\partial h_r}{\partial q} + f_r \right] \circ dW^r , \end{aligned}\]

which is no longer symplectic — a damped oscillator must lose energy — but still arises from a variational principle, the stochastic Lagrange-d'Alembert principle. What a good integrator should preserve is that principle, not the symplectic form.

Strong and weak convergence

Two quite different things can be meant by "the numerical solution is accurate", and a stochastic scheme is designed for one or the other.

A scheme converges with mean-square (strong) order $p$ if

\[\sqrt{ \mathbb{E} \left| q(T) - q_N \right|^2 } \le C \, \Delta t^{\,p} ,\]

that is, it tracks the individual sample path generated by the same noise. This is what is needed to compare against an exact pathwise solution, to study a particular trajectory, or to refine a path.

A scheme converges with weak order $p$ if

\[\left| \mathbb{E}\, \varphi(q(T)) - \mathbb{E}\, \varphi(q_N) \right| \le C \, \Delta t^{\,p}\]

for smooth test functions $\varphi$ — it reproduces expectations, and says nothing about any one path. That is all that is needed for mean energies, ergodic limits or distributions, and it is a much weaker demand, which is why weak schemes can afford much cheaper randomness.

That cost difference is the whole point. A strong scheme must be fed genuine Gaussian increments; a weak scheme is derived only against the moments of the increments, so it can replace them with a discrete random variable that matches those moments and costs one uniform draw. See Noise Processes.

The order 1.0 ceiling

Every method in this package uses only the increments $\Delta W^r = \int_{t_n}^{t_{n+1}} dW^r$, and possibly the iterated integrals $\Delta Z^r = \int_{t_n}^{t_{n+1}} \int_{t_n}^{t} dW^r(\xi) \, dt$. A scheme built from $\Delta W$ alone cannot exceed mean-square order 1.0, however many stages it has — the missing information is in the higher multiple Stratonovich integrals, not in the tableau.

Carrying the $\Delta Z$ terms is therefore necessary to pass that ceiling, but it is not sufficient. Of the methods here only BurrageCL, BurrageE1 and BurrageG5 carry them, and of those three only BurrageCL and BurrageG5 reach order 1.5 — BurrageE1 uses the same terms and remains of order 1.0.

Order 1.0 further requires the noise to be commutative. For a one-dimensional Wiener process commutativity is automatic; for multi-dimensional non-commutative noise the same schemes drop to order 0.5.

Variational integrators

A deterministic variational integrator is derived by discretising the action and applying a discrete variational principle, rather than by discretising the equations of motion. The resulting scheme is symplectic and momentum-preserving by construction.

The same construction works in the stochastic setting. Kraus and Tyranowski derive a discrete stochastic Lagrange-d'Alembert principle and show that a stochastic partitioned Runge-Kutta method

\[\begin{aligned} Q_i &= q_n + \Delta t \sum_j a_{ij} \frac{\partial H}{\partial p}(Q_j,P_j) + \sum_r \Delta W^r \sum_j b_{ij} \frac{\partial h_r}{\partial p}(Q_j,P_j) , \\ P_i &= p_n - \Delta t \sum_j \bar a_{ij} \frac{\partial H}{\partial q}(Q_j,P_j) - \sum_r \Delta W^r \sum_j \bar b_{ij} \frac{\partial h_r}{\partial q}(Q_j,P_j) \\ &\qquad + \Delta t \sum_j \hat a_{ij} F(Q_j,P_j) + \sum_r \Delta W^r \sum_j \hat b_{ij} f_r(Q_j,P_j) , \end{aligned}\]

with final update

\[\begin{aligned} q_{n+1} &= q_n + \Delta t \sum_i \alpha_i \frac{\partial H}{\partial p}(Q_i,P_i) + \sum_r \Delta W^r \sum_i \beta_i \frac{\partial h_r}{\partial p}(Q_i,P_i) , \\ p_{n+1} &= p_n - \Delta t \sum_i \alpha_i \frac{\partial H}{\partial q}(Q_i,P_i) - \sum_r \Delta W^r \sum_i \beta_i \frac{\partial h_r}{\partial q}(Q_i,P_i) \\ &\qquad + \Delta t \sum_i \hat\alpha_i F(Q_i,P_i) + \sum_r \Delta W^r \sum_i \hat\beta_i f_r(Q_i,P_i) , \end{aligned}\]

is such an integrator precisely when its coefficients satisfy the Lagrange-d'Alembert conditions

\[\begin{aligned} \alpha_i \bar a_{ij} + \alpha_j a_{ji} &= \alpha_i \alpha_j , & \beta_i \bar b_{ij} + \beta_j b_{ji} &= \beta_i \beta_j , \\ \beta_i \bar a_{ij} + \alpha_j b_{ji} &= \beta_i \alpha_j , & \alpha_i \bar b_{ij} + \beta_j a_{ji} &= \alpha_i \beta_j , \\ \alpha_i \hat a_{ij} + \hat\alpha_j a_{ji} &= \alpha_i \hat\alpha_j , & \alpha_i \hat b_{ij} + \hat\beta_j a_{ji} &= \alpha_i \hat\beta_j , \\ \beta_i \hat a_{ij} + \hat\alpha_j b_{ji} &= \beta_i \hat\alpha_j , & \beta_i \hat b_{ij} + \hat\beta_j b_{ji} &= \beta_i \hat\beta_j . \end{aligned}\]

The first four are the familiar symplecticity conditions of a partitioned Runge-Kutta method, extended to couple the drift and diffusion tableaus. The last four are new, and they are the reason the $p$ equation is split: they couple the forcing tableaus $\hat a$, $\hat b$ to $a$, $b$. A method that applied one set of coefficients to the whole right-hand side of the $p$ equation could not satisfy them in general, which is why forced systems need an SPSDE and a SISPRK method rather than a PSDE and a SIPRK.

Mean-square order 1.0 is a separate requirement, and a considerably simpler one:

\[\sum_i \alpha_i = \sum_i \hat\alpha_i = \sum_i \beta_i = \sum_i \hat\beta_i = 1 ,\]

\[\sum_{i,j} \beta_i b_{ij} = \sum_{i,j} \beta_i \bar b_{ij} = \sum_{i,j} \beta_i \hat b_{ij} = \sum_{i,j} \hat\beta_i b_{ij} = \sum_{i,j} \hat\beta_i \bar b_{ij} = \sum_{i,j} \hat\beta_i \hat b_{ij} = \tfrac{1}{2} .\]

Note that the drift matrices $a$, $\bar a$, $\hat a$ do not appear at all — only their weights, and the six $\beta$/$\hat\beta$ against $b$/$\bar b$/$\hat b$ products.

The two sets are independent, and a method can satisfy either without the other:

All of this is checked directly on the tableau coefficients by scripts/tableau_conditions.jl, including the two deliberate failures — a run that reported everything satisfied would be reporting a bug in itself.

Verification

Two of the scripts in scripts/ establish the claims made here, and are the place to look before trusting a number:

• tableau_conditions.jl evaluates both sets of conditions above on every tableau in the package and asserts the expected outcome, failures included. It is exact and instant, so it isolates a transcription error in a tableau from an error in the code that consumes it. Run it first.
• convergence_order.jl measures the mean-square order on the Kubo oscillator, whose exact solution is known in closed form, driving each refinement along one common sample path with a GridProcess.