Quarter-Sweep Preconditioned Relaxation Method, Algorithm and Efficiency Analysis for Fractional Mathematical Equation

Abstract

Research into the recent developments for solving fractional mathematical equations requires accurate and efficient numerical methods. Although many numerical methods based on Caputo’s fractional derivative have been proposed to solve fractional mathematical equations, the efficiency of obtaining solutions using these methods when dealing with a large matrix requires further study. The matrix size influences the accuracy of the solution. Therefore, this paper proposes a quarter-sweep finite difference scheme with a preconditioned relaxation-based approximation to efficiently solve a large matrix, which is based on the establishment of a linear system for a fractional mathematical equation. The paper presents the formulation of the quarter-sweep finite difference scheme that is used to approximate the selected fractional mathematical equation. Then, the derivation of a preconditioned relaxation method based on a quarter-sweep scheme is discussed. The design of a C++ algorithm of the proposed quarter-sweep preconditioned relaxation method is shown and, finally, efficiency analysis comparing the proposed method with several tested methods is presented. The contributions of this paper are the presentation of a new preconditioned matrix to restructure the developed linear system, and the derivation of an efficient preconditioned relaxation iterative method for solving a fractional mathematical equation. By simulating the solutions of time-fractional diffusion problems with the proposed numerical method, the study found that computing solutions using the quarter-sweep preconditioned relaxation method is more efficient than using the tested methods. The proposed numerical method is able to solve the selected problems with fewer iterations and a faster execution time than the tested existing methods. The efficiency of the methods was evaluated using different matrix sizes. Thus, the combination of a quarter-sweep finite difference method, Caputo’s time-fractional derivative, and the preconditioned successive over-relaxation method showed good potential for solving different types of fractional mathematical equations, and provides a future direction for this field of research.

1. Introduction

Research into the recent developments for solving fractional mathematical equations requires accurate and efficient numerical methods. The accuracy and efficiency of the numerical method used to solve a fractional mathematical equation determine the precision when interpreting the behaviours of physical phenomena. Accurate and efficient methods enable researchers to explore and solve a greater number of complex fractional mathematical models. Many numerical methods based on Caputo’s fractional derivative have been proposed to approximate the solution of the modelled fractional mathematical equations. For instance, ref. [1] proposed two numerical methods, the non-standard finite difference method and the generalized Euler method, for solving a mathematical model of drug resistance during treatment for human immunodeficiency virus, using the Caputo approach. Subsequently, ref. [2] studied the conversion of an existing tobacco smoking model into a fractional-order Caputo model, and used the combination of the fourth-order Runge–Kutta and Adams–Bashforth–Moulton methods to solve the developed model. In another article, ref. [3] developed a time–space fractional hyperbolic bioheat transfer model for non-Fourier bioheat transfer in living biological tissues during laser irradiation using Caputo’s definition. They applied an L1 finite difference approximation to Caputo’s time-fractional derivative and a central difference approximation to Riesz’s space-fractional derivative for the solution of the model. In addition, ref. [4] investigated Caputo’s fractional Maxwell model of unsteady fluid flow and heat transfer of the natural convection of a viscoelastic non-Newtonian fluid using a finite difference method.

Of the many existing numerical methods, the finite difference method, widely known as FDM, is one of the most popular numerical methods used to solve fractional mathematical equations using the Caputo approach. From our review of several applications of FDM to solve Caputo’s time-fractional mathematical equations, ref. [5] utilized the implicit FDM to solve the time-fractional diffusion equation with a time-invariant type variable. Subsequently, ref. [6] developed a fourth-order FDM to solve a time-fractional diffusion equation after transforming the fractional mathematical equation into a Volterra integro-differential equation via a Laplace transform. In addition, ref. [7] proposed a weighted FDM to effectively solve a system of variable-order time-fractional two-dimensional Burgers’ equations. FDM can also be combined with spline approximation to solve the time-fractional stochastic advection equation [8]. These studies are examples of the numerous works published by global researchers. The abundance of related research in the literature has motivated numerous researchers to propose more accurate and efficient numerical methods for solving fractional mathematical equations.

In our view, the efficiency of numerical methods used to obtain the solution of a fractional mathematical equation when dealing with a large matrix requires further study because the accuracy of the solution is influenced by the matrix size. Therefore, this paper proposes a quarter-sweep finite difference scheme with a preconditioned relaxation-based approximation to efficiently solve a large matrix, which results from establishing a linear system for a fractional mathematical equation. The quarter-sweep finite difference scheme is suitable for solving a fractional mathematical equation because, among other reasons, the derivation of the scheme is similar to that of an implicit FDM, which makes it compatible with Caputo’s time-fractional approximation. The quarter-sweep scheme is derived by modifying an implicit FDM by skipping three grid points for two consecutive unknown grid points. The values of the unknown grid points can be computed using the iterative method until they converge. Then, the remaining grid points resulting from the skipping procedure are computed directly using the approximation function [9,10,11]. Moreover, due to the combination of iterative and direct computation, the quarter-sweep scheme has a lower computational burden compared to the implicit and high-order FDM. This paper introduces a preconditioned relaxation method to improve the convergence rate of the iterations, thereby resulting in a highly efficient numerical method. The contributions of the paper are the presentation of a new preconditioned matrix to restructure the developed linear system and the derivation of an efficient preconditioned relaxation iterative method for solving a fractional mathematical equation.

