Lawson schemes for highly oscillatory stochastic differential equations and conservation of invariants

In this paper, we consider a class of stochastic midpoint and trapezoidal Lawson schemes for the numerical discretization of highly oscillatory stochastic differential equations. These Lawson schemes incorporate both the linear drift and diffusion terms in the exponential operator. We prove that the midpoint Lawson schemes preserve quadratic invariants and discuss this property as well for the trapezoidal Lawson scheme. Numerical experiments demonstrate that the integration error for highly oscillatory problems is smaller than that of some standard methods.

We will also assume that the matrices ∈ R × , = 0, . . . , , are constant, and moreover chosen such that the following commutativity assumption is satisfied: Assumption 1 (Commutativity).
To satisfy this assumption, it is sometimes convenient to split the linear parts of the problem, and let some of it be included in the functions. Applications satisfying Assumption 1 are, e. g., [9] the FitzHugh-Nagumo equation with multiplicative noise, the Lotka-Volterra system and SDEs resulting from spectral spatial discretization of stochastic partial differential equations (SPDEs) with diagonal noise.
Exponential methods for solving such problems have in particular been applied in the SPDE setting, mostly, but not exclusively for problems with additive noise. Cohen [4] proposed an exponential method for stochastic oscillators, Yang et al. [21] suggested one for damped Hamiltonian systems. Recently, Erdoğan and Lord [9] presented a quite general approach for constructing exponential integrators for (2) with multiplicative noise. One of the strategies presented there is an adaptation of a method introduced for ordinary differential equations (ODEs) by Lawson [17] to SDEs of the form (2). The idea is to 1 transform the system by the fundamental solution of the linear part, solve the transformed system by a scheme of preference, and then transform back again. In the ODE literature, this is also referred to as an integrating factor method (Cox and Matthews [6], Maday et al. [18]). This is the procedure which will be applied in this paper, and which is described in detail in Section 2. We are essentially interested in studying highly oscillatory problems, and to show that the Lawson methods can attain good accuracy with larger step sizes than what can be obtained by standard stochastic methods. We will also show that the Lawson methods can, under reasonable assumptions, maintain conservation properties of the underlying scheme. Let us demonstrate the ideas of the paper by an introductory example.

Example 1 (The non-linear Kubo oscillator). As a starting point, consider the linear Kubo oscillator
which models a simple oscillator perturbed by a stochastic term and appears in nuclear magnetic resonance and molecular spectroscopy (Goychuk [10], Kubo [15]). The exact solution of this problem starting at (˜, (˜)) for˜∈ is given by where the fundamental solution˜( ) is a rotation matrix, Clearly, the exact solution ( ) = ( 1 ( ), 2 ( )) of the linear problem is norm-preserving, i. e. satisfies the invariant (4) I ( ( )) = 2 1 ( ) + 2 2 ( ) = constant for all . Cohen [4] proposed an extension to the Kubo oscillator by including a non-linear, skewsymmetric drift term, see also Laurent and Vilmart [16]. We extend this further, with non-linear terms in both the drift and diffusion, in addition to including multidimensional noise: The solutions of these problems all preserve the invariant (4), see Section 3.1 for details. We are interested in studying to which extent a numerical approximation will be able to follow the fast oscillations of the linear parts, as well as how well the invariant (4) is preserved. In this example, the following three methods (described in detail in Section 2) have been applied to the SDE (5): • The standard implicit stochastic midpoint rule ("Midpoint"), which is known to preserve quadratic invariants (Hong et al. [12], Milstein et al. [20]). • The method proposed by Cohen [4] ("TDSL") for highly oscillatory SDEs. This is a drift Lawson scheme based on the trapezoidal rule, but it does not preserve the invariant for the non-linear problem. • A Lawson scheme based on the implicit midpoint rule ("MFSL"). Details of this scheme are given in Section 2, and in Section 3.1 it is proved that this scheme preserves the quadratic invariant I.
In our example we will use = 2 and and the SDE is integrated from 0 to 1, using step size ℎ = 2 −5 . The numerical results are presented in Fig. 1. According to this, the implicit midpoint rule ("Midpoint") preserves the invariant, it is however not able to resolve the fast oscillations. The method proposed by Cohen [4] ("TDSL") resolves the oscillations well, but the solution drifts away from the manifold. The Lawson midpoint rule ("MFSL") both stays on the manifold and resolves the high oscillations. This example is further exploited in Section 4.

