## Abstract

Random ordinary differential equations (RODEs) are ordinary differential equations (ODEs) with a stochastic process in their vector field. They can be analysed pathwise using deterministic calculus, but since the driving stochastic process is usually only Hölder continuous in time, the vector field is not differentiable in the time variable, so traditional numerical schemes for ODEs do not achieve their usual order of convergence when applied to RODEs. Nevertheless deterministic calculus can still be used to derive higher order numerical schemes for RODEs via integral versions of implicit Taylor-like expansions. The theory is developed systematically here and applied to illustrative examples involving Brownian motion and fractional Brownian motion as the driving processes.

## 1. Introduction

Random ordinary differential equations (RODEs) appear to have a shadow existence compared with stochastic differential equations (SDEs), although they have been used and investigated for an even longer time (e.g. Wonham 1970; Soong 1973; Tiwari & Hobbie 1976; Sobczyk 1991; Arnold 1998; Vom Scheidt *et al.* 2002). They enjoy a major advantage over SDEs in that they do not need a special stochastic calculus, but can be analysed pathwise with the usual methods of deterministic calculus. This is also useful for SDEs, since any SDE can be transformed into a RODE (Sussmann 1977; Imkeller & Schmalfuß 2001), which allows stronger pathwise results (Jentzen *et al.* submitted; see also Kloeden & Neuenkirch 2007) to be obtained for SDEs compared with the *L*_{2} results given by the Ito stochastic calculus.

All the same, the theory of RODEs is not just a straightforward repetition of known results for ordinary differential equations (ODEs). This is most apparent in the numerical treatment of RODEs. Since the solution of a RODE is continuously differentiable but not further differentiable in time, conventional numerical methods for ODEs such as the Runge–Kutta schemes do not achieve their traditional order when applied to RODEs (Kloeden *et al.* 1999; Grüne & Kloeden 2001; Carbonell *et al.* 2005). Nevertheless deterministic calculus can still be used to derive higher order numerical schemes for RODEs via integral versions of implicit Taylor-like expansions. These are used in this paper to systematically derive Taylor schemes of arbitrary higher order for RODEs. For notational simplicity this will be done for scalar RODEs with scalar noise, but everything carries over without serious difficulty to the multidimensional setting, i.e. with vector-valued state or noise, or both (Jentzen 2007).

RODEs are defined in §2, where notation and some preliminary results are also presented. Implicit Taylor-like expansions for the solutions of a RODE are then derived in §3. These implicit Taylor-like expansions are approximated recursively in §4, resulting in a class of RODE–Taylor schemes with respect to given multi-indicial sets. Examples of RODE–Taylor schemes are given and compared in §5 for RODEs involving Brownian motion and fractional Brownian motion with Hurst parameter 3/4 as the driving processes. The discretization errors of the RODE–Taylor schemes are then estimated in §6. Finally, some numerical results are presented in §7.

## 2. Preliminaries and notation

Let (*Ω*, , ) be a probability space with the canonical sample space *Ω*≔*C*^{0}(, ), the space of continuous functions from into such that , where *ω*∈*Ω*, is a stochastic process on (*Ω*, , ) with continuous sample paths. In addition, let *U* be an open subset of and let *f* : ×*U*→ be a continuous function. Then,(2.1)is a scalar RODE, that is, an ODEfor each *ω*∈*Ω*. This is a non-autonomous ODE for which the vector field is usually only continuous, but not differentiable in *t*, no matter how smooth the function *f* is in its variables.

For convenience, it will be assumed that *f* is infinitely, often continuously, differentiable in its variables, i.e. *f*∈*C*^{∞}(×*U*,), though *k*-times continuously differentiable with *k* sufficiently large would suffice. In particular, *f* is then locally Lipschitz in *x*, so the initial value problem(2.2)has a unique solution, which will be assumed to exist on a finite time interval [*t*_{0},*T*] under consideration. Since the solution is continuous, there is an *R*(*ω*)=*R*>0 such that |*x*(*t*)|≤*R* for all *t*∈[*t*_{0}*,T*] where the interval [−*R*,*R*]⊂*U*. One notes that *R* depends on *ω*.

Moreover, the stochastic processes *ζ*_{t} will be assumed to have sample paths which are locally Hölder continuous with the same Hölder exponent, i.e. there exists a *γ*∈(0,1] such that for almost all *ω*∈*Ω* and each *T*>0 there exists a *C*=*C*_{ω,T}>0 such that(2.3)Let *θ* be the supremum of all *γ* with this property. Two cases will be distinguished later, case A in which (2.3) also holds for *θ* itself and case B when it does not. Brownian motion with *θ*=1/2 is an example of case B.

