Numerical Picard Iteration Methods for Simulation of Non-Lipschitz Stochastic Differential Equations

: In this paper, we present splitting approaches for stochastic/deterministic coupled differential equations, which play an important role in many applications for modelling stochastic phenomena, e.g., ﬁnance, dynamics in physical applications, population dynamics, biology and mechanics. We are motivated to deal with non-Lipschitz stochastic differential equations, which have functions of growth at inﬁnity and satisfy the one-sided Lipschitz condition. Such problems studied for example in stochastic lubrication equations, while we deal with rational or polynomial functions. Numerically, we propose an approximation, which is based on Picard iterations and applies the Doléans-Dade exponential formula. Such a method allows us to approximate the non-Lipschitzian SDEs with iterative exponential methods. Further, we could apply symmetries with respect to decomposition of the related matrix-operators to reduce the computational time. We discuss the different operator splitting approaches for a nonlinear SDE with multiplicative noise and compare this to standard numerical methods.


Introduction
We are motivated to solve delicate stochastic differential equations, which are non-Lipschitz continuous, e.g., the nonlinear growth is polynomial or exponential, see [1]. Such SDEs are applied in different engineering problems, e.g., fluctuations in hydrodynamics, e.g., liquid jets in spraying and coating apparatuses, see [2][3][4]. For such modelling equation, the standard stochastic scheme have problems and we apply novel so called Picard iteration schemes, which apply Doléans-Dade exponential formula to overcome the local Lipschitzian problem, see [5]. The idea is to overcome the non-Lipschitzian problem with iterative scheme based on exponential functions, see [6]. The novelty in the paper is the additional splitting approaches, which are combined with the Doléans-Dade exponentials and included into the Picard iteration scheme. The splitting use ideas to decompose the underlying operators with respect to reduce the computational time, while we decompose into simpler solvable operators, which can be solved with much more faster solvers, see also [7].
Here, we deal with the following splitting ideas: • Symmetries: We consider symmetries in the operators and decompose to symmetrical operators, such that we could apply fast solver methods for each symmetric operator-part, see [7][8][9]. • Deterministic-Stochastic splitting: We decompose into deterministic and stochastic parts, while we apply fast deterministic and fast stochastic solver, see [10,11].
Such splitting approaches are important to optimize and speedup the solver-processes, see [12]. Further, the symmetries are also important for the specific solver method, while we apply iterative splitting methods, the decomposition into homogeneous and inhomogeneous solver parts are important, see [12]. We decompose the iterative splitting method into a fast solvable exp-matrix operator, which can apply fast symmetric solvers, see [7,13], and a right hand side part, which is a integral operator based on convolution operator. The convolution operator is based on the exp-matrix operators with an additional linear or nonlinear operator, see [14]. Here, we could also apply fast symmetric numerical integration methods and accelerate the solver-process, see [9].
In the paper, we concentrate on the following general form of the non-Lipschitz SDE, which is given as: where we have a(t, X(t)) = a 1 (t, X(t)) + a 2 (t, X(t)) = A(t, X(t))/X(t) ∈ IR\{0}, which is the drift coefficient and b(t, X(t)) = B(t, X(t))/X(t) ∈ IR\{0}, which is the diffusion coefficient. Further we assume a and b are bounded and W t is the Wiener process and we assume W t is independent from X 0 for t ≥ t 0 . For the sake of convenience, we have also assumed, that t 0 = 0, but we can also generalise to t 0 ∈ IR + 0 , which is only a shift in the initial conditions, see [15]. We have the following Assumption 1: Assumption 1. The SDE (1) can be rewritten into the stochastic differential equation: dX(t) = X(t)(a 1 (t, X(t)) + a 2 (t, X(t))) dt + X(t)b(t, X(t))dW t , t ∈ ×[0, T], We assume, that the Equation (1) has a unique strong solution X and it is positive, see also the work [5]. Further, we assume to deal with a decomposition of the full operator A = A 1 + A 2 , while the operator A 1 = a 1 X is fast computable as exp(A 1 ) with the, while A 2 = a 2 X is the more nonlinear operator and applied as a right-hand side, see [16].
In the Section 2, we derive the approximation of stochastic process with the Picard iterations and applied the numerical analysis of the new schemes. The numerical results are given in Section 3. In Section 4, we conclude and summarize our results.

