Symplectic Runge-Kutta Methods for Hamiltonian Systems Driven by Gaussian Rough Paths

TL;DR

This paper develops symplectic Runge-Kutta methods for Hamiltonian systems driven by Gaussian rough paths, demonstrating their structure-preserving properties and convergence, with numerical validation.

Contribution

It introduces symplectic Runge-Kutta methods for rough path driven Hamiltonian systems, establishing their solvability, convergence rates, and structure preservation.

Findings

Phase flow preserves symplectic structure almost surely.

Symplectic Runge-Kutta methods inherit symplecticity.

Numerical experiments confirm theoretical results.

Abstract

We consider Hamiltonian systems driven by multi-dimensional Gaussian processes in rough path sense, which include fractional Brownian motions with Hurst parameter $H\in(1/4,1/2]$. We indicate that the phase flow preserves the symplectic structure almost surely and this property could be inherited by symplectic Runge--Kutta methods, which are implicit methods in general. If the vector fields belong to $Lip^{\gamma}$, we obtain the solvability of Runge--Kutta methods and the pathwise convergence rates. For linear and skew symmetric vector fields, we focus on the midpoint scheme to give corresponding results. Numerical experiments verify our theoretical analysis.

Equations

$$d{Y}_t$$

$$\Vert \text{X}{\Vert }_{p\text{-}var;[s,t]}:=\underset{(t_k)\in \mathcal{D}([s,t])}{\sup }{\left(\sum _{t_k}\mathrm{d}({\text{X}}_{t_k},{\text{X}}_{t_{k+1}}{)}^p\right)}^{1/p},$$

$$\text{X}\in {\mathcal{C}}^{p\text{-}var}([0,T],G^{[p]}({\mathbb{R}}^d)):=\{\text{X}:\Vert \text{X}{\Vert }_{p\text{-}var;[0,T]}<\infty \},$$

$$\Vert \text{X}{\Vert }_{1/p\text{-}\mathrm{H}\ddot{\mathrm{o}}\mathrm{l};[0,T]}:=\underset{0\le s<t\le T}{\sup }\frac{\mathrm{d}({\text{X}}_s,{\text{X}}_t)}{|t-s{|}^{1/p}}<\infty .$$

$$S_{[p]}(x{)}_t:=(1,{\int }_{0\le u_1\le t}d{x}_{u_1},...,{\int }_{0\le u_1<\cdot \cdot \cdot <u_{[p]}\le t}d{x}_{u_1}\otimes \cdot \cdot \cdot \otimes d{x}_{u_{[p]}})\in G^{[p]}({\mathbb{R}}^d).$$

$$d{Y}_t=V(Y_t)d{\text{X}}_t,Y_0=z\in {\mathbb{R}}^m.$$

$$\begin{array}{r}\underset{n\in \mathbb{N}}{\sup }\Vert S_{[p]}(x^n){\Vert }_{p\text{-}var;[0,T]}<\infty ,\\ \underset{n\to \infty }{\lim }\underset{0\le s<t\le T}{\sup }\mathrm{d}(S_{[p]}(x^n{)}_{s,t},{\text{X}}_{s,t})=0.\end{array}$$

$$\begin{array}{rl}& |D^{k}V(y)|\le C,\forall k=0,...,\lfloor \gamma \rfloor ,\forall y\in {\mathbb{R}}^m,\\ & |D^{\lfloor \gamma \rfloor }V(y_1)-D^{\lfloor \gamma \rfloor }V(y_1)|\le C|y_1-y_2{|}^{\gamma -\lfloor \gamma \rfloor },\forall y_1,y_2\in {\mathbb{R}}^m,\end{array}$$

$$d\frac{\partial Y_t}{\partial z}=\sum _{i=1}^{d}D{V}_i(Y_t)\frac{\partial Y_t}{\partial z}d{\text{X}}_t^i,\frac{\partial Y_0}{\partial z}=I_m.$$

$$\underset{t\in [0,T]}{\sup }\left|\frac{\partial Y_t}{\partial z}\right|\le |I_m|+C{\nu }_1\Vert \text{X}{\Vert }_{p\text{-}var;[0,T]}\exp \{C{\nu }_1^p\Vert \text{X}{\Vert }_{p\text{-}var;[0,T]}^p\},$$

$$\underset{t\in [0,T]}{\sup }|Y_t|\le \theta +C(1+\theta ){\nu }_2\Vert \text{X}{\Vert }_{p\text{-}var;[0,T]}\exp \{c{\nu }_2^p\Vert \text{X}{\Vert }_{p\text{-}var;[0,T]}^p\}=:R(\theta ),$$

$$\delta {\int }_0^{T}P\dot{Q}-H(P,Q)dt=0$$

$$H_0(P,Q)+\sum _{i=1}^d{H}_i(P,Q){\dot{\chi }}^i$$

$$\delta S$$

$$={\int }_0^{T}P\delta \dot{Q}+\dot{Q}\delta P-\frac{\delta H_0}{\delta P}\delta P-\frac{\delta H_0}{\delta Q}\delta Q-\sum _{i=1}^d\left(\frac{\delta H_i}{\delta P}{\dot{\chi }}^i\delta P+\frac{\delta H_i}{\delta Q}{\dot{\chi }}^i\delta Q\right)dt\equiv 0.$$

$$\dot{P}=-\frac{\partial H_0(P,Q)}{\partial Q}-\sum _{i=1}^d\frac{\partial H_i(P,Q)}{\partial Q}{\dot{\chi }}^i,\dot{Q}=\frac{\partial H_0(P,Q)}{\partial P}+\sum _{i=1}^d\frac{\partial H_i(P,Q)}{\partial P}{\dot{\chi }}^i,$$

$$\begin{array}{rl}dP=& -\frac{\partial H_0(P,Q)}{\partial Q}dt-\sum _{i=1}^d\frac{\partial H_i(P,Q)}{\partial Q}d{X}^i,P_0=p,\\ dQ=& \frac{\partial H_0(P,Q)}{\partial P}dt+\sum _{i=1}^d\frac{\partial H_i(P,Q)}{\partial P}d{X}^i,Q_0=q.\end{array}$$

$$d{Y}_t=\sum _{i=0}^d{V}_i(Y_t)d{X}_t^i=V(Y_t)d{X}_t,Y_0=z.$$

$$\underset{(t_k),(t_l)\in \mathcal{D}([s,t])}{\sup }{\left(\sum _{t_i,t_l}|\mathbb{E}{X}_{t_k,t_{k+1}}{X}_{t_l,t_{l+1}}{|}^{\rho }\right)}^{1/\rho }\le K|t-s{|}^{1/\rho },$$

$$dP\wedge dQ=dp\wedge dq,a.s.$$

$$dP\wedge dQ=$$

$$=$$

