New preconditioning and half-sweep accelerated overrelaxation solution for fractional differential equation

1 Introduction

Fractional calculus has gained considerable popularity and importance for almost five decades now. It is mainly from various demonstrated applications in biological science, physical science and other branches of sciences. Fractional calculus has significantly contributed to the modelling of transmission of Covid-19 infection (Cui and Liu, 2022), pharmacokinetic compartments (Azizi, 2022), Meningitis with treatment and vaccination (Peter et al., 2022), tumour and immune cells interactions (Tang et al., 2022), mechanical behaviour of asphalt mastic (Lagos-Varas et al., 2022), fluid flow and heat transfer (Turkyilmazoglu, 2022), control behaviour of wearable exoskeletons (Sun et al., 2021) and control behaviour of a knee joint orthosis (Delavari and Jokar, 2021). Many different fractional differential equations (FDE) have arisen from the realistic applications of fractional calculus. FDE is a generalization of differential equations based on the established theory and application of fractional calculus. FDE can also be considered the extended partial differential equations by modifying the integer-order derivative into the fractional-order derivative.

The solutions of FDEs must be obtained to understand the fractional mathematical models. Various solution methods have been proposed to the literature, such as the finite difference method with collocation (Mesgarani et al., 2021;Safdari et al., 2020;Jaleb and Adibi, 2019), finite difference method with preconditioners (Barakitis et al., 2022; Shao and Kang, 2022; Sunarto et al., 2022; Sunarto et al., 2021), finite difference method with Lucas polynomials (Ali et al., 2022), Adomian decomposition method (Turkyilmazoglu, 2022; Ahmad et al., 2022; Turkyilmazoglu, 2021) and variational iteration method (Ibraheem et al., 2022). Following the high interest towards the finite difference method with preconditioners, this paper aims to investigate the approximate solution of a type of FDE, the space-fractional diffusion equation, using a new preconditioning matrix to develop an efficient half-sweep accelerated overrelaxation iterative method. This paper utilizes the Caputo fractional derivative to treat the space-fractional derivative because it allows the inclusion of traditional and conventional initial-boundary conditions in the formulation of the problem (Elsayed and Orlov, 2020). In addition, the Caputo space-fractional derivative's memories affect the dynamics of the considered variables (Sene, 2022). The importance of Caputo space-fractional can be seen in the modelling of biological models (Haghi and Ghanbari, 2022), sediment suspension in ice-covered channels (Wang et al., 2022), drug diffusion through the skin (Caputo and Cametti, 2021) and chaotic processes (Owolabi et al., 2020).

The paper's focus is to assess the improvement in terms of the convergence rate of the solution obtained by the proposed iterative method. Among various iterative methods that can be used to solve the generated system of equations from an FDE (She et al., 2023; AllaHamou et al., 2022; Wen et al., 2022; Sun et al., 2022; Tang and Huang, 2022), the paper proposes a modified accelerated overrelaxation iterative method using a new preconditioning matrix with a half-sweep iteration strategy. The paper's contribution is a new preconditioned iterative method that can solve a space-fractional diffusion equation at a good efficiency level. The following sections of the paper are organized: Section 2 formula tesa discrete approximation to a one-dimensional linear space-fractional diffusion equation using a half-sweep type finite difference method in the Caputo sense. Section 3 derives the proposed iterative method to solve the generated system of equations from the discretized problem. Section 4 illustrates the numerical results of solving several initial-boundary value problems using the proposed numerical method and the comparison analysis against the standard preconditioned accelerated overrelaxation method (Sunarto et al., 2016). The conclusion of the paper is stated in Section 5.

2 Half-sweep type finite difference method in caputosense

This section describes the formulation of a discrete approximation to a one-dimensional linear space-fractional diffusion equation using a half-sweep type finite difference method in the Caputo sense. The paper usesa general space fractional FDE in the formulation, which is given by (Reutskiy and Lin, 2018),

Based on Eq. (1), the variables $c_i,i=1,2$ , and 3 are either constants or functions in terms of $x$ while $g(x,t)$ is a source function.

This paper utilizes half-sweep type implicit finite difference schemes to discretize Eq. (1) for the time derivative, integer-order space derivative and other functions (Ibrahim and Abdullah, 1995; Sunarto et al., 2021; Chew et al., 2021). Meanwhile, Caputo fractional derivative is applied to approximate the fractional-order space derivative. Below is the following established definition of Caputo fractional derivative used in the discretization (Oldham and Spanier, 2006):

Definition 1. Let $x$ be the upper limit of the integral, and a real number $\beta$ be the fractional order, such that $0\le m-1<\beta <m$ where $m$ is a positive integer. Then, $f^{\left(m\right)}\left(\xi \right)$ represents the $m$ -th order derivative of a smooth function $f\left(x\right)$ . Hence, the Caputo fractional derivative of $f\left(x\right)$ is defined as

Combining half-sweep type finite difference schemes and Eq. (3) gives the following discrete approximation to the space-fractional derivative,

