A Numerical Approach of Handling Fractional Stochastic Differential Equations

Abstract

This work proposes a new numerical approach for dealing with fractional stochastic differential equations. In particular, a novel three-point fractional formula for approximating the Riemann–Liouville integrator is established, and then it is applied to generate approximate solutions for fractional stochastic differential equations. Such a formula is derived with the use of the generalized Taylor theorem coupled with a recent definition of the definite fractional integral. Our approach is compared with the approximate solution generated by the Euler–Maruyama method and the exact solution for the purpose of verifying our findings.

1. Introduction

Due to its superior properties for many problems, solutions to fractional stochastic differential equations (FSDEs) driven by the Brownian motion have recently received much attention from scientific researchers. For instance, the authors in [1] looked at a stochastic viscoelastic wave equation with nonlinear damping and logarithmic nonlinear source terms that established a blow-up result. In [2], the existence and uniqueness of mild solutions for neutral delay Hilfer fractional integrodifferential equations were studied with fractional Brownian motion, and consequently certain sufficient conditions for controllability of neutral delay Hilfer fractional differential equations were established with fractional Brownian motion as well. In [3], noninstantaneous impulsive conformable fractional stochastic delay integrodifferential system driven was studied by Rosenblatt process, and accordingly sufficient conditions for approximate controllability and null controllability were established for the considered problem. More recently, the authors in [4] discussed the essential concept behind the multilevel Monte Carlo approach with the exact coupling via performing several numerical implementations. Through this paper, we aim to study the following formula that represents an FSDE formulated in the Caputo sense with a noisy environment:

$$D_{*}^{\alpha }X(t,w)=f(t,X(t,w))dt+g(t,X(t,w))dW(t,w),$$

where f is the drift coefficient, g is the diffusion coefficient, and $W(t,w)$ is the Wiener process, which is also called Brownian motion. For simplicity, we consider $X(t,w)=X\left(t\right)$ and $W(t,w)=W_t$. This gives:

$$D_{*}^{\alpha }X\left(t\right)=f(t,X\left(t\right))dt+g(t,X\left(t\right))d{W}_t,$$

(1)

where $0\le t\le T$. The Wiener process $W\left(t\right)$, which is a stochastic process indexed by nonnegative real number t, satisfies the following three conditions:

• $W_0=1$.
• $W_t–W_s\sim \sqrt{t-s}N(0,1)$ for $0\le s<t$, where $N(0,1)$ indicates a standard normal distribution.
• The two increments $W_t–W_s$ and $W_{\tau }–W_{\upsilon }$ are independent on distinct time intervals for $0\le s<t<\tau <\upsilon$.

Indeed, Equation (1) is an example of an FSDE. Due to the fact that their applications are viewed as stochastic processes, such as those in mechanics, medicine, physics, biology, population dynamics, and finance, these equations are crucial in many areas of business and research [5,6]. In particular, such equations were established based on the fact that the deterministic differential equation can be modified by including a random term. In the same regard, due to the fact that these equations have not frequently exact solutions, numerical methods must be used to approximate their solutions [7,8,9]. The FSDEs are thought of as a natural type of the fractional-order systems that can involve specific random terms because these systems are commonly produced in many real-life models [10,11,12]. However, in order to obtain further insights about some of numerical solutions of stochastic differential equations, the reader may refer to references [13,14,15].

The Euler–Maruyama method is a method for the approximate numerical solution of a stochastic differential equation. It is an extension of the Euler method for ordinary differential equations to stochastic differential equations [16]. It is named after Leonhard Euler and Gisiro Maruyama [13]. As the traditional Euler method, this method is considered unacceptably poor, and requires a too small step size to achieve some serious accuracy [13]. From this point of view, the motivation behind this study is to establish a more efficient novel numerical method than some other existence numerical methods for dealing with FSDEs. This method depends on establishing a new three-point fractional formula for approximating Riemann–Liouville integrator, which is derived with the use of the generalized Taylor theorem coupled with a recent definition of the definite fractional integral.