$$+\sum _{j=1}^m\sum _{k=1}^m\sum _{l=1}^{k-1}\left(P_p^{jk}{Q}_p^{jl}-P_p^{jl}{Q}_p^{jk}\right)d{p}^k\wedge d{q}^l$$

$$+\sum _{j=1}^m\sum _{k=1}^m\sum _{l=1}^{k-1}\left(P_q^{jk}{Q}_q^{jl}-P_q^{jl}{Q}_q^{jk}\right)d{p}^k\wedge d{q}^l.$$

$$\begin{array}{rl}& d{P}_p^{jk}=\sum _{i=0}^d\sum _{r=1}^m\left(-\frac{{\partial }^2{H}_i}{\partial Q^j\partial P^r}{P}_p^{rk}-\frac{{\partial }^2{H}_i}{\partial Q^j\partial Q^r}{Q}_p^{rk}\right)d{X}^i,P_p^{jk}(t_0)={\delta }_{jk},\\ & d{Q}_p^{jk}=\sum _{i=0}^d\sum _{r=1}^m\left(\frac{{\partial }^2{H}_i}{\partial P^j\partial P^r}{P}_p^{rk}+\frac{{\partial }^2{H}_i}{\partial P^j\partial Q^r}{Q}_p^{rk}\right)d{X}^i,Q_p^{jk}(t_0)=0,\\ & d{P}_q^{jk}=\sum _{i=0}^d\sum _{r=1}^m\left(-\frac{{\partial }^2{H}_i}{\partial Q^j\partial P^r}{P}_q^{rk}-\frac{{\partial }^2{H}_i}{\partial Q^j\partial Q^r}{Q}_q^{rk}\right)d{X}^i,P_q^{jk}(t_0)=0,\\ & d{Q}_q^{jk}=\sum _{i=0}^d\sum _{r=1}^m\left(\frac{{\partial }^2{H}_i}{\partial P^j\partial P^r}{P}_q^{rk}+\frac{{\partial }^2{H}_i}{\partial P^j\partial Q^r}{Q}_q^{rk}\right)d{X}^i,Q_q^{jk}(t_0)={\delta }_{jk},\end{array}$$

$$\begin{array}{rl}& d\left(\sum _{j=1}^m\left(P_p^{jk}{Q}_q^{jl}-P_q^{jl}{Q}_p^{jk}\right)\right)=0,\forall k,l,\\ & d\left(\sum _{j=1}^m\left(P_p^{jk}{Q}_p^{jl}-P_p^{jl}{Q}_p^{jk}\right)\right)=0,\forall k\ne l,\\ & d\left(\sum _{j=1}^m\left(P_q^{jk}{Q}_q^{jl}-P_q^{jl}{Q}_q^{jk}\right)\right)=0,\forall k\ne l,\end{array}$$

$$\begin{array}{rl}{Y}_k^h(\alpha )& =Y_k^h+\sum _{\beta =1}^s{a}_{\alpha \beta }V(Y_k^h(\beta ))X_{t_k,t_{k+1}},\\ Y_{k+1}^h& =Y_k^h+\sum _{\alpha =1}^s{b}_{\alpha }V(Y_k^h(\alpha ))X_{t_k,t_{k+1}}\end{array}$$

$$a_{\alpha \beta }{b}_{\alpha }+a_{\beta \alpha }{b}_{\beta }=b_{\alpha }{b}_{\beta },\forall \alpha ,\beta =1,...,s.$$

$$\begin{array}{rl}d{P}_{k+1}^{h,j}& =d{P}_k^{h,j}-\sum _{\alpha =1}^s\sum _{i=0}^d\sum _{r=1}^m{b}_{\alpha }\left(\frac{{\partial }^2{H}_i}{\partial Q^j\partial P^r}d{P}_k^{h,r}(\alpha )+\frac{{\partial }^2{H}_i}{\partial Q^j\partial Q^r}d{Q}_k^{h,r}(\alpha )\right)X_{t_k,t_{k+1}}^i,\\ d{Q}_{k+1}^{h,j}& =d{Q}_k^{h,j}+\sum _{\alpha =1}^s\sum _{i=0}^d\sum _{r=1}^m{b}_{\alpha }\left(\frac{{\partial }^2{H}_i}{\partial P^j\partial P^r}d{P}_k^{h,r}(\alpha )+\frac{{\partial }^2{H}_i}{\partial P^j\partial Q^r}d{Q}_k^{h,r}(\alpha )\right)X_{t_k,t_{k+1}}^i.\end{array}$$

$$\begin{array}{rl}& d{P}_{k+1}^{h,j}\wedge d{Q}_{k+1}^{h,j}=d{P}_k^{h,j}\wedge d{Q}_k^{h,j}\\ & +\sum _{i=0}^d\sum _{r=1}^m\sum _{\alpha =1}^s{b}_{\alpha }d{P}_k^{h,j}\wedge \left(\frac{{\partial }^2{H}_i}{\partial P^j\partial P^r}d{P}_k^{h,r}(\alpha )+\frac{{\partial }^2{H}_i}{\partial P^j\partial Q^r}d{Q}_k^{h,r}(\alpha )\right)X_{t_k,t_{k+1}}^i\\ & -\sum _{i=0}^d\sum _{r=1}^m\sum _{\alpha =1}^s{b}_{\alpha }\left(\frac{{\partial }^2{H}_i}{\partial Q^j\partial P^r}d{P}_k^{h,r}(\alpha )+\frac{{\partial }^2{H}_i}{\partial Q^j\partial Q^r}d{Q}_k^{h,r}(\alpha )\right)\wedge d{Q}_k^{h,j}{X}_{t_k,t_{k+1}}^i\\ & -\left(\sum _{i_1=0}^d\sum _{r_1=1}^m\sum _{\alpha =1}^s{b}_{\alpha }\left(\frac{{\partial }^2{H}_{i_1}}{\partial Q^j\partial P^{r_1}}d{P}_k^{h,r_1}(\alpha )+\frac{{\partial }^2{H}_{i_1}}{\partial Q^j\partial Q^{r_1}}d{Q}_k^{h,r_1}(\alpha )\right)X_{t_k,t_{k+1}}^{i_1}\right)\\ & \wedge \left(\sum _{i_2=0}^d\sum _{r_2=1}^m\sum _{\alpha =1}^s{b}_{\alpha }\left(\frac{{\partial }^2{H}_{i_2}}{\partial P^j\partial P^{r_2}}d{P}_k^{h,r_2}(\alpha )+\frac{{\partial }^2{H}_{i_2}}{\partial P^j\partial Q^{r_2}}d{Q}_k^{h,r_2}(\alpha )\right)X_{t_k,t_{k+1}}^{i_2}\right).\end{array}$$

