 # Exponential integrators for the stochastic Manakov equation

This article presents and analyses an exponential integrator for the stochastic Manakov equation, a system arising in the study of pulse propagation in randomly birefringent optical fibers. We first prove that the strong order of the numerical approximation is 1/2 if the nonlinear term in the system is globally Lipschitz-continuous. Then, we use this fact to prove that the exponential integrator has convergence order 1/2 in probability and almost sure order 1/2, in the case of the cubic nonlinear coupling which is relevant in optical fibers. Finally, we present several numerical experiments in order to support our theoretical findings and to illustrate the efficiency of the exponential integrator as well as a modified version of it.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1. Introduction

Optical fibers play an important role in our modern communication society [agrawal2007nonlinear]. In order to model the light propagation over long distance in randomly varying birefringent optical fibers, the Manakov PMD equation was derived from Maxwell’s equations in [Wai-1996]. As noted in [garnierm06], polarization mode dispersion (PMD) is one of the main limiting effects of high bit rate transmission in optical fiber links. In addition, the work [garnierm06] proves that the asymptotic dynamics of the Manakov PMD equation is given by a stochastic nonlinear evolution equation in the Stratonovich sense: the stochastic Manakov equation, see below. In the present article, we perform a numerical analysis of this stochastic partial differential equation (SPDE).

We now review the literature on the numerical analysis of the stochastic Manakov equation. The work [Gazeau:13], see also [GazeauPhd], numerically studies the impact of noise on Manakov solitons and soliton wave-train propagation by the following time integrators: the nonlinearly implicit Crank–Nicolson scheme, the linearly implicit relaxation scheme, and a Fourier split-step scheme. For instance, it is conjectured that, in the small-noise regime and over short distances, solitons are not strongly destroyed and are stable. Reference [MR3166967], see also [GazeauPhd], proves that the order of convergence in probability of the Crank–Nicolson scheme is . On top of that, it is shown that this numerical integrator preserves the -norm as does the exact solution to the stochastic Manakov equation. Furthermore, it is numerically observed that the almost-sure order of convergence of the relaxation scheme and the split-step scheme is .

The main goal of this article is to present and analyse a linearly implicit exponential integrator for the time discretisation of the stochastic Manakov equation. Exponential integrators for the time integration of deterministic or stochastic (partial) differential equations are nowadays widely used and studied as witnessed by the recent works [MR2413146, MR2536086, MR2578878, MR2652783, MR2995211, MR3047942, MR3033008, MR3353942, MR3463447, MR3484400, MR3736655, MR3771721, acqs18, MR3663004] and references therein. Beside having the same orders of convergence as the nonlinearly implicit Crank–Nicolson scheme from [MR3166967], the proposed exponential integrators offer additional computational advantages as illustrated below.

## 2. Setting and notation

Let be a probability space on which a three-dimensional standard Brownian motion is defined. We endow the probability space with the complete filtration generated by . In the present paper, we consider the nonlinear stochastic Manakov system [MR3166967]

 (1) idX+∂2xXdt+|X|2Xdt+i√γ3∑k=1σk∂xX∘dWk=0,

where is the unknown function with values in with and , the symbol denotes the Stratonovich product, measures the intensity of the noise, is the nonlinear coupling, and , and are the classical Pauli matrices defined by

 σ1=(0110),σ2=(0−ii0),andσ3=(100−1).

The mild form of the stochastic Manakov equation reads

 (2) X(t)=U(t,0)X0+i∫t0U(t,s)F(X(s))ds,

where denotes the initial value of the problem, for with is the random unitary propagator defined as the unique solution to the linear part of (1), and .

Let . We define the Lebesgue spaces of functions with values in . We equip with the real scalar product . Further, for , we denote the space of functions in with their first derivatives in . The norm in is denoted by .

We now recall the local existence and uniqueness result for solutions to (1) obtained in [MR3024974] (see also [GazeauPhd]).

###### Theorem 1 (Theorem 1.2 in [Mr3024974]).