The organization of this paper is arranged as follows: Section 2 aims to recall some basic facts and definitions connected with stochastic differential equations. Section 3 demonstrates the main results of this work so that it contains an established numerical method for solving the FSDEs with the help of using the so-called modified three-point fractional formula for approximating the Riemann–Liouville fractional integrator. Such a formula will be derived first from one of the fractional calculus’s most important results: the generalized Taylor theorem. Then in Section 4, we will apply the derived formula to solve the FSDEs. Section 5 illustrates numerical results that confirm the theoretical findings of this work, followed by the final section that summarizes the conclusion.

2. Preliminaries

In this section, we recall some preliminaries and basic results related to fractional calculus coupled with some needed concepts connected with stochastic differential equations. For more details about fractional calculus and its applications, the reader may refer to references [17,18,19,20].

Definition 1.

Let α be a real nonnegative number. Then the Riemann–Liouville fractional-order integrator $J_a^{\alpha }$ is defined by:

$$J_a^{\alpha }f\left(x\right)=\frac{1}{{\Gamma }\left(\alpha \right)}{\int }_a^x{(x-t)}^{\alpha -1}f\left(t\right)dt,a\le x\le b.$$

(2)

In what follows, we recall certain properties of the Riemann–Liouville fractional-order integral operator for completeness [17]:

$$\begin{array}{cc}\hfill \left(1\right)& J_a^{0}h\left(t\right)=h\left(t\right).\hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill \left(2\right)& J_a^{\mu }{(t-a)}^{\gamma }=\frac{{\Gamma }(\gamma +1)}{{\Gamma }(\mu +\gamma +1)}{(t-a)}^{\mu +\gamma },\gamma \ge -1.\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill \left(3\right)& J_a^{\mu }{J}_a^{\beta }h\left(t\right)=J_a^{\beta }{J}_a^{\mu }h\left(t\right)\mu ,\beta \ge 0.\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill \left(4\right)& J_a^{\mu }{J}_a^{\beta }h\left(t\right)=J_a^{\mu +\beta }h\left(t\right)\mu ,\beta \ge 0.\hfill \end{array}$$

(6)

Definition 2.

Let $\alpha \in {\mathbb{R}}^{+}$ and $m=\lceil \alpha \rceil$ such that $m-1<\alpha \le m$. Then the Caputo fractional-order differentiator of order α is given by:

$$D_{*}^{\alpha }f\left(x\right)=\frac{1}{{\Gamma }(m-\alpha )}{\int }_a^x{(x-t)}^{m-\alpha -1}{f}^{\left(m\right)}\left(t\right)dt,x>a.$$

(7)

In the following content, we list some properties of the Caputo differentiator [17]:

$$\begin{array}{cc}\hfill \left(1\right)& D_{*}^{\mu }c=0,\mathrm{where}c\mathrm{is}\mathrm{constant}.\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill \left(2\right)& D_{*}^{\mu }{(t-a)}^{\rho }=\frac{{\Gamma }(\rho +1)}{{\Gamma }(\rho -\mu +1)}{(t-a)}^{\rho -\mu },\mathrm{where}\rho >\mu -1.\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill \left(3\right)& D_{*}^{\mu }(\mu h\left(t\right)+\omega g\left(t\right))=\mu D_{*}^{\mu }\left(h\left(t\right)\right)+\omega D_{*}^{\mu }\left(g\left(t\right)\right),\hfill \end{array}$$

(10)

where $\mu$ and $\omega$ are constant. In the same regard, we report below some other properties related to the composition between the previous two operators [17]:

$$D_{*}^{\alpha }{J}_0^{\alpha }h\left(t\right)=h\left(t\right),$$

(11)

and

$$J_0^{\alpha }{D}_{*}^{\alpha }h\left(t\right)=h\left(t\right)-\sum _{i=1}^n{h}^i\left(0^{+}\right)\frac{t^i}{i!},$$

(12)

where $t>0$ and $n-1<\alpha \le n$ such that $n\in \mathbb{N}$.

Theorem 1

([21] (generalized Taylor’s theorem)). Suppose that $D_{*}^{k\alpha }f\left(x\right)\in C(0,b]$ for $k=0,1,\dots ,n+1,$ where $0<\alpha \le 1$. Then the function f can be expanded about $x=x_0$ as:

$$\begin{array}{c}\hfill f\left(x\right)=\sum _{i=0}^n\frac{x^{i\alpha }}{{\Gamma }(i\alpha +1)}{D}_{*}^{i\alpha }f\left(x_0\right)+\frac{x^{(n+1)\alpha }}{{\Gamma }\left(\right(n+1)\alpha +1)}{D}_{*}^{(n+1)\alpha }f\left(\xi \right),\end{array}$$

(13)

with $0<\xi <x$, $\forall x\in (0,b]$.

Ito’s formula is a key component in the Ito Calculus, used to determine the derivative of a time-dependent function of a stochastic process. For more overview about such formula, the reader may refer to the reference [22].

Theorem 2

([13] (Itô formula)). Let $d\xi \left(t\right)=adt+bdw\left(t\right)$ and let $f(x,t)$ be a continuous function in $(x,t)\in {\mathbb{R}}^1\times [0,\infty )$ with partial derivatives $f_x,f_{xx},f_t$. Then the process $f\left(\xi \right(t),t)$ has a stochastic differential, given by:

$$\begin{array}{cc}\hfill df\left(\xi \right(t),t)& =[f_t(\xi \left(t\right),t)+f_x(\xi \left(t\right),t)a\left(t\right))+\frac{1}{2}{f}_{xx}(\xi \left(t\right),t)b^2\left(t\right)]dt+f_x(\xi \left(t\right),t)b\left(t\right)dw\left(t\right).\hfill \end{array}$$

Notice that if $w\left(t\right)$ were continuously differentiable in t, then (by the standard calculus formula for total derivatives) the term $\frac{1}{2}{f}_{xx}{b}^{2}dt$ will not appear.

In Figure 1, we illustrate an elementary simulation of the Brownian motion with a step size ${\Delta }t=0.1$. Such simulation can be performed by using a built-in Matlab function (randn) for representing $N(0,1)$ stochastic variable.

Here, we recall two highly significant results that help us in deriving the main results of this work. The first one referred to M. Ortigueira and J. Machado who established a proper formula to find the exact values of given definite fractional integrals [23], while the other one referred to I. Batiha et al. who provide an approach for the Caputo derivative called the modified three-point fractional formula for approximating Caputo derivative [24].

Definition 3

([23] (Definite Fractional Integral)). The definite fractional integral of the function f of order α is given by:

$$J_a^{\alpha }f\left(x\right)={\int }_a^b{f}^{(-\alpha +1)}\left(x\right).dx={\int }_a^b{D}_{*}^{-\alpha +1}f\left(x\right)dx,$$

(14)

where $-\infty <a<b<\infty$ and $\alpha -1<n\le \alpha$ such that $n\in \mathbb{N}$.

Theorem 3

([24]). Suppose $f\in C^3[a,b]$ and $x_0$, $x_1$, $x_2$ are three distinct points in the interval $[a,b]$ such that $a=x_0<x_1=x_0+h<x_2=x_0+2h=b$ with $h>0$. Then the modified three-point fractional formula for approximating the Caputo first derivative is given by:

$$\begin{array}{cc}\hfill D_{*}^{\alpha }f\left(x\right)=& \frac{x^{2-\alpha }}{h^2{\Gamma }(3-\alpha )}\left(f\left(x_0\right)-2f\left(x_1\right)+f\left(x_2\right)\right)\hfill \\ \hfill & -\frac{x^{1-\alpha }}{2h^2{\Gamma }(2-\alpha )}\left(f\left(x_0\right)(x_1+x_2)-2f\left(x_1\right)(x_0+x_2)+f\left(x_2\right)(x_0+x_1)\right)\hfill \\ \hfill & +{\displaystyle \frac{f^{\left(3\right)}\left(\xi \right)}{6}}\left(\frac{6x^{3-\alpha }}{{\Gamma }(4-\alpha )}-\frac{2(x_0+x_1+x_2)x^{2-\alpha }}{{\Gamma }(3-\alpha )}+\frac{(x_0{x}_1+x_0{x}_2+x_1{x}_2)x^{1-\alpha }}{{\Gamma }(2-\alpha )}\right),\hfill \end{array}$$

(15)

for each $x\in [a,b]$, where $\xi \in (a,b)$.