$$\begin{array}{rl}d{P}_k^{h,j}& =d{P}_k^{h,j}(\alpha )+\sum _{\beta =1}^s\sum _{i=0}^d\sum _{r=1}^m{a}_{\alpha \beta }\left(\frac{{\partial }^2{H}_i}{\partial Q^j\partial P^r}d{P}_k^{h,r}(\beta )+\frac{{\partial }^2{H}_i}{\partial Q^j\partial Q^r}d{Q}_k^{h,r}(\beta )\right)X_{t_k,t_{k+1}}^i,\\ d{Q}_k^{h,j}& =d{Q}_k^{h,j}(\alpha )-\sum _{\beta =1}^s\sum _{i=0}^d\sum _{r=1}^m{a}_{\alpha \beta }\left(\frac{{\partial }^2{H}_i}{\partial P^j\partial P^r}d{P}_k^{h,r}(\beta )+\frac{{\partial }^2{H}_i}{\partial P^j\partial Q^r}d{Q}_k^{h,r}(\beta )\right)X_{t_k,t_{k+1}}^i,\end{array}$$

$$a_{\alpha \beta }{b}_{\alpha }+a_{\beta \alpha }{b}_{\beta }=b_{\alpha }{b}_{\beta },\forall \alpha ,\beta =1,...,s.$$

Full text

Authors are supported by National Natural Science Foundation of China (NO. 91530118, NO. 91130003, NO. 11021101 and NO. 11290142).

1. Introduction

We consider a stochastic differential equation (SDE)

[TABLE]

driven by multi-dimensional Gaussian process in rough path sense. Under proper assumptions, $X$ can be lifted to a rough path almost surely and then (1.1) can be interpreted as a rough differential equation (RDE). For example, the lift of the fractional Brownian motion (fBm) with Hurst parameter $H>1/4$ is constructed by piecewise linear approximations and their iterated integrals. We refer to [7, 5] and references therein for more details.

The well-posedness of the RDE is given originally in [13] for the case that the vector fields belong to $Lip^{\gamma}$ (in the sense of E. Stein), i.e., they are bounded and smooth enough with bounded derivatives. When the vector fields are linear, the well-posedness result still holds (see e.g., [7]). In particular, the solution is equivalent to that in Stratonovich sense almost surely when the noise is a semi-martingale. The robustness of the solution allows numerous researches to be developed, such as the density and ergodicity of SDEs driven by non-Markovian noises (see [2, 10] and references therein) and the theory of stochastic partial differential equations (see e.g., [5]).

The author in [3] develops an approximation approach to deal with the equation driven by non-differentiable paths, which are further investigated and called step-$N$ Euler schemes in [7]. When the equation is driven by fBms with $H>1/2$, modified Euler schemes and Taylor schemes are analyzed in [11, 12]. For rougher noises, the step-$N$ Euler schemes are not implementable from the numerical point of view, because they depend on iterated integrals of the noise, which are rather difficult to be simulated. The current implementable methods are the simplified step-$N$ Euler schemes, which are first proposed in [4] for fBm with $H\in(1/3,1/2]$. Authors in [6] give the estimate of the pathwise convergence rate of Wong–Zakai approximations for Gaussian processes under proper assumptions and obtain convergence rates of the simplified step-$N$ Euler schemes. The results are applicable for fBm with $H\in(1/4,1/2]$. Further, the $L^{r}$-convergence rate is obtained in [1] to reduce the complexity of multilevel Monte Carlo.

From the perspective of modeling, Gaussian noises with nontrivial correlations are more general than standard Brownian motions. For a dynamic system with additional nonconservative forces in [8], which could be generally assumed as Gaussian noises, its pathwise movement rule is described by Hamiltonian equations in rough sense. In this paper, we are concerned about a Hamiltonian system driven by multi-dimensional Gaussian process under certain assumption. The fBm with $H\in(1/4,1/2]$ satisfies our assumption with $\rho=\frac{1}{2H}$, where $\rho$ is a parameter related to the regularity of the noise. We prove that the rough Hamiltonian system, as a generalization of the classical case in [9] and the stochastic case driven by standard Brownian motions in [14, 15], also has the characteristic property of the Hamiltonian system that its phase flow preserves the symplectic structure almost surely. The symplecticity is naturally considered to be inherited by numerical methods. However, the simplified step-$N$ Euler schemes, which are explicit, are failed to possess discrete symplectic structure. To our best knowledge, there has been few work about implicit methods for RDEs driven by noises rougher than standard Brownian motions.

We make use of symplectic Runge–Kutta methods to preserve the symplectic structure of the original rough Hamiltonian system. Since these methods are implicit methods in general, the solvability of the numerical methods should be taken into consideration, which is one of differences between explicit methods and implicit methods. For SDEs driven by Brownian motions, authors in [14] truncate each increment of Brownian motions to avoid the possibility that the increments could become unbounded, and give the solvability of implicit methods and convergence rates in mean square sense. However, this truncation skill is not applicable here, as the convergence is considered in pathwise sense. Based on the Brouwer’s theorem, we prove that if the vector fields belong to $Lip^{\gamma}$, then Runge–Kutta methods are solvable for any time step. If the vector fields are linear and skew symmetric, we get the solvability of the $1$-stage symplectic Runge–Kutta method—the midpoint scheme.

Furthermore, we analyze pathwise convergence rates for symplectic Runge–Kutta methods. For the $Lip^{\gamma}$ case, we obtain that the convergence rate of the midpoint scheme is $(\frac{3}{2\rho}-1-\varepsilon)$, for arbitrary small $\varepsilon$ and $\rho\in[1,3/2)$. The convergence rate of another two symplectic Runge-Kutta methods, whose stages are higher, is $(\frac{1}{\rho}-\frac{1}{2}-\varepsilon)$ for $\rho\in[1,2)$. This is limited by the convergence rate of Wong–Zakai approximations, similar to the simplified step-$3$ Euler scheme in [6]. If the vector fields are linear and skew symmetric, we prove the uniform boundedness of numerical solutions for the midpoint scheme and have that its convergence rate is also $(\frac{3}{2\rho}-1-\varepsilon)$ when $\rho\in[1,3/2)$.

This paper is organized as follows. In Section 2, we introduce notations and definitions in rough path theory as well as the result of the well-posedness of RDEs, which contains both the $Lip^{\gamma}$ case and the linear case. In Section 3, we derive Hamiltonian equations with additional nonconservative forces through the variational principle to obtain the rough Hamiltonian system and prove that the phase flow preserves the symplectic structure almost surely. In Section 4, we propose symplectic Runge–Kutta methods for rough Hamiltonian systems and give the solvability of them. In Section 5, we analyze pathwise convergence rates for symplectic Runge–Kutta methods. Numerical experiments in Section 6 verify our theoretical analysis.

2. Preliminaries in rough path theory

