Approximation Solution of the Fractional Parabolic Partial Differential Equation by the Half-Sweep and Preconditioned Relaxation

Abstract

In this study, the numerical solution of a space-fractional parabolic partial differential equation was considered. The investigation of the solution was made by focusing on the space-fractional diffusion equation (SFDE) problem. Note that the symmetry of an efficient approximation to the SFDE based on a numerical method is related to the compatibility of a discretization scheme and a linear system solver. The application of the one-dimensional, linear, unconditionally stable, and implicit finite difference approximation to SFDE was studied. The general differential equation of the SFDE was discretized using the space-fractional derivative of Caputo with a half-sweep finite difference scheme. The implicit approximation to the SFDE was formulated, and the formation of a linear system with a coefficient matrix, which was large and sparse, is shown. The construction of a general preconditioned system of equation is also presented. This study’s contribution is the introduction of a half-sweep preconditioned successive over relaxation (HSPSOR) method for the solution of the SFDE-based system of equation. This work extended the use of the HSPSOR as an efficient numerical method for the time-fractional diffusion equation, which has been presented in the 5th North American International Conference on industrial engineering and operations management in Detroit, Michigan, USA, 10–14 August 2020. The current work proposed several SFDE examples to validate the performance of the HSPSOR iterative method in solving the fractional diffusion equation. The outcome of the numerical investigation illustrated the competence of the HSPSOR to solve the SFDE and proved that the HSPSOR is superior to the standard approximation, which is the full-sweep preconditioned SOR (FSPSOR), in terms of computational complexity.

1. Introduction

In recent years, many effective mathematical physics models have been developed using the theory and applying the partial fractional derivatives. In most of the fractional partial-differential-equation-related works of literature, partial fractional derivatives appear in many anomalous phenomena modelling and in complex systems theory. Examples include a fractional mathematical model of the dynamics of the cancer chemotherapy effect based on both the time instant and the time history [1]; a fractional mathematical model of the population dynamics among cancer stem cells, tumor cells, healthy cells, the effects of excess estrogen and the body’s natural immune response on the cell populations [2]; a fractional mathematical model to investigate the dynamics of tuberculosis between the children and the adults [3]; a fractional model an of energy supply-demand system [4]; and a fractional mathematical model of the COVID-19 pandemic [5].

The success of fractional partial differential equations was attributed to the derivative order’s generalization, from integer-order to arbitrary order. The fractional derivatives in the fractional partial differential equations have an effective memory function that enables many physical phenomena to be described effectively. Fractional partial differential equations can be categorized as a time-fractional type, a space-fractional type, and a time- and space-fractional type. Our research focused on investigating the numerical solution of the space-fractional partial differential equation, particularly the parabolic type equation such as the space-fractional diffusion equation (SFDE). This research focus came from the need to discover an efficient solution for the SFDE problem.

Based on our short review of the existing solution methods, many researchers have suggested the finite-difference method. The method of finite difference can be implicit or explicit, depending on its stability to obtain the solution [6,7,8,9,10]. Since SFDE involves changing the time and space of a continuous variable, numerical treatment such as finite difference discretization is necessary to transform the differential equation into a finite system of linear equations that a computer can solve. Although SFDE can be solved without discretization, this research focused on improving an iterative method that is mostly formulated based on a discretized equation. Furthermore, the iterative method works best in obtaining an accurate approximate solution when a finite system of equations is large and complex. Other than the method of finite difference, several methods have been proposed to solve the space-fractional problems, such as the spectral collocation method [11], the boundary value method [12], the finite volume method [13], and the method of lines and splines [14].

In this study, an unconditionally stable implicit finite difference scheme and the $\beta$ order Caputo fractional partial derivative were executed to discretize the SFDE and to obtain the correct approximation equation. The generated system of equation after the approximation equation was used on the solution domain, which led to a tridiagonal linear system. Since the sparse and large-scale linear system’s coefficient matrix is difficult to be solved analytically, an iterative method is employed as an alternative. In terms of the efficiency of the iterative processes, many authors [8,15,16] have suggested and debated over several numerical iterative methods. Besides that, the block iteration has been presented [17] to show the computation cost efficiency improvement by separating systems of equations. The preconditioned iterative methods have also been widely recognized as the competent methods for solving linear systems. Thus, the symmetry in solving the SFDE problem via a numerical method exists in the use of the finite difference discretization scheme and the iterative method. A finite difference scheme needs to be unconditionally stable so that the solution can be obtained regardless of the time and space step sizes used. On the other hand, an iterative method needs to be efficient, which means the accurate solution of a highly complex linear system can be computed in a short time.