Numerical Analysis of the Splitting Approaches
In the following, we discuss the splitting approaches, which are based on a homogeneous and inhomogeneous part.
For the stochastic differential equation (SDE), we define the homogeneous SDE as an linear or linearised SDE, which can be solved analytically or numerically, see [17][18][19]. Further, there exists also ideas for SDE and stochastic partial differential equations (SPDEs) to decompose into linear and nonlinear parts of the SDE, which can be solved linear and nonlinear stochastic methods, see [20].
We concentrate on linearized SDEs and we use the notation of homogeneous and inhomogeneous parts, which are wel-known for linear SDEs, see [17]. We assume, that the solution of the homogeneous SDE can be used with respect to the variation-of-constants formula, see [21], which is also wel-known from the theory of ODEs, see [12,22]. Such a formula also holds in the infinite dimensional case (under suitable assumptions, see [19]). Therefore, we obtain a solution of the inhomogeneous SDE, which is a mild solution, see [23].
We apply the splitting for the stochastic differential equation in the following homogeneous and inhomogeneous parts:

•
Homogeneous equation: where we have the solution X hom (t) = S m (t), where m are the iterative steps of the Picard iterations, see the Section 2.1.

•
Inhomogeneous equation: where we have the solution where m are the iterative steps of the Picard iterations, see the Section 2.1.
For the approximation of the scheme, we deal with three steps for the homogeneous part:
Picard iterations with Doléans-Dade solutions of the SDE, 3.
Discretisation of the Picard iterations in time.
For the inhomogenous part, we deal with the following parts: 1.
Discretisation of the Picard iterations in time for the inhomogeneous part, 2.
Approximation of the integral-formulation of the inhomogeneous part.

Homogeneous Equation
In the homogeneous part, we approximate the solution of the

Approximate the Diffusion Process
We consider the SDE (1) and we assume Y is the unique solution to the SDE: where c(X, t) = a(X, t) − 1 2 b 2 (X, t) and we assume that for all t ∈ [0, T], we have:

Picard Iterations with Doléans-Dade Solutions of the SDE
In the following, we construct the iterative process, which is based on the idea S m → S for m → ∞.
We have the following iterative scheme: where a 1,m+1 = a 1 (S m (t), t) and b m+1 = b(S m (t), t). Therefore, we obtain: Further, we can rewrite to and we have the following Picard iteration: We have the following convergence result of the Picard iterations: The Picard iteration (9) converge in L 2 with: whereC is a constant, depending on C, C L and S 0 .
Proof. We apply the recursion of the Lemma A1, see in the Appendix A.1. We obtain: we have Y 1 (t) = S 0 and Y 0 (t) = 0, we apply whereC is a constant depending on C and C L .
We have the following Assumption: We have the the stochastic differential equation: and a 1 : [0, T] × IR → IR and b : [0, T] × IR → IR be locally Lipschitz-functions. We assume, that the Equation (11) has a unique strong solution X and it is positive.
A generalization of the Equation (11) is given as: where we define a 1 (t, X(t)) = A 1 (t, X(t))/X(t) and b(t, X(t)) = B(t, X(t))/X(t) on IR\{0} and we assume a and b are bounded. We apply the iterated Doléans-Dade-process for X i , i ∈ IN + such that X i converge to S when i → ∞.
We have the following scheme: where we have . Then, we can prove, that we have the convergence of X i to X for i → ∞, see also [5].

Discretisation of the Picard Iterations in Time
We have the uniform discretization with ∆t = T/n for all ∆t = t j+1 − t j , with j = 0, . . . , n − 1. We have the recursive process with j = 0, . . . , n − 1 with the iterative steps i ≥ 1: where we have