For such sample paths *ω* defineso

### (a) Multi-index notation

Let _{0} denote the non-negative integers. For a multi-index defineIn addition, for a given *γ*∈(0,1], define the weighted magnitude of a multi-index *α* byand for each *K*∈_{+} with *K*≥|*α*|_{γ} let . For *f*∈*C*^{∞}(×*U*,), denotewith and (0,0)!=1. Finally, for brevity, write .

Let *ω* and *R*>0 be fixed, the latter being chosen as an upper bound on the solution of the initial value problem (2.2) corresponding to the sample path *ω* on a fixed interval [*t*_{0},*T*]. Defineand, for brevity, write . Note that the solution of the initial value problem (2.2) is Lipschitz continuous withThe expansion(2.4)will be used later as well.

*Let x and y be real numbers. Then*,*for each p*∈.

▪

## 3. Taylor-like expansions

The solution *x*(*t*) of the initial value problem (2.2) is only once differentiable, so the usual Taylor expansion cannot be continued beyond the linear term. Nevertheless, the special structure of a RODE and the smoothness of *f* in both variables enable one to derive implicit Taylor-like expansions of arbitrary order for the solution.

Fix *k*∈_{+} and writewherefor an arbitrary . Then, the usual Taylor expansion for *f* in both variables gives(3.1)with remainder term(3.2)for some . Substituting this into the integral equation representation of the solution of (2.2),(3.3)givesor, more compactly, as(3.4)whereThe expression (3.4) is implicit in Δ*x*_{s}, so it is not a standard Taylor expansion. Nevertheless, it can be used as the basis for constructing higher order numerical schemes for the RODE (2.1).

## 4. RODE–Taylor schemes

The Taylor-like expansion (3.4) will now be used to derive what will be called RODE–Taylor schemes, a family of explicit one-step schemes for RODEs (2.1), which have the general formwhere *h*_{n}=*t*_{n+1}−*t*_{n}>0 is the step size of the *n*th discretization subinterval.

The simplest case is for *k*=0 as (3.4) reduces towhich leads to the well-known Euler schemeIn order to derive higher order schemes, the Δ*x*_{s} terms inside the integrals must also be approximated, which can be done with a numerical scheme of one order lower than that of the scheme to be derived. The higher order schemes can thus be built up recursively for sets of multi-indices of the formwhere *K*∈_{+} and *θ*∈(0,1] are specified by the noise process in the RODE.

Fix *K*∈_{+} and consider the first step,(4.1)of a numerical approximation at the time instant for a step size *h*∈(0,1] and initial value , where the increments are defined recursively. Specifically, let and define(4.2)where(4.3)with and , i.e. in non-trivial cases, the are evaluated in terms of previously determined with .

This procedure is repeated for each time step to give the _{K}-*RODE–Taylor scheme* (the *K*-RODE–Taylor scheme)(4.4)with discretization subintervals [*t*_{n},*t*_{n+1}] and step size *h*=*t*_{n+1}−*t*_{n}>0. (Note that it is not essential to use the -RODE–Taylor scheme in the recursion here; an arbitrary numerical scheme of order could equally well be used, but would then give another class of schemes.)

## 5. Examples

Some explicit examples of RODE–Taylor schemes will now be presented for scalar RODEs, all but one with scalar noise processes. Two representative noise processes will be considered, Brownian motion and fractional Brownian motion with a Hurst coefficient *H*=3/4.

### (a) RODEs with Brownian motion

The following examples have the Brownian motion or Wiener process as the driving process, which falls into case B with *θ*=1/2. The _{K}-RODE–Taylor schemes, which have pathwise global convergence order *K*−*ϵ*, are presented for *K*=0, 0.5, 1.0, 1.5, 2.0 and 2.5. Since *θ*=1/2for these *K*. In the following schemes, the coefficient functions on the r.h.s. are evaluated at (*t*_{n}, *y*_{n}).

The 0-RODE–Taylor scheme corresponding to _{0=}∅ is *y*_{n}≡*y*_{0}. It is an inconsistent scheme.

The 0.5-RODE–Taylor scheme corresponding to _{0.5}={(0,0)} is the classic Euler scheme(5.1)which has order 0.5−*ϵ* here, cf. Grüne & Kloeden (2001).

The 1.0-RODE–Taylor scheme corresponding to _{1.0}={(0,0), (1,0)} is the ‘improved’ Euler scheme,Its order 1−*ϵ* is comparable to that of the Euler scheme for smooth ODEs.

The 1.5-RODE–Taylor scheme corresponding to _{1.5}={(0,0), (1,0), (2,0), (0,1)} isHere and the last term is obtained from with coming from the Euler scheme (5.1).

