A pseudo-spectral scheme for variable order fractional stochastic Volterra integro-differential equations

1.   Introduction

Fractional calculus ^[1] generalizes the standard calculus and has been applied in many fields of engineering and science ^[2]. Often, finding the exact solutions of fractional order differential equations is not possible, and accurate numerical methods have to be developed.

Stochastic differential equations (SDE) are utilised to simulate a variety of phenomena, including volatile stock prices or systems sensitive to thermal variations. SDEs comprise a variable that corresponds to random white noise, computed as the derivative of a Wiener process or Brownian motion. Other forms of random behaviour, such as jump processes, are also feasible. Stochastic differential equations are conjugate to random differential equations. The solution of stochastic differential equations (SDEs) is a stochastic process, and many SDEs can not be solved analytically. Stochastic integro-differential equations (SIDEs) revealed increasing importance in the last decades, since they can describe accurately many real-world phenomena ^[3,4]. Different numerical techniques for obtaining SIDEs solutions were developed, as those based on finite differences ^[5,6,7], finite elements ^[8] and spectral Galerkin methods ^{[9,10,11,12,13]} and others ^[14,15].

Bernoulli collocation ^[16] and Bernoulli wavelet ^[17] algorithms are used to solve stochastic Ito-Volterra integral equations. Also, Multi-dimensional stochastic integral equations have been solved using bicubic B-spline functions ^[18], Quintic B-spline collocation ^[19], bilinear spline interpolation ^[20], cubic B-spline and bicubic B-spline collocation ^[21], and radial basis functions ^[22]. Moreover, meshless discrete collocation ^[23], Cubic B-spline ^[24], moving least squares-spectral collocation ^[25] and spectral collocation ^[26] techniques have been used to solve stochastic integro-differential equations. The well-posedness of stochastic fractional integro-differential equations is introduced in detail by Dai et al. ^[27,28]. Because the kernel of variable-order operators has a variable exponent, analytical solutions to the associated equations are more difficult to acquire than for constant-order operators. As a result, the numerical solution offers a more realistic method of investigating variable-order evolutional equations.

Spectral algorithms ^{[29,30,31,32,33]} are powerful methods for solving fractional differential equations. There are four different types of spectral methods: Galerkin ^{[34,35,36,37,38]}, collocation ^[39,40,41], tau, and Petrov Galerkin. The idea is to approximate the numerical solution by means of truncated series of basis functions, where the coefficients are chosen so that the error between the exact and approximated solutions is minimized. In the spectral collocation method ^[42,43,44], the approximate solution satisfies the original equation at the collocation points closely, which means that the residuals are zero.

The spectral collocation technique is very useful since it can estimate accurately the solutions of a wide range of equations, including differential, integral, integro-differential and fractional differential equations, and optimal control and variational problems. The collocation technique has been widely applied due to its advantages over other methods, namely in providing highly precise solutions and in yielding exponential rates of convergence. The collocation method, on the other hand, has gained popularity in recent decades for dealing with fractional problems. In this paper we use the shifted fractional order Legendre-Gauss-Radau collocation (SFL-GR-C) method to solve variable order fractional stochastic Volterra integro-differential equations (VOFSV-IDEs). For the independent variable we adopt as interpolation points the shifted fractional order Legendre-Gauss-Radau (SFL-GR) nodes, and we express the solution of the VOFSV-IDE by means of a series of shifted fractional order Legendre orthogonal functions (SFOLOF). Thereon, we estimate the residuals at the SFJ-GR quadrature points. After an adequate mapping, we use the shifted Legendre Gauss-Lobatto quadrature to handle the integral terms. The Brownian motion is treated by means of Lagrange interpolation. Therefore, we obtain an algebraic system that is solved with a suitable method. The precision of the new technique is confirmed with various numerical problems.

The work is structured as follows. Section 2 introduces the preliminary bases and some useful definitions, corollaries and theorems. In Section 3, we apply the developed numerical method to solve FOSV-IDEs with initial condition. In Section 4, we illustrate the benefits of the method by means of numerical examples. Finally, in Section 5, we address the conclusions.

2.   Notations and preliminaries