Inhomogeneous Equation
For the inhomogeneous part, we have the following scheme: where a 2,i (x) is a right hand side based on our decomposition of a(x) = a 1 (x) + a 2 (x).

Discretisation of the Picard Iterations in Time for the Inhomogeneous Part
We have the uniform discretization with ∆t = T/n for all ∆t = t j+1 − t j , with j = 0, . . . , n − 1.
We have the recursively process with j = 0, . . . , n − 1 with the iterative steps i ≥ 1: where we have

Approximation of the Integral-Formulation of the Inhomogeneous Part
The integration of the inhomogeneous part can be done with numerical integration methods, e.g., trapezoidal rule, Romberg's method, see [24].
We have the following convergence results in Lemma 2 Lemma 2. The inhomogeneous part of the Picard iteration (20)-(21) converge in L 2 with: whereC is a constant, depending on C, C L and S 0 .
Proof. We deal with the inhomogeneous part and estimate the integral part, which is given as: andC is a constant. We obtain:S m+1 ≤Ŝ m+1â2 sup 0≤t≤T a 2,i,j (X(t)).
Then, we could deal with the idea in the homogeneous part, we apply: and we have the following Picard iteration: The estimation is done with the same results as in the homogeneous part and we receive: we haveŶ 1 (t) =â 2 andŶ 0 (t) = 0, we apply whereĈ is a constant depending on C and C L .
The convergence results for the homogeneous and inhomogeneous parts can be combined. Therefore, we have a convergence of the Picard iterative scheme in the homogeneous and inhomogeneous case.

Example 1.
For the mid-point rule, we following recursion: Based on the uniform discretization with ∆t = T/n for all ∆t = t j+1 − t j , with j = 0, . . . , n − 1 and the iterative steps i ≥ 1, we have: where we have

Numerical Examples
In the following, we study the different numerical examples, which are based on rational polynomials of the deterministic and stochastic coefficients. Such delicate an only local Lipschitz continuous coefficients are also applied in thin-liquid film models, see [2].
For simpler notation, we apply in the following the entries in the brackets for the labelling of the methods in the tables and figures.

First Example: Nonlinear SDE With Root-Function or Irrational Function
We deal with a test-example based on root-or irrational functions, see also [29]. The function is given as: The analytical solution can be derived by the Ito-Taylor expansion, see [29] and we obtain: We deal with the following numerical methods. Further, we apply the time-intervals ∆t = 1 N with N = 2 9 , 2 10 , 2 11 , 2 12 , 2 13 .

•
Euler-Maruyama-Scheme AB-splitting approach: We initialize t i with i = 0, . . . , N − 1, while t N = T and we have X(0) is the initial condition.
We deal with the 2 steps:

2.
B-part: where we have the solution X(t i+1 ) and we go to the next time-step till i = N − 1.
• ABA-splitting approach: We initialize t i with i = 0, . . . , N − 1, while t N = T and we have X(0) is the initial condition.
We deal with the 2 steps:

3.
A-step (∆t/2): where we have the solution X(t i+1 ) and we go to the next time-step till i = N − 1.
• Iterative Picard approach: we apply an Picard-Iteration where X 0 (t) = x(0), while we apply the implicit method in the drift term and the explicit method in the diffusion term.
The algorithm is given as: We initialize t n with n = 0, . . . , N − 1, while t N = T and we have X(0) is the initial condition.