Consider the initial value , then there exists a maximal stopping time and a unique solution (in the probabilistic sense) to (1) such that -a.s. Furthermore, the norm is almost surely preserved: for . Moreover, the following alternative holds for the maximal existence time of solutions to (1):

 τ∗(X0,ω)=+∞orlimsupt↗τ∗(X0,ω)∥X(t)∥H1=+∞.

Finally, if the initial value belongs to for some , then the corresponding solution also belongs to almost surely.

As seen above, the norm of the solution is preserved just as for the deterministic Manakov equation (i. e. (1) with ). Furthermore, as noted by [MR3166967], the occurrence of blow-up in the stochastic Manakov equation (1) remains an open question.

For the time discretisation of the stochastic Manakov system (1), one has to face two issues. First, the linear part of the equation generates a random unitary propagator which is not easy to compute exactly. In particular, since the Pauli matrices do not commute, it is not the product of the stochastic semi-groups associated to each Brownian motion with the group generated by . Second, the nonlinear coupling term often leads to implicit numerical methods that are costly, see for instance the Crank–Nicolson scheme proposed in [MR3166967].

Therefore, we propose to discretise the stochastic Manakov equation with an exponential integrator, that we now define. Let be a fixed time horizon and consider an integer . We define the step size by and denote discrete times by , for . Discretising the integral present in the mild form (2) (by an explicit Euler step) as well as the random propagator (by a midpoint rule), one gets the following exponential integrator

 (3) Xn+1=Uh,n(Xn+ihF(Xn)),

where , with is the identity operator and . Here, is the , for , are i.i.d. Wiener increments. Since is an approximation of the exponential of the linear random differential operator in (1), we choose to name the scheme (3) exponential integrator.

The exponential integrator (3) thus approximates solutions to the stochastic Manakov equation (1), , at the grid points .

Iterating the recursion (3), one gets the discrete mild form for the exponential scheme

 (4) Xn=Un,0hX0+ihn−1∑ℓ=0Un,ℓhF(Xℓ),

where .

This linearly implicit method is well-defined for all , and one has that for all , (respectively , resp ) provided that (resp. , resp. ). Moreover, if one assumes that is bounded by on , then one has almost surely for all and such that , .

In the following, in the proofs of our results, we denote by a positive constant that may change from one line to the other, but that does not depend on the parameters indicated in the results’ statements.

## 3. Convergence analysis of the exponential integrator

This section presents the main results of the article and gives the corresponding proofs. We start by considering the stochastic Manakov equation (2), where the nonlinearity is assumed to be globally Lipschitz-continuous. We show strong order of convergence for the exponential scheme (3) in that case. Then, we analyse the case of a cubic nonlinearity, i. e. , which is of course not globally Lipschitz-continuous, and we show order of convergence in probability , as well as order of convergence almost surely, for the exponential scheme (3). The main steps of the proofs use similar arguments as in [MR3166967, MR3312594, MR3736655]

as well as other works on the numerical analysis of SPDEs. However, some technical details are handled differently in this paper (see for example the estimation of

and below).

### 3.1. The Lipschitz-continuous case

We present a strong convergence analysis of the exponential integrator (3) when applied to the stochastic Manakov equation (2) when is globally Lipschitz-continuous on . This is the case, for instance, when one introduces a cut-off function for the cubic nonlinearity present in (1): Let and , with , and on . For , we set and define . We thus obtain a bounded globally Lipschitz-continuous function from to , which sends bounded subsets of to bounded subsets from , resp. of to . For ease of presentation, we denote the stochastic processes and solutions to the continuous problem and to the discrete problem, instead of using the notation and that we will use later in the paper, to point out the difference between truncated and untruncated problems and solutions.

###### Theorem 2.

Let , , and . Consider a bounded Lipschitz nonlinearity , defined as above, in the stochastic Manakov equation (2). Then, there exist a constant such that the exponential integrator (3) has strong order of convergence : There exists such that

###### Proof.