Then, putting Eq. (5) together with the half-sweep finite differences for other derivatives such as first-order time derivative, first-order space derivative and source functions, Eq. (1) can be rewritten in the form of finite difference approximation equation in Caputo sense as follows, $$k^{-1}\left(U_{i,n}-U_{i,n-1}\right)=c_1\left(x\right)\rho {\sum }_{j=0,2,4,\cdots }^{i-2}{\sigma }_j^{\beta }\left(U_{i-\left(j-2\right),n}-2U_{i-j,n}+U_{i-\left(j+2\right),n}\right)$$

$$c_2^{\ast }{U}_{i-2,n}+c_3^{\ast }{U}_{i,n}-c_2^{\ast }{U}_{i+2,n}-c_1^{\ast }{\sum }_{j=0,2,4,\cdots }^{i-2}{\sigma }_j^{\beta }\left(U_{i-\left(j-2\right),n}-2U_{i-j,n}+U_{i-\left(j+2\right),n}\right)$$

Based on Eq. (9), one may have the following equations subject to different values of $j$ . For instance, when $j=0$ , Eq. (9) becomes

Hence, when the pattern continues for $j=4,6,\cdots$ , one can easily obtain a general form of the equation that can generate a large-scale system of equations as follows,

When Eq. (12) takes all points bounded by a specified solution domain, the large-scale system of equations can be expressed in the form of a matrix equation,

Noted that the matrix dimensions of matrix $M$ , $\widehat{U}$ and $\widehat{f}$ are $\left(s-2\right)\times \left(s-2\right)$ , $\left(s-2\right)\times 1$ , and $\left(s-2\right)\times 1$ , respectively. This paper suggests that an efficient iterative method needs to be developed to solve a complex matrix equation like Eq. (20). Hence, this paper proposes a new preconditioning matrix that can enhance the convergence rate of the iterated solutions. Moreover, this paper develops a new iterative method called the half-sweep preconditioned accelerated overrelaxation.

3 Derivation of a preconditioned iterative method

This section is devoted to showing the derivation of the proposed preconditioned iterative method to solve the system of equations shown (Eq. (20)). From here, the paper shall use HSPAOR to stand for the proposed method to solve space-fractional diffusion problems. To begin the derivation, let's consider a transformed matrix equation that corresponds to Eq. (20) as follows,

Eq. (24) is obtained using the following matrix transformations with a new preconditioning matrix $P$ ,

Based on the coefficient matrix $A$ that presents in the transformed matrix equation shown in Eq. (24), this paper considers a unique decomposition of matrix $A$ that is given by

Hence, to achieve the desired convergence rate and the objective of the numerical study, which is to investigate the numerical solutions, a manual selection of accelerating parameters is conducted by running the developed simulation program several times until the smallest number of iterations is obtained. The selection procedure can be described as follows. We initially let $\theta =1$ and use different values of $\omega$ within the range $(1,2)$ . When the smallest number of iterations is obtained for some value of $\omega$ , by using the “optimum” value of $\omega$ , we increase the value of $\theta$ gradually until the final smallest number of iterations is obtained. The implementation of the HSPAOR method is programmed using the C++ programming language. The structure of the code and instructions are made thoroughly. Due to the copyright issue, the paper can only provide the following algorithm.

Algorithm 1: HSPAOR method to solve space-fractional diffusion equations

1. Set the initial guess ${\widehat{U}}^{\left(k=0\right)}=0,$ and the tolerance error $∊={10}^{-10}.$

2. For $i=2,4,\cdots ,s-2$ , iterate Eq. (31).

3. For $i=1,3,\cdots ,s-1$ , run linear interpolation module.

4. If $\left|{\widehat{U}}^{\left(k+1\right)}-{\widehat{U}}^{\left(k\right)}\right|\le ∊,$ then go to the next time-step or $n=n+1$ .

5. If the time-step reaches the final step or $n=N$ , display outputs such as numerical solutions, the maximum number of iterations, program execution time, and maximum absolute errors.

4 Numerical experiment and results

Section 4 illustrates the proposed method's results by solving several initial-boundary value problems of space-fractional diffusion. Below are the following test problems considered in this paper.

Based on Eq. (32), the value of $\mathrm{\Gamma }\left(1.5\right)x^{0.5}$ represents the diffusion coefficient, while the function $\left(x^2+1\right)\cos (t+1)-2x\sin (t+1)$ is the source of diffusion. The accuracy of the numerical solution obtained by HSPAOR is compared to the exact solution,

Based on Eq. (32), $\mathrm{\Gamma }\left(1.2\right)x^{\beta }$ represents the diffusion coefficient, while $3x^2\left(2x-1\right)e^{-t}$ is the source function. The accuracy of the numerical solution obtained by HSPAOR is compared to the exact solution,

The results considered take account of numerical solutions, the number of iterations to obtain the final solutions $\left(k_f\right)$ , the final time after completing the C++ program,which is measured in seconds $\left(s\right)$ and the value of absolute errors. Fig. 1 until 6 show the numerical solutions obtained by HSPAOR after solving Examples 1 and 2 using $\beta =1.2,1.5,$ and $1.8$ . The solutions are compared to the exact solutions at various points and time level $T=2.0$ seconds.