Computation
4. i = i + 1, if i = I + 1 then we are done, else we go to Step (c) 5. n = n + 1, if n + 1 = N + 1 then we are done, else we go to Step (b) • Iterative Picard with Doléans-Dade exponential approach: we apply an Picard with Doléans-Dade exponential approach where X 0 (t) = x(0), while we apply the implicit method in the drift term and the explicit method in the diffusion term.
We apply the Picard-iterations with Doléans-Dade exponential approach and the splitting approach: where X 0 (t) = x(0), while we apply the implicit method in the drift term and the explicit method in the diffusion term.
The algorithm is given as: We initialize t n with n = 0, . . . , N − 1, while t N = T and we have X(0) is the initial condition.
We deal with the 2 loops (loop 1 is the computation over the full time-domain and loop 2 is the computation with i = 1, 2, . . .): Computation (we apply the full exp): 4. i = i + 1, if i = I + 1 then we are done, else we go to Step (c) 5. n = n + 1, if n + 1 = N + 1 then we are done, else we go to Step (b) We obtain a critical CFL condition, given as following: where ξ is N(0, 1) distributed. If the criterion is done all is fine and we go on. If we reach such a value, we reduce the time-step ∆t for the recent time step to ∆t = X 2/5 ξ 2 . For the next time-step, we could use the old time step and so on.
For the error analysis, we apply the different errors:

1.
Mean value at t = t n and J-sample paths: we deal with time-step ∆t = 1 N with N = 512, 1024, 2048, 4096 and the time-points are t n , n = 1, . . . , N with end-time-point t N = 1.0. For the sample paths, we apply J = 100 or J = 1000 and for the methods, we have method = {AB, ABA, Em, Mil, Picard, Picard − Doleans, Picard − Splitt − Doleans}.

2.
Local mean square error value at t = t n and J-sample paths:

3.
Global means square error over the full time-scale with different time-steps ∆t and J-sample path: where we apply the Equation (30)

Remark 1.
We could reduce the computational time, while we decompose the exp-operator and apply Pade-approximation. We saw in experiments, that we could accelerate the computation of 2-3 times, such that we are in the same region as the ABA-splitting approach, see the remarks in Appendix A.2.
In the following Table 1, we present the means-values at t = 1.0. In the following Figure 1, we present the exact solution of the first experiment. In the following Figure 2, we present the mean value, while we apply Equation (29), with different time-steps and for all numerical methods. In the following Figure 3, we present the mean square error, while we apply Equation (31), with different time-steps and for all numerical methods. In the following Figure 4, we present the performance of the schemes.   Figures 2 and 3, we present the mean and the mean-square-errors of the different schemes. We see the higher stability and the benefits of the iterative splitting, which applied symmetric techniques like the Doléans-Dade exponential approach. The locally Lipschitz operators are more stable to compute, while we could bound the operators in an exponential-approach, see [12]. In the performance, see Figure 7 in Section 3.2 we see, that the exponential based methods, like the Picard-Doleans or the Picard-Splitt-Doleans, are 2-3 times higher, than the standard schemes. But we obtain much more accurate and stable results and could use much more larger time-steps, such that we accelerate the solver process. Here, the restriction to the CFL-condition for the EM-, Milstein-schemes are higher and therefore the exponential based methods can use larger time-steps or are more accuracy with the same time-steps. For an increasing of the volatility in the stochastic differential equation, we have a stable method based on the Picard-Doleans or Picard-Splitt-Doleans. Here, we deal with an implicit behaviour of the underlying methods and we have to be aware of the numerical diffusion, which appear with larger time-steps, see [30].

Second Example: Linear/nonlinear SDE with Potential Function
We deal with a test-example based on root-or irrational functions, see also [29]. The function is given as: The analytical solution can be derived by the Ito-Taylor expansion, see [29] and we obtain: where W(t) = W t − W 0 = √ t ξ and ξ obeys the Gaussian normal distribution N(0, 1) with E(ξ) = 0 and E(ξ 2 ) = Var(ξ) = 1.
• AB-splitting approach: We initialize t i with i = 0, . . . , N − 1, while t N = T and we have X(0) is the initial condition.
We deal with the 2 steps:
• ABA-splitting approach: We initialize t i with i = 0, . . . , N − 1, while t N = T and we have X(0) is the initial condition.
We deal with the 2 steps:

3.
A-step (∆t/2): where we have the solution X(t i+1 ) and we go to the next time-step till i = N − 1.
The algorithm is given as: We initialize t n with n = 0, . . . , N − 1, while t N = T and we have X(0) is the initial condition.
We apply an Picard with Doléans-Dade exponential approach and right hand side where X 0 (t) = x(0), while we apply the implicit method in the drift term and the explicit method in the diffusion term.
The algorithm is given as: We initialize t n with n = 0, . . . , N − 1, while t N = T and we have X(0) is the initial condition.
We deal with the 2 loops (loop 1 is the computation over the full time-domain and loop 2 is the computation with i = 1, 2, . . .): Computation (we apply Version 1 or Version 2)

4.
Computation (we apply the full exp and integration of the RHS) 5. i = i + 1, if i = I + 1 then we are done, else we go to Step (c) 6. n = n + 1, if n + 1 = N + 1 then we are done, else we go to Step (b) We obtain a critical CFL condition, given as following: where ξ is N(0, 1) distributed. If the criterion is done all is fine and we go on. If we reach such a value, we reduce the time-step ∆t for the recent time step to ∆t = X 4 ξ 2 . For the next time-step, we could use the old time step and so on.

Remark 3.
In the mixed example with additional right hand side, we could reduce the computational time, while we decompose into symmetric operators. We have applied an exp-operator and via Taylor-or Pade-approximation and additional a convolution operator with the right-hand side, see also [12]. We saw the stable behaviour of such a splitting and could additionally accelerate the computation of 2-3 times, such that we are in the same region as the ABA-splitting approach, see the remarks in Appendix A.2.
In the following Figure 5, we present the mean value, which is computed with the Equation (29) with different time-steps and for all the numerical methods. In the following Figure 6, we present the mean square error, which is computed with the Equation (31) with different time-steps and for all the numerical methods. In the following Figure 7, we present the performance of the schemes.

Remark 4.
In the second example, which is a mixed linear and nonlinear stochastic differential equation, we also see the benefit and the efficiency of the iterative splitting approaches with Doléans-Dade exponential approach. Based on the local Lipschitzian of the coefficients, we could also apply larger timesteps for the exponential method with an iterative method. Such a Picard iterations solve the nonlinear parts of the operators and we could stabilize the nonlinear growth-parts with the exponential exponential approach. This allows to neglect the blow-ups in the exponential growth with sufficient large time-steps, see [10,16]. Further the benefit in splitting into symmetric parts of nonlinear and linear parts allows to accelerate the solver processes, see [31]. The benefit of larger time-steps for the iterative Picard-splitting with Doléans-Dade exponential approach can also applied to increasing volatilities in the SDEs. Here, we have to consider the implicit character of the iterative Picard-splitting approaches, see [30]. Such a behaviour can smear or smoothen out the solution of increasing volatility, like increasing diffusion, such that it might be necessary to apply multiscale methods with much more finer time-steps or an restriction of the time-step around the CFL-condition, see [30,32].

Conclusions
We have discussed a new iterative splitting approach, which is based on Picard iterations and Doléans-Dade exponential approach. We proof the convergence of such new schemes, which have embedded the Doléans-Dade exponential approach. Such a combinations allow to circumvent the global non-Lipschitz problems and therefore consider more stable local Lipschitz problems. Such a local problem can be solved with reformulations of exponential based operators, which are more stable in the numerical approach. We could embed such locally Lipschitz problems into a splitting approach which deals with exponential parts and a nonlinear solver given as Picard's method. Such a combination allows to obtain stable and convergent methods. Due to symmetric approaches of the operators, we could split into more appropriate linear, nonlinear and stochastic operators, which allows faster computations. The first numerical results are tested with linear and nonlinear stochastic differential equations with rational functions of the drift and diffusion parts. We present the benefits of the iterative splitting approach with the Doléans-Dade exponential, while we obtain an implicit method with stable results. Therefore, we are independent of the CFL conditions and accelerate the solver process. In future, we will test the splitting approach to more inhomogeneous problems and in the case of increasing volatility, e.g., real-life problems with stochastic lubrication models. By integration, we receive the results, that proves our lemma.