Midpoint
The outline of this paper is as follows: In Section 2 we derive the midpoint and trapezoidal stochastic Lawson schemes. In Section 3 it is proved that, under reasonable assumptions, the Lawson transformation preserves linear and quadratic invariants, and consequently, if the underlying scheme preserves such invariants, so will the corresponding Lawson scheme. Section 4 is devoted to numerical experiments.

S L
In this section we will shortly outline how to derive the stochastic Lawson (SL) schemes, with particular emphasis on the implicit midpoint and the trapezoidal Lawson scheme.
In this paper, only Lawson schemes based on the implicit trapezoidal and midpoint rules will be discussed. These two methods, applied to (10), are given by where the stochastic increments are Δ = ( +1 ) − ( ), = 0, . . . , , and we used that ( ) = 0 ∈ R +1 and ( +1 ) = Δ . Under appropriate conditions on the smoothness and boundedness of theˆ(e. g. for mean square convergence of order up to one thatˆ0 + 1 2 =1ˆ ˆandˆ1, . . . ,ˆsatisfy a global Lipschitz condition and are, together with all the associated elementary differentials up to order three, of linear growth), these approximations will be of mean square order = 0.5 ( = 1 for commutative noise) and weak order˜= 1, see e. g. Debrabant and Kvaernø [7], Hong et al. [12], Kloeden and Platen [14], Milstein and Tretyakov [19]. Note that for methods that conserve a bounded manifold, smooth functions can be replaced by smooth bounded functions being zero outside a suitable ball, see e.g. Cohen et al. [5].
Applying the back-transformation results in the two Lawson methods of interest in this paper: Trapezoidal Lawson rule: Midpoint Lawson rule: where Δ = ( +1 ) = =0 Δ . We would like to emphasize that there is some freedom in how to choose the matrices in the SDE (2) and thus the Lawson methods (15) and (16), as for any ∈ R × an integrand˜+˜( ) can be written as see also Debrabant et al. [8]. The matrices have to be chosen such that Assumption 1 is satisfied. Different splittings lead to different schemes. The underlying numerical scheme is restored by setting all the linear parts to 0 ( = 0, 1, . . . , ). In a drift stochastic Lawson scheme (DSL), only the linear drift term is included in the exponential, that is 0 ≠ 0 and = 0 for = 1, . . . , . The exponential Δ = ℎ +1 0 is computed only once, and the DSL schemes can be quite efficient. The scheme suggested by Cohen [4] (denoted by TDSL in the present article) is an example of a DSL scheme based on the trapezoidal rule (15). If ≠ 0 for at least one ≠ 0, the scheme is denoted as a full stochastic Lawson scheme (FSL). In this case the exponentials depend on the stochastic increments and have to be calculated for each step, so the performance of the FSL schemes depends on how efficient this can be done.
For a continuous splitting between the linear and the non-linear part, see Erdoğan and Lord [9].

P
Hong et al. [12]) studied preservation of linear and quadratic invariants for SDEs. In particular, they showed that the stochastic midpoint rule (as a representative for the Gauss methods), preserves quadratic invariants. The aim of this section is to extend these results to the stochastic Lawson (SL) schemes. It is also well known that the standard trapezoidal rule almost preserves quadratic invariants for deterministic differential equations, see e.g. Hairer et al. [11,Example V.4.2]. This property does unfortunately not extend to the stochastic trapezoidal rule in general; we will study under which conditions it does hold. The commutativity condition Assumption 1 is always assumed to be satisfied throughout this section.
In the following, we first investigate when stochastic Lawson schemes preserve quadratic invariants before briefly commenting on the preservation of linear invariants.