Let us denote the difference by . Using the definitions of the numerical and exact solutions, we thus obtain

 ∥en∥1 =∥∥ ∥∥Un,0hX0+ihn−1∑ℓ=0Un,ℓhF(Xℓ)−U(tn,0)X0−i∫tn0U(tn,s)F(X(s))ds∥∥ ∥∥1 ≤∥∥(Un,0h−U(tn,0))X0∥∥1+∥∥ ∥∥n−1∑ℓ=0∫tℓ+1tℓ(Un,ℓhF(Xℓ)−U(tn,s)F(X(s)))ds∥∥ ∥∥1 =:In1+In2.

We begin by estimating the term using

 In2 −U(tn,tℓ)F(X(s))+U(tn,tℓ)F(X(s))−U(tn,s)F(X(s)))ds∥∥1 +∥∥ ∥∥n−1∑ℓ=0∫tℓ+1tℓ(Un,ℓh−U(tn,tℓ))F(X(s))ds∥∥ ∥∥1+∥∥ ∥∥n−1∑ℓ=0∫tℓ+1tℓ(U(tn,tℓ)−U(tn,s))F(X(s))ds∥∥ ∥∥1 =:Jn1+Jn2+Jn3+Jn4.

We next bound the expectation of each of the four terms above to the power . Using the fact that the nonlinearity is globally Lipschitz-continuous from to , and that is an isometry on all (see Appendix 5), the first term can be estimated as follows

 E[maxn=0,1,…,N(Jn1)2p] ≤E⎡⎣maxn=0,1,…,N(n−1∑ℓ=0∫tℓ+1tℓds∥∥Un,ℓh(F(Xℓ)−F(X(tℓ)))∥∥1)2p⎤⎦ ≤CT2pE[maxℓ=0,1,…,N∥∥Xℓ−X(tℓ)∥∥2p1]=CT2pE[maxℓ=0,1,…,N∥∥eℓ∥∥2p1].

Similarly, for the second term we obtain

 E[maxn=0,1,…,N(Jn2)2p] ≤CE⎡⎣maxn=0,1,…,N(n−1∑ℓ=0∫tℓ+1tℓ∥X(tℓ)−X(s)∥1ds)2p⎤⎦ ≤CE⎡⎣(N−1∑ℓ=0suptℓ≤s≤tℓ+1∥X(tℓ)−X(s)∥1h)2p⎤⎦.

Using Hölder’s inequality, one then gets

 E[maxn=0,1,…,N(Jn2)2p] ≤Ch2pE⎡⎣(N−1∑ℓ=012p/(2p−1))2p−1N−1∑ℓ=0(suptℓ≤s≤tℓ+1∥X(tℓ)−X(s)∥1)2p⎤⎦ ≤Ch2pN2p−1N−1∑ℓ=0E[suptℓ≤s≤tℓ+1∥X(tℓ)−X(s)∥2p1]≤Ch2pN2php≤Chp,

where we have used the estimate from Lemma 5.4 (temporal regularity of the mild solution) in [MR3166967].

Using Lemma 5.3 (uniform boundedness of the mild solution in ) in [MR3166967], as well as the fact that sends bounded sets of to bounded sets of , we infer, using Proposition 2.2 (strong convergence for linear problems, i. e. when ) in [MR3166967] that one has

 ∀s∈[0,T],E(maxn∈{0,…,N}maxℓ∈{0,…,n}∥∥(Un,ℓh−U(tn,tℓ))F(X(s))∥∥2p1)≤Chp.

Therefore, we estimate the third term, using Hölder’s inequality, as follows

 E[maxn=0,1,…,N(Jn3)2p] ≤E⎡⎣maxn=0,1,…,N(n−1∑ℓ=0∫tℓ+1tℓ∥∥(Un,ℓh−U(tn,tℓ))F(X(s))∥∥1ds)2p⎤⎦ ≤E⎡⎣(∫T0maxn=0,1,…,Nmaxℓ=0,…,n∥∥(Un,ℓh−U(tn,tℓ))F(X(s))∥∥1ds)2p⎤⎦ ≤CT2p−1E[∫T0maxn=0,1,…,Nmaxℓ=0,…,n∥∥(Un,ℓh−U(tn,tℓ))F(X(s))∥∥2p1ds] ≤CT2p−1∫T0E[maxn=0,1,…,Nmaxℓ=0,1,…,n∥∥(Un,ℓh−U(tn,tℓ))F(X(s))∥∥2p1]ds≤Chp.

