# Improvement of Split-Step Forward Milstein Schemes for SODEs Arising in Mathematical Physics

Authors: Hassan Ranjbar; Kazem Nouri; Leila Torkzadeh
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

In the present investigation, new explicit approaches by the Milstein method and increment function of the Jacobian derivative of the drift coefficient are designed. Several numerical tests such as Cox–Ingersoll–Ross process, stochastic Brusselator, and Davis-Skodje system are presented to illustrate the accuracy and the efficiency of our schemes. Furthermore, we show that the strong convergence rate of our procedures is approximately one.


---

## Body

## 1. Introduction

Many applied problems in different fields of science, engineering, and technology are required to find mathematical models [1–3]. So, the design of numerical techniques for solving these equations is a challenging task in the numerical analysis [4–6]. Some robust and efficient schemes for solving problems are brought forward [7, 8].Stochastic ordinary differential equations (SODEs) play a pivotal role in explaining some physical phenomena such as chemical reactions [9], financial mathematics [10], mathematical ecology [11], epidemiology [12], medicine [13], and population dynamics [14]. Generally, SODEs cannot be solved analytical, but many numerical solutions can be found, for instance, the split-step theta Milstein method [15], the least-squares method [16], the discrete Temimi–Ansari method [17], the improved Euler-Maruyamamethod [18], the five-stage Milstein method [19], the split-step Milstein method [20], the split-step Adams–Moulton Milstein method [21], the split-step forward Milstein method [22], and the Runge-Kutta method [23, 24].In the last two decades, the explicit type of numerical schemes for stiff SODEs has attracted many researchers’ attention see [25–27], for example. ODE solvers are the key to the success of this kind of numerical scheme. Several studies have recently introduced new explicit Milstein-type schemes using suitable ODE solvers. Yin and Gan [28] constructed an error-corrected Milstein scheme by arming the Milstein method with an error correction term. Nouri et al. [29] developed and applied a new explicit split-step Milstein method for stiff SODEs such as the stochastic Davis-Skodje system. These explicit numerical methods have a broad stability region, so they have been proposed to solve stiff SODEs.There are several attempts to construct numerical methods based on split-step forward Milstein schemes [30] to improve the numerical results for stiff SODEs.


## 2. Schemes Formulation