Preservation of quadratic invariants.
In this section we study when the SDE (2) for a ∈ R × satisfies the quadratic invariant (17) I ( ( )) := ( ) ( ) = I ( ( 0 )), and the ability of the numerical schemes to preserve it. Without loss of generality we will assume that = . For the rest of this section we will use a combination of the following assumptions:
The necessity of these two assumptions can be seen by applying the chain rule for Stratonovich integrals and using the symmetry of to obtain To preserve I ( ( )) we need both terms to be zero, the first is zero if Assumption 2 holds and the second is zero if Assumption 3 holds. In summary, it holds (see e.g. Hong et al. [12]): (17).
We are now ready to state the main result of this section: Theorem 2. Let Assumptions 1 to 3 hold and let denote the discrete time approximation of the SDE (2) at time point using an SL scheme with an underlying numerical onestep method that preserves quadratic invariants. Then the SL scheme preserves the same quadratic invariant as the SDE (2), i.e. (20) I ( ) = I ( ( )) = I ( ( 0 )).
Requirements for stochastic Runge-Kutta methods to preserve quadratic invariants are given in Hong et al. [12]. Strictly speaking, those results only apply to autonomous systems, it is however straightforward to extend them to the nonautonomous transformed SDE (9). We can then conclude with the following corollary: Under the assumptions of Theorem 2 the midpoint stochastic Lawson scheme preserves the quadratic invariants (17).
It is well known that the standard trapezoidal rule almost preserves quadratic invariants for deterministic differential equations, see e.g. Hairer et al. [11,Example V.4.2]. This property does unfortunately not extend to the stochastic trapezoidal rule in general. But if all the drift and diffusion terms are linear and included in the operator, a similar result can be obtained, which is also a good argument for using full stochastic Lawson schemes instead of the drift ones.

Theorem 3. Let Assumptions 1 to 3 hold, let
≡ 0 for = 1, . . . , in the SDE (2) (all diffusion terms are linear and commute) and let denote its discrete time approximation using the trapezoidal FSL scheme with equidistant step sizes ℎ. Then Proof. The trapezoidal Lawson rule (15) can be considered as composed by two half steps, one with an exponential forward Euler step, and one with an exponential backward Euler step, as follows:ˆ= Doing one more half step with the forward Euler method giveŝ which together with (23) could be considered as one step of a modified version of the midpoint rule starting fromˆ, and we would like to see if it conserves quadratic invariants: By Assumption 2, Δ is orthogonal and commutes with , so where we also used Assumption 1 and that by (23) it holds that The second term of (24) vanishes since ( ) = 0 for all ∈ R by Assumption 3. Due to the last term, the quadratic invariant is not preserved in general. Since Δ 0 = ℎ, if = 0 for = 1, . . . the last sum is 0, andˆ +1ˆ+ 1 =ˆ ˆ, which by (17)

Direct computations show that
and in conclusion Notice that a similar argument can also be used as a direct proof of Corollary 1: One step of the stochastic midpoint Lawson scheme is composed by one half step of the exponential backward Euler scheme followed by one half step of the exponential forward Euler scheme, both using the same stochastic increments Δ /2, so the last sum of (24) is 0.
Notice also that if the quadratic invariant is associated to a bounded manifold, i. e., the eigenvalues of are all positive, then the proof of Theorem 3 implies that under the conditions given there, there exists a constant such that |I ( ) − I ( 0 )| ≤ ℎ 2 for all = 0, 1, . . . , and ℎ sufficiently small.

Preservation of linear invariants.
For completeness we also study when the SDE (2) preserves linear invariants (25)Ĩ ( ( )) = ( ) =Ĩ ( ( 0 )) for an ∈ R and the ability of the numerical schemes to preserve it as well.
We now note that by Assumption 4 we have that Using this property, Lemmas 2 and 3 and Theorem 2 extend directly to linear invariants.