We first recall some notations in rough path theory. For $p\in[1,\infty)$, we are interested in a continuous map $\textbf{X}:[0,T]\rightarrow G^{[p]}(\mathbb{R}^{d})$, where $[p]$ is an integer satisfying $p-1<[p]\leq p$ and $G^{[p]}(\mathbb{R}^{d})$ is the free step-$[p]$ nilpotent Lie group of $\mathbb{R}^{d}$, equipped with the Carnot-Carathéodory metric $\rm d$.

Since the group $G^{[p]}(\mathbb{R}^{d})$ is embedded in the truncated step-$[p]$ tensor algebra, i.e., $G^{[p]}(\mathbb{R}^{d})\subset\bigoplus_{n=0}^{[p]}(\mathbb{R}^{d})^{\otimes n}$ with $(\mathbb{R}^{d})^{\otimes 0}=\mathbb{R}$, the increment of X is defined by $\textbf{X}_{s,t}:=\textbf{X}_{s}^{-1}\otimes\textbf{X}_{t}$. Denoting by $\pi_{1}(G^{[p]}(\mathbb{R}^{d}))$ the projection of $G^{[p]}(\mathbb{R}^{d})$ into $\mathbb{R}^{d}$, we have $\pi_{1}(\textbf{X}_{s,t})=\pi_{1}(\textbf{X}_{t})-\pi_{1}(\textbf{X}_{s})$.

Define the $p$-variation of X by

[TABLE]

where $\mathcal{D}([s,t])$ is the set of dissections of $[s,t]$. We call X a weak geometric $p$-rough path if

[TABLE]

where the set $\mathcal{C}^{p\text{-}var}([0,T],G^{[p]}(\mathbb{R}^{d}))$ contains drivers of RDEs. In addition, we say X is of H$\rm\ddot{o}$lder-type if

[TABLE]

This implies $|\pi_{1}(\textbf{X}_{s,t})|\leq\|\textbf{X}\|_{1/p\text{-}{\rm H\ddot{o}l};[0,T]}|t-s|^{1/p}$ for any $0\leq s<t\leq T$. Similar to the $1/p$-H$\rm\ddot{o}$lder continuity in classical case, a larger $p$ implies lower regularity of X.

The solution of an RDE is constructed by means of a sequence of bounded variation functions on $\mathbb{R}^{d}$, which are appropriate approximations of X. Let $x:[0,T]\rightarrow\mathbb{R}^{d}$ be a bounded variation function, we consider its canonical lift to a weak geometric $p$-rough path defined by $S_{[p]}(x)$ with

[TABLE]

One can observe that $S_{[p]}(x)$ contains information about iterated integrals of $x$ up to $[p]$-level, which is corresponding to the regularity of X. If $S_{[p]}(x)$ is close to X in $p$-variation sense, then $x$ is a proper approximation of X. In other words, the rougher X is the more information is required.

We now introduce the definition of the solution of an RDE

[TABLE]

Definition 2.1*.*

([7, Definition 10.17]) Let $p\in[1,\infty)$ and $\textbf{X}\in\mathcal{C}^{p\text{-}var}([0,T],G^{[p]}(\mathbb{R}^{d}))$. Suppose there exists a sequence of bounded variation functions $\{x^{n}\}$ on $\mathbb{R}^{d}$, such that

[TABLE]

Suppose in addition $\{y^{n}\}$ are solutions of equations $dy^{n}_{t}=V(y^{n}_{t})dx^{n}_{t}$, $y^{n}_{0}=z$, in Riemann-Stieltjes integral sense. If $y^{n}_{t}$ converges to $Y_{t}$ in $L^{\infty}([0,T])$-norm, then we call $Y_{t}$ a solution of (2.1).

To ensure the well-posedness of an RDE, proper assumptions are given for $V$, which will be described by the notation $Lip^{\gamma}$. Throughout this paper, $|\cdot|$ is the Euclidean norm, and we will use $C$ as generic constants, which may be different from line to line.

Definition 2.2*.*

Let $\gamma>0$ and $\lfloor\gamma\rfloor$ be the largest integer strictly smaller than $\gamma$, i.e., $\gamma-1\leq\lfloor\gamma\rfloor<\gamma$. We say $V\in Lip^{\gamma}$, if $V$ is $\lfloor\gamma\rfloor$-Fr$\rm\acute{e}$chet differentiable and the $k$th-derivative of $V$, $D^{k}V$, satisfies that

[TABLE]

for some constant $C$. The smallest constant $C$ satisfying these two inequalities is denoted by $|V|_{Lip^{\gamma}}$.

In the sequel, we give the theorem of the well-posedness of an RDE.

Theorem 2.1*.*

Let $\textbf{X}=(\textbf{X}^{i})_{1\leq i\leq d}\in\mathcal{C}^{p\text{-}var}([0,T],G^{[p]}(\mathbb{R}^{d}))$. If $V=(V_{i})_{1\leq i\leq d}$ is a collection of vector fields in $Lip^{\gamma}$ with $\gamma>p$, or a collection of linear vector fields of the form $V_{i}(Y)=A_{i}Y$, then (2.1) has a unique solution on $[0,T]$. Moreover, the Jacobian of the flow, $\frac{\partial Y_{t}}{\partial z}$, exists and satisfies the linear RDE

[TABLE]

Proof.

For the case of $V\in Lip^{\gamma}$, the existence and uniqueness of the solution follow from Theorem 10.14 and Theorem 10.26 in [7]. From Theorem 10.26 in [7] and Propsition 1 in [2], we know that $\frac{\partial Y_{t}}{\partial z}$ exists and satisfies (2.2) with a bound

[TABLE]

where $\nu_{1}\geq|V|_{Lip^{\gamma}}$ and $C$ depends only on $p$.

In the linear case, for any fixed initial value $z_{0}$, the existence and uniqueness of the solution follow from Theorem 10.53 in [7]. It implies that $Y$ will not blow up on $[0,T]$ and

[TABLE]

where $\nu_{2}\geq\max\{|A_{i}|:i=0,...,d\}$, $\theta\geq|z_{0}|$ and $C$ depends only on $p$. This guarantees that we can localize the problem as Theorem 10.21 in [7]. Considering the set $\{z:|z|\leq 2\theta\}$, which is a neighborhood of $z_{0}$, and replacing $V$ by a compactly supported $\tilde{V}\in Lip^{\gamma}$ which coincides with $V$ in the ball $B_{2R}:=\{y:|y|\leq 2R(\theta)\}$, we turn the problem into the first case and obtain that $\frac{\partial Y_{t}}{\partial z}\mid_{z=z_{0}}$ exists and satisfies (2.2). ∎

3. Symplectic structure of rough Hamiltonian system