2.1.   Fractional calculus

There are various definitions of fractional integration and differentiation, but the most reliable and extensively used are the famous Caputo and Riemann-Liouville formulations. There exist fractional derivatives of fixed, variable, distributed, and tempered orders for both definitions. In this paper, we will look at VOFSV-IDEs with the Caputo formulation.

Definition 2.1. The Caputo fractional derivative of variable order $\rho(x)$ is:

where $s = \lceil\rho(x)\rceil$ and the Gamma function $\Gamma(.)$ is given by

The Legendre polynomials $\Omega_{m}(t),$ $m = 0, 1\ldots,$ obey to the Rodrigues formula:

Moreover, $\Omega_{m}(t)$ corresponds to a Legendre polynomial ^[45,46] of degree $m$ and we can get the $p$th derivative of $\Omega_{m}(t)$ as:

where

We obtain the orthogonality by:

where $\omega(t) = 1$, $h_{m} = \frac{2}{2m+1}.$

The integrals given above were evaluated efficiently using the Legendre-Gauss-Lobatto quadrature. For $\psi\in S_{2N-1}$ $[-1, 1]$, we have:

Let the discrete inner product be given as:

In the case of smooth solutions, the employment of classical polynomials in spectral techniques, such as Legendre and Jacobi, suffices to find approximate solutions with great accuracy. On the contrary, in the case of non-smooth solutions, this precision degrades, prompting us to utilize polynomials of fractional orders to prevent this issue, see ^[26,47,48] for more details.

Definition 2.2. Let SFOLOF, the inferred function from the Legendre polynomial, be provided by ^[26,47,48]

Theorem 2.1. For $\varpi_{\mathcal{L}, f}^{(\varepsilon)}(x) = x^{\varepsilon-1}$, a complete $L_{\varpi_{\mathcal{L}, f}^{(\varepsilon)}}^{2}[0, \mathcal{L}]$-orthogonal system is obtained ^[26,47,48]

where $h^{(\varepsilon)}_{\mathcal{L}, k} = \frac{\mathcal{L}^{\varepsilon}}{2(2k+1)\varepsilon}$.

Corollary 2.2. Let $\Omega_{\mathcal{M}} = span\{\Omega_{\mathcal{L}, r}^{(\varepsilon)}:0\leq r\leq \mathcal{M}\}$, be the fractional-polynomial space of finite dimensions. Along with the orthogonal characteristic $(2.8)$, the function $\Theta(\xi)\in L_{\mathcal{W}_f^{(\varepsilon)}}^{2}[0, \mathcal{L}]$ may be extracted as

Theorem 2.3. The relevant nodes and Christoffel numbers of the shifted fractional Legendre-Gauss (Gauss-Radau or Gauss-Lobatto) interpolation in the interval $[0, \mathcal{L}]$ may be obtained from ^[26,47,48]

where $x_{\mathcal{K}, s}, \mathit{{and}}\, \varpi_{\mathcal{K}, s}, \, 0\leqslant s\leqslant \mathcal{M},$are the nodes and Christoffel numbers of the standard Legendre-Gauss (Gauss-Radau or Gauss-Lobatto) interpolation in the interval $[-1, 1]$.

Theorem 2.4. Let us consider that $D^{k\varepsilon}\mathcal{X}(x)\in C[0, \mathcal{L}]$ $\forall k = 0, 1, ..., \mathcal{M}$. Given the best approximation $\mathcal{X}_{\mathcal{M}, \varepsilon}(x)$ to $\mathcal{X}(x)$ from $\Omega_{\mathcal{M}}$, the error bound is ^[26]

In addition, the exponential error convergence of fractional Legendre (as a special case of fractional Jacobi) has already been gained, see ^[47,48]. In the case of a smooth solution, the errors drop exponentially as $N\rightarrow \infty$.

2.2.   Stochastic calculus

Scalar Brownian motion is a stochastic process $\omega(t)$ that varies continuously on $t \in [0, \mathcal{L}]$ and fulfills the conditions:

1) $\omega(t) = 0$, with probability 1.