3. Modified Three-Point Fractional Formula

This parts aims to develop a novel fractional-order version of the classical three-point formula that might be used to approximate integrals. Such formula will be called the modified three-point fractional formula for approximating the Riemann–Liouville fractional integral operator. For this purpose and based on Theorem 3, we can easily deduce the next result that establishes an approximation for the Caputo derivative operator of order $2\alpha$.

Corollary 1.

Under the same assumptions of Theorem 3, the modified three-point fractional formula for approximating the Caputo derivative operator of order $2\alpha$ is given by:

$$\begin{array}{cc}\hfill D_{*}^{2\alpha }f\left(x\right)=& \frac{x^{2-2\alpha }}{h^2{\Gamma }(3-2\alpha )}(f\left(x_0\right)-2f\left(x_1\right)+f\left(x_2\right))\hfill \\ & +\frac{f^{\left(3\right)}\left(\xi \right)}{6}\left(\frac{6x^{3-2\alpha }}{{\Gamma }(4-2\alpha )}-\frac{2(x_0+x_1+x_2)x^{2-2\alpha }}{{\Gamma }(3-2\alpha )}\right).\hfill \end{array}$$

(16)

In light of the previous discussion and based on the generalized Taylor Theorem 1 coupled with the Definite Fractional Integral Definition 3, we establish the next so-called modified three-point fractional formula for approximating Riemann–Liouville integrator.

Theorem 4.

Let $D_{*}^{5\alpha }f\in C^4[a,b]$, where $\alpha =\frac{n}{m}$ such that $n\le m$ with $n,m\in {\mathbb{Z}}^{+}$ and $m=2k-1$ for $k\in {\mathbb{Z}}^{+}$. Suppose $x_0$, $x_1$ and $x_2$ are three distinct points in the interval $[a,b]$ such that $a=x_0<x_1=x_0+h<x_2=x_0+2h=b$ with $h>0$. Then the three-point central fractional formula for approximating Riemann–Liouville integrator is given by:

$$\begin{array}{cc}\hfill J_a^{\alpha }f\left(x\right)=2hf\left(x_1\right)& +\frac{2h^{3\alpha }{x}_1^{2-2\alpha }}{{\Gamma }(3\alpha +1){\Gamma }(3-2\alpha )}(f\left(x_0\right)-2f\left(x_1\right)+f\left(x_2\right))\hfill \\ & +\frac{2h^{3\alpha }}{6{\Gamma }(3\alpha +1)}{f}^{\left(3\right)}\left(\xi \right)\left(\frac{6x_1^{(3-2\alpha )}}{{\Gamma }(4-2\alpha )}-\frac{2(x_0+x_1+x_2)}{{\Gamma }(3-2\alpha )}{x}_1^{(2-2\alpha )}\right)+\frac{2h^{5\alpha }}{{\Gamma }(5\alpha +1)}{D}_{*}^{4\alpha }f\left(\xi \right).\hfill \end{array}$$

(17)

Proof. 

In order to prove this result, we expand first the generalized Taylor Theorem 1 on f about $x=x_1$ to obtain:

$$\begin{array}{cc}\hfill f\left(x\right)=f\left(x_1\right)& +D_{*}^{\alpha }f\left(x_1\right)\frac{{(x-x_1)}^{\alpha }}{{\Gamma }(\alpha +1)}+D_{*}^{2\alpha }f\left(x_1\right)\frac{{(x-x_1)}^{2\alpha }}{{\Gamma }(2\alpha +1)}\hfill \\ & +D_{*}^{3\alpha }f\left(x_1\right)\frac{{(x-x_1)}^{3\alpha }}{{\Gamma }(3\alpha +1)}+D_{*}^{4\alpha }f\left(\xi \right)\frac{{(x-x_1)}^{4\alpha }}{{\Gamma }(4\alpha +1)}.\hfill \end{array}$$

(18)

By applying a $J_a^{\alpha }$ to both sides of the above equality, coupled with using Equation (14) and Definition 3, we obtain:

$$\begin{array}{cc}\hfill J_a^{\alpha }f\left(x\right)=2hf\left(x_1\right)& +\frac{D_{*}^{\alpha }f\left(x_1\right)}{{\Gamma }\left(2\alpha \right)}{\int }_a^b{(x-x_1)}^{2\alpha -1}dx\hfill \\ & +\frac{D_{*}^{2\alpha }f\left(x_1\right)}{{\Gamma }\left(3\alpha \right)}{\int }_a^b{(x-x_1)}^{3\alpha -1}dx\hfill \\ & +\frac{D_{*}^{3\alpha }f\left(x_1\right)}{{\Gamma }\left(4\alpha \right)}{\int }_a^b{(x-x_1)}^{4\alpha -1}dx\hfill \\ & +\frac{D_{*}^{4\alpha }f\left(\xi \right)}{{\Gamma }\left(5\alpha \right)}{\int }_a^b{(x-x_1)}^{5\alpha -1}dx.\hfill \end{array}$$

(19)

This immediately gives:

$$\begin{array}{cc}\hfill J_a^{\alpha }f\left(x\right)=2hf\left(x_1\right)& +\frac{D_{*}^{\alpha }f\left(x_1\right)}{{\Gamma }(2\alpha +1)}\left({\left(h\right)}^{2\alpha }-{(-h)}^{2\alpha }\right)\hfill \\ & +\frac{D_{*}^{2\alpha }f\left(x_1\right)}{{\Gamma }(3\alpha +1)}\left({\left(h\right)}^{3\alpha }-{(-h)}^{3\alpha }\right)\hfill \\ & +\frac{D_{*}^{3\alpha }f\left(x_1\right)}{{\Gamma }(4\alpha +1)}\left({\left(h\right)}^{4\alpha }-{(-h)}^{4\alpha }\right)\hfill \\ & +\frac{D_{*}^{4\alpha }f\left(\xi \right)}{{\Gamma }(5\alpha +1)}\left({\left(h\right)}^{5\alpha }-{(-h)}^{5\alpha }\right).\hfill \end{array}$$

(20)

Observe that we can clearly assert that ${(-h)}^{2\alpha }={\left(h\right)}^{2\alpha }$ and ${(-h)}^{4\alpha }={\left(h\right)}^{4\alpha }$. However, for the other similar terms, we can find where $\alpha$ will be defined by taking ${(-h)}^{3\alpha }=-{\left(h\right)}^{3\alpha }$ and ${(-h)}^{5\alpha }=-{\left(h\right)}^{5\alpha }$. Actually, this is valid for $\alpha =\frac{n}{m}$ such that $n\le m$, where $n,m\in {\mathbb{Z}}^{+}$ with $m=2k-1$ for $k\in {\mathbb{Z}}^{+}$. In other words and without loss of generality, if $\alpha =1$, then we can have ${(-h)}^{2\alpha }={\left(h\right)}^{2\alpha }$, ${(-h)}^{3\alpha }=-{\left(h\right)}^{3\alpha }$, ${(-h)}^{4\alpha }={\left(h\right)}^{4\alpha }$, and ${(-h)}^{5\alpha }=-{\left(h\right)}^{5\alpha }$. Based on this discussion, we can obtain the following equation:

$$\begin{array}{cc}\hfill J_a^{\alpha }f\left(x\right)& =2hf\left(x_1\right)+\frac{2h^{3\alpha }}{{\Gamma }(3\alpha +1)}{D}_{*}^{2\alpha }f\left(x_1\right)+\frac{2h^{5\alpha }}{{\Gamma }(5\alpha +1)}{D}_{*}^{4\alpha }f\left(\xi \right).\hfill \end{array}$$

(21)

Now, by substituting (16) into (21), we obtain:

$$\begin{array}{cc}\hfill J_a^{\alpha }& f\left(x\right)=2hf\left(x_1\right)\hfill \\ & +\frac{2h^{3\alpha }}{{\Gamma }(3\alpha +1)}\left[\frac{x_1^{2-2\alpha }}{{\Gamma }(3-2\alpha )}(f\left(x_0\right)-2f\left(x_1\right)+f\left(x_2\right))+\frac{f^{\left(3\right)}\left(\xi \right)}{6}\left(\frac{6x_1^{(3-2\alpha )}}{{\Gamma }(4-2\alpha )}-\frac{2(x_0+x_1+x_2)}{{\Gamma }(3-2\alpha )}{x}_1^{(2-2\alpha )}\right)\right]\hfill \\ & +\frac{2h^{5\alpha }}{{\Gamma }(5\alpha +1)}{D}_{*}^{4\alpha }f\left(\xi \right).\hfill \end{array}$$