In this section, we derive autonomous Hamiltonian equations with additional nonconservative forces via a modified Hamilton’s principle. Under Assumption 3.1, we show that these equations form a rough Hamiltonian system when the nonconservative forces are characterized by rough signals. Further, this rough Hamiltonian system possesses the symplectic structure similar to the deterministic case.

The classical Hamilton’s principle indicates that the motion $Q(t)$ extremizes the action functional $S(Q):=\int_{0}^{T}L(Q,\dot{Q})dt$ with $\delta Q(0)=\delta Q(T)=0$ under variation, where $L$ is the Lagrangian of a deterministic Hamiltonian system. The author in [8] gives the modified Hamilton’s principle

[TABLE]

based on the Legendre transform $H(P,Q)=P\dot{Q}-L(Q,\dot{Q})$. However, for a Hamiltonian system influenced by additional nonconservative forces, its Hamiltonian energy turns to be

[TABLE]

instead of $H(P,Q)$. The second term is the total work done by the additional forces and $\dot{\chi}^{i}=\frac{dX^{i}_{t}}{dt}$ represents a formal time derivative of $X^{i}$. In this case, the modified Hamilton’s principle reads

[TABLE]

Integration by parts eventually yields the Hamiltonian equations of the motion with additional forces:

[TABLE]

which can be rewritten as

[TABLE]

Suppose $P=(P^{1},...,P^{m})^{\top},Q=(Q^{1},...,Q^{m})^{\top}$, $p=(p^{1},...,p^{m})^{\top},q=(q^{1},...,q^{m})^{\top}$ $\in\mathbb{R}^{m}$. In order to get a compact form of (3.1), we denote $Y=(P^{\top},Q^{\top})^{\top}$, $z=(p^{\top},q^{\top})^{\top}$ $\in\mathbb{R}^{2m}$, $X_{t}^{0}=t$, $X_{t}=(X^{0}_{t},X^{1}_{t},X^{2}_{t},...,X^{d}_{t})^{\top}\in\mathbb{R}^{d+1}$ and $V_{i}=(V_{i}^{1},...,V_{i}^{2m})=(-\frac{\partial H_{i}}{\partial Q^{1}},...,-\frac{\partial H_{i}}{\partial Q^{m}},\frac{\partial H_{i}}{\partial P^{1}},...,\frac{\partial H_{i}}{\partial P^{m}})$, $i=0,...,d$. Then (3.1) is equivalent to

[TABLE]

In the sequel, we consider $X^{i}_{t}=X^{i}_{t}(\omega)$, $i=1,...,d$, as Gaussian noises under the following assumption.

Assumption 3.1*.*

Let $X^{i}_{t}$, $i=1,...,d$, be independent centered Gaussian processes with continuous sample path on $[0,T]$. There exist some $\rho\in[1,2)$ and $K\in(0,+\infty)$ such that for any $0\leq s<t\leq T$, the covariance of $X$ satisfies

[TABLE]

where $X_{t_{k},t_{k+1}}=X_{t_{k+1}}-X_{t_{k}}$.

For any $p>2\rho$, by piecewise linear approximations, $X_{t}$ can be naturally lifted to a H$\rm\ddot{o}$lder-type weak geometric $p$-rough path $\textbf{X}_{t}$ almost surely, which takes values in $G^{[p]}(\mathbb{R}^{d+1})$ and $\pi_{1}(\textbf{X}_{s,t})=X_{s,t}$ (see e.g., [7, Theorem 15.33]). Therefore, (3.1) or (3.2) can be transformed into an RDE almost surely and interpreted in the framework of rough path theory. In this sense, we call (3.1) or (3.2) a rough Hamiltonian system.

The fBm with Hurst parameter $H\in(1/4,1/2]$ satisfies Assumption 3.1 with $\rho=\frac{1}{2H}$, which is frequently used in practical applications. When $X_{t}^{i}$, $i=1,...,d$, are Brownian motions, i.e., $H=1/2$, $Y_{t}$ equals to the solution of a stochastic Hamiltonian system driven by standard Brownian motions (see [15, 14]) in Stratonovich sense almost surely.

As a function of time $t$ and initial value $z$, $Y$ is a phase flow for almost every $\omega$. In the deterministic case and the stochastic case, we know that the phase flow preserves the symplectic structure, i.e., the differential $2$-form $dP\wedge dQ$ is invariant. Note that the differential here is made with respect to the initial value, which is different from the formal time derivative in (3.1). The geometric interpretation is that the sum of the oriented areas of a two-dimensional surface, obtained by projecting the phase flow onto the coordinate planes ($p^{1},q^{1}$),…,($p^{m},q^{m}$), is an integral invariant. The rough Hamiltonian systems, a generalization for the deterministic case and the stochastic case but allowing rougher noises, share this characteristic property of Hamiltonian systems, which is proved in the next Theorem.

Theorem 3.1*.*

The phase flow of the rough Hamiltonian system (3.1) preserves the symplectic structure:

[TABLE]

Proof.

From Theorem 2.1, we know that $P$ and $Q$ are differentiable with respect to $p$ and $q$. Denote $P^{jk}_{p}=\frac{\partial P^{j}}{\partial p^{k}}$, $Q^{jk}_{p}=\frac{\partial Q^{j}}{\partial p^{k}}$, $P^{jk}_{q}=\frac{\partial P^{j}}{\partial q^{k}}$, $Q^{jk}_{q}=\frac{\partial Q^{j}}{\partial q^{k}}$. Since $dP^{j}=\sum\limits^{m}_{k=1}P^{jk}_{p}dp^{k}+\sum\limits^{m}_{l=1}P^{jl}_{q}dq^{l}$ and $dQ^{j}=\sum\limits^{m}_{k=1}Q^{jk}_{p}dp^{k}+\sum\limits^{m}_{l=1}Q^{jl}_{q}dq^{l}$, we have

[TABLE]

Meanwhile, the time derivatives of $P^{jk}_{p}$, $Q^{jk}_{p}$, $P^{jk}_{q}$, $Q^{jk}_{q}$ yield

[TABLE]

where all coefficients are calculated at $(P,Q)$. Then one can check

[TABLE]

to get $\sum\limits^{m}_{j=1}dP^{j}\wedge dQ^{j}=\sum\limits^{m}_{j=1}dp^{j}\wedge dq^{j}$. ∎

Because of the symplecticity of the phase flow of a rough Hamiltonian system, it is natural to construct numerical methods to inherit this property. However, the simplified step-$N$ Euler schemes cannot inherit this property. Consequently, in the next three sections we will propose and analyze Runge–Kutta methods for rough Hamiltonian systems.

4. Runge–Kutta methods

Given a time step $h$, we construct $s$-stage Runge–Kutta method by

[TABLE]