In the next case _{2.0}={(0,0), (1,0), (2,0), (3,0), (0,1), (1,1)} and , so the terms and corresponding to the 0.5-, 1.0-RODE–Taylor schemes are required in the r.h.s. of the new scheme. The resulting 2.0-RODE–Taylor scheme is then

The 2.5-RODE–Taylor scheme iswith the multi-index set _{2.5}={(0,0), (1,0), (2,0), (3,0), (4,0), (0,1), (1,1), (2,1), (0,2)}.

Finally, consider a two-dimensional Brownian motion . Here the 1.5-RODE–Taylor scheme corresponds to the indicial set _{1.5}={(0,0,0), (1,0,0), (0,1,0), (1,1,0), (2,0,0), (0,2,0), (0,0,1)} and is given by

### (b) RODEs with fractional Brownian motion

Fractional Brownian motion with Hurst coefficient *H*=3/4 also falls into case B with *θ*=3/4. The RODE–Taylor schemes generally contain fewer terms than the schemes of the same order with Brownian motion or attain a higher order when they contain the same terms.

_{K}={(0,0)} for *K*∈(0,3/4]: the RODE–Taylor scheme is the classic Euler scheme in example 5.2, but now the order is (3/4)−*ϵ*.

_{K}={(0,0), (1,0)} for *K*∈(3/4,1]: the RODE–Taylor scheme is the same as that in example 5.3 and also has order 1−*ϵ*.

_{K}={(0,0), (1,0), (0,1)} for *K*∈(1,3/2]: the RODE–Taylor scheme,which has order 1.5−*ϵ*, omits one of the terms in the RODE–Taylor scheme of the same order for Brownian motion given in example 5.4.

_{K}={(0,0), (1,0), (0,1), (2,0)} for *K* is in (3/2,7/4]: the Taylor scheme is the same as that in the Brownian motion case in example 5.4, but now has order (7/4)−*ϵ* instead of order (3/2)−*ϵ*.

## 6. Discretization error

The increment function of a RODE–Taylor schemeis continuous in its variables as well as continuously differentiable, hence locally Lipschitz, in the variable (see Jentzen 2007) withso the RODE–Taylor schemes are consistent and hence convergent for each *K*>0 (e.g. Deuflhard & Bornemann 2002). Moreover, the classic theorem for ODEs on the loss of a power from the local to global discretization errors also holds for the RODE–Taylor schemes. Consequently, it suffices to estimate the local discretization error, which is defined here aswhere is the value of the solution of the RODE at time with initial value and is the first step of the numerical scheme with step size *h* for the same initial value.

It is necessary to distinguish two cases, case A in which the Hölder estimate (2.3) also holds for the supremum *θ* of the admissible exponents itself and case B when it does not.

Define and, for *K*>0,In addition, letand define(6.1)

*The local discretization error for a RODE–Taylor scheme in case A satisfies**for each* 0*≤h≤*1, *where**In case B it satisfies**for ϵ>0 arbitrarily small, where*

Case A will be proved in detail with some remarks provided on the essential differences which arise in case B.

To begin, note that the *T*_{α} terms in the Taylor-like expansion of (3.4) of the RODE solution are easily estimated to give(6.2)where and . Then the estimate(6.3)of the remainder term follows from

Case A is proved by mathematical induction on *K*≥0. Here estimates such as (6.2) and (6.3) hold for all *γ*≤*θ*, including *θ* itself, which will be used below. For *K*=0 the assertion is trivial, becauseFor *K*∈(0,*θ*), the index *k*=*k*_{K}=0 and the multi-index set _{K}={(0,0)}. HenceThis completes the first step of the induction argument.

Suppose now that the assertion is true for *K*<*L* for a fixed *L*=*lθ*, where *l*∈. The aim is to show that the assertion is then true for any *K*∈[*L*,*L*+*θ*). Then *k*=*k*_{K}=1 and

Estimates of the first and third errors were given, respectively, in (6.2) and (6.3). The second is given in the following lemma.

*For α*∈_{K} *with α*_{2}>0, *we have*(6.4)

*with the convention that C*_{r}*≔*1 *for r*<0.

For brevity, omit the variables and in and write

ThenBy lemma 2.1 the difference in the integrand can be estimated by(6.5)where, by the definition of *R*_{K},(6.6)Hence

By the induction assumption, , sowhere *C*_{r}≔1 for *r*<0.

This completes the proof of lemma 6.2. ▪

Returning to the proof of theorem 6.1, the three error estimates (6.2)–(6.4) combine to giveFinally, using the double exponential expansion (2.4), one obtainsThis completes the proof of case A of theorem 6.1. ▪