(22)

Simplifying the above equation yields:

$$\begin{array}{cc}\hfill J_a^{\alpha }f\left(x\right)=2hf\left(x_1\right)& +\frac{2h^{3\alpha }{x}_1^{2-2\alpha }}{{\Gamma }(3\alpha +1){\Gamma }(3-2\alpha )}(f\left(x_0\right)-2f\left(x_1\right)+f\left(x_2\right))\hfill \\ & +\frac{2h^{3\alpha }}{6{\Gamma }(3\alpha +1)}{f}^{\left(3\right)}\left(\xi \right)\times \left(\frac{6x_1^{(3-2\alpha )}}{{\Gamma }(4-2\alpha )}-\frac{2(x_0+x_1+x_2)}{{\Gamma }(3-2\alpha )}{x}_1^{(2-2\alpha )}\right)\hfill \\ & +\frac{2h^{5\alpha }}{{\Gamma }(5\alpha +1)}{D}_{*}^{4\alpha }f\left(\xi \right).\hfill \end{array}$$

(23)

Hence, the three-point central fractional formula for approximating Riemann–Liouville integrator is given by:

$$\begin{array}{cc}\hfill J_a^{\alpha }f\left(x\right)& \approx 2hf\left(x_1\right)+\frac{2h^{3\alpha }{x}_1^{2-2\alpha }}{{\Gamma }(3\alpha +1){\Gamma }(3-2\alpha )}(f\left(x_0\right)-2f\left(x_1\right)+f\left(x_2\right)),\hfill \end{array}$$

(24)

which completes the desired result. □

4. Handling FSDE Using the Modified Three-Point Fractional Formula

In this section, we introduce a novel numerical solution for the FSDE by using the modified three-point fractional formula for approximating Riemann–Liouville integrator. To this aim, we reconsider again the FSDE (1) again as follows:

$$D_{*}^{\alpha }X\left(t\right)=f(t,X(t,w))dt+g(t,X\left(t\right))dW\left(t\right),0\le t\le T,$$

(25)

with the initial condition $X\left(0\right)=\eta$. Now, by taking $J_0^{\alpha }$ to the both sides of (25), we obtain:

$$X\left(t\right)=\eta +J_0^{\alpha }f(t,X(t,w))dt+J_0^{\alpha }g(t,X\left(t\right))dW\left(t\right).$$

(26)

Herein, we suppose that $t_0,t_1,t_2,\dots ,t_n$ are distinct points in the interval $[0,T]$ such that $0=t_0<t_1=t_0+h<t_2=t_0+2h<\dots <t_n=t_0+nh=T$, where $h>0$ is the step size of the discretization. Now, by applying (24) into (26), we have then an approximate numerical solution for FSDE (25), which would be:

$$\begin{array}{cc}\hfill X\left(t_i\right)\approx \eta & +2hf(t_i,X\left(t_i\right))+\frac{2h^{3\alpha }{t}_i^{2-2\alpha }}{{\Gamma }(3\alpha +1){\Gamma }(3-2\alpha )}\times [f(t_{i-1},X\left(t_{i-1}\right))-2f(t_i,X\left(t_i\right))+f(t_{i+1},X\left(t_{i+1}\right))]\hfill \\ & +2hg(t_i,X\left(t_i\right))+\frac{2h^{3\alpha }{t}_i^{2-2\alpha }}{{\Gamma }(3\alpha +1){\Gamma }(3-2\alpha )}\times [g(t_{i-1},X\left(t_{i-1}\right))-2g(t_i,X\left(t_i\right))+g(t_{i+1},X\left(t_{i+1}\right))],\hfill \end{array}$$

(27)

for $i=1,2,3,\dots ,n$, where $n=\frac{T}{h}$.

5. Applications