N
In this section four different stochastic Lawson schemes are tested numerically and compared with the standard midpoint rule: The trapezoidal DSL (TDSL), the trapezoidal FSL (TFSL), the midpoint DSL (MDSL) and the midpoint FSL (MFSL). All schemes are of weak order 1, and of strong order 1 for SDEs with commutative noise, otherwise of strong order 0.5. The methods are tested on two oscillatory problems preserving quadratic invariants. The error vs. step size for different degrees of oscillatory behaviour is measured, as well as the ability of the schemes to preserve the invariants. Finally, the methods are tested on a stochastic Fermi-Pasta-Ulam-Tsingou (FPUT) problem. This is highly oscillatory, but does not preserve quadratic invariants exactly. The main findings of the experiments are: • The SL schemes resolve the oscillatory behaviour better than the midpoint rule.
For small linear diffusion terms, there are no significant differences between the DSL and the FSL schemes. For high noise problems, the FSL schemes are superior. • The midpoint SL schemes preserve quadratic invariants.
• The trapezoidal FSL scheme almost preserves quadratic invariants for linear diffusions, but not for non-linear diffusions. • The overhead caused by calculating the exponentials is insignificant, even for the FSL-schemes. In all cases, the underlying non-linear algebraic equations of the implicit methods are solved by Newton's method with tolerance 10 −12 .
4.1. Stochastic rigid body problem. The first example is the classical rigid body problem perturbed with linear skew-symmetric drift and diffusion terms, see e.g. Abdulle et al. [1], Anmarkrud and Kvaernø [2], Cohen [4]. The SDE is given by with, following the above references,  Fig. 2. The 95%-confidence intervals are in all cases found to be less than 12% of the corresponding error value.

TDSL
TFSL MDSL MFSL Midpoint reference line with slope 1 2. The stochastic rigid body problem: Strong error vs. step size for different values of and .
The results of experiment (a) are shown in Fig. 2a. In this experiment, there is a significant contribution to the solution from the noise term. All schemes demonstrate an acceptable performance when = = 1, the case with slow oscillations. For higher values of and the FSL schemes are superior. In the case of = = 10, only the FSL-schemes are able to produce reliable solutions for ℎ > 2 −7 . This behaviour is as expected, since the FSL schemes solve the oscillatory parts of the problem exactly.
The results of experiment (b) are given in Fig. 2b. Here, the oscillatory noise is quite small, and the oscillations are dominated by the contribution from the linear drift in all cases but = 0. For = 0, we observe no significant differences between the methods. The dominating part here is the nonlinear drift term, which is handled similarly by the schemes. For higher values of , all the SL schemes yield similar results, and all are superior to the midpoint rule. In particular the SL schemes are able to resolve the high oscillations for larger step sizes, up to ℎ = 2 −6 even when = 100. The midpoint rule fails to give reliable results for step sizes above approximately ℎ = 2 −9 in this case.
All schemes exhibit, as expected, strong order 1 for sufficiently small step sizes. In Fig. 3 TDSL we depict the computational efforts, measured as wall-clock time per batch of 25 paths vs. the strong error averaged over all batches. Somewhat surprisingly, the FSL schemes are slightly cheaper in terms of computational work per step. All the methods are implicit, thus the efforts of solving the nonlinear equations dominate. The computational overhead required for calculating the exponentials for the FSL methods seems in this example to be outweighted by the fact that there is no diffusion term to calculate in the right hand side of the equation, as it is absorbed in the exponential.
TDSL TFSL reference line with slope 1 reference line with slope 2 4. The stochastic rigid body problem: Weak error of vs. step size for different values of and .
We will next see how well the methods preserve invariants, also over long time integration. In the case of the stochastic rigid body problem (26), Assumptions 2 and 3 are satisfied, and by Lemma 1 the exact solution of (26) preserves the invariant (28) I ( ( )) = ( ) ( ) and thus, when using the given initial values, stays on the unit sphere. The matrices 0 and 1 commute, so the requirements of Corollary 1 are satisfied and the MDSL and MFSL schemes should preserve the invariant. As 1 ( ) = 0 the requirements for Theorem 3 are satisfied, and the TFSL scheme nearly preserves the invariant, and the weak order with respect to is 2. This is clearly demonstrated in Fig. 4, the weak order is one for the TDSL scheme, and two for the TFSL scheme. For the small noise problem, the error from the drift term dominates, and the order, in particular for larger values of ℎ is then close to 2. The midpoint schemes preserve the invariant exact, and are not included in the figure.
Lastly, as 0 ( ) ≠ 0, Theorem 3.1 in Cohen [4] does not apply, and we expect the TDSL scheme to drift off. To confirm this, two extreme cases from the experiments above have been chosen, that is = = 10, and = 100 and = 0.3. The SDE is solved on the interval [0, 100], using the step size ℎ = 25/2 10 . Figure 5 shows the value of the invariant I ( ). For visibility, only each 128th step is plotted. In Fig. 6, we zoom in onto the interval [60, 70], for better visibility neglecting the TFSL method and for the midpoint methods plotting only every 6th value. Notice that the drift of the TDSL schemes is severe in the experiments with much noise, while in the case of small noise, it is reasonably small. To this aim, the stochastic rigid problem is solved by the TFSL scheme over the interval [0, 100] for different step sizes, with 1000 independent simulations. Fig. 7 shows the maximum average deviation over the interval, clearly demonstrating the second order of the deviation predicted by Theorem 3. The parameters used in this experiment are = 10 and = 0.3.
The stochastic rigid body problem with = 10 and = 0.3: Maximum average deviation on the interval [0,100] vs. step size when solved by the TFSL method. denotes the numerically determined order of convergence.