2) For $0\leq s < t\leq \mathcal{L}$ the random variable given by the increment $\omega(t)-\omega(s)$ is normally distributed with zero mean and variance $t-s$, which means $\omega(t)-\omega(s)\sim \sqrt{t-s} \quad \mathcal{N}(0, 1),$ where $\mathcal{N}(0, 1)$ denotes a normally distributed random variable with zero mean and unit variance.

3) For $0\leq s < t < u < v\leq \mathcal{L},$ $\omega(t)-\omega(s)$ and $\omega(v)-\omega(u)$ are independent.

A stochastic integral is the integral of some function $\psi(t)$ over some interval $[0, \mathcal{L}]$, but with respect to a Brownian motion $\omega(t)$ as $\int_{0}^{\mathcal{L}} \psi(t) d\omega_t.$

Definition 2.3. (The Itô integral) Let $\psi \in \upsilon(s_1, s_2)$, then the Itô integral of $\psi$ is defined by

where $\psi_{\mathcal{M}}$ is a sequence of elementary functions such that

Property 2.5. (Integration by parts) Suppose that $\psi(s, \tau) = \psi(s)$ depends only on $s$, and $\psi$ is continuous and bounded in $[0, t]$. Then,

Property 2.6. (The Itô isometry) Let $\psi\in \nu(s_1, s_2)$, then

3.   VOFSV-IDEs with initial condition

We solve the VOFSV-IDEs:

with initial condition:

where $\mathcal{X}(x), \, \mathcal{H}(x, \mathcal{X}), \, \mathcal{Q}_{1}(\zeta, x)$ and $\mathcal{Q}_{2}(\zeta, x)$, with $x\in[0, \mathcal{L}]$, denote stochastic processes defined on the probability space $(\Omega, \digamma, P)$, and $\mathcal{X}$ is unknown. Moreover, $\int\limits_{0}^{x}\mathcal{Q}_{2}(\zeta, x)\mathcal{X}(\zeta)d\omega(\zeta)$ is the Itô integral. Therefore, using integration by parts, Eq (3.1) yields:

with

where

The solution of Eq (3.2) can be approximated by:

Despite the location of the nodes is optional, we take $x^{\varepsilon}_{\mathcal{L}, \mathcal{Q}, j}$ as SFL-GR nodes.

The $\Omega^{(\varepsilon)}_{\mathcal{L}, j}(x)$ is given analytically by:

where

By means of Eq (2.1), we find:

and, thus, we have:

Therefore, we get:

Adopting $\zeta = \frac{x}{\mathcal{L}}\lambda$, the integrals in (3.3) are converted into:

Using the quadrature, we have:

where $\lambda^{\mathcal{L}, \mathcal{R}}_{r}$ are the shifted Legendre-Gauss-Lobatto nodes. As a result, the integral terms can be written as:

The residual of (3.2) can thus be computed as:

which is complied to be zero at $\mathcal{M}$ points:

where $r = 1, 2, \ldots, \mathcal{M}$. Hence, for $\mathcal{M}+1$ unknowns, we have $\mathcal{M}$ algebraic equations. One additional equation may be constructed from the beginning condition (3.2) as follows:

Finally, from Eqs (3.11) and (3.12), we have a system of algebraic equations.

Consider the discretized Brownian motion for computational purposes, in which $\omega(x)$ is given at $x$ discrete values and $\omega(x)$ is evaluated using Lagrange interpolation. We have:

where $\omega_0 = 0$ with probability 1, $\omega_i = \omega(x^{\varepsilon}_{\mathcal{L}, \mathcal{Q}, i})$, $\triangle x^{\varepsilon}_{\mathcal{L}, \mathcal{Q}, i} = x^{\varepsilon}_{\mathcal{L}, \mathcal{Q}, i}-x^{\varepsilon}_{\mathcal{L}, \mathcal{Q}, i-1}$, and $\triangle \omega_i$ is a random variable such that $\sqrt{\triangle x^{\varepsilon}_{\mathcal{L}, \mathcal{Q}, i}} \mathcal{N}(0, 1)$.

4.   Numerical results and assessment

We present numerical results obtained with the new algorithm, illustrating its effectiveness and high accuracy. Therefore, three problems are solved. The algorthim code is run via MATHEMATICA version 12.2.