Consider the Itô SODE of the form [23](1)dYt=PYt,tdt+∑k=1nQkYt,tdWkt,with Yt0=Y0. In (1), P:ℝd×0,T⟶ℝd is the drift coefficient, Qk:ℝd×0,T⟶ℝd×n is the diffusion coefficient, and Wktk=1n are a one-dimensional standard Brownian motion.In [23, 24], Milstein’s scheme showed in the following manner:(2)Ym+1=Ym+hPYm,tm+∑k=1nQkYm,tmℐktm,tm+1+∑k1,k2=1nℒk1Qk2Ym,tmℐk1,k2tm,tm+1,where tm=t0+hm, h=T−t0/N, and(3)ℒk1=∑i=1dgk1,iYm,tm∂∂Yim,ℐktm,tm+1=∫tmtm+1dWkς≅ΔWk,m=hξk,ℐk1,k2tm,tm+1=∫tmtm+1∫tmς1dWk1ς2dWk2ς1,where ξk=N0,1 is a normally distributed random variable [23]. To construct new numerical methods, we divide scheme (2) as follows:(4)X1m·=hP·,tm,X2m·=−12h∑k=1nℒkQk·,tm,X3m·=Q·,tmℐktm,tm+1+∑k1,k2=1nℒk1Qk2·,tmJk1,k2tm,tm+1,where(5)ℐk1,k2tm,tm+1=Jk1,k2tm,tm+1,ifk1≠k2,Jk1,k2tm,tm+1−12h,ifk1=k2.Remark 1.
(see [20]) In the above equations, ℐktm,tm+1, ℐk1,k2tm,tm+1, and Jk1,k2tm,tm+1 denote stochastic Itô and Stratonovich integrals, respectively.
Now usingφYm=ehJP−1/hJP, we present the following fully explicit schemes based on the relations (4), the first improved three-stage Milstein (ITSM-i) scheme(6)Y^1m=Ym+φYmX1mYm,Y^2m=Y^1m+X2mY^1m,Ym+1=Y^2m+X3mY^2m.
The second improved three-stage Milstein (ITSM-ii) scheme is as follows:(7)Y^1m=Ym+φYmX1mYm,Y^2m=Y^1m+X3mY^1m,Ym+1=Y^2m+X2mY^2m.
The third improved three-stage Milstein (ITSM-iii) scheme is as follows:(8)Y^1m=Ym+X2mYm,Y^2m=Y^1m+φYmX1mY^1m,Ym+1=Y^2m+X3mY^2m.
The fourth improved three-stage Milstein (ITSM-iv) scheme is as follows:(9)Y^1m=Ym+X2mYm,Y^2m=Y^1m+X3mY^1m,Ym+1=Y^2m+φYmX1mY^2m.
The fifth improved three-stage Milstein (ITSM-v) scheme is as follows:(10)Y^1m=Ym+X3mYm,Y^2m=Y^1m+φYmX1mY^1m,Ym+1=Y^2m+X2mY^2m.
The sixth improved three-stage Milstein (ITSM-vi) scheme is as follows:(11)Y^1m=Ym+X3mYm,Y^2m=Y^1m+X2mY^1m,Ym+1=Y^2m+φYmX1mY^2m.
By denoting in approximation methods (6)–(11), JP is a Jacobian Matrix of drift coefficient P.


## 3. Numerical Experiments