To bound the last term, we first use the isometry property of the continuous random propagator and Hölder’s inequality to get

 E[maxn=0,1,…,N(Jn4)2p] ≤E⎡⎣maxn=0,1,…,N(n−1∑ℓ=0∫tℓ+1tℓ∥(U(tn,tℓ)−U(tn,s))F(X(s))∥1ds)2p⎤⎦ ≤E⎡⎣maxn=0,1,…,N(n−1∑ℓ=0∫tℓ+1tℓ1×∥(Id−U(s,tℓ))F(X(s))∥1ds)2p⎤⎦ ≤Ch2pE⎡⎢ ⎢⎣⎛⎜ ⎜⎝(N−1∑ℓ=012p2p−1)2p−12p(N−1∑ℓ=0suptℓ≤s≤tℓ+1∥(Id−U(s,tℓ))F(X(s))∥2p1)12p⎞⎟ ⎟⎠2p⎤⎥ ⎥⎦ (5) ≤Ch2pN2p−1N−1∑ℓ=0E[suptℓ≤s≤tℓ+1∥(Id−U(s,tℓ))F(X(s))∥2p1].

In order to estimate the expectation above, we write this term as

 E [suptℓ≤s≤tℓ+1∥(Id−U(s,tℓ))(F(X(tℓ))−F(X(tℓ))+F(X(s)))∥2p1] (6)

The first term in the equation above is the exact solution to the linear SPDE with initial value at initial time which has the mild Ito form

 Z(t)−F(X(tℓ))=(S(t−tℓ)−Id)F(X(tℓ))+i√γ3∑k=1∫ttℓS(t−u)σk∂xZ(u)dWk(u),

where is the group solution to the free Schrödinger equation. Owing at the regularity property of the group (see for instance the first inequality in the proof of [GazeauPhd, Lemma 4.2.1]), the fact that the exact solution is almost surely bounded in , that sends bounded sets of to bounded sets of , and Burkholder–Davis–Gundy’s inequality (for the second term), one obtains the following bound for the first term in (3.1)

 E[suptℓ≤s≤tℓ+1∥(Id−U(s,tℓ))F(X(tℓ))∥2p1]≤Chp.

Using the fact that the random propagator is an isometry, that is globally Lipschitz-continuous, and the regularity property of the exact solution , one gets the estimate

for the second term in (3.1).

Combining the estimates above, one finally arrives at the bound

 E[maxn=0,1,…,N(Jn4)2p]≤Ch2pN2p−1N−1∑ℓ=0hp≤Chp.

Altogether we thus obtain

 ≤CE[maxn=0,1,…,N(In1)2p]+CT2pE[maxn=0,1,…,N∥en∥2p1]+Chp

using once again [MR3166967, Proposition 2.2].

For small enough, i. e. such that , the inequality above gives

 ≤C1−CT2p1hp,

on . In order to iterate this procedure, we impose, if necessary, that is small enough (or, equivalently, that is big enough), to ensure that can be chosen as before and as some integer multiple of (say for some positive integer ), while is some multiple integer of (say for some positive integer ). To obtain a bound for the error on the longer time interval , we iterate the procedure above by choosing and estimate the error on the interval . We repeat this procedure, times, up to final time . This can be done since the above error estimates are uniform on the intervals for (with a slight abuse of notation for the time interval):

 E[max[Tm,Tm+1]∥Xn−Xm(tn)∥2p1] ≤CEhp,

where is the error constant obtained above, are discrete times in , is the exact solution with initial value , denotes the exact solution with initial value at time , and corresponds to numerical solutions at time for . For the total error, we thus obtain (details are only written for the first two intervals)

 +E[max[T1,T2]∥Xn−X(tn)∥2p1]+…+E[max[TM−1,TM]∥Xn−X(tn)∥2p1] ≲CEhp+CEhp+CLE[∥∥Y1−X(T1)∥∥2p1]+…≲CEhp+CLCEhp+… ≲CEhp+CLCEhp+C2LCEhp+…+CM−1LCEhp≲Chp,