To investigate the efficacy of this new numerical method, a well-known fractional mathematical equation, namely, the fractional diffusion equation (FDE), was selected as a mathematical problem. FDE is a fractional mathematical equation that has been successfully used in the mathematical modelling of physical phenomena. Due to the limitation of the classic diffusion equation in the modelling of the anomalous diffusion process, the theory and application of fractional diffusion have evolved. Hence, FDE has been used to successfully model shipping water events [12], image denoising [13], dye-sensitized solar cells [14], diffusion of soluble substances [15], groundwater pollution [16], option pricing and risk calculation [17], and signal smoothing performance analysis [18]. A better understanding of FDE models can be achieved using an efficient numerical method called the quarter-sweep preconditioned relaxation method. It should be noted that this paper extends the works from [19] and [20], which implemented the standard implicit FDM with preconditioned relaxation and the half-sweep difference scheme with preconditioned relaxation, respectively, for solving the time FDE. The work in [19] showed the improvement in the number of iterations and execution time after implementing a preconditioned successive over-relaxation with implicit FDM using the Caputo approach to solve the time FDE. Another article [20] applied a computational complexity reduction technique called the half-sweep iteration to successfully reduce the computational cost using preconditioned successive over-relaxation, and eventually improve the efficiency of Caputo’s finite difference scheme. This paper extends these works by employing the quarter-sweep FDM to reduce the computational cost of solving the FDE using a preconditioned successive over-relaxation when dealing with a large matrix.

This paper is organized as follows: Section 2 presents Caputo’s fractional operator, one of the important definitions in the fractional derivative theory. In Section 3, an approximation to a general time FDE is formulated using Caputo’s fractional operator and the quarter-sweep numerical discretization procedure. Section 4 explains the stability of the approximation to the considered mathematical equation. Section 5 discusses the convergence of the quarter-sweep finite difference approximation in the Caputo sense. In Section 6, the concept of preconditioned successive over-relaxation (PSOR) used to formulate the quarter-sweep PSOR (QSPSOR) iterative method is discussed. Section 7 shows the implementation and application of an algorithm written in C++ programming language for the numerical experiment and simulation with QSPSOR, and provides a comparison with selected existing methods. Finally, a conclusion is given in Section 8.

2. Preliminary

Formulating a quarter-sweep approximation to a time FDE requires good knowledge of fractional calculus theory and several useful operators. In the preliminary stage of this paper, Caputo’s fractional operator is used to approximate the time-fractional derivative term. The definition of Caputo’s derivative with a fractional order is as follows:

Definition 1

([21]).Let$m-1<\alpha <m$, where$m$is an element of the set of natural numbers, and$\alpha$is a real number, and a function$f\left(x\right)$such that$D_0^{\alpha }f\left(x\right)$exists. Caputo’s fractional operator can be defined in the form of:

$$D_0^{\alpha }f\left(x\right)={{\displaystyle \int }}_0^x\frac{f^{\left(m\right)}\left(r\right)}{{\left(x-r\right)}^{\alpha -\left(m-1\right)}} dr.$$

(1)

Using Definition 1, a time-fractional derivative term in a differential equation can be approximated by:

$$D_t^{\alpha }z\left(x,t\right)=\frac{1}{\Gamma \left(1-\alpha \right)}{{\displaystyle \int }}_0^t\frac{\partial z\left(x,\tau \right)}{\partial t}{\left(t-\tau \right)}^{-\alpha } d\tau ,$$

(2)

where $0<\alpha <1$ and $z\left(x,t\right)>0$.

The next section discusses the formulation of a quarter-sweep finite difference approximation to the time FDE, which uses Equation (2) to approximate the time-fractional derivative.

3. Mixed Caputo Fractional Operator and Quarter-Sweep Discretization Scheme

To demonstrate the formulation of the mixed Caputo fractional operator and quarter-sweep discretization scheme, let us consider a well-known model of the time FDE, namely:

$$\frac{{\partial }^{\alpha }z\left(x,t\right)}{\partial t^{\alpha }}=a\left(x\right)\frac{{\partial }^{2}z\left(x,t\right)}{\partial x^2}+b\left(x\right)\frac{\partial z\left(x,t\right)}{\partial x}+c\left(x\right)z\left(x,t\right),$$

(3)

where $a\left(x\right), b\left(x\right),$ and $c\left(x\right)$ are either predefined functions, constants, or mixed.

A suitable boundary condition must be imposed on Equation (3) before the solution can be derived using an FDM. In this paper, a Dirichlet-type boundary condition is used to construct a solution domain. Furthermore, the solution domain is restricted to the usual finite domain $0\le x\le \gamma$ because the aim of this paper is to investigate of the efficacy of a new iterative method developed from the use of Equation (2) and a quarter-sweep FDM. Based on Equation (3), Caputo’s fractional approximation to the fractional order in time is illustrated as follows:

$$\frac{{\partial }^{\alpha }z\left(x_i,t_n\right)}{\partial t^{\alpha }}=\frac{1}{\Gamma \left(1-\alpha \right)}{{\displaystyle \int }}_0^{t_n}\frac{\partial z\left(x_i,\tau \right)}{\partial t}{\left(t_n-\tau \right)}^{-\alpha } d\tau .$$

(4)

Equation (4) can be derived further and simplified to:

$$\begin{gathered} \frac{{\partial }^{\alpha }z\left(x_i,t_n\right)}{\partial t^{\alpha }}=\frac{1}{\Gamma \left(1-\alpha \right)\left(1-\alpha \right)k^{\alpha }}{\displaystyle \sum }_{j=1}^n\left(Z_{i,j}-Z_{i,j-1}\right)\left[{(n-\left(j-1\right))}^{1-\alpha }-{(n-j)}^{1-\alpha }\right] \\ +\frac{1}{\Gamma \left(1-\alpha \right)\left(1-\alpha \right)}{\displaystyle \sum }_{j=1}^n\left[{(n-\left(j-1\right))}^{1-\alpha }-{(n-j)}^{1-\alpha }\right]O\left(k^{2-\alpha }\right). \end{gathered}$$

(5)

Next, by defining:

$${\sigma }_{\alpha ,k}=\frac{1}{\Gamma \left(1-\alpha \right)\left(1-\alpha \right)k^{\alpha }},$$

(6)

where $k$ is a real number, and:

$${\omega }_{\alpha ,j}=j^{1-\alpha }-{(j-1)}^{1-\alpha },$$

(7)

a discrete approximation to the fractional order in time becomes:

$$\frac{{\partial }^{\alpha }z\left(x_i,t_n\right)}{\partial t^{\alpha }}={\sigma }_{\alpha ,k}{\displaystyle \sum }_{j=0}^n{\omega }_{\alpha ,j}\left(Z_{i,n-\left(j-1\right)}-Z_{i,n-j}\right)+\frac{1}{\Gamma \left(1-\alpha \right)\left(1-\alpha \right)}{n}^{1-\alpha }O\left(k^{2-\alpha }\right),$$

(8)

and, finally:

$$\frac{{\partial }^{\alpha }z\left(x_i,t_n\right)}{\partial t^{\alpha }}={\sigma }_{\alpha ,k}{\displaystyle \sum }_{j=0}^n{\omega }_{\alpha ,j}\left(Z_{i,n-\left(j-1\right)}-Z_{i,n-j}\right)+O\left(k\right).$$

(9)

To discretize the right-hand side part of Equation (3) using a quarter-sweep implicit scheme of FDM, let the solution domain of the problem be partitioned uniformly. Let some integers $M>0$ and $N>0$ denote the grid sizes in space and time, respectively. Similar to the standard finite difference framework, we use $h=\gamma /M$ and $k=T/N$, respectively. Then, a uniform mesh of the solution can be formed using the distributed points in the space $\left[0, \gamma \right]$ as $x_i=ih,i=1,2,\dots ,M-1$, and the distributed points in the time $\left[0, T\right]$ are labelled $t_n=nk,n=1,2,\dots ,N$. All of the unknown values of the solution function $z\left(x_i,t_n\right)$ located at these grid points are represented by $Z_{i,n}$.

Using Equation (9), in combination with a quarter-sweep implicit scheme, the quarter-sweep Caputo fractional approximation of Equation (3) to the point centered at $z\left(x_i,t_n\right)=z\left(ih,nk\right)$ can be formulated and written as:

$${\sigma }_{\alpha ,k}{\displaystyle \sum }_{j=0}^n {\omega }_{\alpha ,j}\left(Z_{i,n-\left(j-1\right)}-Z_{i,n-j}\right)+O\left(k\right)=\frac{a_i}{16h^2}\left(Z_{i-4,n}-2Z_{i,n}+Z_{i+4,n}\right)+\frac{b_i}{8h}\left(Z_{i+4,n}-Z_{i-4,n}\right)+c_i{Z}_{i,n}+O\left(h^2\right),$$

(10)

where $i=4,6\dots ,m-4$.

The approximation structure presented in Equation (10) is similar to a standard implicit approximation equation, with the exception that the distance between any two points is quadrupled. Equation (10) is also consistent with the first-order accuracy in time and second-order accuracy in space. When Equation (10) is written in a particular form, for instance, at time level $n\ge 2$, one can obtain:

$${\sigma }_{\alpha ,k}{\displaystyle \sum }_{j=2}^n {\omega }_{\alpha ,j}\left(Z_{i,n-\left(j-1\right)}-Z_{i,n-j}\right)=a_i^{\ast }{Z}_{i-4,n}+b_i^{\ast }{Z}_{i,n}+c_i^{\ast }{Z}_{i+4,n},$$

(11)

where:

$$a_i^{\ast }=\frac{a_i}{16h^2}-\frac{b_i}{8h}, b_i^{\ast }=c_i-\frac{a_i}{8h^2}, c_i^{\ast }=\frac{a_i}{16h^2}+\frac{b_i}{8h}.$$

(12)

Moreover, for $n=1$, one can obtain:

$$F_{i,1}=-a_i^{\ast }{Z}_{i-4,1}+b_i^{\ast \ast }{Z}_{i,1}-c_i^{\ast }{Z}_{i+4,1},$$

(13)

where:

$$b_i^{\ast \ast }={\sigma }_{\alpha ,k}-b_i^{\ast },F_{i,1}={\sigma }_{\alpha ,k}{Z}_{i,0}.$$

(14)

Furthermore, using Equation (14), a corresponding linear system can be constructed properly in the form of a matrix as:

$$A\tilde{Z}=\tilde{F},$$

(15)

where:

$$A=\left[\begin{array}{ccc}{b}_4^{\ast \ast }& -c_4^{\ast }& \\ -a_8^{\ast }& b_8^{\ast \ast }& -c_8^{\ast }\\ \begin{array}{c}\\ \\ \begin{array}{c}\\ \end{array}\end{array}& \begin{array}{c}-a_{12}^{\ast }\\ \\ \begin{array}{c}\\ \end{array}\end{array}& \begin{array}{c}{b}_{12}^{\ast \ast }\\ \ddots \\ \begin{array}{c}\\ \end{array}\end{array}\end{array} \begin{array}{ccc}& & \\ & & \\ \begin{array}{c}-c_{12}^{\ast }\\ \ddots \\ \begin{array}{c}-a_{M-8}^{\ast }\\ \end{array}\end{array}& \begin{array}{c}\\ \ddots \\ \begin{array}{c}{b}_{M-8}^{\ast \ast }\\ -a_{M-4}^{\ast }\end{array}\end{array}& \begin{array}{c}\\ \\ \begin{array}{c}-c_{M-8}^{\ast }\\ b_{M-4}^{\ast \ast }\end{array}\end{array}\end{array}\right],$$

(16)

$$\tilde{Z}={\left[Z_{4,1}, Z_{8,1}, Z_{12,1},\dots ,Z_{M-8,1},Z_{M-4,1}\right]}^T,$$

(17)

and:

$$\tilde{F}={\left[F_{4,1}+a_4^{\ast }{F}_{0,1}, F_{8,1}, F_{12,1},\dots ,F_{M-8,1},F_{M-4,1}+c_{M-4}^{\ast }{F}_{M,1}\right]}^T.$$

(18)

4. Stability of Caputo’s Fractional Approximation with a Quarter-Sweep Scheme

This section discusses the stability analysis of the quarter-sweep Caputo fractional approximation to the time FDE as shown in Equation (3). In this paper, two common stability verification approaches are applied: Von Neumann’s method and the Lax equivalence theorem. The developed theorem for the stability of the quarter-sweep Caputo fractional approximation is presented below.

Theorem 1.

The quarter-sweep implicit approximation with Caputo’s fractional operator for$0<\alpha <1$, on the finite space$0\le x\le 1$and all$t>0$, is unconditionally stable.

Proof of Theorem 1.

Suppose that a solution function of Equation (3) has the form:

$$Z_{i,n}={\xi }_n{e}^{I\left(\lambda ih\right)},$$

(19)

where $I=\sqrt{-1}$ and $\lambda$ is a real number. Using Equation (19), the quarter-sweep Caputo fractional approximation shown by Equation (11) can be rewritten as:

$${\sigma }_{\alpha ,k}{\displaystyle \sum }_{j=2}^n{\omega }_{\alpha ,j}\left({\xi }_{n-\left(j-1\right)}{e}^{I\left(\lambda ih\right)}-{\xi }_{n-j}{e}^{I\left(\lambda ih\right)}\right)=a_i^{\ast }{\xi }_n{e}^{I\left(\lambda \left(i-4\right)h\right)}+b_i^{\ast }{\xi }_n{e}^{I\left(\lambda ih\right)}+c_i^{\ast }{\xi }_n{e}^{I\left(\lambda \left(i+4\right)h\right)}.$$

(20)

Equation (20) can be transformed into:

$$\begin{gathered} -{\sigma }_{\alpha ,k}{\xi }_{n-1}{e}^{I\left(\lambda ih\right)}+{\sigma }_{\alpha ,k}{\displaystyle \sum }_{j=2}^n {\omega }_{\alpha ,j}\left({\xi }_{n-\left(j-1\right)}{e}^{I\left(\lambda ih\right)}-{\xi }_{n-j}{e}^{I\left(\lambda ih\right)}\right)= \\ {a}_i^{\ast }{\xi }_n{e}^{I\left(\lambda \left(i-4\right)h\right)}+\left(b_i^{\ast }-{\sigma }_{\alpha ,k}\right){\xi }_n{e}^{I\left(\lambda ih\right)}+c_i^{\ast }{\xi }_n{e}^{I\left(\lambda \left(i+4\right)h\right)}. \end{gathered}$$

(21)

By simplifying and rearranging the distinct terms in Equation (21), one can obtain:

$$-{\sigma }_{\alpha ,k}{\xi }_{n-1}+{\sigma }_{\alpha ,k}{\displaystyle \sum }_{j=2}^n{\omega }_{\alpha ,j}\left({\xi }_{n-\left(j-1\right)}-{\xi }_{n-j}\right)={\xi }_n\left(\left(\left(a_i^{\ast }+c_i^{\ast }\right)\cos 4\lambda h\right)+\left(b_i^{\ast }-{\sigma }_{\alpha ,k}\right)\right),$$

(22)

which can be reduced to:

$${\xi }_n=\frac{{\xi }_{n-1}+{\sum }_{j=2}^n {\omega }_{\alpha ,j}\left({\xi }_{n-j}-{\xi }_{n-\left(j-1\right)}\right)}{\left(1-\frac{\left(a_i^{\ast }+c_i^{\ast }\right)\cos 4\lambda h+b_i^{\ast }}{{\sigma }_{\alpha ,k}}\right)}.$$

(23)

From Equation (23), for all values of $j, n, \omega ,\alpha , \lambda ,h,\sigma$, and $k$, it follows that:

$${\xi }_1\le {\xi }_0,$$

(24)

and when $n=2$, we have:

$${\xi }_2\le {\xi }_1+{\omega }_{\alpha ,2}\left({\xi }_0-{\xi }_1\right),$$

(25)

which can generally be stated as:

$${\xi }_n\le {\xi }_{n-1}+{\displaystyle \sum }_{j=2}^n {\omega }_{\alpha ,j}\left({\xi }_{n-j}-{\xi }_{n-\left(j-1\right)}\right),n\ge 2.$$

(26)

Since:

$${\xi }_n\le {\xi }_{n-1}+{\displaystyle \sum }_{j=2}^n{\omega }_{\alpha ,j}\left({\xi }_{n-j}-{\xi }_{n-\left(j-1\right)}\right)\le {\xi }_{n-j},$$

(27)

the inequalities in Equations (24)–(26) imply:

$${\xi }_n\le {\xi }_{n-1}\le {\xi }_{n-2}\le \cdots \le {\xi }_1\le {\xi }_0.$$

(28)

Thus, the stability of the quarter-sweep scheme via Caputo’s time-fractional operator is established as:

$${\xi }_n=\left|Z_{i,n}\right|\le {\xi }_0=\left|Z_{i,0}\right|.$$

(29)

 □

5. Convergence of Caputo’s Fractional Approximation with a Quarter-Sweep Scheme

This section discusses the convergence analysis of the quarter-sweep Caputo fractional approximation to the time FDE in Equation (3). To begin the analysis, we define $e_{i,n}=z\left(x_i,t_n\right)-Z_{i,n},i=1,2,\dots ,M-1,n=1,2,\dots ,N$, where $z_{i,n}$ and $Z_{i,n}$ are the exact and approximate solutions, respectively. Substitution of $e_{i,n}$ into Equations (11) and (13) gives:

$${\sigma }_{\alpha ,k}{\displaystyle \sum }_{j=2}^n {\omega }_{\alpha ,j}\left(e_{i,n-\left(j-1\right)}-e_{i,n-j}\right)=a_i^{\ast }{e}_{i-4,n}+b_i^{\ast }{e}_{i,n}+c_i^{\ast }{e}_{i+4,n}+R_{i,n},n\ge 2,$$

(30)

where:

$$R_{i,n}={\sigma }_{\alpha ,k}{\displaystyle \sum }_{j=2}^n {\omega }_{\alpha ,j}\left(z\left(x_i,t_{n-\left(j-1\right)}\right)-z\left(x_i,t_{n-j}\right)\right)-a_i^{\ast }z\left(x_{i-4},t_n\right)-b_i^{\ast }z\left(x_i,t_n\right)-c_i^{\ast }z\left(x_{i+4},t_n\right),$$

(31)

and:

$${\sigma }_{\alpha ,k}{e}_{i,0}=-a_i^{\ast }{e}_{i-4,1}+\left({\sigma }_{\alpha ,k}-b_i^{\ast }\right)e_{i,1}-c_i^{\ast }{e}_{i+4,1}+R_{i,1}.$$

(32)

Then, we have:

$$\frac{{\partial }^{\alpha }z\left(x_i,t_n\right)}{\partial t^{\alpha }}+C_{1}k={\sigma }_{\alpha ,k}{\displaystyle \sum }_{j=0}^n{\omega }_{\alpha ,j}\left(z\left(x_i,t_{n-\left(j-1\right)}\right)-z\left(x_i,t_{n-j}\right)\right),$$

(33)

and the quarter-sweep first and second-order central difference operators as follows:

$$\frac{{\partial }^{2}z\left(x_i,t_n\right)}{\partial x^2}+C_2{h}^2=\frac{z\left(x_{i-4},t_n\right)-2z\left(x_i,t_n\right)+z\left(x_{i+4},t_n\right)}{16h^2},$$

(34)

and:

$$\frac{\partial z\left(x_i,t_n\right)}{\partial x}+C_{3}h=\frac{z\left(x_{i-4},t_n\right)-z\left(x_{i+4},t_n\right)}{8h}.$$

(35)

Hence:

$$R_{i,n}=\frac{1}{{\sigma }_{\alpha ,k}}\left[\frac{{\partial }^{\alpha }z\left(x_i,t_n\right)}{\partial t^{\alpha }}-a\left(x\right)\frac{{\partial }^{2}z\left(x_i,t_n\right)}{\partial x^2}-b\left(x\right)\frac{\partial z\left(x_i,t_n\right)}{\partial x}-c\left(x\right)z\left(x_i,t_n\right)\right]+C_1{k}^{1+\alpha }+k^{\alpha }\left(C_2{h}^2+C_{3}h\right),$$

(36)

and also:

$$\left|R_{i,n}\right|\le C\left(k^{1+\alpha }+k^{\alpha }\left(h^2+h\right)\right),i=1,2,\dots M-1,n=1,2,\dots ,N,$$

(37)

where $C$ is an arbitrary constant.

Because $kn\le N$ is finite, we obtain the following theorem.

Theorem 2.

Let$Z_{i,n}$be the approximation of the exact value$z\left(x_i,t_n\right)$, computed by the quarter-sweep implicit approximation with Caputo’s fractional operator. Then, there is a positive constant$C$such that:

$$\left|z\left(x_i,t_n\right)-Z_{i,n}\right|\le C\left(k+h^2\right),i=1,2,\dots M-1,n=1,2,\dots ,N.$$

(38)