with coefficients $a_{\alpha\beta}$, $b_{\alpha}$, $\alpha,\beta=1,...,s$, $t_{k}=kh$, $k=1,...,N$ and $Y^{h}_{0}=z$. Here, we suppose $N=T/h\in\mathbb{N}$ for simplicity. The numerical solution $Y^{h}_{k}$ if exists is an approximation for $Y_{t_{k}}$.

Theorem 4.1*.*

The $s$-stage Runge–Kutta method (4.1) inherits the symplectic structure of a rough Hamiltonian system, if the coefficients satisfy

[TABLE]

Proof.

Denoting the $j$th-component of $(P^{h}_{k},Q^{h}_{k})$ by $(P^{h,j}_{k},Q^{h,j}_{k})$, $j=1,...,m$, then the exterior derivatives of $P^{h,j}_{k+1}$ and $Q^{h,j}_{k+1}$ show that

[TABLE]

Exterior product performed between above equations yield

[TABLE]

Replacing $dP^{h,j}_{k}$ and $dQ^{h,j}_{k}$ in the above equation by the following expressions

[TABLE]

we obtain the symplectic condition

[TABLE]

This ensures the discrete symplectic structure $dP^{h}_{k+1}\wedge dQ^{h}_{k+1}=dP^{h}_{k}\wedge dQ^{h}_{k}$. ∎

Indeed, symplectic Runge–Kutta methods are implicit in general. As we said before, the solvability of symplectic Runge–Kutta methods should be given. Based on the Brouwer’s theorem, we consider the bounded case and the linear case for the vector field $V$ respectively.

Proposition 4.1*.*

If $V=(V_{i})_{0\leq i\leq d}$ is a collection of vector fields in $Lip^{\gamma}$ for some $\gamma>0$, then for arbitrary time step $h>0$, initial value and coefficients $\{a_{\alpha\beta}:\alpha,\beta=1,...,s\}$, the $s$-stage Runge–Kutta method (4.1) has at least one solution for any $\omega$.

Proof.

Fix $h$ and $Y_{k}^{h}$. Denote $Z_{1},...,Z_{s}\in\mathbb{R}^{2m}$ and $Z=(Z_{1}^{\top},...,Z_{s}^{\top})^{\top}\in\mathbb{R}^{2ms}$. We define a map $\phi:\mathbb{R}^{2ms}\rightarrow\mathbb{R}^{2ms}$ with

[TABLE]

Let $c\geq\max\{|a_{\alpha\beta}|:\alpha,\beta=1,...,s\}$, $\nu\geq|V|_{Lip^{\gamma}}$ and $R=\sqrt{s}|Y_{k}^{h}|+s^{2}c\nu|X_{t_{k},t_{k+1}}(\omega)|+1$, we have that for any $|Z|=R$,

[TABLE]

We aim to show that $\phi(Z)=0$ has a solution in the ball $B_{R}:=\{Z:|Z|\leq R\}$. Assume by contradiction that $\phi(Z)\neq 0$ for any $|Z|\leq R$. We define a continuous map $\psi$ by $\psi(Z)=-\frac{R\phi(Z)}{|\phi(Z)|}$. Since $\psi:B_{R}\rightarrow B_{R}$, $\psi$ has at least one fixed point $Z^{*}$ such that $|Z^{*}|=|\psi(Z^{*})|=R$. This leads to a contradiction that $|Z^{*}|^{2}=\psi(Z^{*})^{\top}Z^{*}=-\frac{RZ^{\top}\phi(Z)}{|\phi(Z)|}<0$. Therefore, $\phi$ has at least one solution and then (4.1) is valid. ∎

Things in the linear case are more complicated, then additional conditions should be considered. We focus on the $1$-stage symplectic Runge–Kutta method ($a_{11}=\frac{1}{2},b_{1}=1$), i.e., the midpoint scheme:

[TABLE]

Proposition 4.2*.*

Suppose $V=(V_{i})_{0\leq i\leq d}$ is a collection of linear vector fields of the form $V_{i}(Y)=A_{i}Y$. If $A_{i}$, $i=0,...,d$ are all skew symmetric, i.e., $A_{i}=-A_{i}^{\top}$, then for arbitrary time step $h>0$ and initial value, the the midpoint scheme is solvable for any $\omega$.

Proof.

Fix $h$ and $Y_{k}^{h}$. We define a map $\phi:\mathbb{R}^{2m}\rightarrow\mathbb{R}^{2m}$ with

[TABLE]

Since $A_{i}$, $i=0,...,d$, are all skew symmetric, we have $Z^{\top}A_{i}Z=0$, which leads to

[TABLE]

From Proposition 4.1, we know that $\phi$ has at least a solution in the ball $\{Z:|Z|\leq|Y_{k}^{h}|+1\}$. This implies that the midpoint scheme is solvable. ∎

Remark 4.1*.*

For a linear rough Hamiltonian system, i.e., (3.2) can be rewritten equivalently into the form (3.1), if in addition $A_{i}\in\mathbb{R}^{2m\times 2m}$ is skew symmetric, then $A_{i}$ is in the form of \left(\begin{array}[]{cc}A_{i}^{1}&-A_{i}^{2}\\ A_{i}^{2}&A_{i}^{1}\end{array}\right), where $A_{i}^{1}=-(A_{i}^{1})^{\top}\in\mathbb{R}^{m\times m}$ and $A_{i}^{2}=(A_{i}^{2})^{\top}\in\mathbb{R}^{m\times m}$.

5. Convergence analysis

In Definition 2.1, the solution of an RDE depends on the information about iterated integrals of the noise, which is very difficult if possible to be simulated. The Runge–Kutta methods are implementable because they only make use of the increments of the noise and omit the antisymmetric part of iterated integrals, which is called L$\acute{e}$vy’s area in the Brownian motion setting. To analyze the error caused by this, we introduce piecewise linear approximations, which is a special case of Wong–Zakai approximations. We define $x^{h}_{t}:=(x^{h,0}_{t},...,x^{h,d}_{t})$ with

[TABLE]

and consider the following differential equation in Riemann-Stieltjes integral sense:

[TABLE]

Theorem 5.1*.*

Suppose $X$ and $x^{h}$ are as in Assumption 3.1 and (5.1) respectively. If $V=(V_{i})_{0\leq i\leq d}$ is a collection of vector fields in $Lip^{\gamma}$ for some $\gamma>2\rho$, or a collection of linear vector fields of the form $V_{i}(Y)=A_{i}Y$, then both (3.2) and (5.2) have unique solutions $Y$ and $Y^{h}$ almost surely.

Further, let $\theta\geq|z|$, $\nu\geq\max\{|V|_{Lip^{\gamma}},|A_{i}|,i=0,,,,d\}$, then for any $0\leq\eta<\min\{\frac{1}{\rho}-\frac{1}{2},\frac{1}{2\rho}-\frac{1}{\gamma}\}$, there exits a finite random variable $C(\omega)$ and a null set $M$ such that