where is the Lipschitz constant of the exact flow of (1) from to itself and the last constant is independent of and with for big enough. This concludes the proof of the theorem. ∎

### 3.2. Convergence in the non-Lipschitz case

Using the above result as well as ideas from [MR1873517, MR2268663, MR3166967, MR3312594, MR3736655], one can show convergence in probability of order and almost sure convergence of order for the exponential integrator (3) when applied to the stochastic Manakov equation (1).

###### Proposition 3.

Let and . Denote by the maximum stopping time for the existence of a strong adapted solution, denoted by , of the stochastic Manakov equation (1). For all stopping time a.s. there exists such that we have

 ∀h∈(0,h0),limC→∞P(max0≤n≤Nτ∥Xn−X(tn)∥H1≥Ch1/2)=0,

where denotes the numerical solution given by the exponential integrator (3) with time step and .

###### Proof.

For , let us denote by , resp. , the exact, resp. numerical, solutions to the stochastic Manakov equation (2) with a truncated nonlinearity . We denote by a positive constant such that for all , .

Fix , , . Let be a stopping time such that a.s. . By Theorem 1, there exists an such that . We have the inclusion

 {max0≤n≤Nτ∥Xn−X(tn)∥1≥ε} ⊂{max0≤n≤Nτ∥X(tn)∥1≥R0−1} ∪({max0≤n≤Nτ∥Xn−X(tn)∥1≥ε}∩{max0≤n≤Nτ∥X(tn)∥1

Taking probabilities, we obtain

 P({max0≤n≤Nτ∥Xn−X(tn)∥1≥ε})≤ε/2+P({max0≤n≤Nτ∥Xn−X(tn)∥1≥ε}∩{max0≤n≤Nτ∥X(tn)∥1

In order to estimate the terms on the right-hand side, we define the random variable

, with the convention that if the set is empty. If then we have by triangle inequality

 max0≤n≤nε−1∥Xn∥1=max0≤n≤nε−1∥Xn−X(tn)+X(tn)∥1≤ε+R0−1≤R0.

By definition of the exponential integrator (3), for , we have

 ∥Xnε∥1≤∥∥Xnε−1∥∥1+h∥∥|Xnε−1|2Xnε−1∥∥1≤R0+κh∥∥Xnε−1∥∥31≤R0+κhR30≤4R0,

in this case and thus for .

If , then thanks to the definition of . Therefore we get . Furthermore, by definition of , we have . We then deduce that

 {max0≤n≤Nτ∥Xn−X(tn)∥1≥ε}={max0≤n≤Nτ∥Xn−X(tn)∥1≥ε}∩{nε≤Nτ}.

Combining the above, using Markov’s inequality as well as the strong error estimates from Theorem 2, since a.s., there exists such that

 P(max0≤n≤Nτ∥Xn−X(tn)∥1≥ε,nε≤Nτ,max0≤n≤Nτ∥X(tn)∥1

This last term is smaller than for small enough. All together we obtain

and thus convergence in probability.

To get the order of convergence in probability, we choose such that for all small enough, . As above, for all positive real number , we have

 {max0≤n≤Nτ∥Xn−X(tn)∥1≥Ch1/2} ⊂{max0≤n≤Nτ∥X(tn)∥1≥R1} ∪({max0≤n≤Nτ∥Xn−X(tn)∥1≥Ch1/2}∩{max0≤n≤Nτ∥X(tn)∥1

Taking probabilities and using Markov’s inequality as well as the strong error estimate from Theorem 2, we obtain

 P({max0≤n≤Nτ∥Xn−X(tn)∥1≥Ch1/2}) ≤ε2+P({max0≤n≤Nτ∥∥Xn4R1−X4R1(tn)∥∥1≥Ch1/2}) ≤ε2+K(4R1,γ,T,p,∥X0∥6)C2p,

since almost surely. For large enough, we infer

 P({max0≤n≤Nτ∥Xn−X(tn)∥1≥Ch1/2})≤ε2+ε2=ε,

uniformly for