6. Concept and Formulation of Quarter-Sweep Preconditioned Relaxation Method

The presented PSOR iterative method is a continuation of the derivation of the quarter-sweep finite difference approximation via Caputo’s fractional derivative. The formulation of the PSOR method begins with the linear system, as shown by Equation (15). The coefficient matrix considered in Equation (15) has a large scale and is sparse. Thus, the proposed PSOR iterative method solves Equation (15) using the following matrix restructure. Firstly, a preconditioned linear system can be defined as:

$$A^{\prime }\tilde{Z}={\tilde{F}}^{\prime },$$

(39)

where the new matrices $A^{\prime }=PA{P}^T, {\tilde{F}}^{\prime }=P\tilde{F}$, and $\tilde{Z}=P^T\psi$. The preconditioned matrix labelled $P$ that is used in this paper is given by:

$$P=I+S,$$

(40)

where:

$$S=\left[\begin{array}{cccccc}0& -r_1& 0& 0& 0& 0\\ 0& 0& -r_2& 0& 0& 0\\ 0& 0& 0& -r_3& 0& 0\\ 0& 0& \ddots & \ddots & \ddots & 0\\ 0& 0& 0& 0& 0& -r_{m-1}\\ 0& 0& 0& 0& 0& 0\end{array}\right].$$

(41)

Secondly, let us consider that the new coefficient matrix $A^{\prime }$ in Equation (39) can be defined as the summation of the matrix components of:

$$A^{\prime }=Di-Lo-Va,$$

(42)

where $Di$, $Lo$, and $Va$ are diagonal, lower triangular, and upper triangular matrices, respectively. By combining Equations (39) and (42), the resultant QSPSOR iteration method can be derived as:

$${\psi }^{\left(\eta +1\right)}={\left(Di-\omega Lo\right)}^{-1}\left(\left(1-\omega \right)Di+\omega Va\right){\psi }^{\left(\eta \right)}+{\left(Di-\omega Lo\right)}^{-1}{\tilde{F}}^{\prime },$$

(43)

where ${\psi }^{\left(\eta +1\right)}$ represents an unknown vector at the $\left(\eta +1\right)$th iteration. The algorithm of the QSPSOR method used to compute the solution of the time FDE is described in the following

Table 1. QSPSOR computational algorithm.

	$\mathbf{Set}{\mathit{\psi }}^{\left(0\right)}=0$ and the Tolerance Error $\mathit{\epsilon }={10}^{-10}.$
(i)	$\mathrm{For}n=1,2,\dots ,N$ and for $i=4, 8,\dots ,M-4$, iterate the formula shown in Equation (43),
(ii)	$\mathrm{Compute}{\tilde{Z}}^{\left(\eta +1\right)}=P^T{\psi }^{\left(\eta +1\right)},$
(iii)	$\mathrm{Convergence}\mathrm{criterion}{\Vert {\tilde{Z}}^{\left(\eta +1\right)}-{\tilde{Z}}^{\left(\eta \right)}\Vert }_{\infty }<\epsilon ,$
(iv)	If the criterion is achieved, display approximate solutions.

.

7. Implementation and Application of C++ for Numerical Experiment

In this work, the C++ programming language was applied to implement the designed QSPSOR computational algorithm. C++ is a high-level programming language used in the field of numerical analysis. The C++ language also facilitates low-level coding because it is an extension of the medium-level programming language, C. Because C++ has all of the features and advantages of C, C++ enables the manipulation of low-level memory and the development of numerical iterative formulas to solve mathematical problems. In this work, the C++ simulation code was written to prioritize performance, speed, efficiency, and flexibility of use. Hence, implementing the quarter-sweep approximation based on Caputo’s fractional operator via the QSPSOR C++ code can yield important data, such as the number of iterations, execution time, and absolute errors.

For the numerical experiment, two problems of the time-fractional diffusion equation were selected to test the performance of the QSPSOR. For the efficiency comparison and analysis, two existing methods from our previous work were used, abbreviated as FSPSOR [19] and HSPSOR [20]. To compare the performance of these methods, three criteria were considered, namely, $\eta$—representing the number of iterations, $sec.$—representing the execution time of the C++ simulation code, and $\widehat{\epsilon }$—representing the magnitude of the absolute error. The three criteria for comparison and efficiency analysis were observed and taken based on three different values of fractional order, namely, $\alpha =0.25$, $\alpha =0.50$, and $\alpha =0.75$. Moreover, the simulation considered five different grid points, $M=128, 256, 512, 1024,$ and $2048$. The problems are considered below.

Example 1

([22]).Let the time-fractional initial boundary value problem be:

$$\frac{{\partial }^{\alpha }z\left(x,t\right)}{\partial t^{\alpha }}=\frac{{\partial }^{2}z\left(x,t\right)}{\partial x^2},0<\alpha \le 1,0\le x\le \gamma ,t>0,$$

(44)

where the boundary conditions are stated in fractional terms:

$$z\left(0,t\right)=\frac{2k{t}^{\alpha }}{\mathsf{\Gamma }\left(\alpha +1\right)},z\left(\gamma ,t\right)={\gamma }^2+\frac{2k{t}^{\alpha }}{\mathsf{\Gamma }\left(\alpha +1\right)},$$

(45)

and the initial condition is:

$$z\left(x,0\right)=x^2.$$

(46)

The analytical solution of Equation (44) is given by:

$$z\left(x,t\right)=x^2+\frac{2k{t}^{\alpha }}{\mathsf{\Gamma }\left(\alpha +1\right)}.$$

(47)

Example 2

([22]).Let the next time-fractional initial boundary value problem be given as:

$$\frac{{\partial }^{\alpha }z\left(x,t\right)}{\partial t^{\alpha }}=\frac{1}{2}{x}^2\frac{{\partial }^{2}z\left(x,t\right)}{\partial x^2},0<\alpha \le 1,0\le x\le \gamma ,t>0,$$

(48)

where the boundary conditions are given in fractional terms:

$$z\left(0,t\right)=0,z\left(1,t\right)=e^t,$$

(49)

and the initial condition is given by:

$$z\left(x,0\right)=x^2.$$

(50)

The exact solution of Equation (38) is:

$$z\left(x,t\right)=x^2{e}^t.$$

(51)

The collected numerical results of the implementation of FSPSOR, HSPSOR, and QSPSOR to solve Examples 1 and 2 are recorded in

Table 2. Comparison between FSPSOR, HSPSOR, and QSPSOR in solving Example 1.

$\mathit{M}$	Method	$\mathit{\alpha }=0.25$			$\mathit{\alpha }=0.50$			$\mathit{\alpha }=0.75$
		$\mathit{\eta }$	$\mathit{s}\mathit{e}\mathit{c}.$	$\widehat{\mathit{\epsilon }}$	$\mathit{\eta }$	$\mathit{s}\mathit{e}\mathit{c}.$	$\widehat{\mathit{\epsilon }}$	$\mathit{\eta }$	$\mathit{s}\mathit{e}\mathit{c}.$	$\widehat{\mathit{\epsilon }}$
128	FSPSOR	28	0.84	2.36 × 102	80	1.90	6.20 × 104	246	5.76	3.99 × 102
	HSPSOR	16	0.18	2.36 × 102	37	0.54	6.99 × 104	94	2.36	3.99 × 102
	QSPSOR	8	0.05	2.37 × 102	14	0.29	6.19 × 104	32	0.09	4.21 × 102
256	FSPSOR	53	5.33	2.43 × 102	211	17.84	5.69 × 104	806	67.75	3.97 × 102
	HSPSOR	34	2.20	2.43 × 102	94	6.90	6.21 × 104	303	34.65	3.97 × 102
	QSPSOR	15	0.27	2.44 × 102	39	2.37	6.99 × 104	101	12.84	4.03 × 102
512	FSPSOR	120	41.43	2.46 × 102	566	182.83	5.36 × 104	2635	843.91	3.96 × 102
	HSPSOR	67	21.65	2.46 × 102	246	86.09	5.36 × 104	988	421.58	3.96 × 102
	QSPSOR	31	5.04	2.47 × 102	100	40.61	6.21 × 104	337	198.20	3.96 × 102
1024	FSPSOR	250	372.35	2.48 × 102	1514	726.29	5.13 × 104	6012	1699.87	3.95 × 102
	HSPSOR	141	189.58	2.48 × 102	655	434.72	5.13 × 104	2413	1003.78	3.95 × 102
	QSPSOR	66	47.29	2.49 × 102	266	214.51	5.69 × 104	1095	501.76	3.95 × 102
2048	FSPSOR	808	901.76	2.49 × 102	4052	3469.73	5.13 × 104	35,289	7052.28	3.93 × 102
	HSPSOR	305	308.80	2.49 × 102	1788	1956.43	5.13 × 104	18,143	4025.90	3.93 × 102
	QSPSOR	143	95.25	2.50 × 102	709	365.43	5.35 × 104	3574	2000.83	3.93 × 102

and

Table 3. Comparison between FSPSOR, HSPSOR, and QSPSOR in solving Example 2.

$\mathit{M}$	Method	$\mathit{\alpha }=0.25$			$\mathit{\alpha }=0.50$			$\mathit{\alpha }=0.75$
		$\mathit{\eta }$	$\mathit{s}\mathit{e}\mathit{c}.$	$\widehat{\mathit{\epsilon }}$	$\mathit{\eta }$	$\mathit{s}\mathit{e}\mathit{c}.$	$\widehat{\mathit{\epsilon }}$	$\mathit{\eta }$	$\mathit{s}\mathit{e}\mathit{c}.$	$\widehat{\mathit{\epsilon }}$
128	FSPSOR	28	0.84	2.36 × 102	80	1.90	6.20 × 104	246	5.76	3.99 × 102
	HSPSOR	16	0.18	2.36 × 102	37	0.54	6.99 × 104	94	2.36	3.99 × 102
	QSPSOR	8	0.05	2.37 × 102	14	0.29	6.19 × 104	32	0.09	4.21 × 102
256	FSPSOR	53	5.33	2.43 × 102	211	17.84	5.69 × 104	806	67.75	3.97 × 102
	HSPSOR	34	2.20	2.43 × 102	94	6.90	6.21 × 104	303	34.65	3.97 × 102
	QSPSOR	15	0.27	2.44 × 102	39	2.37	6.99 × 104	101	12.84	4.03 × 102
512	FSPSOR	120	41.43	2.46 × 102	566	182.83	5.36 × 104	2635	843.91	3.96 × 102
	HSPSOR	67	21.65	2.46 × 102	246	86.09	5.36 × 104	988	421.58	3.96 × 102
	QSPSOR	31	5.04	2.47 × 102	100	40.61	6.21 × 104	337	198.20	3.96 × 102