In the present segment, some numerical examples are reported to demonstrate the efficiency and applicability of the proposed schemes. We compare the numerical methods (6)–(11) with the procedures TSM 1a-TSM 1f [30]. For the computations, we have used Matlab R2010a.Example 1.
Consider the following nonlinear stiff SDE [23],(12)dYt=−1+Yt1−Y2tdt+1−Y2tdWt,Y0=0,t∈0,1.
The exact solution is given by(13)Yt=1+Y0exp−2t+2Wt+Y0−11+Y0exp−2t+2Wt−Y0+1.
In Figure1, we compare the means of absolute errors (MAEs) and roots of mean-square errors (RMSEs) of our methods with slopes of dashed line 1.0. We use discrete Brownian paths over 0,1 with δ=10−14. Denoting Yi,N as the numerical approximation to YtN at endpoint T=tN in the i-th simulation of all M simulations with step-size h, we compare the analytic solution (13) with the numerical approximation using h=2rδ, r=0,1,…,7 over M=2000 different sample paths. Here, the MAEs and RMSEs are denoted as follows [18, 23]:(14)MAEs=1M∑i=1MYi,N−YitN,and(15)RMSEs=1M∑i=1MYi,N−YitN21/2.Figure 1
MAEs and RMSEs of methods (6)–(11) applied to the (12).Example 2.
Consider the Cox–Ingersoll–Ross (CIR) process [31, 32](16)dYt=ab−Ytdt+cYtdWt.
We know that the CIR process (16) is strictly positive if c2/2≤ab [33]. For initial value Y0=1, constants a=5, b=4, and c=4, Table 1 presents the percentage of negative values produced by the Milstein, (6)–(11) and TSM 1a-TSM 1f [30] schemes, for three different intervals T=1,2,4 and different values of the step-sizes h=1/2,1/4,1/8, along 105 paths. The results in Table 1 clearly show that our schemes other than the ITSM-ii method (7) can indeed guarantee positivity. Table 1 also shows that for step-size h=1/8, TSM 1a-TSM 1f schemes other than TSM 1b method ensure the positivity of the numerical solution.Table 1
Numerical results.
T=1T=2T=4Schemesh=1/2h=1/4h=1/8h=1/2h=1/4h=1/8h=1/2h=1/4h=1/8Milstein65.35%41.59%13.78%90.50%72.38%29.03%99.31%93.93%51.72%TSM1 a49.80%2.44%0%81.65%5.24%0%98.22%10.78%0%TSM1 b71.93%52.98%32.69%94.77%79.04%52.06%99.85%95.92%75.70%TSM1 c56.73%1.13%0%82.42%2.34%0%97.94%4.74%0%TSM1 d54.94%1.25%0%84.27%2.69%0%98.40%5.53%0%TSM1 e60.66%4.52%0%88.75%9.49%0%99.07%19.01%0%TSM1 f58.56%1.96%0%86.41%4.19%0%98.66%8.62%0%ITSM-i0%0%0%0%0%0%0%0%0%ITSM-ii55.58%53.26%41.65%79.52%75.62%61.63%95.54%93.33%83.31%ITSM-iii0%0%0%0%0%0%0%0%0%ITSM-iv0%0%0%0%0%0%0%0%0%ITSM-v0%0%0%0%0%0%0%0%0%ITSM-vi0%0%0%0%0%0%0%0%0%Example 3.
Let us consider the following stochastic Brusselator [23]:(17)dY1t=α−1Y1t+αY12t+(1+Y1t)2Y2tdt+σY1t1+Y1tdWt,dY2t=−αY1t−αY12t−(1+Y1t)2Y2tdt−σY1t1+Y1tdWt,with α=1.9 and σ=0.1. Figure 2 displays the RMSEs with step-sizes h=2−r, r=7,…,14 at T=1, initial data (Y10,Y20)=−0.1,0, and M=2000 sample paths. The numerical methods with step-size 10−14 are used to generate the reference solutions. In Figure 2, besides our results, also the results of the TSM 1a-TSM 1f and Milstein methods are presented. This figure shows that the convergence rate of the methods for the stochastic Brusselator system (17) is approximately 1.0.Figure 2
RMSEs of methods applied to the stochastic Brusselator system (17).Example 4.
We investigate a Davis-Skodje system [29, 34, 35].(18)dY1t=−Y1tdt+μY1tdW1t,dY2t=−γY2t+γY1t1+Y1t−Y1t(1+Y1t)2dt+σγY2tdW2t,with initial (Y10,Y20)=6,0.853. For analysis of the Davis-Skodje system (18), we choose the parameters as follows:(i)
Case I:μ=σ=0.1,γ=100,h=1/50Case II:μ=σ=0.1,γ=1000,h=1/500
The mean of 1000 sample paths of the nonlinear system (18) for two cases with T=10 is shown in Figure 3. We can see that our schemes converge quickly towards the asymptotic solution 0,0.Figure 3
Numerical results of methods applied to the Davis-Skodje system (18).
(a)(b)


## 4. Conclusions

New explicit schemes have proved helpful for efficiently solving stochastic problems with both stiff and nonstiff components. The present paper proposes new methods for the numerical approximation of SODEs driven by multidimensional standard Brownian motion. For construction methods (6)–(11), explicit Milstein schemes and increment function of the Jacobian derivative of drift coefficient are used. The efficiency and applicability of the suggested methods are confirmed by some numerical examples such as Cox–Ingersoll–Ross process, stochastic Brusselator, and Davis-Skodje system. In Table 1, we observe the positive solution of the CIR model using numerical methods Milstein (6)–(11) and TSM 1a-TSM 1f [30]. From the results of Table 1, we find that our methods preserve the positivity of the CIR model for all step sizes. For two set parameter values of the Davis-Skodje system, Figure 3 reveals that our schemes converge quickly towards the asymptotic solution. Also, our numerical tests show that the strong convergence rate of the schemes (6)–(11) is approximately 1.0 (see Figures 1 and 2).



---

*Source: 7498865-2022-06-29.xml*