Close

Fig. 1

Numerical solutions by HSPAOR against exact solutions of Example 1 at $\beta =1.2$ .

Export to PPT

Based on Figs. 1 through 3, the effectiveness of HSPAOR in computing numerical solutions of Example 1 at various orders of space-fractional is illustrated. The numerical solutions are sufficiently close to the provided exact solutions at $\beta =1.2$ and well-fitted to the exact solutions at both $\beta =1.5$ and $1.8$ . However, HSPAOR shows some disadvantages in computing the numerical solutions of Example 2 at $\beta =1.2$ and $1.5$ compared to the exact solutions, see Figs. 4 and 5. The accuracy of the solutions by HSPAOR is better when the value of space-fractional order is set to be greater than 1.5 or $\beta =1.8;$ for instance, see Fig. 6.

Close

Fig. 2

Numerical solutions by HSPAOR against exact solutions of Example 1 at $\beta =1.5$ .

Export to PPT

Close

Fig. 3

Numerical solutions by HSPAOR against exact solutions of Example 1 at $\beta =1.8$ .

Export to PPT

Close

Fig. 4

Numerical solutions by HSPAOR against exact solutions of Example 2 at $\beta =1.2$ .

Export to PPT

Close

Fig. 5

Numerical solutions by HSPAOR against exact solutions of Example 2 at $\beta =1.5$ .

Export to PPT

Close

Fig. 6

Numerical solutions by HSPAOR against exact solutions of Example 2 at $\beta =1.8$ .

Export to PPT

Next, comparison in terms of the number of iterations, program completion time and maximum absolute error between HSPAOR and two tested methods, such as the standard or full-sweep preconditioned accelerated overrelaxation (FSPAOR) (Sunarto et al., 2016) and implicit Euler (Meerschaert and Tadjeran, 2006) is shown in Tables 1 until 6. The comparison is conducted using three different values of space-fractional order, $\beta =1.2,1.5$ , and $1.8$ , and five different numbers of domain points for the consistency of the solutions.

Based on Tables 1 until 6, the comparison results show that the HSPAOR method is more efficient than he FSPAOR and implicit Euler methods in solving Examples 1 and 2. The number of iterations and program completion time required by the HSPAOR method to obtain the final numerical solutions at all different points are significantly lesser than the other two tested methods. However, the absolute errors produced by the HSPAOR method are slightly larger than the FSPAOR and implicit Euler methods for Example 1 using $\beta =1.5$ and $1.8$ and Example 2 using $\beta =1.8$ . Furthermore, by observing the consistency of the numerical solutions with the increasing number of points in computation, this paper found that the absolute errors show some sign of gradual growth for Example 1 at $\beta =1.2$ and Example 2 at all values of $\beta$ .

To complete the numerical experiment, this paper compares the maximum absolute errors produced by the proposed HSPAOR method (with a time-step 0.2) with some numerical methods, including the methods that utilize the Chebyshev polynomial of degree $n_p$ . The error comparison is made using a similar setting of Example 2 that has been done (Khader, 2011; Saadatmandi and Dehghan, 2011; Azizi and Loghmani, 2013). Table 7 shows the comparison in terms of maximum absolute errors against the selected three methods.

Based on the findings through the numerical experiment, HSPAOR possesses the advantage in terms of computational efficiency, especially when a large system of equations is considered. The reason is that the iteration procedure by the preconditioned accelerated overrelaxation is highly efficient in computing the generated system of equations. Besides that, using a half-sweep strategy in formulating the finite difference approximation in the Caputo sense has successfully reduced the computational complexity in the developed program. However, to achieve a greater efficiency level, the accuracy of the solution becomes the trade-off. The disadvantage of the HSPAOR method is revealed when it is used to solve Example 2 using $\beta =1.2$ and $1.5$ . Since the development of HSPAOR is based on implicit finite difference schemes, the accuracy of HSPAOR is limited by the properties of implicit finite difference schemes, which are second-order accurate in space. This paper hypothesized that the magnitude of absolute errors could be reduced using higher-order finite difference schemes and different fractional definitions.

5 Conclusion

This paper successfully developed an efficient half-sweep accelerated overrelaxation iterative method using a new preconditioning matrix to solve several space-fractional diffusion problems. The Caputo fractional derivative is compatible with formulating a discrete approximation equation via implicit finite difference schemes. The numerical results showed the superiority of the proposed method in terms of solution efficiency against the standard preconditioned accelerated overrelaxation and implicit Euler methods. When the absolute errors by the proposed method are compared against several existing numerical methods, the errors are slightly larger than all considered methods. The magnitude of errors can be reduced by using higher-order finite difference schemes and different fractional definitions. Based on the performance of the proposed method in terms of efficiency, it has the potential to solve a variety of space-fractional diffusion models efficiently. Future investigation will improve the solutions' absolute errors so that the proposed method's reliability can be increased.