The contribution of this study was to build and examine the efficacy of the half-sweep preconditioned SOR (HSPSOR) method, which was formulated from the use of implicit finite difference and the Caputo fractional derivative, for resolving the SFDE implicitly. We also used the full-sweep preconditioned SOR (FSPSOR) iterative techniques as a control method to analyze the efficacy of the HSPSOR method. To begin the formulation and the investigation of the HSPSOR method, we considered the following general form of the SFDE as follows:

$$\frac{\partial U\left(x,t\right)}{\partial t}=A_1\frac{{\partial }^{\beta }U\left(x,t\right)}{\partial x^{\beta }}+A_2\frac{\partial U\left(x,t\right)}{\partial x}+A_{3}U\left(x,t\right)+B\left(x,t\right),$$

(1)

where $A_i,i=1, 2,$ and $3$ were arbitrary constants, and $B\left(x,t\right)$ was a known function. Notice that when the value of $\beta$ in Equation (1) equals 2, one can get a usual diffusion equation with a second-order derivative in space. The derivative of order-$\beta$ for an SFDE usually lies between $1<\beta \le 2$ to provide flexibility in the study of the effect of the medium and the space interaction with the fluid [18]. The solution of Equation (1), which is $U\left(x,t\right)$, can be approximated using the method of finite difference subject to the initial and boundary conditions as given by:

$$U\left(x,0\right)=U_0, 0\le x\le l,$$

(2)

and

$$U\left(0,t\right)=U_L,U\left(l,t\right)=U_R,0<t<T.$$

(3)

Before the space-fractional term in Equation (1) is discretized by the finite difference mean, the following established definitions from the theory of fractional derivatives must be defined as follows [19]:

Definition 1.

Let a real number$\alpha >0$, and let the function$f$be continuous on$F^{\prime }=\left(0,\infty \right)$and integrable on any finite subinterval of$F=\left[0,\infty \right)$. Then, for$x>0$, the Riemann–Liouville fractional integral,$f$with the order$\alpha$is defined as:

$$J^{\alpha }=\frac{1}{\Gamma \left( \alpha \right)}{{\displaystyle \int }}_0^x{\left(x-\xi \right)}^{ \alpha -1}f\left(\xi \right)d\xi .$$

(4)

Definition 2.

Let a real number$\beta >0$, such that$m-1<\mathsf{\beta }\le \mathrm{m}$, where$m$is a natural number element. Let the function$f$be continuous on$F^{\prime }=\left(0,\infty \right)$and integrable on any finite subinterval of$F=\left[0,\infty \right)$. Then, for$x>0$, the Caputo fractional,$f$with the order$\beta$is defined as:

$$D_x^{\beta }f\left(x\right)=\frac{1}{\Gamma \left(m-\beta \right)}{{\displaystyle \int }}_0^x\frac{f^m\left(\xi \right)}{{\left(x-\xi \right)}^{\beta -m+1}}d\xi .$$

(5)

From Definition 2, the following properties hold:

$$D_x^{\beta }C=0,$$

(6)

and

