Runge-Kutta methods for rough differential equations
Martin Redmann Sebastian Riedel
math.NA cs.NA math.PR 60H10, 60H35, 60L20, 60L70, 65C30

Runge-Kutta for RDEs]Runge-Kutta methods for rough differential equations


Martin Redmann, Martin Luther University Halle-Wittenberg, Institute of Mathematics, Theodor-Lieser-Str. 5, 06120 Halle (Saale), Germany
martin.redmann@mathematik.uni-halle.de


Sebastian Riedel 

Institut für Mathematik, Technische Universität Berlin, Germany and Weierstraß-Institut, Berlin, Germany
riedel@math.tu-berlin.de

 
2020 Mathematics Subject Classification. 60H10, 60H35, 60L20, 60L70, 65C30.     
 










  We study Runge-Kutta methods for rough differential equations which can be used to calculate solutions to stochastic differential equations driven by processes that are rougher than a Brownian motion. We use a Taylor series 
  representation (B-series) for both the numerical scheme and the solution of the rough differential equation in order to determine conditions that guarantee the desired order of the local error for the 
  underlying Runge-Kutta method. Subsequently, we prove the order of the global error given the local rate. In addition, we simplify the numerical approximation by introducing a Runge-Kutta scheme 
  that is based on the increments of the driver of the rough differential equation. This simplified method can be easily implemented and is computational cheap since it is derivative-free. 
  We provide a full characterization of this implementable Runge-Kutta method meaning that we provide necessary and sufficient algebraic conditions for an optimal order of convergence in case that the driver, e.g., is 
  a fractional Brownian motion with Hurst index 1/4 < H ≤1/2. We conclude this paper by conducting numerical experiments verifying the theoretical rate of convergence.