[TABLE]

holds for any $h\geq 0$. $C(\omega)$ depends on $\rho$, $\eta$, $\gamma$, $\nu$, $\theta$, $K$ and $T$.

Proof.

For the first case of $V$, see Theorem 6 and Corollary 8 in [6]. $C(\omega)$ depends on $\rho$, $\eta$, $\gamma$, $\nu$, $K$ and $T$.

For the linear case, we know if $p>2\rho$, then for any $\omega\in\Omega\setminus M$, $\|S_{[p]}(x^{h})(\omega)\|_{p\text{-}var;[0,T]}$ has an upper bound for all $h\geq 0$ (see [7, Theorem 15.33]). According to the estimate in (2.3), there exists some constant $R>0$ depending on $\omega$, $\rho$, $\theta$, $\eta$, $\nu$ and $T$, such that it is a bound for $Y_{t}$ and $y^{h}_{t}$, for all $h\geq 0$ on $[0,T]$. Using again localization, we replace $V$ by a compactly supported $\tilde{V}\in Lip^{\gamma}$ which coincides with $V$ in the ball $B_{R}:=\{y:|y|\leq R\}$ and then the result can be generalized from the first case to this linear case. ∎

We are in position to analyze the convergence rate. If $V\in Lip^{\gamma}$, we analyze the midpoint scheme for $\rho\in[1,3/2)$ in Theorem 5.2 and two $s$-stage Runge–Kutta methods ($s>1$) for $\rho\in[1,2)$ in Corollary 5.1. If $V$ is linear, we focus on the convergence rate of the midpoint scheme for $\rho\in[1,3/2)$ in Theorem 5.3, using the uniform boundedness of numerical solutions in Proposition 5.1.

Theorem 5.2*.*

Suppose $\rho\in[1,3/2)$ and $Y_{k}^{h}$ is the numerical solution to the midpoint scheme (4.2). If $|V|_{Lip^{\gamma}}\leq\nu<\infty$ with $\gamma>2\rho$, then for any $0<\eta<\min\{\frac{3}{2\rho}-1,\frac{1}{2\rho}-\frac{1}{\gamma}\}$, there exists a finite random variable $C(\omega)$ and a null set $M$ such that

[TABLE]

holds for any $0<h\leq 1$. $C(\omega)$ depends on $\rho$, $\eta$, $\gamma$, $\nu$, $K$ and $T$.

Proof.

We divide the error into two parts

[TABLE]

The second part has been estimated in Theorem 5.1, then it suffices to estimate the first part.

$\it{Step\ 1.\ Local\ error.}$ According to the definition of $x^{h}_{t}$, for $k=1,...,N$, $n\geq 1$, $i_{1},...,i_{n}\in\{0,...,d\}$, we have

[TABLE]

which concludes $\int_{t_{k-1}\leq u_{1}<\cdot\cdot\cdot<u_{n}\leq t_{k}}dx^{h}_{u_{1}}\otimes\cdot\cdot\cdot\otimes dx^{h}_{u_{n}}=X_{t_{k-1},t_{k}}^{\otimes n}/n!$. Taylor expansion yields

[TABLE]

with

[TABLE]

As for $V(Y^{h}_{1/2})$ in (4.2), we have

[TABLE]

where $\tilde{Y}^{h}_{0}$ is between $Y^{h}_{0}$ and $Y^{h}_{1/2}$. Then,

[TABLE]

with

[TABLE]

Noticing $|Y^{h}_{1}-Y^{h}_{0}|\leq|V(Y^{h}_{1/2})\|X_{0,t_{1}}|$ and $\gamma>2$, we obtain the difference between $Y^{h}_{1}$ and $y^{h}_{t_{1}}$ is

[TABLE]

For any $0<\eta<\frac{3}{2\rho}-1$, there exists a $p>2\rho$, such that $\eta=\frac{3}{p}-1$. Since for any $\omega\in\Omega\setminus M$, $|X_{0,t_{1}}(\omega)|=|\pi_{1}(\textbf{X}_{0,t_{1}}(\omega))|\leq\|\textbf{X}(\omega)\|_{1/p\text{-}{\rm H\ddot{o}l};[0,T]}h^{1/p}$, the local error is

[TABLE]

$\it{Step\ 2.\ Global\ error.}$ We use the notation $\pi(t_{0},y_{0},x)_{t}$, $t\geq t_{0}$, to represent the solution of (5.2) with the driven noise $x$ and the initial value $y_{0}$ at time $t_{0}$. For any $\omega\in\Omega\setminus M$, we know that $\|S_{[p]}(x^{h})(\omega)\|^{p}_{p\text{-}var;[0,T]}$ has an upper bound for all $h>0$ (see [7], Theorem 15.33). Since $\gamma>2\rho$, similar to Theorem 10.30 in [7], there exists some constant $C=C(\gamma,p)$ such that for any $k=1,...,N$,

[TABLE]

Therefore, for any $\omega\in\Omega\setminus M$,

[TABLE]

The fact that $\frac{3}{2\rho}-1\leq\frac{1}{\rho}-\frac{1}{2}$ for $\rho\in[1,3/2)$ and Theorem 5.1 complete the proof. ∎

When the components of the noise are independent standard Brownian motions, i.e., $\rho=1$, the convergence order is (almost) consistent with the mean square convergence rate of SDEs in Stratonovich sense with multiplicative noise. Moreover, note that the estimate in Theorem 5.2 is valid for the case when $\rho\in[1,3/2)$, which is similar to the simplified step-2 Euler scheme in [6]. To construct a numerical scheme for $\rho\in[1,2)$ and fill the gap between the convergence rates of two parts in (5.3), we use a higher stage symplectic Runge–Kutta method with local order $\tau\geq 4$ when applied to classical ordinary differential equations. If $\gamma>\max\{2\rho,\tau-1\}$, then the estimate in (5.4) for the global error will hold for any $0<\eta<\frac{\tau}{2\rho}-1$. Two symplectic Runge–Kutta methods with $\tau=5$ and $\tau=4$ are given in Corollary 5.1 to illustrate this point.

Corollary 5.1*.*

Suppose $\rho\in[1,2)$ and $Y_{k}^{h}$ represents the solution to the Runge–Kutta methods with coefficients expressed in the Butcher tableaus below.

Method I:

[TABLE]

Method II:

[TABLE]

where $a=1.351207$ is the real root of $6x^{3}-12x^{2}+6x-1=0$.

If $|V|_{Lip^{\gamma}}\leq\nu<\infty$ with $\gamma>4$ for Method I ($\gamma>\max\{2\rho,3\}$ for Method II), then for any $0<\eta<\min\{\frac{1}{\rho}-\frac{1}{2},\frac{1}{2\rho}-\frac{1}{\gamma}\}$, there exists a finite random variable $C(\omega)$ and a null set $M$ such that