1024	FSPSOR	250	372.35	2.48 × 102	1514	726.29	5.13 × 104	6012	1699.87	3.95 × 102
	HSPSOR	141	189.58	2.48 × 102	655	434.72	5.13 × 104	2413	1003.78	3.95 × 102
	QSPSOR	66	47.29	2.49 × 102	266	214.51	5.69 × 104	1095	501.76	3.95 × 102
2048	FSPSOR	808	901.76	2.49 × 102	4052	3469.73	5.13 × 104	35,289	7052.28	3.93 × 102
	HSPSOR	305	308.80	2.49 × 102	1788	1956.43	5.13 × 104	18,143	4025.90	3.93 × 102
	QSPSOR	143	95.25	2.50 × 102	709	365.43	5.35 × 104	3574	2000.83	3.93 × 102

, respectively. Based on Table 2 and Table 3, it can be observed that the QSPSOR method requires fewer iterations to obtain satisfactorily accurate approximate solutions of Examples 1 and 2. This significant improvement in the number of iterations results in a shorter execution time, which indicates QSPSOR can more efficiently solve Examples 1 and 2 than the two tested methods. The accuracy of the solution obtained by the QSPSOR method is comparable and almost equivalent to those of the FSPSOR and HSPSOR methods. A comparison of the different fractional orders of the selected examples shows that the accuracy of the solutions of the three numerical methods, which are based on the implicit finite difference scheme and Caputo’s time-fractional approximation, is greatest at $\alpha =0.50$, followed by $\alpha =0.75$ and $\alpha =0.25$. Although the QSPSOR method can solve the selected problems efficiently, the magnitude of the absolute errors of the solutions is slightly bigger than those of the FSPSOR and HSPSOR methods. Thus, the use of the quarter-sweep scheme as the complexity reduction approach for solving time-fractional behavior is a disadvantage due to the less effective direct computation of the remaining points. Future work will further investigate the appropriate treatment to reduce the absolute errors of the quarter-sweep difference scheme in the Caputo sense.

In addition, Figure 1, Figure 2, Figure 3 and Figure 4 show the performance graphs of the number of iterations and the execution time of the three methods. These graphs provide illustrations of the differences in terms of the efficiency of the combination of Caputo’s time-fractional approximation and preconditioned iterations via three different schemes: the standard implicit, half-sweep, and quarter-sweep. The percentage improvement in terms of the number of iterations for the proposed QSPSOR method, at $\alpha =0.25$, indicated that the number of iterations declined significantly, by 74.64% and 53.18%, compared to the FSPSOR and HSPSOR methods, respectively. Similarly, for execution time, the QSPSOR method was significantly faster, by 90.71% and 76.18%, than the FSPSOR and HSPSOR methods, respectively. At $\alpha =0.50$, the number of iterations needed by QSPSOR was 82.26% and 59.95% less than that required by the FSPSOR and HSPSOR methods, respectively. As a result, the C++ program execution time was 81.83% and 59.35% faster than that of the FSPSOR and HSPSOR methods, respectively. Finally, for the fractional-order $\alpha =0.75$, QSPSOR required 86.67% and 66.69% fewer iterations than the FSPSOR and HSPSOR methods, respectively, and the C++ program execution was 79.62% and 62.49% quicker than that of the FSPSOR and HSPSOR methods, respectively.

8. Conclusions

This paper discusses and presents the numerical solution of a fractional mathematical equation, namely, the time FDE, using a quarter-sweep mixed Caputo time-fractional approximation, and a C++ program to implement the QSPSOR iterative method. A preconditioned matrix was successfully implemented to modify the matrix structure to improve the convergence rate, and the QSPSOR method was formulated to efficiently solve the matrix. A numerical experiment using two alternative methods, namely FSPSOR and HSPSOR, showed the superiority of the proposed QSPSOR method in terms of the efficiency of solving two time-fractional initial boundary value problems of the diffusion equation. The findings of this paper can be highlighted as follows:

• The QSPSOR method significantly reduced the number of iterations and execution time compared to the existing FSPSOR and HSPSOR methods. On average, QSPSOR reduced the number of iterations and execution time by 81.19% and 84.06%, respectively, compared to FSPSOR. Moreover, QSPSOR reduced the number of iterations and execution time compared to HSPSOR by 59.94% and 66.00%, respectively. The result also showed that using the quarter-sweep scheme and PSOR iteration can reduce the computational complexity for solving the time FDE with a large matrix size.
• Observations regarding the accuracy of all of the implemented numerical methods indicated that their numerical solutions were in good agreement. Furthermore, the accuracy of the solutions of the time-FDE problems obtained by each of the three numerical methods was greater at $\alpha =0.50$, followed by $\alpha =0.75$ and $\alpha =0.25$. The combination of the quarter-sweep implicit finite difference scheme and Caputo’s time-fractional derivative enabled an accurate solution for the time FDE to be computed.
• However, a disadvantage of the quarter-sweep difference scheme is that the magnitude of the absolute errors is slightly larger than that of the two previous methods. The accuracy of the quarter-sweep scheme can be improved by applying a suitable treatment.

Thus, the combination of a quarter-sweep finite difference method, Caputo’s time-fractional derivative, and the preconditioned successive over-relaxation method showed good potential to solve different types of fractional mathematical equations, and provides a future direction for this field of research.