4.1.   Problem I

We solve the VOFSV-IDEs:

where, for $\sigma = 0$, we get $\mathcal{X}(x) = x^3$.

Table 1 shows the maximum absolute error between the approximate and exact solutions, given $\sigma = 0$ and $\rho(x) = x^2$. We verify that the exact solution is obtained when $\varepsilon\mathcal{M} = p$, with $p$ denoting the power exponent in the exact solution. Table 2 lists the absolute error for $\sigma = \frac{1}{2}$ and $\rho(x) = x^2$, while Figure 1 depicts the numerical and exact solutions for $\sigma = 0$, $\varepsilon = \frac{1}{5}$ and $\mathcal{M} = \mathcal{Q} = 12$. The variation of the absolute errors are shown in Figure 2. Additionally, we can see the error degradation in Figure 3. When adopting $\sigma = \frac{1}{2}, \, \rho(x) = x^2$, we get the approximate solution of Problem I as:

4.2.   Problem II

The following VOFSV-IDEs are solved:

where the solution is $\mathcal{X}(x) = x^{5}$, for $\sigma = 0$ and $\rho(x) = \frac{1}{10} \left(1-x^2\right)$.

Table 3 summarizes the maximum absolute error for $\sigma = 0$ and $\rho(x) = \frac{1}{10} \left(1-x^2\right)$, where we can see that for $\varepsilon\mathcal{M} = p$ the exact solution is obtained. Table 4 shows that accurate results are gotten when $\sigma = \frac{1}{5}$ and $\rho(x) = \frac{1}{10} \left(1-x^2\right)$. Figure 4 illustrates the estimated and exact solutions for $\sigma = 0$, $\varepsilon = \frac{1}{5}$ and $\mathcal{M} = \mathcal{Q} = 20$, showing an excellent matching between them. The variation of the absolute error is shown in Figure 5, while the numerical solution is plotted in Figure 6 for $\sigma = \frac{1}{5}$, $\varepsilon = 1$ and $\mathcal{M} = \mathcal{Q} = 25$. Taking $\sigma = \frac{1}{5}, \, \rho(x) = \frac{1}{10} \left(1-x^2\right)$, the numerical solution of Problem II is:

4.3.   Problem III

Here, we treat the VOFSV-IDEs:

With $\sigma = 0$ and $\rho(x) = \frac{1}{20} \left(1-x^2\right)$, we have $\mathcal{X}(x) = \sqrt{x}$.

Table 5 compiles the maximum absolute error for $\sigma = 0$ and $\rho(x) = \frac{1}{20} \left(1-x^2\right)$, while Table 6 lists the errors for $\sigma = \frac{1}{5}$ and $\rho(x) = \frac{1}{20} \left(1-x^2\right)$. Figures 7 and 8 depict the exact and numerical solutions and the variation of the absolute error, respectively, for $\sigma = 0$, $\varepsilon = \frac{1}{5}$ and $\mathcal{M} = \mathcal{Q} = 20$. Furthermore, in Figure 9, the numerical solution of Problem III is depicted for $\sigma = \frac{1}{5}$, $\varepsilon = 1$ and $\mathcal{M} = \mathcal{Q} = 20$. With $\sigma = \frac{1}{5}, \, \rho(x) = \frac{1}{20} \left(1-x^2\right),$ the numerical solution of Problem III is:

5.   Conclusions

A collocation strategy was proposed to solve VOFSV-IDEs. Numerical examples proved the new method's accuracy and applicability. Indeed, the results showed that good precision is achieved even with a small number of points. Other methods need a larger number of points to achieve identical results and, thus, involve a higher computational burden. For example, the Euler-Maruyama approximation arrives at $10^{-4}$ as the best solution with step size $\frac{1}{512}$. It is important to note that the proposed approach also thrives in situations with non-smooth solutions, where the accuracy of many other methods is affected.

Acknowledgments

The authors extend their appreciation to the Deputyship for Research & Innovation, Ministry of Education in Saudi Arabia for funding this research work through project number IFP-IMSIU202109.

Conflict of interest

The authors declare no conflicts of interest.