[TABLE]

In addition, both of the two methods inherit the symplectic structure of a rough Hamiltonian system.

Proof.

According to [9], we know that the corresponding $\eta$ in (5.4) belongs to $(0,\frac{5}{2\rho}-1)$ for Method I and $(0,\frac{2}{\rho}-1)$ for Method II. This result together with Theorem 5.1 finally leads to the global error estimate. Moreover, the coefficients satisfy the symplectic condition in Theorem 4.1. ∎

Remark 5.1*.*

This convergence rate equals to the simplified step-3 Euler scheme for $\gamma$ large enough. Actually, using Runge–Kutta methods with larger $\tau$ will not improve the convergence rate since the error of piecewise linear approximations persists.

The next proposition is essential to get the convergence rate of the midpoint scheme in the linear case.

Proposition 5.1*.*

If $V$ in (3.2) is a collection of skew symmetric linear vector fields of the form $V_{i}(Y)=A_{i}Y$, then numerical solutions from the midpoint scheme (4.2) are uniformly bounded. More precisely, $|Y^{h}_{k}|=|z|$, $k=1,...,N$.

Proof.

For any $k=0,...,N-1$, since $Y^{h}_{k}=Y^{h}_{k+1/2}-\frac{1}{2}\sum^{d}_{i=0}A_{i}Y^{h}_{k+1/2}X_{t_{k},t_{k+1}}^{i},$ we have

[TABLE]

Then, $|Y^{h}_{k+1}|^{2}=|Y^{h}_{k}|^{2}=|z|$ from $(Y^{h}_{k+1/2})^{\top}A_{i}Y^{h}_{k+1/2}=0$. ∎

Theorem 5.3*.*

Suppose $\rho\in[1,3/2)$. Let $V$ be as in Proposition 5.1 and $Y_{k}^{h}$ represents the numerical solution to the midpoint scheme. If $\theta\geq|z|$, $\nu\geq\max\{|A_{i}|:i=0,...,d\}$, then for any $\gamma>2\rho$, $0<\eta<\min\{\frac{3}{2\rho}-1,\frac{1}{2\rho}-\frac{1}{\gamma}\}$, there exists a finite random variable $C(\omega)$ and a null set $M$ such that

[TABLE]

holds for any $0<h\leq 1$. $C(\omega)$ depends on $\rho$, $\eta$, $\gamma$, $\nu$, $\theta$, $K$ and $T$.

Proof.

Proposition 5.1 implies that there exists some constant $R>0$ depending on $\omega$, $\rho$, $\theta$, $\eta$, $\nu$ and $T$, such that $|\pi(t_{k-1},Y^{h}_{k-1},x^{h})_{t}|\leq R$ for all $h$, $k$, $t$. Based on this result, we know that the localization skill is applicable. Combining Theorem 5.1 and Theorem 5.2, we complete the proof. ∎

6. Numerical experiments

In this section, we illustrate our theoretical results by two examples. In the first example, we verify the theoretical convergence rate given in Corollary 5.1. In the second one, we consider the Kubo oscillator and compare the performance of the midpoint scheme, the simplified step-$2$ Euler scheme and the simplified step-$3$ Euler scheme.

6.1. Example 1

We consider a rough Hamiltonian system in $\mathbb{R}^{2}$ driven by a two-dimensional fBm with independent components and suppose

[TABLE]

More precisely, the corresponding equations are

[TABLE]

Convergence rate. Since the vector fields are bounded with bounded derivatives, the theoretical convergence rates are almost $0.3,0.2,0.1$ for both Method I and Method II in Corollary 5.1 when the Hurst parameter $H$ are $0.4,0.35,0.3$ respectively. Figure 1 shows the maximum error in the discretisation points, $\max_{k=1,...,N}|Y^{h}_{k}(\omega)-Y_{t_{k}}(\omega)|$, which we call the pathwise maximum error for short, for Method I. Pictures in the same row are from three sample paths. Numerical results are in acceptable accordance with theoretical analysis. In addition, the fact that the error becomes larger as $H$ increases indicates the influence of the roughness of the noise.

6.2. Example 2

For the linear case, we consider the Kubo oscillator driven by a three-dimensional fBm with independent components and $H=0.4$:

[TABLE]

We have the expression of the exact solution

[TABLE]

[TABLE]

Denoting $Y=(P,Q)^{\top}\in\mathbb{R}^{2}$ and $X^{0}:=t$, we have that (6.1) is equivalent to $dY=\sum^{d}_{i=0}A_{i}YdX^{i}$, where $A_{0}$=\left(\begin{array}[]{cc}0&-1\\ 1&0\end{array}\right) and $A_{i}$=\left(\begin{array}[]{cc}0&-\epsilon\\ \epsilon&0\end{array}\right), $i=1,2,3$. Since $A_{i}$ are all skew symmetric, the midpoint scheme is solvable. We compare it with the simplified step-$2$ and step-$3$ Euler schemes in the following experiments.

Evolution of domains. Figure 2 shows the evolution of domains in the phase plane for one sample path. The initial domain is a square with four corners at $(1,1)$, $(2,1)$, $(2,2)$ and $(1,2)$. Images at $t=0.4,1.6,8$, are presented under the exact mapping and the three numerical methods. One can observe that the images are still squares. The exact mapping is area preserving, which is equivalent to that it preserves the symplectic structure of the system. This property is inherited by the midpoint scheme, since the area of red squares are the same. However, the area of images under the simplified step-$2$ and step-$3$ Euler schemes increase and decrease respectively, which proves that they are not symplectic methods.

Conservation of quadratic invariant. For the initial value $(1,1)$, the exact solution is on the circle with center at the origin and radius $r=\sqrt{2}$, which implies another invariant of this system (see also [1, Example 1]). We present the phase trajectory for three sample paths in Figure 3. The midpoint scheme preserves this property such that its solutions are always on that circle, which is proved in Proposition 5.1. The phase trajectory of the simplified step-$2$ Euler scheme deviates from the exact one a lot and that of the simplified step-$3$ Euler scheme shrinks gradually. Figure 2 and Figure 3 verify that the midpoint scheme, as a symplectic Runge–Kutta method, preserves the invariants of the system (6.1) and show its long time behavior is better than another two methods.

Stability for $\epsilon$ and $T$. We investigate the influence of the size $\epsilon$ of the noise and the total time $T$ in Figure 4. When $T=10$, the simplified step-$2$ scheme is unstable. The pathwise maximum error of the midpoint scheme is the smallest among the three methods. As $\epsilon$ and $T$ become larger, this advantage is more clearer. This result implies the stability of an implicit method.