Non-linear Kubo oscillator.
We next study the methods applied to a problem with multiple noise terms. To this aim, the non-linear Kubo oscillator, already introduced in Section 1, is chosen. The SDE is In the current experiments, = 2 and the nonlinear terms are thus there are two independent diffusion terms, one linear and one non-linear. These are scaled to ensure that the highly oscillatory parts come from the linear terms. As in the rigid body case, and vary to investigate how the schemes respond to highly oscillatory parts in the drift and / or diffusion terms. The experiment setup is as for the stochastic rigid problem:  Fig. 8. The 95%-confidence intervals have been calculated and span in all cases less than 15% of the corresponding error values.
The case in which = is presented in Fig. 8a. Again, there are no significant differences between the methods for = = 1 and the strong order is 0.5, as expected.
For higher values of and , the advantage of the FSL methods for larger step sizes is evident. In the cases of the midpoint and the DSL schemes, the error is dominated by the linear stochastic diffusion term, thus causing the first order behaviour. This term is incorporated in the exponentials of the FSL schemes, and will there thus not contribute to the errors. For problems with small linear noise contributions in Fig. 8, all the SL schemes outperform the midpoint scheme for the oscillatory cases, when ∈ {10, 50}. For step sizes larger than ℎ = 2 −8 , the Newton solver frequently failed to solve the underlying nonlinear algebraic equations, in particular for the TDSL scheme. The computational efforts depicted in Fig. 9 again demonstrate that there are only small differences in execution time between the methods, thus the more accurate methods are also the most efficient in terms of wall-clock terms vs. accuracy.
We next investigate how well the numerical solution preserves quadratic invariants. Assumptions 2 and 3 are satisfied, and by Lemma 1 the exact solution of (29) is normpreserving, thus is constant, and the solution will stay on a circle. The matrices 0 , 1 and 2 trivially commute; thus Corollary 1 applies and all the midpoint based methods are expected to preserve the invariant. As 2 ( ) ≠ 0, Theorem 3 no longer applies. We thus expect to see a weak order one for the two trapezoidal SL schemes, which is exactly what is observed in Fig. 10. Finally, we want to study how well the quadratic invariants are preserved over time. Using each of the schemes, one solution path is calculated on the interval [0, 50], using step size ℎ = 2 −5 . The following two sets of parameters are chosen: = = 10 and = 100 and = 0.3. The values of I ( ) are shown in Fig. 11. Again, for visibility, only every 40th step is plotted.
In both cases, the mid-point based schemes preserve the invariant, while there is a significant drift-off for the two trapezoidal based schemes. The case = = 10 is the example depicted in Fig. 1, although there, the integration interval is only [0, 1]. In Fig. 11, where the integration is done over longer time, the drift-off leaves the TDSL scheme basically useless. The midpoint based methods preserve the invariant. In the low noise case, the drift-off is almost the same for the two trapezoidal schemes, and again, there is none or little drift-off in the midpoint schemes. The 95%-confidence intervals have been calculated and in all cases they span less than 14% of the shown mean values.