The main difference in the proof of case B of theorem 6.1 is in the estimate of the second error given in lemma 6.2. The required estimate is given by the following lemma, where the convention that *C*_{r}≔1 for *r*<0 is also used.

*For α*∈_{k} *with α*_{2}>0, *h*∈(0,1] *and ϵ*>0 *arbitrarily small*,(6.7)*where*

As in lemma 6.2, write *E*_{2} for the expression to be estimated. Nowsince Δ*s*∈(0,1] andfor *α*_{1}≤*K*/*θ*≤*k*_{k}+1=*k*+1. Thusso, using the inequalities (6.5) and (6.6) to estimate the difference in the integrand,where by the induction assumption. Hencesince .

This completes the proof of lemma 6.3. ▪

## 7. Numerical results

The above theoretical results are illustrated here with two numerical examples. The first one involves the initial value problem for the scalar ODE (figure 1)(7.1)which has the explicit solution(7.2)where *ω*(*t*) denotes a sample path of a standard Brownian motion.

The second example is the scalar ODE(7.3)which has the explicit solution(7.4)In this situation, the noise process *ω*(*t*) is given bywhere *W*_{t} and *V*_{t} denote two sample paths of independent standard Brownian motions (figure 2).

In the first example, we use the *K*-RODE–Taylor schemes for *K*=1, 1.5, 2, 2.5 and the Heun scheme, i.e. the Runge–Kutta scheme, for ODEs defined bywhich has order 2 for ODEs with a smooth vector field (see Deuflhard & Bornemann 2002). In the second example, the *K*-RODE–Taylor schemes for *K*=0.5, 1.5, 2.5 are calculated. In both figures the error of the approximation of the solution in the interval [0,1] is calculated and shown in the log–log plot.

In the second example the three order lines correspond to the orders 0.5, 1.5 and 2.5, so the three RODE–Taylor schemes attain their predicted order.

In the first example, the three order lines correspond to the orders 1, 2 and 3. One sees that the Heun and the 1-RODE–Taylor schemes have order 1, the 1.5 and the 2-RODE–Taylor schemes have order 2, and finally the 2.5-RODE–Taylor scheme has order 3 here. The Heun, the 1.5-RODE–Taylor and the 2.5-RODE–Taylor schemes actually attain an order 0.5 higher than predicted. This can be explained as follows. The RODE here is the same as the SDEand the _{1.5}-RODE–Taylor scheme corresponds to the strong Taylor scheme of strong order *γ*=2.0 applied to this SDE (see Kloeden & Platen 1992; Jentzen 2007), which is shown in (Jentzen *et al.* (submitted), see also Kloeden & Neuenkirch (2007)) to have pathwise order 2.0−*ϵ* for arbitrarily small *ϵ*>0. For the second example, which does not have an equivalent SDE form, the _{1.5}-RODE–Taylor schemes order attains the predicted order 1.5−*ϵ*. (A similar explanation holds for the other schemes.)

### (a) Approximation of the integrals

The *K*-RODE–Taylor schemes for *K*≥1 require the approximation of integrals such as . How this can be done will be explained in the context of the 2.5-RODE–Taylor scheme. For *i*, *j*∈_{0}, denoteThe 2.5-RODE–Taylor scheme can thus be rewritten asThe integrals *J*_{(i)} for *i*=1, …, 4, and *J*_{(1,0)}, *J*_{(2,0)} can then be approximated by Riemann sumswith *δ*=(*h*/*m*) and *m* sufficient large so that the desired order is obtained, i.e. *m*=⌈(1/*h*)^{4}⌉ in this case. The other integrals in the 2.5-RODE–Taylor schemes are then evaluated from these using formulae such as *J*_{(i,j)}+*J*_{(j,i)}=*J*_{(i)}.*J*_{(j)}, which can be verified by integration by parts.

Note that the number of evaluations of the Brownian motion does not change for different RODE–Taylor schemes for the same specified accuracy. Why then use high-order RODE–Taylor schemes? Larger step sizes can be used to achieve a given desired accuracy, so fewer evaluations of the vector field *f* and its derivatives are required. This is particularly advantageous when the spatial dimension is high and the noise dimension is small (see Grüne & Kloeden 2001).

## Acknowledgments

This work has been partially supported by the DFG project ‘Pathwise numerical analysis of stochastic evolution equations’ and the Ministerio de Educación y Ciencia (Spain) and FEDER (European Community) grant MTM2005-01412.

## Footnotes

- Received May 29, 2007.
- Accepted July 25, 2007.

- © 2007 The Royal Society