Herein, we intend to test the validation of the approximate numerical solution (27) proposed for the FSDE (25). For this purpose, we list the following examples.

Example 1.

Assume that $f\left(X\right(t\left)\right)=-10X\left(t\right)$ and $g\left(X\right(t\left)\right)=1$ with the initial condition $X\left(0\right)=1$, i.e., we have the following nonlinear FSDE [25]:

$$\begin{array}{c}\hfill D_{*}^{\alpha }X\left(t\right)=-10X\left(t\right)dt+dW\left(t\right),0\le t\le 1,\end{array}$$

(28)

subject to the initial condition $X\left(0\right)=1$. Accordingly, by applying the approximate numerical solution given in (27), we obtain:

$$\begin{array}{cc}\hfill X\left(t_i\right)\approx 1& +2h\left(-10X\left(t_i\right)\right)+\frac{2h^{3\alpha }{t}_i^{2-2\alpha }}{{\Gamma }(3\alpha +1){\Gamma }(3-2\alpha )}\times \left\{\left(-10X\left(t_{i-1}\right)\right)-2\left(-10X\left(t_i\right)\right)+\left(-10X\left(t_{i+1}\right)\right)\right\}\hfill \\ & +2h+\frac{2h^{3\alpha }{t}_i^{2-2\alpha }}{{\Gamma }(3\alpha +1){\Gamma }(3-2\alpha )}\times \left\{\left(1\right)-2\left(1\right)+\left(1\right)\right\},\hfill \end{array}$$

(29)

for $i=1,2,3,\dots ,10$.

To see how the numerical solution (29) looks according to different values of α, we plot Figure 2. In addition, in order to validate our proposed numerical scheme, we plot once again our numerical solution (29) in Figure 3 and compare it with the exact solution and with another numerical solution generated by Euler–Maruyama method. We also plot the absolute error gained from such a comparison in Figure 4.

We can notice that our proposed approximate solution (29) generated by our numerical method is closer to the exact solution of the FSDE (28) than that of Euler–Maruyama’s solution.

Example 2.

Assume that $f\left(X\left(t\right)\right)=\frac{3}{10}X\left(t\right)$ and $g\left(X\left(t\right)\right)=X^{\frac{2}{10}}\left(t\right)$ with the initial condition $X\left(0\right)=1$, i.e., we have the following nonlinear FSDE [25]:

$$\begin{array}{c}\hfill D_{*}^{\alpha }X\left(t\right)=\frac{3}{10}X\left(t\right)dt+X^{\frac{2}{10}}\left(t\right)dW\left(t\right),0\le t\le 1,\end{array}$$

(30)

subject to the initial condition $X\left(0\right)=1$. With the use of the approximate numerical solution given in (27), we obtain:

$$\begin{array}{cc}\hfill X\left(t_i\right)\approx 1& +2h\left(\frac{3}{10}X\left(t_i\right)\right)+\frac{2h^{3\alpha }{t}_i^{2-2\alpha }}{{\Gamma }(3\alpha +1){\Gamma }(3-2\alpha )}\times \left\{\left(\frac{3}{10}X\left(t_{i-1}\right)\right)-2\left(\frac{3}{10}X\left(t_i\right)\right)+\left(\frac{3}{10}X\left(t_{i+1}\right)\right)\right\}\hfill \\ & +2h\left(X^{\frac{2}{10}}\left(t_i\right)\right)+\frac{2h^{3\alpha }{t}_i^{2-2\alpha }}{{\Gamma }(3\alpha +1){\Gamma }(3-2\alpha )}\times \left\{\left(X^{\frac{2}{10}}\left(t_{i-1}\right)\right)-2\left(X^{\frac{2}{10}}\left(t_i\right)\right)+\left(X^{\frac{2}{10}}\left(t_{i+1}\right)\right)\right\},\hfill \end{array}$$

(31)

for $i=1,2,3,\dots ,10$.

In this regard, the numerical solution (31) according to different values of α can be seen in Figure 5. Moreover, we plot once again our numerical solution (31) in Figure 6 and compare it with the exact solution and with another numerical solution generated by Euler–Maruyama method. Furthermore, we also plot the absolute error gained from such a comparison in Figure 7.