TDSL
TFSL MDSL MFSL Midpoint reference line with slope 1 12. Convergence of the TDSL, TFSL, MDSL, MFSL and implicit midpoint rule for the stochastic Fermi-Pasta-Ulam-Tsingou problem In all cases the Lawson schemes perform significantly better than the implicit midpoint rule. For small values of the two FSL schemes are only slightly better than the two DSL schemes. For larger values of , the FSL schemes outperform the DSL schemes. This is in line with what was observed for the Kubo oscillator in Fig. 8. But even in the small noise case, the feasible step sizes of the Lawson schemes are larger than the feasible step sizes for the implicit midpoint rule -the convergence deteriorates for the implicit midpoint rule approximately when ℎ = 2 −6 , whereas the convergence of the Lawson schemes only deteriorates around ℎ = 2 −3 , the SL schemes thus allowing to use much larger step sizes than the implicit midpoint rule. For larger values of , the midpoint rule is in practice useless for the step sizes applied here. The work-precision diagram Fig. 13 again demonstrates that even for this slighly larger problem the higher accuracy of the SL schemes, in particular the FSL schemes, compensates for the disadvantage of the additional computational work required for the calculation of the matrix exponentials. In the current example, the TFSL method is slightly the most efficient.

TDSL
TFSL MDSL MFSL Midpoint reference line with slope 1 reference line with slope 2 by = 1 2 ( 2 1, + 2 2 1, ) for the FPUT problem, and the result is presented in Fig. 14. In the small noise case, when = 0.02, the error is clearly dominated by the diffusion term, and the order of the method is close to the deterministic order two. Still, the FSL schemes are significantly more accurate than the midpoint rule, although the difference between those again is minuscule. When the noise level increases to = 0.2, the midpoint and the DSL schemes perform almost the same, while the errors of the FSL schemes are significantly smaller. We also observe the order two behaviour of the trapezoidal FSL scheme, although, as already mentioned, no quadratic invariants are conserved in this case. However, in the case where = 0 for = 1, . . . , , the transformed system (10) is in fact an SDE with additive noise, for which the weak order of the trapezoidal rule is two, see e. g. Milstein and Tretyakov [19]. In Fig. 15  In this paper, we proved that stochastic Lawson schemes under suitable conditions preserve linear and quadratic invariants. We proved that the trapezoidal stochastic Lawson scheme nearly preserves quadratic invariants if the diffusion terms are linear and fully included in the exponential. These results have been verified by numerical experiments.
For stochastic differential equations with highly oscillatory drift and diffusion we numerically demonstrated that full stochastic Lawson schemes allow for larger step-sizes than standard schemes, in the sense of being able to resolve high oscillations. All methods are implicit, and in our implementations, the cpu-time used for solving the nonlinear systems was dominant compared to the time used to evaluate the matrix exponentials, and in terms of accuracy versus computational work, the more accurate methods were the most efficient. This, of course, depends heavily on the problem at hand, for larger problems with more demanding matrix exponential evaluations, the drift Lawson methods might turn out to be preferable.
A Nicky Cordua Mattsson would like to thank the SDU e-Science centre for partially funding his PhD and the Department of Mathematics at the Norwegian University of Science and Technology for kindly hosting him during his visit. The authors would like to thank the unknown referees for their very helpful comments.