$$D_x^{\beta }{x}^{\eta }=\left\{\begin{array}{c}0,\\ \frac{\Gamma \left(\eta +1\right)}{\Gamma \left(\eta +1-\beta \right)}{x}^{\eta -\beta },\end{array}\right.$$

(7)

where $\eta$ is an element of natural number and $\Gamma \left(.\right)$ is a gamma function.

2. Approximation to a Space-Fractional Diffusion Equation

Equation (1) approximates the space-fractional term using the Caputo definition and the second-order half-sweep finite-difference. For more details about half-sweep see [17,20,21,22]. Suppose that $h=\frac{𝓵}{P},$ with $P$ as any positive integer. Then, the space-fractional term can be equated as:

$$\frac{{\partial }^{\beta }U\left(x,t\right)}{\partial x^{\beta }}=\frac{1}{\Gamma \left(2-\beta \right)}{{\displaystyle \int }}_0^x\frac{{\partial }^{2}U\left(\xi ,t\right)}{\partial x^2}{\left(x-\xi \right)}^{1-\beta }d\xi ,$$

(8)

which is also equivalence to:

$$\frac{1}{\Gamma \left(2-\beta \right)}{\displaystyle \sum }_{j=0, 2,4,\dots }^{i-2}{{\displaystyle \int }}_{jh}^{\left(j+1\right)h}\frac{U_{i+2-j,n}-2U_{i-j,n}+U_{i-2-j,n}}{2h^2}{\left(Ph-\xi \right)}^{1-\beta }d\xi .$$

(9)

Let $\phi =\frac{{2\mathrm{h}}^{-2}}{\Gamma \left(3-\beta \right)}$ and $g_j^{\beta }={\left(\frac{j}{2}+1\right)}^{2-\beta }-{\left(\frac{j}{2}\right)}^{2-\beta }$. Then, from Equation (9), we obtain:

$$\frac{{\partial }^{\beta }U\left( x_i, t_n\right)}{\partial x^{\beta }}=\phi {\displaystyle \sum }_{j=0,2,4,\dots }^{i-2}{g}_j^{\beta }\left(U_{i+2-j,n}-2U_{i-j,n}+U_{i-2-j,n}\right).$$

(10)

Using Equation (10) and the implicit finite difference scheme for the remaining derivatives in Equation (1), the simplified approximation to the SFDE via the finite difference and Caputo partial derivative can be formulated into:

$$\begin{gathered} \gamma U_{i,n}-A_1\phi {\displaystyle \sum }_{j=0,2,4,\dots }^{i-2}{g}_j^{\beta }\left(U_{i+2-j,n}-2U_{i-j,n}+U_{i-2-j,n}\right)-\frac{A_2}{4h}\left(U_{i+2,n}-U_{i-2,n}\right) \\ -A_3{U}_{i,n}-B_{i,n}=\gamma U_{i,n-2}, \end{gathered}$$

(11)

for $i=2,4,\dots ,P-2$. Equation (11) can be rearranged neatly into:

$$b_i{U}_{i-2,n}+\left({\gamma }_i-c_i\right)U_{i,n}-b_i{U}_{i+2,n}-a_i^{*}=f_i,$$

(12)

where $b_i=\frac{A_2}{4h},{\gamma }_i=\Delta t,c_i=A_3,a_i^{*}=a_i{\displaystyle \sum }_{j=0,2,4,\dots }^{i-2}{g}_j^{\beta }\left(U_{i+2-j,n}-2U_{i-j,n}+U_{i-2-j,n}\right),$ $a_i=A_1\phi ,$ and $f_i=\gamma U_{i,n-2}+B_{i,n}$. Based on Equation (12), we can simply state an equation for $n>3$ as:

$$p_i{U}_{i-6,n}+q_i{U}_{i-4,n}+r_i{U}_{i-2,n}+s_i{U}_{i,n}+v_i{U}_{i+2,n}-R_i=f_i,$$

(13)

where $R_i=a_i{\displaystyle \sum }_{j=0,2,4,\dots }^{i-2}{g}_j^{\beta }\left(U_{i+2-j,n}-2U_{i-j,n}+U_{i-2-j,n}\right),p_i=-a_i{g}_2^{\beta },q_i=-a_i{g}_1^{\beta }+2a_i{g}_2^{\beta },r_i=b_i-a_i{g}_2^{\beta }+2a_i{g}_1^{\beta }-a_i,s_i=-a_i{g}_1^{\beta }+2a_i+\gamma -c_i,$ and $v_i=-a_i-b_i.$ From Equation (13), the system of the linear equation can be expressed as follows:

$$A\underset{˜}{U}=\underset{˜}{f},$$

(14)

where

$$A={\left[\begin{array}{cccccccc}{s}_2& v_2& & & & & & \\ r_4& s_4& v_4& & & & & \\ q_6& r_6& s_6& v_6& & & & \\ p_8& q_8& r_8& s_8& v_8& & & \\ & p_{10}& q_{10}& r_{10}& s_{10}& v_{10}& & \\ & & \ddots & \ddots & \ddots & \ddots & \ddots & \\ & & & p_{P-4}& q_{P-4}& r_{P-4}& s_{P-4}& v_{P-4}\\ & & & & p_{P-2}& q_{P-2}& r_{P-2}& s_{P-2}\end{array}\right]}_{\left(P-2\right)\times \left(P-2\right)},$$

$$\underset{˜}{U}={\left[\begin{array}{cccccc}{U}_2& U_4& U_6& \cdots & U_{P-4}& U_{P-2}\end{array}\right]}^T,$$

and

$$\underset{˜}{f}={\left[\begin{array}{cccccc}{f}_2-r_2{U}_2& f_4+q_4{U}_4& f_6+p_6{U}_6& \cdots & f_{P-4}+R_{P-4}& f_{P-2}-r_{P-2}{U}_P\end{array}+R_{p-2}\right]}^T$$

3. Half-Sweep Preconditioned SOR Formulation

The system (Equation (14)) is transformed in the form of the linear equation’s preconditioned system as follows:

$$A^{*}\underset{˜}{x}=\underset{~}{f^{*}},$$

(15)

where $A^{*}=\rho A{\rho }^T,\underset{~}{f^{*}}=\rho \underset{~}{f,}$ and $\underset{˜}{U}={\rho }^T\underset{˜}{x}$.

Further, to discover the numerical tridiagonal linear system solutions of Equation (15), we considered the half-sweep preconditioned SOR (HSPSOR) iterative method. The $\rho$ was the preconditioned matrix and was given as [23]:

$$\rho =I+S,$$

(16)

where

$$S={\left[\begin{array}{cccccc}0& -v_1& 0& 0& 0& 0\\ 0& 0& -v_2& 0& 0& 0\\ 0& 0& 0& -v_3& 0& 0\\ 0& 0& \ddots & \ddots & \ddots & 0\\ 0& 0& 0& 0& 0& -v_{P-1}\\ 0& 0& 0& 0& 0& 0\end{array}\right]}_{\left(P-1\right)\times \left(P-1\right)}.$$

Meanwhile, $I$ was an identity matrix. The generated coefficient matrix $A^{*},$ as in Equation (15), can be split as follows:

$$A^{*}=D-L-V,$$

(17)

where $D, L,$ and $V$ were the diagonal matrices, the lower triangular, and the upper triangular, respectively. Considering Equation (17), the HSPSOR iterative method can be rewritten in the form of:

$${\underset{˜}{x}}^{\left(k+1\right)}={\left(D-\omega L\right)}^{-1}\left[\omega V+\left(1-\omega \right)D\right]{\underset{˜}{x}}^{\left(k\right)}+{\left(D-\omega L\right)}^{-1}\underset{~}{f^{*}},$$

(18)

where the unknown vector, ${\underset{˜}{x}}^{\left(k+1\right)}$, was given the solution at the $\left(k+1\right)\mathrm{th}$ iteration, and the relaxation parameters were chosen within the range, $1<\omega <2$. The computational algorithm of the HSPSOR iterative method was as follows.

Algorithm 1 The HSPSOR iterative method.

i. $\mathrm{Initialize} \underset{˜}{U}\leftarrow 0 \mathrm{and} \epsilon \leftarrow {10}^{-10}.$

ii. $\mathrm{For} j=0, 1,\dots ,n, \mathrm{and} \mathrm{for} i=0, 1, 2,\dots ,P, \mathrm{calculate} \mathrm{Equation} \left(18\right), \mathrm{then} \mathrm{approximate} \mathrm{solutions}:$ $${\underset{˜}{U}}^{\left(k+1\right)}={\rho }^T{\underset{˜}{x}}^{\left(k+1\right)}.$$ $\mathrm{If} \mathrm{the} \mathrm{convergence} \mathrm{criterion} \mathrm{is} \mathrm{satisfied}, \mathrm{that} \mathrm{is}:$ $$⎥⎥{\underset{˜}{U}}^{\left(k+1\right)}-{\underset{˜}{U}}^{\left(k\right)}⎥⎥\le \epsilon ,$$ $\mathrm{go} \mathrm{to} \mathrm{Step} \left(\mathrm{iii}\right). \mathrm{Otherwise}, \mathrm{go} \mathrm{back} \mathrm{to} \mathrm{Step} \left(\mathrm{ii}\right).$

iii Stop.

4. Numerical Evaluation via C++

In this section, we implemented two examples of SFDE, using C++ programming to verify the HSPSOR method’s effectiveness. The C++ programming language was used to code Algorithm 1 because of a better coding organization and comprehension, which enabled numerical computation to be conducted accurately. Furthermore, C++ programming was preferred for this work compared to other mathematics applications such as Mathematica and Maple because the number of iterations and the execution time could be recorded and optimized according to the code arrangement.

For the comparison purpose, we considered a standard or full-sweep preconditioned SOR by taking into account the number of iterations, the execution time (seconds), and the maximum error at the values of $\beta =1.2,\hspace{0.33em}\beta =1.5,\hspace{0.33em}\mathrm{and} \beta =1.8$. The tolerance error considered in Algorithm 1 and the C++ simulation was set at $\epsilon ={10}^{-10}$. This stopping condition was fixed for the different mesh sizes used, M; that is, 128, 256, 512, 1024, and 2048. There was no best value of $\epsilon$ because, when a large value of $\epsilon$ was chosen such as $\epsilon ={10}^{-2}$, the number of iterations became small, and the result became inaccurate. Moreover, when a small value of $\epsilon$ was used, such as $\epsilon ={10}^{-100}$, the number of iterations became large without improving the accuracy. Therefore, $\epsilon ={10}^{-10}$ was selected arbitrarily to check the convergence of the solution. The solution can only be obtained through the simulation of Algorithm 1 after the iteration process was completed successfully.

Example 1.

Let us consider the SFDE’s initial boundary value problem:

$$\frac{\partial U\left(x,t\right)}{\partial t}=A_1\frac{\partial U^{\beta }\left(x,t\right)}{\partial x^{\beta }}+B\left(x,t\right),$$

(19)

Example 2.

Let us consider the SFDE’s initial boundary value problem:

$$\frac{\partial U\left(x,t\right)}{\partial t}=\Gamma \left(1.2\right)x^{\beta }\frac{\partial U^{\beta }\left(x,t\right)}{\partial x^{\beta }}+3x^2\left(2x-1\right)e^{-t},$$

(20)

Based on

Table 1. Computing the results with $\beta =1.2, 1.5,$ and $1.8$.

M	Method	$\mathit{\beta }=1.2$			$\mathit{\beta }=1.5$			$\mathit{\beta }=1.8$
		K	Time	Max Error	K	Time	Max Error	K	Time	Max Error
128	FSPSOR	34	0.84	2.37 × 10−2	80	1.90	6.20 × 10−4	246	5.76	3.99 × 10−2
	HSPSOR	23	0.38	2.37 × 10−2	37	0.54	6.99 × 10−4	94	2.36	3.99 × 10−2
256	FSPSOR	67	5.33	2.44 × 10−2	211	17.84	5.69 × 10−4	806	67.75	3.97 × 10−2
	HSPSOR	34	2.73	2.44 × 10−2	94	6.90	6.21 × 10−4	303	34.65	3.97 × 10−2
512	FSPSOR	129	41.43	2.47 × 10−2	566	182.83	5.36 × 10−4	2635	843.91	3.96 × 10−2
	HSPSOR	67	22.65	2.47 × 10−2	246	86.09	5.69 × 10−4	988	421.58	3.96 × 10−2
1024	FSPSOR	278	472.35	2.49 × 10−2	1514	898.29	5.13 × 10−4	11,829	2099.87	3.95 × 10−2
	HSPSOR	141	206.58	2.49 × 10−2	655	434.72	5.36 × 10−4	5413	1033.78	3.95 × 10−2
2048	FSPSOR	608	1219.76	2.50 × 10−2	4052	4299.73	5.02 × 10−4	47,289	8852.28	3.93 × 10−2
	HSPSOR	305	608.80	2.50 × 10−2	2188	2133.43	5.13 × 10−4	23,143	4425.90	3.93 × 10−2

and

Table 2. Computing the results with $\beta =1.2, 1.5,$ and $1.8$.

M	Method	$\mathit{\beta }=1.2$			$\mathit{\beta }=1.5$			$\mathit{\beta }=1.8$
		K	Time	Max Error	K	Time	Max Error	K	Time	Max Error
128	FSPSOR	34	0.84	2.37 × 10−2	80	1.90	6.20 × 10−4	246	5.76	3.99 × 10−2
	HSPSOR	23	0.38	2.37 × 10−2	37	0.54	6.99 × 10−4	94	2.36	3.99 × 10−2
256	FSPSOR	67	5.33	2.44 × 10−2	211	17.84	5.69 × 10−4	806	67.75	3.97 × 10−2
	HSPSOR	34	2.73	2.44 × 10−2	94	6.90	6.21 × 10−4	303	34.65	3.97 × 10−2
512	FSPSOR	129	41.43	2.47 × 10−2	566	182.83	5.36 × 10−4	2635	843.91	3.96 × 10−2
	HSPSOR	67	22.65	2.47 × 10−2	246	86.09	5.69 × 10−4	988	421.58	3.96 × 10−2
1024	FSPSOR	278	472.35	2.49 × 10−2	1514	898.29	5.13 × 10−4	11,829	2099.87	3.95 × 10−2
	HSPSOR	141	206.58	2.49 × 10−2	655	434.72	5.36 × 10−4	5413	1033.78	3.95 × 10−2
2048	FSPSOR	608	1219.76	2.50 × 10−2	4052	4299.73	5.02 × 10−4	47,289	8852.28	3.93 × 10−2
	HSPSOR	305	608.80	2.50 × 10−2	2188	2133.43	5.13 × 10−4	23,143	4425.90	3.93 × 10−2

, it can be observed that for the five different mesh sizes and the three different $\beta$ used for the simulation of FSPSOR and HSPSOR to solve Example 1 and 2, respectively, the number of iterations and the execution time required by HSPSOR are always smaller than FSPSOR. These results illustrated the success of the Half-Sweep finite difference scheme with Caputo’s space fractional to approximate the solution of SFDE with lower computational complexity. The maximum errors produced by PSOR iteration to solve both Example 1 and 2 decreases with the increasing mesh sizes for $\beta =1.8$ and $1.5$. This means that with these values of $\beta$, the numerical solution of SFDE via HSPSOR can be computed nearly to its exact solution with a sufficiently large mesh size is used. However, for $\beta =1.2$, the maximum error is getting larger when the mesh is made to be narrower. This shows the sign of limitation from the use of PSOR iteration, and a modification on PSOR iteration needs to be conducted to handle this problem.

5. Conclusions

This paper describes the mathematical derivation of the implicit finite difference in Caputo’s approximation equations in which this approximation equation leads to a linear system. By imposing the iterative methods of the FSPSOR and the HSPSOR, based on observation of all experimental effects, it was evident that the number of iterations of the HSPSOR decreased by approximately 31.30–85.45 per cent compared with the iterative methods of the FSPSOR. Meanwhile, the execution time was much quicker, by around 41.18–95.33% more than the FSPSOR method. This implies that, relative to the FSPSOR iterative methods, the HSPSOR method needs the minimum number of iterations and computational time. It can be inferred from the precision of both iterative approaches that their numerical resolutions were in the acceptable domain. When the tabulated numerical results of the HSPSOR to solve the SFDE was compared against the results from the application of the HSPSOR on the time-fractional diffusion equation (TFDE) [24], it was found that the overall maximum errors produced by the HSPSOR in solving the SFDE were slightly greater than when the HSPSOR was used to solve the TFDE. A thorough study of the error of the method used will be conducted in the future. From this work, the capability of the half-sweep finite difference scheme to reduce the computational complexity for solving space- and time-fractional diffusion equations and the compatibility of the scheme with the PSOR iteration was shown.