We can see here that our proposed approximate solution (31) generated by our numerical method is closer to the exact solution of the FSDE (30) than that of Euler–Maruyama’s solution.

Example 3.

Assume here $f\left(X\left(t\right)\right)=\frac{2}{5}{X}^{\frac{3}{5}}\left(t\right)+5X^{\frac{4}{5}}\left(t\right)$ and $g\left(X\left(t\right)\right)=X^{\frac{4}{5}}\left(t\right)$ with the initial condition $X\left(0\right)=10$, i.e., we have the following nonlinear FSDE [26]:

$$\begin{array}{cc}\hfill D_{*}^{\alpha }X\left(t\right)& =\left(\frac{2}{5}{X}^{\frac{3}{5}}\left(t\right)+5X^{\frac{4}{5}}\left(t\right)\right)dt+X^{\frac{4}{5}}\left(t\right)dW\left(t\right),0\le t\le 1,\hfill \end{array}$$

(32)

subject to the initial condition $X\left(0\right)=10$. Applying the approximate numerical solution given in (27) yields:

$$\begin{array}{cc}\hfill X\left(t_i\right)\approx 10& +2h\left(\frac{2}{5}{X}^{\frac{3}{5}}\left(t_i\right)+5X^{\frac{4}{5}}\left(t_i\right)\right)+\frac{2h^{3\alpha }{t}_i^{2-2\alpha }}{{\Gamma }(3\alpha +1){\Gamma }(3-2\alpha )}\hfill \\ \hfill & \times \left\{\left(\frac{2}{5}{X}^{\frac{3}{5}}\left(t_{i-1}\right)+5X^{\frac{4}{5}}\left(t_{i-1}\right)\right)-2\left(\frac{2}{5}{X}^{\frac{3}{5}}\left(t_i\right)+5X^{\frac{4}{5}}\left(t_i\right)\right)+\left(\frac{2}{5}{X}^{\frac{3}{5}}\left(t_{i+1}\right)+5X^{\frac{4}{5}}\left(t_{i+1}\right)\right)\right\}\hfill \\ & +2h\left(X^{\frac{4}{5}}\left(t_i\right)\right)+\frac{2h^{3\alpha }{t}_i^{2-2\alpha }}{{\Gamma }(3\alpha +1){\Gamma }(3-2\alpha )}\times \left\{\left(X^{\frac{4}{5}}\left(t_{i-1}\right)\right)-2\left(X^{\frac{4}{5}}\left(t_i\right)\right)+\left(X^{\frac{4}{5}}\left(t_{i+1}\right)\right)\right\},\hfill \end{array}$$

(33)

for $i=1,2,3,\dots ,10$.

Herein, Figure 8 depicts the numerical solution (33) according to different values of α. Moreover, Figure 9 demonstrates the validity of our proposed numerical scheme by making a numerical comparison between our numerical solution (33) and the exact solution coupled with Euler–Maruyama’s solution. In the same regard, we also plot the absolute error gained from such a comparison in Figure 10 for completeness.

Obviously, one might undoubtedly observe that our proposed approximate solution (33) generated by our numerical method is closer to the exact solution of the FSDE (32) than that of Euler–Maruyama’s solution.

It should be noted that because the stochastic differential equation consists of a deterministic differential equation coupled with a random term, then the approximate solution for such an equation will be slightly changed from time to time. This is because the random term in that equation is usually programmed by using a built-in Matlab function (called randn) for representing the $N(0,1)$ stochastic variable. This would generate the Brownian motion with a step size of $h=0.1$, which would affect the approximate solution each time the prepared code is run. So, the approximate solution generated by both methods (our method and the Euler–Maruyama method) are not consistent. Despite these changes, our proposed scheme stills presents better results than that of the Euler–Maruyama scheme.

6. Conclusions

In this paper, a numerical solution of the FSDE has been proposed by using a new modification of the classical three-point formula called the modified three-point fractional formula for approximating Riemann–Liouville integrator. The generalized Taylor theorem with the recent definition where the definite fractional integral have been used to derive this formula. Based on several numerical experiments, one can clearly observe that the numerical solution generated by our scheme is closer to the exact solution of the FSDE than that of Euler–Maruyama’s solution.