[
    S. Riedel
    
=============






§ INTRODUCTION


Ordinary differential equations (ODEs) have many real life applications. They, e.g., describe chemical, physiological and ecological processes or they appear as spatially discretized 
partial differential equations like the heat equation. Often analytic solutions to ODEs do not exist which requires numerical approximations in order to solve these equations. 
An important class of such schemes are Runge-Kutta methods <cit.> which can be of arbitrary order of convergence.
These are often preferred in practice since they are derivative-free in contrast to Taylor methods. Computing derivatives of the right hand side function 
f_0 of an ODE can either be very costly or closed form expressions might not be available.

However, in many applications uncertainties need to be taken into account.
Therefore, for a more accurate modeling in such cases, a noise term can be added to an ODE leading to stochastic differential equations (SDEs). Runge-Kutta schemes for SDEs driven by a Brownian motion have already 
been established, see, e.g., <cit.>.

Lyons' rough paths theory provides an alternative way to SDEs which goes far beyond the scope of usual Itō equations. In this paper, we are interested in studying numerical schemes to solve 
rough differential equations (RDEs) of the form

    dy(t) = f_0(y(t))  dt + f(y(t))  d𝐗(t),   y(t_0) = y_0 ∈^m,

where 𝐗 is a suitable rough path above some α-Hölder path X = (X^1,…,X^d)  [0,T] →^d, 

f = (f_1,…,f_d) and f_i ^m →^m are vector fields for every i = 0,…,d. Such equations represent SDEs driven by stochastic processes that are potentially rougher than a
Brownian motion if 𝐗 is a random rough path, i.e. a stochastic process with sample paths lying in a rough paths space. One benefit of rough paths theory compared to Itō's theory
is that one is not restricted to the martingale framework. In fact, there is a large class of stochastic processes which have “natural extensions” to rough paths valued processes, cf. <cit.>.
For instance, many Gaussian processes possess such a “natural lift” including the fractional Brownian motion with Hurst parameter H > 1/4, but even more general processes like the bifractional Brownian motion,
Volterra processes or processes which can be represented by random Fourier series, cf. <cit.> for a discussion.

In the context of rough paths theory, numerical schemes are indispensable when simulating the solution to an RDE driven by a random rough path or when discretizing rough stochastic partial differential
equations <cit.>. In fact, numerical schemes played a fundamental role in rough paths theory from the very beginning. This is probably most visible in the work of Davie <cit.>,
where the Milstein scheme is used to solve rough differential equations theoretically. This approach was generalized to higher order Taylor-type schemes by Friz and Victoir <cit.>. However,
in a stochastic context, these schemes are of little use in practice since they contain iterated stochastic integrals whose distribution is unknown in general. 
To overcome this difficulty, Deya, Neuenkirch and Tindel introduced so-called simplified schemes in <cit.> in which the iterated stochastic integrals are replaced by products of increments of the driving process.
These schemes were successfully used in different contexts, cf. e.g. <cit.>. 


However, as Taylor methods, these numerical approximations suffer from the need to calculate or simulate derivatives of the vector fields f_k. As mentioned above, even if the derivatives are available,
this can be very expensive especially in a large scale setting (e.g. spatially discretized rough partial differential equation). Moreover, the simplified scheme is difficult to implement in general. 
Therefore, we see the need of studying Runge-Kutta methods for rough differential equations that can easily be implemented and are derivative-free. 

Our approach to establish Runge-Kutta methods is classical, both in the deterministic and the stochastic context: First, we define a class of equations which can be expanded in a B-series. 

Second, we have to find a B-series representation of the equation (<ref>). Comparing both series, we can, in principle, deduce the order conditions of the Runge-Kutta method by matching their coefficients.
A B-series representation of an ODE contains combinations of products and derivatives of the defining vector field which can be described in the language of trees. For SDEs, integrated products of iterated stochastic integrals
have to be considered in addition which can be described in the same language. We call such objects tree-iterated (stochastic) integrals in the sequel. 
A rough path in the sense of Lyons <cit.> is a collection of objects which “mimic” the iterated integrals of the underlying path. 
Lyons' theory is able to solve differential equations driven by geometric rough paths, i.e. those which obey the usual chain rule. Gubinelli realized in <cit.> that one can even solve rough differential equations
driven by non-geometric rough paths if one additionally assumes that all tree-iterated integrals are known. He calls such objects branched rough paths. Thinking of B-series representations of SDEs,
this is a very natural approach to solve equations of the form (<ref>). For us, it is therefore reasonable to use his theory and to interpret the equation (<ref>) as a
rough differential equation driven by a branched rough path. Doing this, we are able to deduce the B-series expansion of (<ref>) in Theorem <ref>.
Comparing both B-series and matching their coefficients up to a given order for an arbitrary multidimensional driving process X and its tree-iterated integrals can be very hard,
cf.  <cit.> for a 2-dimensional example. However, we already pointed out that in practice, one is not able to simulate the tree-iterated integrals anyway. 
We therefore make the same ansatz as in <cit.> and replace the tree-iterated integrals by products of increments. This simplifies the task of matching the coefficients a lot, and one is able to deduce
the order conditions, in principle, up to any order, cf. Theorem <ref> and the subsequent remark. We call such schemes simplified Runge-Kutta methods.
As in <cit.>, the Wong-Zakai error plays a fundamental role in their convergence analysis. Loosely speaking, our main result is the following:


 Let 𝐗 be an α-Hölder rough path (branched or geometric) and assume that f_0 and f are sufficiently smooth and bounded with bounded derivatives.
 If the Wong-Zakai error to approximate (<ref>) is of order r_0, a simplified Runge-Kutta method (<ref>) of order p converges with rate . 


As an application, we can study the scheme when the driving process is a fractional Brownian motion with Hurst parameter H ∈ (1/4,1/2]. In this case, the Wong-Zakai error is arbitrarily close
to 2H -1/2, cf. <cit.>. We therefore obtain:


 For a fractional Brownian motion with Hurst parameter H ∈ (1/4,1/2], a simplified Runge-Kutta scheme of order 3 converges with rate arbitrarily close to 2H - 1/2.


We already pointed out that numerical schemes studied in the context of rough paths theory are mostly of Taylor-type. To our knowledge, the only exception is the article by
Hong, Huang and Wang <cit.> where a class of symplectic Runge-Kutta methods is considered to solve Hamiltonian equations driven by Gaussian processes.
Our article differs from  <cit.> in several regards. On the technical level, no B-series are used in  <cit.>, the authors have to prove all necessary
estimates “by hand” in the framework of geometric rough paths. Consequently, they do not provide general order conditions. For instance, no explicit Runge-Kutta methods are deduced in <cit.>.
Moreover, their approach is probably hard to generalize to schemes of arbitrary order, whereas our approach does not have any limitations in this regard. 

The article is structured as follows. In Section <ref>, we define the equations which can be expanded to obtain the desired B-series.
Section <ref> explains the concept of branched rough paths, deduces the B-series representation for equation (<ref>) and discusses the local error of full Runge-Kutta methods.
Simplified Runge-Kutta methods are defined in Section <ref>, where the necessary order conditions are derived to obtain the local error of the numerical scheme.
In Section <ref>, we deduce the global error for our methods. The article closes with numerical experiments presented in Section <ref>. 

Let us finally mention that in the whole article, we will discard the drift in (<ref>) and consider equations of the form

    dy(t) = f(y(t))  d𝐗(t),   y(t_0) = y_0,

only which simplifies the exposition a lot. Furthermore, this is not a real limitation if we assume that the first component of X is just the path t ↦ t.  
















































§ NOTATION AND BASIC DEFINITIONS




 §.§ General notation
 Let I be an interval in  and V be a linear space. We call a function X  I → V a path and X_s,t := X(t) - X(s) an increment. For a general two-parameter function X  I × I → V, we will often write X_s,t instead of X(s,t). If |·| is a norm on V and X  I × I → V, we set

    X_α := X_α;I := sup_s,t ∈ Is ≠ t|X_s,t|/|t-s|^α

for α∈ (0,1] and call it the α-Hölder (semi-)norm of X. For x ∈, we use the notation ⌊ x ⌋ := max{k ∈ :  k ≤ x }.  
Let γ > 0 and γ = [γ] + {γ} where [γ] is an integer and {γ}∈ (0,1]. We will say that a vector field f ^m →^m belongs to the class Lip^γ if f is [γ]-times continuously differentiable and the [γ]-th derivative is locally {γ}-Hölder continuous. f is of class Lip^γ_b if, in addition, f and all its derivatives are bounded and if the [γ]-th derivative is globally {γ}-Hölder continuous. More generally, a collection of vector fields f = (f_1,…,f_d) is of class Lip^γ resp. Lip^γ_b if every f_i, i = 1, …, d, is of class Lip^γ resp. Lip^γ_b.





 §.§ Trees and the Connes-Kreimer Hopf algebra
 Let 𝒯 be the set of all rooted, labeled trees with vertex decorations from the set {1,…,d}. We will use a recursive definition to construct trees. The empty tree will be denoted by 1. We use the convention that 1 ∉𝒯 and set 𝒯^0 := 𝒯∪{1}. If τ_1,…τ_m ∈𝒯^0, [τ_1 ⋯τ_m]_a denotes the tree obtained by attaching all trees τ_1,…τ_m to a new vertex which we label by a  ∈{1,…,d}. We use the notation ∙ _a = [1]_a for the single vertex tree with label a. The order of the branches of the tree does not matter, i.e. [τ_σ(1)⋯τ_σ(m)]_a = [τ_1⋯τ_m]_a holds for every permutation σ. If τ∈𝒯^0 is a tree, |τ| denotes the number of vertices. The set 𝒯_N consists of all trees τ∈𝒯 such that |τ| ≤ N and we set 𝒯_N^0 := 𝒯_N ∪{1}. We let

    ℱ := {τ_1 ⋯τ_m  : τ_i ∈𝒯, m ∈}

denote the set of unordered forests and define ℱ^0 := ℱ∪{1}. The map | · | is extended to ℱ^0 by setting

    |τ_1 ⋯τ_m| := |τ_1| + … + |τ_m|.

As before, ℱ_N contains all 𝔥∈ℱ with |𝔥| ≤ N and ℱ_N^0 := ℱ_N ∪{1}.



We define (ℋ,·) to be the commutative polynomial algebra generated by the variables 𝒯. Alternatively, we can view ℋ as the real vector space spanned by the elements in ℱ^0. A coproduct Δℋ→ℋ⊗ℋ is defined recursively by setting Δ 1 := 1 ⊗ 1 and 

    Δ [τ_1 ⋯τ_m]_a := [τ_1 ⋯τ_m]_a ⊗ 1 + (id⊗ B_+^a) (Δτ_1 ⋯Δτ_m)

for a tree [τ_1 ⋯τ_m]_a where B_+^a is the operator defined by B_+^a(τ_1 ⋯τ_n) := [τ_1 ⋯τ_n]_a on the forest  τ_1 ⋯τ_n. We then extend the definition to forests by setting Δ (τ_1 ⋯τ_m) := Δτ_1 ⋯Δτ_m and eventually define Δ on ℋ by linear extension. We will use Sweedler's notation 

    Δ𝔥 = ∑_(𝔥)𝔥^(1)⊗𝔥^(2).

One can also construct an antipode S ℋ→ℋ, i.e. a map which satisfies

    M(id⊗ S) Δ x =  M(S  ⊗id) Δ x = x

for every x ∈ℋ where M(x ⊗ y) := xy. Then, (ℋ,·,Δ,S) is called the Connes-Kreimer Hopf algebra <cit.>, cf. also <cit.>. The dual Hopf algebra will be denoted by (ℋ^*,⋆,δ,S^*). 

For a general account on Hopf algebras, cf. <cit.> or <cit.>. The Connes-Kreimer Hopf algebra is also discussed in <cit.>.













§ THE FULL RUNGE-KUTTA METHOD


In this section, we will define s-stage Runge-Kutta methods. We follow the approach developed by Burrage and Burrage in <cit.>. Let Z^(1), …, Z^(d) be given s× s-matrices and z^(1),…, z^(d) vectors in ^s.
For given y_n ∈^m, consider the equations

    Y_i    = y_n + ∑_k = 1^d ∑_j=1^s Z^(k)_ij f_k(Y_j) 
    
          y_n+1   = y_n + ∑_k = 1^d ∑_i=1^s z^(k)_i f_k(Y_i).

In applications, Z and z can (and will) be random. Moreover, both values will depend on the step size h>0 of the numerical scheme (<ref>). Note that the first equation can be implicit in which case the existence of a solution is not guaranteed. In fact, this question will depend on the properties of the vector fields f_i. For instance, it can be shown that solutions exist in case that all vector fields are bounded, cf. <cit.>. However, we will not address this question here and just assume that solutions exist.

Set 
Φ(1)(h):= 1_s := (1,⋯,1)^T ∈^s for i = 1,…,d and for a tree τ = [τ_1 ⋯τ_n]_i,

    Φ(τ)(h) := Π_j=1^n (Z^(i)Φ(τ_j)(h)),   
     a(τ)(h) := ⟨ z^(i), Π_k=1^nΦ(τ_k)(h) ⟩.


Notice that above, the product of two vectors has to be understood component-wise. 


 Let f = (f_1,…,f_d) be sufficiently smooth such that all derivatives below exist. For a tree τ∈𝒯^0, we define the elementary differentials F(τ) ^m →^m recursively by setting 
 
  
  (i) F(1)(y) := y,
  
  (ii) F(∙_i)(y) := f_i(y)  and
  
  (iii) F(τ)(y) := f_i^(n)(y) (F(τ_1)(y),…,F(τ_n)(y)) for a tree τ = [τ_1 ⋯τ_n]_i where f_i^(n) denotes the n-th total derivative of f_i
 
 for y ∈^m.










Next,we define some combinatoric quantities. For unlabeled trees, we set

    γ(1) = 0,  γ(∙) = 1,  γ([τ_1,…,τ_k]) = |[τ_1,…,τ_k]| ∏_i=1^k γ(τ_i),

and

    β(1) = 1,   β(∙) = 1,   β(τ ) := |τ| - 1|τ_1|,…, |τ_k|1/r_1! ⋯ r_q!∏_j = 1^k β(τ_j)

where τ = [τ_1,…,τ_k] = [(τ_1)^r_1,…,(τ_q)^r_q], τ_1,…,τ_q being pairwise distinct trees. For labeled trees, we use the same definition. The main result in <cit.> we are going to use is the following:

 The Taylor series expansion of (<ref>) is
 
    y_1 = y_0 + ∑_τ∈𝒯γ(τ)/|τ|!β(τ)  a(τ)(h) F(τ)(y_0).



 <cit.>.

The coefficients in (<ref>) are sometimes noted in a different form which we recall now. Following <cit.>, we define the symmetry factor σ for unlabeled trees as

    σ(1) = 1,   σ(∙) = 1,   σ(τ) = r_1! ⋯ r_q! ∏_j = 1^k σ(τ_j)

where τ = [τ_1,…,τ_k] = [(τ_1)^r_1,…,(τ_q)^r_q] with τ_1,…,τ_q being pairwise distinct trees. For labeled trees, the same definition is used.


  For every tree τ∈𝒯^0,
 
    1/σ(τ) = γ(τ)/|τ|!β(τ).





We prove this lemma by induction on the height of τ for unlabeled trees. The equality is true for τ = 1 and τ = ∙. Let us assume that the claim is true for each sub-tree of 
τ = [τ_1,…,τ_k] = [(τ_1)^r_1,…,(τ_q)^r_q]. Then, we have

    β(τ)    = (|τ| - 1)!/|τ_1|! ⋯ |τ_k|!1/r_1! ⋯ r_q!∏_j = 1^k β(τ_j) 
       = (|τ| - 1)!/|τ_1|! ⋯ |τ_k|!1/r_1! ⋯ r_q!∏_j = 1^k |τ_j|!/γ(τ_j)1/σ(τ_j)
       = (|τ| - 1)!/r_1! ⋯ r_q!|τ|/γ(τ)∏_j = 1^k 1/σ(τ_j)
       = |τ|!/γ(τ)1/σ(τ).



 




§ BRANCHED ROUGH PATHS AND B-SERIES EXPANSION FOR ROUGH DIFFERENTIAL EQUATIONS


In this section, we recall the concept of a branched rough paths introduced by Gubinelli in <cit.>. We use a similar approach and notation as Hairer and Kelly in <cit.>. Our main goal is to deduce the B-series expansion of (<ref>) which we will eventually achieve in Theorem <ref>. Note that Hairer and Kelly state a similar result in <cit.>. However, we can not use their result for two reasons: First, it is not quantitative, i.e. the order of the truncation error is not specified. Second, it is not entirely correct since the expansion in <cit.> lacks the symmetry factor. The proof of <cit.> was corrected in the Master thesis of Rosa Preiß, and we are grateful to her for providing us with the corrected version.




 Let α∈ (0,1]. A α-branched rough path is a map 𝐗 [0,T] × [0,T] →ℋ^* such that 
 
  
  * for all s, t ∈ [0,T] and all 𝔥_1,𝔥_2 ∈ℋ,
  
    ⟨𝐗_s,t, 𝔥_1 ⟩⟨𝐗_s,t,𝔥_2 ⟩ = ⟨𝐗_s,t, 𝔥_1 ·𝔥_2 ⟩

  
  
  * for all s,u,t ∈ [0,T],
  
    𝐗_s,t = 𝐗_s,u⋆𝐗_u,t

  
  
  * for every τ∈𝒯,
  
    sup_s ≠ t | ⟨𝐗_s,t, τ⟩ |/|t-s|^α |τ| < ∞.




 
 
 The space of α-branched rough paths will be denoted by 𝒞^α([0,T],^d). It is a complete metric space with metric
 
    ϱ_α(𝐗,𝐘) := ∑_τ∈𝒯_Nsup_s ≠ t | ⟨𝐗_s,t - 𝐘_s,t,τ⟩ |/|t-s|^α |τ|

 where N = ⌊ 1/α⌋.






 Let X = (X^1 ,…, X^d)  [0,T] →^d be a piece-wise C^1-path. We define 𝐗 by setting
 
    ⟨𝐗_s,t, ∙_i ⟩ = X^i(t) - X^i(s)   and  ⟨𝐗_s,t, [τ_1 ⋯τ_n]_i ⟩ = ∫_s^t ⟨𝐗_s,u, τ_1 ⟩⋯⟨𝐗_s,u, τ_n ⟩  dX^i(u)

 for any tree τ = [τ_1 ⋯τ_n]_i ∈𝒯 and
 
    ⟨𝐗_s,t, τ_1 ⋯τ_n ⟩ = ⟨𝐗_s,t, τ_1 ⟩⋯⟨𝐗_s,t, τ_n ⟩

 for any forest τ_1 ⋯τ_n. We can now extend 𝐗 linearly to a map on ℋ and therefore obtain a map 𝐗 [0,T] × [0,T] →ℋ^*. It can be shown that this map defines a α-branched rough path for every α∈ (0,1]. If X is α-Hölder continuous for some α∈ (1/2,1], we can use the Young integral to define 𝐗 as above. In this case, it can be shown that 𝐗 defines a α'-branched rough path for every α' ∈ (0,α]. In this case, 𝐗 is called the natural lift of the (smooth) path X to a branched rough path.




Next, we define a class of paths which we can integrate against a branched rough path.



 Let 𝐗 be a α-branched rough path and N = ⌊ 1/α⌋. A path 𝐙 [0,T] →ℋ_N-1 satisfying
 
    ⟨𝔥, 𝐙(t) ⟩ = ⟨𝐗_s,t⋆𝔥, 𝐙(s) ⟩ + R^𝔥_s,t

 for each 𝔥∈ℱ^0_N-1 where |R^𝔥_s,t| ≤ C |t-s|^(N - |𝔥|)α is called controlled by 𝐗. We will also say that 𝐙 is a controlled path above the path t ↦ Z(t) := ⟨ 1, 𝐙(t) ⟩. More generally, we call a path 𝐙 [0,T] → (ℋ_N-1)^m controlled by 𝐗 if (<ref>) holds, understood as an equation in ^m. The space of controlled paths 𝐙 [0,T] → (ℋ_N-1)^m will be denoted by 𝒬_𝐗(^m) which is a Banach space with the norm
 
    𝐙_𝒬_𝐗(^m) := |𝐙(0)| + ∑_𝔥∈ℱ^0_N-1R^𝔥_(N - |𝔥|)α.





The following lemma is given as an exercise in <cit.>. We provide a full proof here for the reader's convenience.



 Let 𝐗 be a α-branched rough path and 𝐙 be controlled by 𝐗. Set 
 
    Z̃_s,t := ∑_𝔥∈ℱ_N-1^0⟨𝔥, 𝐙(s) ⟩⟨𝐗_s,t,[𝔥]_i⟩.

 Then,
 
    Z̃_s,t - Z̃_s,u - Z̃_u,t = - ∑_𝔥∈ℱ^0_N-1⟨𝐗_u,t, [𝔥]_i ⟩ R^𝔥_s,u.






 Let 𝔥̃∈ℱ^0_N-1. Then,
 
    ⟨𝐗_s,t⋆𝔥̃, 𝐙(s) ⟩   = ∑_𝔥∈ℱ_N-1^0⟨𝐗_s,t⋆𝔥̃, 𝔥⟩⟨𝔥, 𝐙(s) ⟩
       = ∑_𝔥∈ℱ_N-1^0⟨𝐗_s,t⊗𝔥̃, Δ𝔥⟩⟨𝔥, 𝐙(s) ⟩
       = ∑_𝔥∈ℱ_N-1^0∑_(𝔥)⟨𝐗_s,t⊗𝔥̃, 𝔥^(1)⊗𝔥^(2)⟩⟨𝔥, 𝐙(s) ⟩
       = ∑_𝔥∈ℱ_N-1^0∑_(𝔥)1_{𝔥^(2) = 𝔥̃}⟨𝐗_s,t, 𝔥^(1)⟩⟨𝔥, 𝐙(s) ⟩.

    For 𝔥∈ℱ^0_N-1,
    
    ⟨𝐗_s,t, [𝔥]_i ⟩   = ⟨𝐗_s,u⋆𝐗_u,t, [𝔥]_i ⟩
       = ⟨𝐗_s,u⊗𝐗_u,t, Δ [𝔥]_i ⟩
       = ⟨𝐗_s,u⊗𝐗_u,t, [𝔥]_i ⊗ 1 ⟩ + ⟨𝐗_s,u⊗𝐗_u,t, (id⊗ B_+^i) Δ𝔥⟩
       = ⟨𝐗_s,u , [𝔥]_i  ⟩ + ∑_(𝔥)⟨𝐗_s,u, 𝔥^(1)⟩⟨𝐗_u,t, [𝔥^(2)]_i ⟩.

    
    It follows that
    
    Z̃_s,t - Z̃_s,u - Z̃_u,t   = ∑_𝔥∈ℱ^0_N-1⟨𝔥, 𝐙(s) ⟩( ⟨𝐗_s,t, [𝔥]_i ⟩ - ⟨𝐗_s,u, [𝔥]_i ⟩) - ⟨𝔥, 𝐙(u) ⟩⟨𝐗_u,t,[𝔥]_i⟩
       = ∑_𝔥∈ℱ^0_N-1∑_(𝔥)⟨𝔥, 𝐙(s) ⟩⟨𝐗_s,u, 𝔥^(1)⟩⟨𝐗_u,t, [𝔥^(2)]_i ⟩ - ⟨𝔥, 𝐙(u) ⟩⟨𝐗_u,t,[𝔥]_i⟩
       = ∑_𝔥∈ℱ^0_N-1∑_(𝔥)∑_𝔥̃∈ℱ^0_N-11_{𝔥^(2) = 𝔥̃}⟨𝔥, 𝐙(s) ⟩⟨𝐗_s,u, 𝔥^(1)⟩⟨𝐗_u,t, [𝔥̃]_i ⟩ - ⟨𝔥, 𝐙(u) ⟩⟨𝐗_u,t,[𝔥]_i⟩
       = ∑_𝔥̃∈ℱ^0_N-1∑_𝔥∈ℱ^0_N-1∑_(𝔥)1_{𝔥^(2) = 𝔥̃}⟨𝔥, 𝐙(s) ⟩⟨𝐗_s,u, 𝔥^(1)⟩⟨𝐗_u,t, [𝔥̃]_i ⟩ - ⟨𝔥̃, 𝐙(u) ⟩⟨𝐗_u,t,[𝔥̃]_i⟩
       = ∑_𝔥̃∈ℱ^0_N-1⟨𝐗_u,t, [𝔥̃]_i ⟩( ⟨𝐗_s,u⋆𝔥̃, 𝐙(s) ⟩ - ⟨𝔥̃, 𝐙(u) ⟩) 
       = - ∑_𝔥̃∈ℱ^0_N-1⟨𝐗_u,t, [𝔥̃]_i ⟩ R^𝔥̃_s,u.






 Let T > 0, 𝐗 be a α-branched rough path and 𝐙 be controlled by 𝐗. Then,
 
    ∫_s^t Z(r)  d𝐗^i(r) := lim_|𝒫| → 0∑_[u,v] ∈𝒫Z̃_u,v

 exists for every i ∈{1,…, d} and [s,t] ⊆ [0,T] where
 
    Z̃_u,v := ∑_𝔥∈ℱ_N-1^0⟨𝔥, 𝐙(u) ⟩⟨𝐗_u,v,[𝔥]_i⟩

 and 𝒫 denotes a partition of [s,t] with mesh size |𝒫|. Moreover,there exists a constant C depending only on α and T such that
  
    |∫_s^t Z(r)  d𝐗^i(r) - Z̃_s,t| ≤ C |t-s|^(N+1)α∑_𝔥∈ℱ^0_N-1⟨𝐗_·,·, [𝔥]_i ⟩_(|𝔥| + 1)α;[s,t] R^𝔥_(N - |𝔥|)α; [s,t].






 This is a consequence of the sewing lemma <cit.> and Lemma <ref>.




 The above theorem defines a map which sends a controlled path 𝐙 to a path t ↦∫_0^t Z(r)  d𝐗^i(r) ∈^m. In fact, this map can be naturally extended to a map 𝐙↦∫_0^·𝐙(r)  d𝐗^i(r) where t ↦∫_0^t𝐙(r)  d𝐗^i(r) is a controlled path above t ↦∫_0^t Z(r)  d𝐗^i(r). To do this, we have to specify ⟨𝔥, ∫_0^t𝐙(r)  d𝐗^i(r) ⟩ for every dual element 𝔥∈ℱ^*_N-1∪{1}. We set 
 
    ⟨ 1, ∫_0^t𝐙(r)  d𝐗^i(r) ⟩ := ∫_0^t Z(r)  d𝐗^i(r)

 and 
 
    ⟨ [τ_1 ⋯τ_n]_i, ∫_0^t𝐙(r)  d𝐗^i(r) ⟩ := ⟨τ_1 ⋯τ_n, 𝐙(t) ⟩.

 In all other cases, we define
 
    ⟨τ_1 ⋯τ_n, ∫_0^t𝐙(r)  d𝐗^i(r) ⟩ := 0.

 More generally, if 𝐙 = (𝐙^1,…,𝐙^d) and every 𝐙^i is controlled by 𝐗, we define a controlled path t ↦∫_0^t 𝐙(r) · d 𝐗(r) by setting 
 
    ⟨ 1, ∫_0^t 𝐙(r) · d 𝐗(r) ⟩ := ∑_i = 1^d ⟨ 1, ∫_0^t 𝐙^i(r)  d𝐗^i(r) ⟩,

 
    ⟨ [τ_1 ⋯τ_n]_i, ∫_0^t 𝐙(r) · d 𝐗(r) ⟩ := ⟨τ_1 ⋯τ_n, 𝐙^i(t) ⟩

 for a tree [τ_1 ⋯τ_n]_i ∈𝒯_N-1, i ∈{1,…,d} and 

    ⟨τ_1 ⋯τ_n, ∫_0^t𝐙(r) · d𝐗(r) ⟩ := 0

 otherwise.
 
 
  The map
  
    I 𝒬_𝐗(^m)^d    →𝒬_𝐗(^m) 
    𝐙   ↦∫_0^·𝐙(r) · d𝐗(r)

  is well-defined and continuous.

 
 
 
  <cit.>.
 
 
 The next lemma shows that controlled paths composed with sufficiently smooth functions are again controlled.
 
 
  Let ϕ^m →^m be sufficiently smooth such that all derivatives below exist. For 𝐙∈𝒬_𝐗(^m), we define ⟨ 1,   ϕ(𝐙(t)) ⟩ := ϕ(Z(t)) and
  
    ⟨𝔥, ϕ(𝐙(t)) ⟩ := ∑_n = 1^N-1∑_𝔥_1 ⋯𝔥_n = 𝔥1/n!ϕ^(n)(Z(t)) ( ⟨𝔥_1, 𝐙(t) ⟩, …, ⟨𝔥_n, 𝐙(t) ⟩)

  for 𝔥∈ℱ_N-1^*. Then, ϕ(𝐙) ∈𝒬_𝐗(^m). 

 
 
 
  <cit.>.
 
 
 We are now able to say what a solution to (<ref>) actually means.
 
 
  A path y  [0,T] →^m is a solution to (<ref>) if y(t_0) = y_0 and if there exists a controlled path 𝐘∈𝒬_𝐗(^m) above y such that
  
    𝐘(t) - 𝐘(s) = ∫_s^t f(𝐘(r)) · d 𝐗(r)

  holds for every s ≤ t, s,t ∈ [t_0,T], where we set f(𝐘(r)) := (f_1(𝐘(r)), … f_d(𝐘(r))).

  
  
  
 
 
 Proving that (<ref>) admits a (unique) solution is done by a standard fixed-point argument <cit.>. If f is of class Lip^γ-1 for some γ > 1/α, a local solution to (<ref>) exists. If Lip^γ-1_b, the solution exists on every time interval. For f being of class Lip^γ resp. Lip^γ_b, the local resp. global solution is unique. Moreover, in the second case, the solution map is continuous.
 
 Recall the definition of the elementary differential F(τ) for τ∈𝒯^0 given in Definition <ref>. 
 








 Let 𝐘 [0,T] →ℋ_N-1 with y(t) = ⟨ 1, 𝐘(t) ⟩ be a solution to (<ref>). Then, the coefficients of 𝐘 are given by 
  
    ⟨τ, 𝐘(t) ⟩ = 1/σ(τ) F(τ)(y(t))

  for τ∈𝒯^*_N-1∪{1} and ⟨τ_1 ⋯τ_n, 𝐘(t) ⟩ = 0 for τ_1 ⋯τ_n ∈ℱ^*_N-1∖𝒯^*_N-1.




 Being a solution to (<ref>) means that
 
    y(t) - y(s) = ⟨ 1, ∫_s^t f(𝐘(r)) · d𝐗(r) ⟩

 and
 
    ⟨ [τ_1 ⋯τ_n]_i , 𝐘(t) ⟩ = ⟨ [τ_1 ⋯τ_n]_i , ∫_0^t f(𝐘(r)) · d𝐗(r) ⟩

 for all [τ_1 ⋯τ_n]_i ∈𝒯^*_N-1, i ∈{1,…,d}, and
 
    ⟨τ_1 ⋯τ_n , 𝐘(t) ⟩ = 0

 for all τ_1 ⋯τ_n ∈ℱ^*_N-1∖𝒯^*_N-1. We prove the assertion for all trees τ∈𝒯^*_N-1 by induction on the height of τ. For τ = 1, the claim follows by definition. Now let τ = [τ_1 ⋯τ_n]_i = [(τ_1)^r_1⋯ (τ_q)^r_q]_i ∈𝒯^*_N-1 for some i ∈{1,…,d} and pairwise distinct trees τ_1,…,τ_q. From (<ref>), we have
  
    ⟨τ , 𝐘(t) ⟩   = ⟨ [τ_1 ⋯τ_n]_i , ∫_0^t f(𝐘(r)) · d𝐗(r) ⟩
       = ⟨τ_1 ⋯τ_n, f_i(𝐘(t)) ⟩
       =  ∑_λ_1, ⋯, λ_n ∈𝒯^*_N-2
    λ_1 ⋯λ_n = τ_1 ⋯τ_n1/n! f_i^(n)(y(t)) ( ⟨λ_1, 𝐘(t) ⟩, …,  ⟨λ_n, 𝐘(t) ⟩) 
       = 1/r_1! ⋯ r_q!1/n!∑_σ∈sym(n) f_i^(n)(y(t)) ( ⟨τ_σ(1), 𝐘(t) ⟩, …,  ⟨τ_σ(n), 𝐘(t) ⟩) 
       = 1/r_1! ⋯ r_q! f_i^(n)(y(t)) ( ⟨τ_1, 𝐘(t) ⟩, …,  ⟨τ_n, 𝐘(t) ⟩) 
       = 1/r_1! ⋯ r_q!1/σ(τ_1) ⋯σ(τ_n) f_i^(n)(y(t)) ( F(τ_1)(y(t)), …,  F(τ_n)(y(t))) 
       = 1/σ(τ) F(τ)(y(t))

  by induction hypothesis.














































 
 
 Let h>0. 
 Then, (<ref>) has the expansion
  
    y(t_0 + h) 
       =y_0 + ∑_τ∈𝒯_p1/σ(τ) F(τ)(y_0) ⟨𝐗_t_0, t_0 + h,τ⟩ + 𝒪(h^(p+1)α)

  for every p ≥⌊ 1/α⌋.

 
 
 
  Let s < t. Note that
  
    y(t) - y(s) = ∑_i = 1^d ∫_s^t f_i(y(r))  d𝐗^i(r).

  For i ∈{1,…, d}, set 
  
    Z̃_s,t^i := ∑_𝔥∈ℱ_p-1^0⟨𝔥, f_i(𝐘(s)) ⟩⟨𝐗_s,t,[𝔥]_i ⟩  = ∑_𝔥∈ℱ_p-1^01/σ([𝔥]_i) F([𝔥]_i)(y(s)) ⟨𝐗_s,t,[𝔥]_i ⟩

  where we use Lemma <ref> for the equality. We therefore obtain
  
    y(t) - y(s) =  ∑_i = 1^d Z̃_s,t^i + R_s,t = ∑_τ∈𝒯_p1/σ(τ) F(τ)(y(s)) ⟨𝐗_s,t,τ⟩ + R_s,t

  where 
  
    R_s,t = ∑_i = 1^d ∫_s^t f_i(y(r))  d𝐗^i(r) - Z̃_s,t^i.

  Using Lemma <ref> and the sewing Lemma <cit.>, we conclude that R_s,t is of order 𝒪((t-s)^(p+1)α).




 


































































































We introduce the local error by le(t_0, y_0; h) : = y(t_0 + h) - y_1 which is the error of one step with the iterative scheme (<ref>) starting in the exact value y_0.
Comparing Theorems <ref> and <ref> and exploiting Lemma <ref>, we see that 




the local error is 
    |le(t_0, y_0; h)| = 𝒪(h^(p+1)α)

for sufficiently small h>0 if and only if 

    ⟨𝐗_t_0, t_0+h,τ⟩  = a(τ)(h)  ∀τ∈𝒯 with  |τ|≤ p.







§ SIMPLIFIED RUNGE KUTTA METHODS


 


In the following, 𝐗 denotes an α-branched rough path for some α∈ (0,1]. Assume that there is a smooth path X^h such that its natural lift 𝐗^h to a branched rough path, cf. Example <ref>, approximates 𝐗, i.e., ϱ_α(𝐗^h, 𝐗)→ 0 for h→ 0. This implies that 𝐗 is a geometric rough path <cit.> and that ϱ^g_α(𝐗^h, 𝐗)→ 0 where ϱ^g_α denotes the inhomogeneous rough paths metric for geometric rough paths <cit.>. We introduce the equation associated to the smooth driver by 

    dy^h(t) = f(y^h(t))  d𝐗^h(t),   y^h(t_0) = y_0.

This equation can be solved by considering 𝐗 as a branched rough path or a geometric rough path, the solution is the same in both cases. It also coincides with the solution to the corresponding Riemann-Stieltjes equation which is well-defined since X^h is smooth by assumption. Since the solution to (<ref>) is a locally Lipschitz continuous function of 𝐗, cf. <cit.> in the case of geometric rough paths or <cit.> for branched rough paths, i.e., 

    sup_t∈ [t_0, T]|y(t)-y^h(t)| ≲ϱ^g_α(𝐗^h, 𝐗),

we find that y^h is close to y for sufficiently small h. In this section, we restrict ourselves to branched rough paths 𝐗 that are the limit of a lifted piece-wise linear approximation X^h of X. This, e.g., includes semi-martingales, fractional Brownian motions with Hurst index H>1/4 and other Gaussian processes <cit.>. This piece-wise linear approximation to X on some grid t_0< t_1 <… <t_N =T is constructed as follows:

    X^h(t) = X(t_k) + t-t_k/h_k[X(t_k+1) - X(t_k) ],   t∈(t_k, t_k+1],

where h_k=t_k+1-t_k and k= 0, 1, …, N - 1. We assume that this piece-wise linear 
approximation converges with rate r_0>0, meaning that 

    ϱ^g_α(𝐗^h, 𝐗) = 𝒪(h^r_0)

for sufficiently small h, where h=max_k=0,…, N-1| t_k+1-t_k|.


 In <cit.>, the almost sure convergence rate of 𝐗^h to 𝐗 is calculated for the natural lift 𝐗 (in the sense of <cit.>) of a large class of Gaussian processes X in the metric ϱ_α^g. In particular, for the lift of a fractional Brownian motion with Hurst parameter H ∈ (1/4,1,2], one can show that the rate r_0 is arbitrarily close to 2H - 1/2 provided one chooses α sufficiently small.




Below, we analyze the order of the local error of some simplified Runge-Kutta scheme if the underlying driver is 𝐗^h, considered as α-branched rough path. This scheme is obtained by setting 
    Z^(k)_ij = a_ij X^k_t_n, t_n+1  and   z^(k)_i =b_i X^k_t_n, t_n+1

in (<ref>), where X^k_t_n, t_n+1 denotes the increment of the kth component of X on [t_n, t_n+1].
Method (<ref>) then becomes 

    Y_i^h    = y_n^h + ∑_k = 1^d ∑_j=1^s a_ij f_k(Y^h_j)  X^k_t_n, t_n+1 = y_n^h + ∑_j=1^s a_ij f(Y^h_j)  X_t_n, t_n+1
    
          y_n+1^h    = y_n^h + ∑_k = 1^d ∑_i=1^s b_i f_k(Y^h_i)  X^k_t_n, t_n+1 = y_n^h + ∑_i=1^s b_i f(Y^h_i)  X_t_n, t_n+1,

where 𝒜 = (a_ij) is a deterministic matrix and b=(b_i) a deterministic vector. This Runge-Kutta method based on the increments of X was considered in <cit.> in the context 
of implicit schemes for equations driven by a certain class of Gaussian processes. We aim to find general conditions on the coefficients b and 𝒜 that guarantee the desired order of the
local error when approximating (<ref>). We begin with a result characterizing the branched rough path if the underlying path is given by (<ref>). 
 

Let τ∈𝒯 be a tree of order p, i.e, |τ| = p. Then, for the branched rough path associated to the piece-wise linear approximation in (<ref>), we have 
    ⟨𝐗^h_t_k, t_k+1, τ⟩ = 1/γ(τ) X^i_1_t_k, t_k+1   X^i_2_t_k, t_k+1… X^i_p_t_k, t_k+1

where the i_ℓ∈{1, 2, …, d} are the labels of the tree τ and X^i_t_k, t_k+1 is the increment of the ith component of X on [t_k, t_k+1].
                                                                                  

We prove by induction on the height of τ that 
    ⟨𝐗^h_t_k, t, τ⟩ = 1/γ(τ)(t-t_k/h_k)^p X^i_1_t_k, t_k+1   X^i_2_t_k, t_k+1… X^i_p_t_k, t_k+1

for t∈ (t_k, t_k+1]. Setting t=t_k+1 then yields the claim. The identity is true for τ = 1 and τ = ∙_i_1.
Let us assume that (<ref>) is true for all sub-trees of τ = [τ_1, …, τ_n]_i_p. Then, according to 
Example <ref>, we have 
    ⟨𝐗^h_t_k, t, τ⟩  = ∫_t_k^t ⟨𝐗^h_t_k,u, τ_1 ⟩⋯⟨𝐗^h_t_k,u, τ_n ⟩  dX^h, i_p(u) 
    = ∫_t_k^t ⟨𝐗^h_t_k,u, τ_1 ⟩⋯⟨𝐗^h_t_k,u, τ_n ⟩X^i_p_t_k, t_k+1/h_k  du
       = ∫_t_k^t 1/γ(τ_1)(u-t_k/h_k)^p_1 X^ĩ_1_t_k, t_k+1⋯ X^ĩ_p_1_t_k, t_k+1⋯1/γ(τ_n)(u-t_k/h_k)^p_n X^i̅_1_t_k, t_k+1⋯ X^i̅_p_n_t_k, t_k+1X^i_p_t_k, t_k+1/h_k  du,

where p_i:= |τ_i| (i=1, …, n). Since ∑_i=1^n p_i= p-1 and Π_i=1^n 1/γ(τ_i) = p/γ(τ), we obtain

    ⟨𝐗^h_t_k, t, τ⟩  = 
    ∫_t_k^t p/γ(τ)(u-t_k/h_k)^p-11/h_k X^i_1_t_k, t_k+1⋯ X^i_p_t_k, t_k+1  du 
    = 1/γ(τ)(t-t_k/h_k)^p X^i_1_t_k, t_k+1⋯ X^i_p_t_k, t_k+1

which concludes the proof of this proposition.

The local error of the simplified Runge-Kutta scheme applied to (<ref>) is defined as  le^h(t_0, y_0; h) : = y^h(t_0 + h) - y^h_1. We can now rewrite (<ref>) and (<ref>)
using Proposition <ref>. Moreover, within the series representation given in Theorem <ref>, a(τ)(h) is replaced by a^h(τ)(h) if the simplifying ansatz
(<ref>) is used. Now, the order of the local error of (<ref>) is
 
    |le^h(t_0, y_0; h)| = 𝒪(h^(p+1)α)

for sufficiently small h>0 if and only if 

    1/γ(τ) X^i_1_t_0, t_0+h   X^i_2_t_0, t_0+h… X^i_|τ|_t_0, t_0+h   = a^h(τ)(h)  ∀τ∈𝒯 with  |τ|≤ p,

where α is the Hölder regularity of X. Based on (<ref>), we aim to find proper choices of 𝒜 and b in (<ref>) that provide
the desired local rate in (<ref>). In order to simplify the notation in the result below, we introduce c_i:= ∑_j=1^s a_ij. We now formulated conditions for the order of the local 
error associated to the simplified Runge-Kutta scheme.

The simplified Runge-Kutta method (<ref>) approximating (<ref>) has a local error of order (p+1)α, i.e., 

    |le^h(t_0, y_0; h)| = 𝒪(h^(p+1)α)

if and only if the following conditions are satisfied for all ℓ=1, …, p:



Let i_1, i_2, i_3∈{1, …, d}. We start analyzing (<ref>) for all trees of order one, i.e.,  τ= ∙_i_1. Using the definition of a^h(τ)(h), i.e., we plug in 
(<ref>) in the definition of a(τ)(h), we obtain 
    a^h(∙_i_1)(h) := ⟨ z^(i_1), Φ(1)(h) ⟩ = ∑_i=1^s z^(i_1)_i = ∑_i=1^s b_i X^i_1_t_0, t_0+h.

Inserting this into (<ref>), we find 
    1/γ(∙_i_1) X^i_1_t_0, t_0+h =   ∑_i=1^s b_i X^i_1_t_0, t_0+h

which is equivalent to ∑_i=1^s b_i=1. We continue with the trees of order two. These are of the form τ = [∙_i_2]_i_1. Again, we determine a^h(τ)(h) which is 

    a^h(τ)(h) := ⟨ z^(i_1), Φ(∙_i_2)(h) ⟩ =  ⟨ z^(i_1), Z^(i_2)Φ(1)(h)⟩ =  ⟨ b, 𝒜1_s⟩ X^i_1_t_0, t_0+hX^i_2_t_0, t_0+h,

using the representations in (<ref>). With this expression for a^h(τ)(h), (<ref>) becomes 
    1/2 X^i_1_t_0, t_0+h   X^i_2_t_0, t_0+h   =  ⟨ b, 𝒜1_s⟩ X^i_1_t_0, t_0+hX^i_2_t_0, t_0+h

exploiting that γ(τ) = 2. This is equivalent to ∑_i=1^s b_i c_i = 1/2. We conclude this proof by considering the order three trees. We start with trees of the form
τ = [[∙_i_3]_i_2]_i_1. Then, a^h(τ)(h) is

    a^h(τ)(h)    = ⟨ z^(i_1), Φ([∙_i_3]_i_2)(h) ⟩ = ⟨ z^(i_1), Z^(i_2)Φ(∙_i_3)(h) ⟩ = ⟨ z^(i_1), Z^(i_2) Z^(i_3)1_s⟩
       =  ⟨ b, 𝒜 (𝒜1_s)⟩ X^i_1_t_0, t_0+hX^i_2_t_0, t_0+h X^i_3_t_0, t_0+h.

Moreover, we see that γ(τ) = 6. Using the above, (<ref>) for τ = [[∙_i_3]_i_2]_i_1 is equivalent to 

    1/6 = ⟨ b, 𝒜 (𝒜1_s)⟩ = ∑_i=1^s∑_j=1^s b_i a_ij c_j.

Now, the only type of tree left is the branched tree τ = [∙_i_2, ∙_i_3]_i_1. The corresponding a^h(τ)(h) is 

    a^h(τ)(h)    = ⟨ z^(i_1), Φ(∙_i_2)(h) Φ(∙_i_3)(h)⟩  = ⟨ z^(i_1), Z^(i_2)1_s Z^(i_3)1_s⟩
       =  ⟨ b, (𝒜1_s) (𝒜1_s)⟩ X^i_1_t_0, t_0+hX^i_2_t_0, t_0+h X^i_3_t_0, t_0+h.

Notice that the product of two vectors is meant component-wise. For this tree, (<ref>) therefore is equivalent to 
    1/3 = 1/γ(τ) = ⟨ b, (𝒜1_s) (𝒜1_s)⟩ = ∑_i=1^s b_i c_i^2

which finally proves the claim.
                                                                                                                                   
                                                                                                                                   
                                                                                                                                   
In fact, we can easily find algebraic conditions for any ℓ>3 in Table <ref> by considering trees τ∈𝒯 with |τ|>3 in (<ref>). This means that we can achieve
a local rate of (p+1)α for the simplified Runge-Kutta method for arbitrary p∈ℕ.

The conditions given in Table <ref> are nothing but the consistency conditions known for the case of f≡ 0 in (<ref>), see, e.g., <cit.>.
The consistency order is the order in the step size h 
of the expression le(t_0, x_0; h)/h. If all conditions in Table <ref> are fulfilled, then one has a scheme of consistency order 3 assuming f≡ 0 in (<ref>). Such 3rd order schemes 
are well-studied in the ordinary differential equation scenario. Below, we provide just a few examples that satisfy these conditions.

We introduce the general Butcher-scheme:
    ℬ𝒮 := [  c_1 a_11 a_12    … a_s1;  c_2 a_21 a_22    … a_2s;    ⋮    ⋮    ⋮    ⋱    ⋮;  c_s a_s1 a_s2    ⋯ a_ss;       b_1  b_2    ⋯  b_s ].


 
  (i) An explicit Runge-Kutta scheme satisfying all conditions in Table <ref> is Heun's third-order method:
 
    ℬ𝒮 = [   0   0   0   0; 1/3 1/3   0   0; 2/3   0 2/3   0;     1/4   0 3/4 ].

Hence, the iterative scheme (<ref>) is 
    y_n+1^h = y_n^h + 1/4[f(y_n^h ) + 3f(Y^h_3)]  X_t_n, t_n+1

where Y_3^h is given by 
    Y_3^h     = y_n^h +  2/3 f(Y^h_2)  X_t_n, t_n+1  with       Y_2^h  = y_n^h +  1/3 f(y_n^h)  X_t_n, t_n+1.

 
  (ii) Another explicit method fulfilling the conditions in Table <ref> is Kutta's third order scheme:
  
    ℬ𝒮 = [   0   0   0   0; 1/2 1/2   0   0;   1  -1   2   0;     1/6 2/3 1/6 ].

 Consequently, the simplified Runge-Kutta method is 
 
    y_n+1^h    = y_n^h + 1/6[f(y_n^h) + 4 f(Y^h_2) +f(Y^h_3) ] X_t_n, t_n+1,

 where Y_2^h and Y_3^h are computed by 
    Y_2^h = y_n^h + 1/2 f(y_n^h)  X_t_n, t_n+1  and   Y_3^h = y_n^h + [- f(y_n^h) + 2 f(Y^h_2) ]  X_t_n, t_n+1.

 
 
Notice that there is much more schemes satisfying the above conditions, e.g., <cit.> provide two implicit Runge-Kutta methods (for stochastic differential equations 
driven by a certain class of Gaussian processes) that satisfy the requirements in Table <ref>.
 
 


§ GLOBAL RATES



























































 §.§ Global rate of the full Runge-Kutta scheme



Let y_t^s, y_0 denote the solution to (<ref>) at time t starting in y_0 at s, i.e., y_s^s, y_0= y_0.
Let a numerical scheme be given as the following one step method: 
  
    y_n+1 = y_n + Φ(y_n, 𝐗_t_n,t_n+1),

 where t_0 < t_1 < … < t_N = T is a partition of [t_0,T].
Below, we analyze the order of convergence of the numerical method (<ref>). The next proposition shows that we loose one order from 
the local to the global error.

 If there is a constant C_1>0 such that 
    | y_t^s, y_0 - y_0 - Φ(y_0,𝐗_s,t)|≤ C_1 | t-s|^1+r

for | t-s| being sufficiently small and if

    | y_t^s, y_0 - y_t^s, ỹ_0|≤ C_2 | y_0-ỹ_0|

for some constant C_2>0, where y_0,ỹ_0∈ℝ^m and 0≤ s≤ t≤ T. Then, there is some C>0 such that

    max_k=0, …, N| y(t_k) - y_k |≤ C h^r

for r>0, where h=max_k=0,…, N-1| t_k+1-t_k|.


We write the global error as follows 
  
    y_k - y(t_k) =∑_j=0^k-1(y_t_k^t_j+1, y_j+1- y_t_k^t_j, y_j)
                                     
using that y_t_k^t_0, y_0=y(t_k) and y_t_k^t_k, y_k= y_k. We combine 
    y_t_k^t_j, y_j= y_t_k^t_j+1, y_t_j+1^t_j, y_j  and   y_j+1  = y_j + Φ(y_j,𝐗_t_j,t_j+1)

with (<ref>) which yields 
    | y(t_k) - y_k|   =∑_j=0^k-1| y_t_k^t_j+1, y_j + Φ(y_j,𝐗_t_j,t_j+1) - y_t_k^t_j+1, y_t_j+1^t_j, y_j|≤ C_2 ∑_j=0^k-1| y_j + Φ(y_j,𝐗_t_j,t_j+1) -  y_t_j+1^t_j, y_j|
       ≤ C_1 C_2 ∑_j=0^k-1| t_j+1-t_j |^r+1≤ C_1 C_2 h^r ∑_j=0^k-1(t_j+1-t_j)≤ C_1 C_2 (T-t_0) h^r
 
exploiting assumptions (<ref>) and (<ref>). This concludes the proof of this proposition.




 §.§ Global rate of the simplified Runge-Kutta scheme


In this section, we study a particular case of 𝐗 being an α-Hölder geometric rough path, 0 < α≤ 1, that can be approximated by the lift of its piece-wise linear approximated underlying path X. For such driver, the simplified Runge-Kutta method (<ref>) converges. Its order is shown in the following theorem.


Let 𝐗 be an α-Hölder geometric rough path in (<ref>), 0 < α≤ 1,
and let its piece-wise linear approximation X^h be given by (<ref>). We assume that  the Wong-Zakai approximation converges with rate r_0>0, meaning that

    sup_t∈ [t_0, T]|y(t)-y^h(t)| =  𝒪(h^r_0)

for sufficiently small h, where h=max_k=0,…, N-1| t_k+1-t_k| is the maximal step size of the underlying grid, y and y^h are the solutions to (<ref>) and (<ref>), respectively. 
If all conditions in Table <ref> are satisfied and the right hand side f is of class Lip^γ_b for some γ > 1/α,

then the simplified Runge-Kutta method (<ref> converges 
with rate η=min{r_0, 4α-1} to the solution of (<ref>), i.e., 
there is a constant C>0 such that 
    max_k=0, …, N| y(t_k) - y_k^h |≤ C h^η

for sufficiently small h.





It holds that  
    | y(t_k) - y_k^h |≤sup_t∈ [t_0, T]|y(t)-y^h(t)| + max_k=0, …, N| y^h(t_k) - y_k^h |.

Theorem <ref> gives us a rate of 4α for the local error of simplified Runge-Kutta method. Proposition <ref> now provides that 

    max_k=0, …, N| y^h(t_k) - y_k^h | = 𝒪(h^4 α - 1)

if assumption (<ref>) holds true. Let y^h, s, y_0_t denote the solution to (<ref>) with initial time s and initial state y_0.
Since 𝐗 is α Hölder and since 𝐗^h is convergent and hence bounded in h, there is a constant K>0 independent of h such that

    | y_t^h, s, y_0 - y_t^h, s, ỹ_0|≤ K ^K| t-s|^α| y_0-ỹ_0|,

 cf. <cit.>. This implies (<ref>) and concludes the proof.




 
  (i) From (<ref>), a sufficient condition for (<ref>) is that ϱ_α'^g(𝐗^h,𝐗) = 𝒪(h^r_0) for some 0 < α' ≤α in which case one has to assume f ∈Lip^γ_b for some γ > 1/α'.
 
  (ii) Theorem <ref> is formulated for any roughness parameter α > 0. For a fractional Brownian motion with Hurst parameter H ∈ (1/4,1), 
 it gives an optimal rate of convergence in the case when H ∈ (1/4,1/2]. Indeed, from <cit.>, we know that r_0 can be chosen arbitrarily close to 2H - 1/2. Since 2H - 1/2 < 4H -1,
 the convergence rate of the simplified Runge-Kutta scheme is arbitrarily close to 2H - 1/2. This rate is the same as for the simplified Milstein scheme introduced in <cit.>, cf. <cit.>,
 which is believed to be optimal due to the results obtained in <cit.>. 



  






















































§ NUMERICAL EXPERIMENTS


We illustrate the rate of convergence of a scheme presented in Example <ref>. In particular, we apply Heun's method to (<ref>) with f_1(y) = cos(y), f_2(y)=sin(y) and y(t)∈ℝ.
Then, we have 
    dy(t) = cos(y(t))  d𝐗^1(t)+sin(y(t))  d𝐗^2(t),   y(0) = 1,    t ∈ [0, T].

We assume that X is a path of a two-dimensional fractional Brownian motion with independent components and Hurst index H=0.4. Moreover, 𝐗 denotes its geometric lift and T=0.25. This example was considered by Deya, Neuenkirch and Tindel in <cit.> in the context of rates of convergence for a Milstein scheme. We use equidistant grid points, i.e., h=T/N. We determine the maximal discretization error 

    ℰ(h):= max_k=0, …, N| y(t_k) - y_k^h |

for different h, where y_k^h is the kth iterate of the simplified scheme in (<ref>) with coefficients defined in Example <ref> (i).
There is no explicit representation for the solution to 
(<ref>). Therefore, we create a reference solution based on the numerical method for a very small step size. In Figure <ref>, the red circles show ℰ(h) in dependence of h 
for three paths X of a fractional Brownian motion. The theoretical rate of convergence for the Heun method is 2H -0.5 = 0.3. The slopes of the regression lines in blue confirm this rate up to an acceptable deviation. 
This deviation can be explained by the fluctuations that can be expected due to the underlying small rate of convergence.
                                                          

































 §.§ Acknowledgements



SR is supported by the MATH+ project AA4-2 Optimal control in energy markets using rough analysis and deep networks.
Work on this paper was started while MR and SR were supported by the DFG via Research Unit FOR 2402. Both authors would like to thank Rosa Preiß for
related discussions and for providing us with her Master thesis in which she corrected the proof of <cit.>.





alpha

'
FGGR16

[Abe80]Abe80
Eiichi Abe.
Hopf algebras, volume 74 of Cambridge Tracts in
  Mathematics.
Cambridge University Press, Cambridge-New York, 1980.
Translated from the Japanese by Hisae Kinoshita and Hiroko Tanaka.

[BB00]BB00
Kevin Burrage and Pamela M. Burrage.
Order Conditions of Stochastic Runge–Kutta Methods by B-Series.
SIAM Journal on Numerical Analysis, 38(5):1626–1646, 2000.

[BBR+18]BBRRS18
Christian Bayer, Denis Belomestny, Martin Redmann, Sebastian Riedel, and John
  Schoenmakers.
Solving linear parabolic rough partial differential equations.
arXiv:1803.09488, 2018.

[BFRS16]BFRS16
Christian Bayer, Peter K. Friz, Sebastian Riedel, and John Schoenmakers.
From rough path estimates to multilevel Monte Carlo.
SIAM J. Numer. Anal., 54(3):1449–1483, 2016.

[But87]butcher
John C. Butcher.
The Numerical Analysis of Ordinary Differential Equations:
  Runge-Kutta and General Linear Methods.
Wiley-Interscience, USA, 1987.

[CK98]CK98
Alain Connes and Dirk Kreimer.
Hopf algebras, renormalization and noncommutative geometry.
Comm. Math. Phys., 199(1):203–242, 1998.

[Dav07]Dav07
Alexander M. Davie.
Differential equations driven by rough paths: an approach via
  discrete approximation.
Appl. Math. Res. Express. AMRX, (2):Art. ID abm009, 40, 2007.

[DK09]debrabant08bao
Kristian Debrabant and Anne Kværnø.
B-series analysis of stochastic Runge-Kutta methods that use an
  iterative scheme to compute their internal stage values.
SIAM J. Numer. Anal., 47(1):181–203, 2008/09.

[DNT12]DNT12
Aurélien Deya, Andreas Neuenkirch, and Samy Tindel.
A Milstein-type scheme without Lévy area terms for SDEs
  driven by fractional Brownian motion.
Ann. Inst. Henri Poincaré Probab. Stat., 48(2):518–550,
  2012.

[FGGR16]FGGR16
Peter K. Friz, Benjamin Gess, Archil Gulisashvili, and Sebastian Riedel.
The Jain-Monrad criterion for rough paths and applications to
  random Fourier series and non-Markovian Hörmander theory.
Ann. Probab., 44(1):684–738, 2016.

[FH14]FH14
Peter K. Friz and Martin Hairer.
A Course on Rough Paths with an introduction to regularity
  structures, volume XIV of Universitext.
Springer, Berlin, 2014.

[FR14]FR14
Peter K. Friz and Sebastian Riedel.
Convergence rates for the full Gaussian rough paths.
Ann. Inst. Henri Poincaré Probab. Stat., 50(1):154–194,
  2014.

[FV10a]FV10-2
Peter K. Friz and Nicolas B. Victoir.
Differential equations driven by Gaussian signals.
Ann. Inst. Henri Poincaré Probab. Stat., 46(2):369–413,
  2010.

[FV10b]FV10
Peter K. Friz and Nicolas B. Victoir.
Multidimensional stochastic processes as rough paths, volume
  120 of Cambridge Studies in Advanced Mathematics.
Cambridge University Press, Cambridge, 2010.
Theory and applications.

[Gub10]Gub10
Massimiliano Gubinelli.
Ramification of rough paths.
J. Differential Equations, 248(4):693–721, 2010.

[HHW18]HHW18
Jialin Hong, Chuying Huang, and Xu Wang.
Symplectic Runge-Kutta methods for Hamiltonian systems driven
  by Gaussian rough paths.
Appl. Numer. Math., 129:120–136, 2018.

[HK15]HK15
Martin Hairer and David Kelly.
Geometric versus non-geometric rough paths.
Ann. Inst. Henri Poincaré Probab. Stat., 51(1):207–251,
  2015.

[HLW06]HairerLubichWanner
Ernst Hairer, Christian Lubich, and Gerhard Wanner.
Geometric Numerical Integration: Structure-Preserving
  Algorithms for Ordinary Differential Equations; 2nd ed., volume 31.
Springer, 2006.

[HNW10]hairer10sodI
Ernst Hairer, Syvert Paul Nørsett, and Gerhard Wanner.
Solving ordinary differential equations. I: Nonstiff problems.
  2nd revised ed., 3rd corrected printing.
Springer, Berlin, 2010.

[HW10]hairer10sodII
Ernst Hairer and Gerhard Wanner.
Solving ordinary differential equations. II: Stiff and
  differential-algebraic problems. Reprint of the 1996 2nd revised ed.
Springer, Berlin, 2010.

[KP99]kloeden99nso
Peter Eris Kloeden and Eckhard Platen.
Numerical solution of stochastic differential equations,
  volume 21 of Applications of Mathematics.
Springer-Verlag, Berlin, 2nd edition, 1999.

[LCL07]LCL07
Terry J. Lyons, Michael Caruana, and Thierry Lévy.
Differential equations driven by rough paths, volume 1908 of
  Lecture Notes in Mathematics.
Springer, Berlin, 2007.
Lectures from the 34th Summer School on Probability Theory held in
  Saint-Flour, July 6–24, 2004, With an introduction concerning the Summer
  School by Jean Picard.

[LQ02]LQ02
Terry J. Lyons and Zhongmin Qian.
System control and rough paths.
Oxford Mathematical Monographs. Oxford University Press, Oxford,
  2002.
Oxford Science Publications.

[Lyo98]Lyo98
Terry J. Lyons.
Differential equations driven by rough signals.
Rev. Mat. Iberoamericana, 14(2):215–310, 1998.

[Man06]Man06
Dominique Manchon.
Hopf algebras, from basics to applications to renormalization.
arXiv:math/0408405v2, 2006.

[MT04]milstein04snf
Grigori N. Milstein and Michael V. Tretyakov.
Stochastic numerics for mathematical physics.
Scientific Computation. Berlin: Springer. ixx, 594 p., 2004.

[NTU10]NTU10
Andreas Neuenkirch, Samy Tindel, and Jérémie Unterberger.
Discretizing the fractional Lévy area.
Stochastic Process. Appl., 120(2):223–254, 2010.

[Röß10]roessler10rkm
Andreas Rößler.
Runge-Kutta methods for the strong approximation of solutions of
  stochastic differential equations.
SIAM J. Numer. Anal., 48(3):922–952, 2010.

[Swe69]Swe69
Moss E. Sweedler.
Hopf algebras.
Mathematics Lecture Note Series. W. A. Benjamin, Inc., New York,
  1969.






