
In this paper, a new modified hybrid explicit group (MHEG) iterative method is presented for the efficient and accurate numerical solution of a time-fractional diffusion equation in two space dimensions. The time fractional derivative is defined in the Caputo sense. In the proposed method, a Laplace transformation is used in the temporal domain, and, for the spatial discretization, a new finite difference scheme based on grouping strategy is considered. The unique solvability, unconditional stability and convergence are thoroughly proved by the matrix analysis method. Comparison of numerical results with analytical and other approximate solutions indicates the viability and efficiency of the proposed algorithm.
Citation: Fouad Mohammad Salama, Nur Nadiah Abd Hamid, Norhashidah Hj. Mohd Ali, Umair Ali. An efficient modified hybrid explicit group iterative method for the time-fractional diffusion equation in two space dimensions[J]. AIMS Mathematics, 2022, 7(2): 2370-2392. doi: 10.3934/math.2022134
[1] | Fouad Mohammad Salama, Nur Nadiah Abd Hamid, Umair Ali, Norhashidah Hj. Mohd Ali . Fast hybrid explicit group methods for solving 2D fractional advection-diffusion equation. AIMS Mathematics, 2022, 7(9): 15854-15880. doi: 10.3934/math.2022868 |
[2] | Muhammad Asim Khan, Norma Alias, Umair Ali . A new fourth-order grouping iterative method for the time fractional sub-diffusion equation having a weak singularity at initial time. AIMS Mathematics, 2023, 8(6): 13725-13746. doi: 10.3934/math.2023697 |
[3] | Mubashara Wali, Sadia Arshad, Sayed M Eldin, Imran Siddique . Numerical approximation of Atangana-Baleanu Caputo derivative for space-time fractional diffusion equations. AIMS Mathematics, 2023, 8(7): 15129-15147. doi: 10.3934/math.2023772 |
[4] | Fouad Mohammad Salama, Faisal Fairag . On numerical solution of two-dimensional variable-order fractional diffusion equation arising in transport phenomena. AIMS Mathematics, 2024, 9(1): 340-370. doi: 10.3934/math.2024020 |
[5] | Zihan Yue, Wei Jiang, Boying Wu, Biao Zhang . A meshless method based on the Laplace transform for multi-term time-space fractional diffusion equation. AIMS Mathematics, 2024, 9(3): 7040-7062. doi: 10.3934/math.2024343 |
[6] | Zunyuan Hu, Can Li, Shimin Guo . Fast finite difference/Legendre spectral collocation approximations for a tempered time-fractional diffusion equation. AIMS Mathematics, 2024, 9(12): 34647-34673. doi: 10.3934/math.20241650 |
[7] | Ailing Zhu, Yixin Wang, Qiang Xu . A weak Galerkin finite element approximation of two-dimensional sub-diffusion equation with time-fractional derivative. AIMS Mathematics, 2020, 5(5): 4297-4310. doi: 10.3934/math.2020274 |
[8] | Lei Ren . High order compact difference scheme for solving the time multi-term fractional sub-diffusion equations. AIMS Mathematics, 2022, 7(5): 9172-9188. doi: 10.3934/math.2022508 |
[9] | Bin Fan . Efficient numerical method for multi-term time-fractional diffusion equations with Caputo-Fabrizio derivatives. AIMS Mathematics, 2024, 9(3): 7293-7320. doi: 10.3934/math.2024354 |
[10] | Farman Ali Shah, Kamran, Zareen A Khan, Fatima Azmi, Nabil Mlaiki . A hybrid collocation method for the approximation of 2D time fractional diffusion-wave equation. AIMS Mathematics, 2024, 9(10): 27122-27149. doi: 10.3934/math.20241319 |
In this paper, a new modified hybrid explicit group (MHEG) iterative method is presented for the efficient and accurate numerical solution of a time-fractional diffusion equation in two space dimensions. The time fractional derivative is defined in the Caputo sense. In the proposed method, a Laplace transformation is used in the temporal domain, and, for the spatial discretization, a new finite difference scheme based on grouping strategy is considered. The unique solvability, unconditional stability and convergence are thoroughly proved by the matrix analysis method. Comparison of numerical results with analytical and other approximate solutions indicates the viability and efficiency of the proposed algorithm.
In recent years, fractional calculus has gained a great deal of attention from the research community due to its widespread applications in physics, chemistry, fluid mechanics, viscoelasticity, finance, control systems and other areas of science and engineering. The phenomena in the aforesaid fields can be modelled very successfully by means of equations containing fractional derivatives and fractional integrals; and therefore, the investigation of fractional differential equations has become a hot topic for many researchers. Fractional kinetic equations including Fokker-Planck equation, fractional diffusion equation, fractional cable equation and fractional advection-diffusion equation are proved to be powerful instruments for modelling transport dynamics in complex systems such as electrodiffusion of ions in spiny dendrites, transport of proteins molecules, movement of a solution in an aquifer, pollution of underground water, movement under the effect of optical tweezers and more [1,2,3,4,5]. To get more insight into the applications of fractional calculus, the interested reader can refer to the classical books in [6,7], besides the recent books in [8,9].
In this work, we study the following two-dimensional time-fractional diffusion equation with a non-homogeneous source term:
{C0Dαtu(x,y,t)=ax∂2u(x,y,t)∂x2+ay∂2u(x,y,t)∂y2+f(x,y,t),(x,y)∈Ω,0<t≤T,(1.1)u(x,y,t)=p(x,y,t),(x,y)∈∂Ω,0<t≤T,(1.2)u(x,y,0)=g(x,y),(x,y)∈Ω∪∂Ω,(1.3) |
where f(x,y,t), p(x,y,t) and g(x,y) are known functions, Ω=(0,L)×(0,L) is a bounded domain in R2, ∂Ω is the boundary, ax>0 and ay>0 are the diffusion coefficients, and α(0<α≤1) is the anomalous diffusion exponent. The term C0Dαtu(x,y,t) stands for the Caputo fractional derivative of the function u(x,y,t) that is defined as,
C0Dαtu(x,y,t)={1Γ(1−α)∫t0(t−ξ)−α∂u(x,y,ξ)∂ξdξ, 0<α<1,∂u(x,y,t)∂t, α=1. |
From the above definition, it can be observed that the fractional diffusion problems (1.1)–(1.3) corresponds to the classical diffusion model as α=1. Fractional diffusion equation is one of the most fundamental equations in the literature, and its capability to model many problems in science and engineering is a well established fact. Generally speaking, it has been widely used for describing random walks and modelling phenomena that are governed by anomalous diffusion in various fields such as porous systems, nuclear magnetic resonance and transport in fractal geometries [10,11,12]. Since most fractional differential equations are difficult to handle analytically, many researchers have resorted to numerical techniques, especially the finite difference method for solving fractional diffusion problems [4,10,11,13,14,15,16,17,18,19,20,21,22,23,24,25]. The complexity of fractional differential equations stems from the presence of non-local fractional derivatives that have the property of global dependence on time or space. This forms the principal obstacle to the development of efficient simulation algorithms in terms of CPU time and memory usage, especially for multi-dimensional problems [26,27,28]. To surmount this issue, techniques including fast Poisson solver [14], parallel implementation [15,29], preconditioning [15] and multigrid method [30] were suggested in the literature. Therefore, the question of how we can remove the fractional derivative from our computations to reduce the complexity and establish efficient solution algorithms is considered a big question arising in the numerical simulations of fractional differential equations. In this regard, Ren et al. [31] have introduced a relatively new approach for solving Caputo-type fractional differential equations. In this approach, a Laplace transform technique is proposed to approximate the fractional differential equation by its corresponding integer-order differential equation, which can be solved with less effort by using a suitable numerical or analytical method (see [32,33]). Recently in [34], a new hybrid standard point (HSP) iterative method based on a combination of the Laplace transform technique and implicit difference scheme has been proposed to solve the time-fractional diffusion problems (1.1)–(1.3). The authors showed that their method produces accurate numerical solutions with an optimal computational complexity of O(N) and optimal memory requirement of O(Ms), where N and Ms are the total number of time levels and spatial unknowns, respectively.
It is well known that the explicit group iterative methods [35,36,37,38,39] based on finite difference approximations are proposed for the numerical solutions of integer-order differential equations. From one side, explicit group methods are unconditionally stable as the conventional implicit schemes. On the other side, they reduce the number of spatial unknowns taken in the iterative process and thus lead to computationally efficient algorithms. Although the grouping methods have been successfully employed to abroad spectrum of classical differential models, very little work has been reported to deal with fractional differential equations using grouping techniques (see [40,41,42,43]). In addition, the handling of nonlinear and variable order fractional differential equations by explicit group methods, besides their parallel implementation, is another topic that is still at its infancy. Therefore, this subject needs a major development.
The main purpose of this paper is to combine the Laplace transform technique with an explicit group scheme to solve the two-dimensional time-fractional diffusion Eqs (1.1)–(1.3). The unconditional stability, convergence and solvability of the resulting method, namely the modified hybrid explicit group (MHEG) method are rigorously proved. To illustrate the efficiency and feasibility of the proposed method, the fast HSP iterative method based on the recent work in [34] is also presented.
The rest of this paper is organized as follows. In Section 2, we briefly review the existing HSP iterative method [34] for solving the problems (1.1)–(1.3). In Section 3, the MHEG iterative method is proposed, and its unconditional stability, convergence and solvability are rigorously proved in Section 4. In Section 5, numerical simulations are performed to indicate the efficiency of the proposed method and support our theoretical analysis. Finally, this work ends with a brief summary in Section 6.
Since the Caputo fractional derivative in (1.1) is non-local and has the character of history dependence, conventional difference schemes (such as an implicit or explicit scheme with a particular discretiztion formula for the fractional derivative) for the problems (1.1)–(1.3) result in costly simulations in terms of CPU time and memory consumption. In order to surmount the computational challenge, and utilizing the Laplace transform technique, the Caputo's time fractional derivative is approximated as follows [34]:
C0Dαtu(x,y,t)≈α∂u(x,y,t)∂t+(1−α)[u(x,y,t)−u(x,y,0)]. | (2.1) |
By substituting (2.1) into (1.1), the original time-fractional diffusion Eq (1.1) is reduced to its corresponding integer-order partial differential equation (PDE) with the following initial and boundary conditions [34]:
{∂u(x,y,t)∂t=Ax∂2u(x,y,t)∂x2+Ay∂2u(x,y,t)∂y2−(r−1)u(x,y,t)+(r−1)g(x,y)+rf(x,y,t),(x,y)∈Ω,0<t≤T,(2.2)u(x,y,t)=p(x,y,t),(x,y)∈∂Ω,0<t≤T,(2.3)u(x,y,0)=g(x,y),(x,y)∈Ω∪∂Ω,(2.4) |
where Ax=ax/α, Ay=ay/α and r=1/α are positive constants. Next, we let
xi=ih, i=0,1,…,M,yj=jh, j=0,1,…,M,tk=kτ, k=0,1,…,N, |
where h=L/M and τ=T/N are the uniform space and time step sizes, respectively. We also introduce the grid functions given by
uki,j=u(xi,yj,tk),u0i,j=g(xi,yj),fki,j=f(xi,yj,tk),0≤i,j≤M, 0≤k≤N. |
Here, the notations δtUk+1i,j, δ2xUk+1i,j and δ2yUk+1i,j are defined as follows:
δtUk+1i,j=1τ(Uk+1i,j−Uki,j)+O(τ),δ2xUk+1i,j=1h2(Uk+1i+1,j−2Uk+1i,j+Uk+1i−1,j)+O(h2),δ2yUk+1i,j=1h2(Uk+1i,j+1−2Uk+1i,j+Uk+1i,j−1)+O(h2). | (2.5) |
Now, applying the difference operators in (2.5) to the approximating PDEs (2.2)–(2.4) and disregarding the truncation errors, we obtain the implicit difference scheme in the following form:
{uk+1i,j=1(1+(r−1)τ+2q1+2q2)[q1(uk+1i+1,j+uk+1i−1,j)+q2(uk+1i,j+1+uk+1i,j−1),+uki,j+(r−1)τu0i,j+rτfk+1i,j],1≤i,j≤M−1,0≤k≤N−1,(2.6)uki,j|∂Ω=p(xi,yj,tk),0≤k≤N,(2.7)u0i,j=g(xi,yj),0≤i,j≤M,(2.8) |
where q1=Axτ/h2, q2=Ayτ/h2 and uki,j is the numerical approximation of the function Uki,j after neglecting the truncation errors.
Considering the fact that the spatial and temporal discretizations usually lead to large and sparse systems of linear equations to be solved, iterative solvers, such as Gauss-Seidel method, are more practical and economic than the direct solvers. Here, and at each time level, the HSP method functions through the iterative evaluation of solutions at all mesh points ⧫ shown in Figure 1 by using Eq (2.7) until certain convergence is achieved. Once the converged solution values uk+1=(uk+111,uk+11,2,…,uk+11,M−1,…,uk+1M−1,1,uk+1M−1,2,…,uk+1M−1,M−1)T are attained, they are utilized as an initial guess for the next time level uk+2=(uk+211,uk+21,2,…,uk+21,M−1,…,uk+2M−1,1,uk+2M−1,2,…,uk+2M−1,M−1)T. Such process is repeated until it reaches the targeted time level. The HSP iterative method has been shown to reduce the memory requirement, computational complexity as well as the CPU time greatly when compared to the fractional standard point (FSP) iterative methods derived based on the conventional implicit difference schemes. For extra details, please refer to [34]. Seeking for more efficient solution algorithm, we propose our method in the next section.
Consider the standard mesh of points with the spatial step size 2h=2L/M illustrated in Figure 2. We begin this section by discretizing the derivatives of the space variables in (2.2) around these 2h spaced points. Recalling the δt definition, the spatial difference operators δ2x and δ2y are defined as follows:
δtUk+1i,j=1τ(Uk+1i,j−Uki,j)+O(τ),δ2xUk+1i,j=1h2(Uk+1i+2,j−2Uk+1i,j+Uk+1i−2,j)+O(h2),δ2yUk+1i,j=1h2(Uk+1i,j+2−2Uk+1i,j+Uk+1i,j−2)+O(h2). | (3.1) |
Applying the difference operators in (3.1) to the approximating PDE (2.2), we obtain
uk+1i,j−uki,jτ=Ax(uk+1i+2,j−2uk+1i,j+uk+1i−2,j4h2)+Ay(uk+1i,j+2−2uk+1i,j+uk+1i,j−24h2)−(r−1)uk+1i,j+(r−1)u0i,j+rfk+1i,j+O(τ+h2). | (3.2) |
After simplification and imposing the initial and boundary conditions, the following implicit scheme is attained:
{uk+1i,j=1(1+(r−1)τ+q1/2+q2/2)[q14(uk+1i+2,j+uk+1i−2,j)+q24(uk+1i,j+2+uk+1i,j−2),+uki,j+(r−1)τu0i,j+rτfk+1i,j],2≤i,j≤M−2,0≤k≤N−1,(3.3)uki,j|∂Ω=p(xi,yj,tk),0≤k≤N,(3.4)u0i,j=g(xi,yj),0≤i,j≤M,(3.5) |
where uki,j is the approximate numerical solution of Uki,j after omitting the truncation error terms. The MHEG method can now be constructed by applying Eq (3.3) to a group of four mesh points of the solution domain. This will result in the following (4×4) system of equations:
(D−q1/40−q2/4−q1/4D−q2/400−q2/4D−q1/4−q2/40−q1/4D)(uk+1i,juk+1i+2,juk+1i+2,j+2uk+1i,j+2)=(rhsi,jrhsi+2,jrhsi+2,j+2rhsi,j+2), | (3.6) |
where
D=1+(r−1)τ+q1/2+q2/2,rhsi,j=q14uk+1i−2,j+q24uk+1i,j−2+uki,j+(r−1)τu0i,j+rτfk+1i,j,rhsi+2,j=q14uk+1i+4,j+q24uk+1i+2,j−2+uki+2,j+(r−1)τu0i+2,j+rτfk+1i+2,j,rhsi+2,j+2=q14uk+1i+4,j+2+q24uk+1i+2,j+4+uki+2,j+2+(r−1)τu0i+2,j+2+rτfk+1i+2,j+2,rhsi,j+2=q14uk+1i−2,j+2+q24uk+1i,j+4+uki,j+2+(r−1)τu0i,j+2+rτfk+1i,j+2. |
The coefficients matrix in (3.6) can be inverted to get the MHEG equation given as below:
(uk+1i,juk+1i+2,juk+1i+2,j+2uk+1i,j+2)=1η(η1η2η3η4η2η1η4η3η3η4η1η2η4η3η2η1)(rhsi,jrhsi+2,jrhsi+2,j+2rhsi,j+2), | (3.7) |
where
η=1256(4+4(r−1)τ+q1+q2)(4+4(r−1)τ+3q1+q2)(4+4(r−1)τ+q1+3q2)×(4+4(r−1)τ+3q1+3q2),η1=132(2+2(r−1)τ+q1+q2)(16+16q1+16q2+8q1q2+3q21+3q22+16q1(r−1)τ+16q2(r−1)τ+32(r−1)τ+16(r−1)2τ2),η2=164q1(16+16q1+16q2+8q1q2+3q21+5q22+16q1(r−1)τ+16q2(r−1)τ+32(r−1)τ+16(r−1)2τ2),η3=116q1q2(2+2(r−1)τ+q1+q2),η4=164q2(16+16q1+16q2+8q1q2+5q21+3q22+16q1(r−1)τ+16q2(r−1)τ+32(r−1)τ+16(r−1)2τ2). |
Back to Figure 2, it can be seen that the discretized solution domain comprises three different types of mesh points (⧫, ◻, ○). Here, the MHEG method functions through the iterative evaluation of solutions at ⧫ mesh points by utilizing Eq (3.7) until a certain convergence criterion is met. Afterwards, the solution values at the remaining mesh points of types ◻ and ○ are computed directly once in a particular sequence. The computational cost of the iterative method still depends on the number of unknowns involved in the iteration process. For the MHEG method, it is obvious that the number of mesh points taken in the iteration process is lesser than that of the prespecified HSP method, which is expected to speed up the convergence rate of the proposed method. For convenience, the MHEG method in combination with the Gauss-Seidel iterative solver is elaborated in Algorithm 1.
Algorithm 1: MHEG method utilizing the Gauss-Seidel iterative solver |
1: Branch the mesh points of the discretized solution domain into the following categories:
● Iterative points of type ⧫. ● Direct points of types ◻ and ○. 2: Arrange all the ⧫ points into groups of four mesh points as depicted in Figure 2. 3: Set an initial guess for the solution at the current time level. 4: For each four-point group at the current time level, iterate the intermediate solutions at ⧫ mesh points utilizing (ˆuk+1,n+1i,jˆuk+1,n+1i+2,jˆuk+1,n+1i+2,j+2ˆuk+1,n+1i,j+2)=1η(η1η2η3η4η2η1η4η3η3η4η1η2η4η3η2η1)(rhsi,jrhsi+2,jrhsi+2,j+2rhsi,j+2), where rhsi,j=q14uk+1i−2,j+q24uk+1i,j−2+uki,j+(r−1)τu0i,j+rτfk+1i,j,rhsi+2,j=q14uk+1i+4,j+q24uk+1i+2,j−2+uki+2,j+(r−1)τu0i+2,j+rτfk+1i+2,j,rhsi+2,j+2=q14uk+1i+4,j+2+q24uk+1i+2,j+4+uki+2,j+2+(r−1)τu0i+2,j+2+rτfk+1i+2,j+2,rhsi,j+2=q14uk+1i−2,j+2+q24uk+1i,j+4+uki,j+2+(r−1)τu0i,j+2+rτfk+1i,j+2, and perform the Gauss-Seidel solver (uk+1,n+1i,juk+1,n+1i+2,juk+1,n+1i+2,j+2uk+1,n+1i,j+2)=ω(ˆuk+1,n+1i,jˆuk+1,n+1i+2,jˆuk+1,n+1i+2,j+2ˆuk+1,n+1i,j+2)+(1−ω)(uk+1,ni,juk+1,ni+2,juk+1,ni+2,j+2uk+1,ni,j+2), where n is the iteration number and ω=1 is the relaxation parameter of the Gauss-Seidel solver. η, η1, η2, η3 and η4 are as defined before in this section. 5: Test the convergence. If convergence is achieved, go to step 6. Otherwise, step 4 is repeated until convergence is attained. 6: The solutions on the rest of mesh points of types ◻ and ○ are evaluated directly once in the following manner: ● For ◻ points, the approximating PDE (2.2) is discretized by using skewed finite difference operators. Such difference operators are derived by rotating the standard mesh an angle 45∘ clockwise and applying the Taylor series expansion afterwards [45]. This will lead to the following skewed difference scheme for (2.2): uk+1i,j=1(1+(r−1)τ+q1+q2)[q12(uk+1i+1,j−1+uk+1i−1,j+1)+q22(uk+1i+1,j+1+uk+1i−1,j−1)+uki,j+(r−1)τu0i,j+rτfk+1i,j]. ● For ○ points, the difference scheme formula defined by (2.6) is employed. 7: If the final time level is reached, go to step 8. Otherwise, adopt the solution values of the last time level as an initial guess for the next time level and go to step 4. 8: Display the numerical solutions. |
In this section, we turn to consider the stability and convergence of the MHEG scheme defined by (3.7) by using the matrix analysis method. For the convenience of our analysis, three remarks are given at first.
Remark 4.1. Let z∈N. The symbol ‖Az×z‖ denotes the infinity norm of matrix A=[ai,j] which is defined as,
‖Az×z‖∞=max1≤i≤z{z∑j=1|ai,j|}. |
Remark 4.2. Az×z is called a stricatly diagonally dominant (SDD) matrix if |ai,i|>ri(A), where ri(A) is the i-th deleted absolute row sum (i.e., ri(A)=∑zj≠i,j=1|ai,j|, 1≤i≤z).
Remark 4.3. Given that Az×z is a SDD matrix, then Az×z is non-singular and the following norm's bound does hold [44]:
‖A−1‖∞≤1min1≤i≤M{|ai,i|−ri(A)}. |
Here, we analyze the stability of the MHEG scheme (3.7). For simplicity and without loss of generality, we assume that q1=q2=q=τ/h2. Indeed, the MHEG scheme (3.7) can be represented in the following matrix form:
Auk+1=Buk+Cu0+b,0≤k≤N−1, | (4.1) |
where ukis an (M−2)24-dimensional vector expressed in the block form
uk=(uk1,uk2,…,uk(M−2)216)T, ukl=(uki,j,uki+2,j,uki+2,j+2,uki,j+2)T, 1≤l≤(M−2)216, |
and
A=(J1J2J3J1J2⋱J3J1J2J3J1),B=(H1H1⋱H1H1),b=(S1S1⋮S1S1), |
C=(P1P1⋱P1P1),J1=(Q1Q3Q2Q1Q3⋱Q2Q1Q3Q2Q1),J2=(Q5Q5⋱Q5Q5),J3=(Q4Q4⋱Q4Q4),H1=(I4I4⋱I4I4),P1=(T1T1⋱T1T1),S1=(G1G1⋮G1G1),Q1=(1+(r−1)τ+q−q/40−q/4−q41+(r−1)τ+q−q/400−q/41+(r−1)τ+q−q/4−q/40−q/41+(r−1)τ+q),Q2=(000−q/400−q/4000000000),Q3=(000000000−q/400−q/4000),Q4=(0−q/4000000000000−q/40),Q5=(0000−q/4000000−q/40000), |
T1=((r−1)τ0000(r−1)τ0000(r−1)τ0000(r−1)τ),I4=(1000010000100001),G1=rτ(fi,jfi+2,jfi+2,j+2fi,j+2). |
Taking note of the above block matrices, the matrix representation in (4.1) can be rewritten in the more elaborated general form
[A(M−2)24×(M−2)24]uk+1=[B(M−2)24×(M−2)24]uk+[C(M−2)24×(M−2)24]u0+b. | (4.2) |
Next, we prove the following result:
Theorem 4.1. The MHEG scheme defined by (4.2) is unconditionally stable.
Proof. Suppose uk and ˜uk are respectively the exact and approximate solutions of (4.2). Then the error at time level k is expressed as ek=uk−˜uk. From Remarks 4.2 and 4.3, it follows that A is invertible and therefore Eq (4.2) can be rewritten as
uk+1=A−1Buk+A−1Cu0+A−1b,0≤k≤N−1. | (4.3) |
The error function ek satisfies (4.3) and leads to the following round-off error equation:
ek+1=A−1Bek+A−1Ce0,0≤k≤N−1, | (4.4) |
where
ek+1=(ek+10ek+10⋮ek+10),ek+10=(ψk+11ψk+12⋮ψk+1(M−2)216),ψk+1=(ψk+1i,jψk+1i+2,jψk+1i+2,j+2ψk+1i,j+2), |
and ψk+1i,j=uk+1i,j−˜uk+1.
In order to demonstrate the stability, we shall prove that ‖ek+1‖≤‖e0‖ for all 0≤k≤N−1. To this end, the mathematical induction will be used. For k=0, we obtain
e1=A−1Be0+A−1Ce0. |
Using the fact that the matrix and the vector infinity norms are consistent, yields
‖e1‖=‖A−1Be0+A−1Ce0‖≤‖A−1B‖‖e0‖+‖A−1C‖‖e0‖≤‖A−1‖‖B‖‖e0‖+‖A−1‖‖C‖‖e0‖=(‖B‖+‖C‖)‖A−1‖‖e0‖. |
Define ri(A) as in Remark 4.2. Then using Remark 4.3, we get
‖e1‖≤(‖B‖+‖C‖)min1≤i≤M{|ai,i|−ri(A)}‖e0‖1+(r−1)τ|1+(r−1)τ+q|−(|−q/4|+|−q/4|+|−q/4|+|−q/4|)=1+(r−1)τ1+(r−1)τ‖e0‖=‖e0‖. |
Next, suppose that
‖es+1‖≤‖e0‖,s=1,2,…,k−1. | (4.5) |
We will show the above inequality is true for s=k. Noticing (4.4) and (4.5), we get
‖ek+1‖=‖A−1Bek+A−1Ce0‖≤‖A−1‖‖B‖‖ek‖+‖A−1‖‖C‖‖e0‖≤‖A−1‖‖B‖‖e0‖+‖A−1‖‖C‖‖e0‖≤(‖B‖+‖C‖)min1≤i≤M{|ai,i|−ri(A)}‖e0‖=1+(r−1)τ|1+(r−1)τ+q|−(|−q/4|+|−q/4|+|−q/4|+|−q/4|)=1+(r−1)τ1+(r−1)τ‖e0‖=‖e0‖. |
Thus, the proof is completed.
We now analyze the convergence of the MHEG scheme (3.7) by following a similar approach to that in the previous subsection. As a start, we suppose that the truncation errors on each four-point group of mesh points at any time level are represented by the block vector of the following form:
Rk+1=(Rk+11,Rk+12,…,Rk+1(M−2)216)T,Rn+1l=(Rk+1i,j,Rk+1i+2,j,Rk+1i+2,j+2,Rk+1i,j+2)T, 1≤l≤(M−2)216. |
Noticing Eq (3.2) and since i, j and k are finite, there is a positive constant C∗ such that
‖Rk+1‖≤C∗(τ+h2),0≤k≤N−1. | (4.6) |
Subtracting (4.2) from the following equation:
AUk+1=BUk+CU0+b+Rk+1, |
will result in error equation given by
AEk+1=BEk+CE0+Rk+1, | (4.7) |
where
Ek+1=(Ek+10Ek+10⋮Ek+10),Ek+10=(ϕk+11ϕk+12⋮ϕk+1(M−2)216),ϕk+1=(ϕk+1i,jϕk+1i+2,jϕk+1i+2,j+2ϕk+1i,j+2), |
and ϕk+1i,j=Uk+1i,j−uk+1i,j.
The convergence property is given in the next theorem.
Theorem 4.2. Assume problem (1.1) has a smooth solution u(x,y,t). Then, the MHEG scheme defined by (3.7) is convergent and the convergence order is O(τ+h2).
Proof. We will use mathematical induction for the proof. For k=0 and using that E0=0, we obtain
E1=A−1R1. |
In view of of Remark 4.3 and utilizing (4.6), yields
‖E1‖=‖A−1R1‖≤‖A−1‖‖R1‖≤1min1≤i≤M{|ai,i|−ri(A)}C∗(τ+h2)=11+(r−1)τC∗(τ+h2)=C0(τ+h2), |
where C0=C∗/(1+(r−1)τ).
∴‖E1‖≤C0(τ+h2). |
Now, assume that
Es+1≤Cs(τ+h2),s=1,2,…,k−1. | (4.8) |
We prove the above result is true for s=k. Noticing (4.7) and (4.8), we get
‖Ek+1‖=‖A−1BEk+A−1Rk+1‖≤‖A−1‖‖B‖‖Ek‖+‖A−1‖‖Rk+1‖≤1min1≤i≤M{|ai,i|−ri(A)}[Ck−1(τ+h2)+C∗(τ+h2)]=11+(r−1)τ(Ck−1+C∗)(τ+h2)=Ck(τ+h2), |
where Ck=Ck−1+C∗ as limn→∞τ=0.
∴‖En+1‖≤Ck(τ+h2),0≤n≤N−1. |
Hence, the proof is completed by induction.
Theorem 4.3 The MHEG scheme defined by (3.7) is uniquely solvable.
Proof. In view of (4.1), The coefficients matrix A of the MHEG scheme (3.7) is an SDD matrix. Consequently, and based on Remark 4.3, it follows that A is a non-singular matrix. This proves the existence and uniqueness of the solution of the MHEG scheme.
In this section, we report on numerical simulations for (1.1)–(1.3). In this endeavor, two test problems with known exact solutions are given to demonstrate the accuracy and efficiency of the newly developed MHEG method. For comparison purposes, the fast HSP method proposed recently in [34] is applied to solve the two test problems. All simulations using the both methods are performed in Mathematica 11.3 software and run on a computer with an Intel (R) Core (TM) i7-8550U CPU and 8.00 GB of RAM. In these simulations, the Gauss-Seidel iterative solver with a stopping criterion selected to be 10−5 are employed to generate the corresponding numerical results. To get more insight into these results, and based on the total number of arithmetic operations to be implemented before and after the convergence process, an analysis on the computational complexity of the MHEG and HSP iterative methods is established and presented in Table 1.
Method | Per iteration | After convergence | Total operations |
HSP | 13σ2∗Ite | - | 13σ2∗Ite |
MHEG | 4.5(σ−1)2∗Ite | 3.25(3σ2+2σ−1) | 4.5(σ−1)2∗Ite+3.25(3σ2+2σ−1) |
Test problem 5.1. We consider the following diffusion equation of fractional order:
C0Dαtu(x,y,t)=∂2u(x,y,t)∂x2+∂2u(x,y,t)∂y2+(2t2−αΓ(3−α)+2t2)sin(x)sin(y), |
subject to the boundary and initial conditions
u(x,y,t)|∂Ω=t2sin(x)sin(y),u(x,y,0)=0, |
in the spatial domain (0,L)×(0,L)=(0,1)2, T=1, and the exact analytic solution is u(x,y,t)=t2sin(x)sin(y).
In this test, numerical outputs including maximum absolute error Emax=max1≤i,j≤M−1|u(xi,yj,tN)−uNi,j|, CPU time (in seconds), number of iterations (Ite) as well as total number of arithmetic operations are computed for different values of τ, h and α and listed in Table 2. It can be observed that the MHEG iterative method results in significantly faster simulations without compromising too much of the accuracy when compared to the HSP iterative method. To illustrate this further, Figures 3 and 4 show respectively the graphs of CPU time and total operations of both methods against various mesh sizes. Form these figures, it is apparent that the shape of the computing time plots is consistent with the shape of the total computing effort plots. This indicates that the experimental results are in good agreement with our theoretical analysis. The comparison of the exact and the numerical solutions besides the graphical error representation when τ=h=1/18 and α=0.55 are depicted in Figures 5 and 6, respectively. A comparison of the numerical errors of the test problems for different fractional orders shows that the accuracy of the proposed method is greatest at α=0.55, followed by α=0.75 and α=0.95. It is clear that the numerical solution matches well with the exact solution.
α | τ=h | Method | CPU time (sec) | Ite | Total operations | Emax |
0.55 | ||||||
1/6 | HSP | 0.062 | 27 | 8,775 | 1.6737E-03 | |
MHEG | 0.015 | 2 | 417 | 1.7983E-03 | ||
1/10 | HSP | 0.968 | 56 | 58,968 | 9.5405E-04 | |
MHEG | 0.140 | 11 | 4,013 | 1.0764E-03 | ||
1/14 | HSP | 3.937 | 86 | 188,942 | 6.0571E-04 | |
MHEG | 0.359 | 16 | 12,097 | 7.6144E-04 | ||
1/18 | HSP | 10.062 | 117 | 439,569 | 3.8681E-04 | |
MHEG | 0.734 | 22 | 28,269 | 5.7584E-04 | ||
0.75 | ||||||
1/6 | HSP | 0.062 | 26 | 8,450 | 2.3916E-03 | |
MHEG | 0.015 | 2 | 417 | 2.4750E-03 | ||
1/10 | HSP | 0.937 | 53 | 55,809 | 1.4477E-03 | |
MHEG | 0.109 | 10 | 3,725 | 1.5593E-03 | ||
1/14 | HSP | 3.109 | 80 | 175,760 | 9.9823E-04 | |
MHEG | 0.250 | 15 | 11,449 | 1.1488E-03 | ||
1/18 | HSP | 8.062 | 107 | 401,999 | 6.9546E-04 | |
MHEG | 0.640 | 21 | 27,117 | 9.0926E-04 | ||
0.95 | ||||||
1/6 | HSP | 0.062 | 25 | 8,125 | 2.9204E-03 | |
MHEG | 0.015 | 2 | 417 | 2.9789E-03 | ||
1/10 | HSP | 0.656 | 50 | 52,650 | 1.7502E-03 | |
MHEG | 0.062 | 10 | 3,725 | 1.8611E-03 | ||
1/14 | HSP | 2.562 | 75 | 164,775 | 1.1935E-03 | |
MHEG | 0.203 | 14 | 10,801 | 1.3410E-03 | ||
1/18 | HSP | 7.937 | 100 | 375,700 | 8.2842E-04 | |
MHEG | 0.515 | 19 | 24,813 | 1.0461E-03 |
Test problem 5.2. We consider another diffusion equation of fractional order as follows:
C0Dαtu(x,y,t)=∂2u(x,y,t)∂x2+∂2u(x,y,t)∂y2+(2t2−αΓ(3−α)−2t2)e(x+y), |
with the boundary and initial conditions
u(x,y,t)|∂Ω=t2e(x+y),u(x,y,0)=0, |
in the spatial domain Ω=(0,1)2, T=1, and the exact analytic solution u(x,y,t)=t2e(x+y).
The numerical results of solving this problem for various values of τ, h and α are presented in Table 3. Figure 7 highlights the CPU time plots versus several mesh sizes. In view of Table 3 and Figure 7, it can be seen that the MHEG iterative method reduces the CPU time greatly without jeopardizing the accuracy of numerical solutions compared to the HSP iterative method. This can be attributed to the reduction in the number of executed arithmetic operations that leads to a lower computational complexity as shown in Figure 8. Figure 9 introduces the comparison between the exact and the numerical solutions, whereas Figure 10 presents the 3D plot of the maximum errors when τ=h=1/30 and α=0.75. Again, these demonstrate the effectiveness of the proposed MHEG method. Table 4 illustrates the experimental convergence order in the temporal direction using the formula ρ(τ,h)=log2(Emax(2τ,h)/Emax(τ,h)). It can be observed that the performance is in good agreement with the theoretical analysis.
α | τ=h | Method | CPU time (sec) | Ite | Total operations | Emax |
0.55 | ||||||
1/18 | HSP | 12.671 | 170 | 638,690 | 6.3352E-03 | |
MHEG | 0.859 | 29 | 36,333 | 6.8017E-03 | ||
1/22 | HSP | 31.906 | 221 | 1,266,993 | 5.0040E-03 | |
MHEG | 2.359 | 38 | 72,833 | 5.5331E-03 | ||
1/26 | HSP | 64.234 | 271 | 2,201,875 | 4.0351E-03 | |
MHEG | 3.390 | 46 | 125,485 | 4.6340E-03 | ||
1/30 | HSP | 118.625 | 320 | 3,498,560 | 3.2629E-03 | |
MHEG | 7.937 | 55 | 202,425 | 3.8992E-03 | ||
0.75 | ||||||
1/18 | HSP | 12.640 | 156 | 586,092 | 1.0239E-02 | |
MHEG | 0.625 | 27 | 34,029 | 1.0665E-02 | ||
1/22 | HSP | 28.828 | 199 | 1,140,867 | 8.4919E-03 | |
MHEG | 1.343 | 34 | 65,633 | 9.0001E-03 | ||
1/26 | HSP | 57.593 | 242 | 1,966,250 | 7.1948E-03 | |
MHEG | 3.359 | 42 | 115,117 | 7.8080E-03 | ||
1/30 | HSP | 102.688 | 284 | 3,104,972 | 6.1807E-03 | |
MHEG | 7.328 | 49 | 181,257 | 6.9414E-03 | ||
0.95 | ||||||
1/18 | HSP | 12.359 | 143 | 537,251 | 1.1950E-02 | |
MHEG | 0.562 | 25 | 31,725 | 1.2403E-02 | ||
1/22 | HSP | 25.968 | 182 | 1,043,406 | 9.7841E-03 | |
MHEG | 1.328 | 31 | 60,233 | 1.0285E-02 | ||
1/26 | HSP | 51.796 | 219 | 1,779,375 | 8.1651E-03 | |
MHEG | 3.203 | 38 | 104,749 | 8.8048E-03 | ||
1/30 | HSP | 94.796 | 255 | 2,787,915 | 6.9147E-03 | |
MHEG | 6.890 | 44 | 163,617 | 7.6882E-03 |
Test problem | τ | Emax | Computational order |
1 | 1/10 | 9.9154E-04 | 0.9835 |
1/20 | 5.0149E-04 | ||
2 | 1/10 | 1.1829E-02 | 0.9658 |
1/20 | 6.0570E-03 |
In this paper, a new MHEG iterative method for solving the two-dimensional diffusion equation with time fractional derivative has been developed. The unique solvability, unconditional stability and convergence are proved by the matrix analysis method. Moreover, the feasibility and efficiency of the proposed numerical method are confirmed through the computational outputs drawn from the conducted numerical simulations. The advantages of being uncomplicated, easy to implement and diminishing the amount of computational complexity indicate that the proposed method and numerical analysis in this article can be extended to solve other types of non-linear and variable order fractional differential equations. In addition, the parallel implementation of the proposed MHEG method is an interesting line to study which can be considered as a part of future work.
This research was supported by Ministry of Higher Education Malaysia, Fundamental Research Grant Scheme with Project Code: FRGS/1/2019/STG06/USM/01/1. We would like to express our gratitude to the anonymous reviewers for their valuable comments.
The authors have no conflict of interest.
[1] |
H. R. Ghehsareh, A. Zaghian, S. M. Zabetzadeh, The use of local radial point interpolation method for solving two-dimensional linear fractional cable equation, Neural Comput. Appl., 29 (2018), 745–754. doi: 10.1007/s00521-016-2595-y. doi: 10.1007/s00521-016-2595-y
![]() |
[2] |
Y. L. Zhao, T. Z. Huang, X. M. Gu, W. H. Luo, A fast second-order implicit difference method for time-space fractional advection-diffusion equation, Numer. Func. Anal. Opt., 41 (2020), 257–293. doi: 10.1080/01630563.2019.1627369. doi: 10.1080/01630563.2019.1627369
![]() |
[3] |
M. Hussain, S. Haq, Weighted meshless spectral method for the solutions of multi-term time fractional advection-diffusion problems arising in heat and mass transfer, Int. J. Heat Mass Tran., 129 (2019), 1305–1316. doi: 10.1016/j.ijheatmasstransfer.2018.10.039. doi: 10.1016/j.ijheatmasstransfer.2018.10.039
![]() |
[4] |
M. Abbaszadeh, A. Mohebbi, A fourth-order compact solution of the two dimensional modified anomalous fractional sub-diffusion equation with a nonlinear source term, Comput. Math. Appl., 66 (2013), 1345–1359. doi: 10.1016/j.camwa.2013.08.010. doi: 10.1016/j.camwa.2013.08.010
![]() |
[5] |
G. M. Zaslavsky, Chaos, fractional kinetics, and anomalous transport, Phys. Rep., 371 (2002), 461–580. doi: 10.1016/S0370-1573(02)00331-9. doi: 10.1016/S0370-1573(02)00331-9
![]() |
[6] | S. G. Samko, A. A. Kilbas, O. I. Marichev, Fractional integrals and derivatives: Theory and applications, Yverdon: Gordon and Breach, 1993. |
[7] | I. Podlubny, Fractional differential equations, mathematics in science and engineering, New York: Academic Press, 1999. |
[8] | B. Guo, X. Pu, F. Huang, Fractional partial differential equations and their numerical solutions, Singapore: World Scientific, 2015. |
[9] | A. Kochubei, Y. Luchko, V. E. Tarasov, I. Petra, Handbook of fractional calculus with applications, Berlin: De Gruyter Grand Forks, 2019. |
[10] |
J. Y. Shen, Z. Z. Sun, R. Du, Fast finite difference schemes for time fractional diffusion equations with a weak singularity at initial time, E. Asian J. Appl. Math., 8 (2018), 834–858. doi: 10.4208/eajam.010418.020718. doi: 10.4208/eajam.010418.020718
![]() |
[11] |
A. Chen, C. Li, A novel compact adi scheme for the time-fractional subdiffusion equation in two space dimensions, Int. J. Comput. Math., 93 (2016), 889–914. doi: 10.1080/00207160.2015.1009905. doi: 10.1080/00207160.2015.1009905
![]() |
[12] |
G. H. Gao, Z. Z. Sun, Y. N. Zhang, A finite difference scheme for fractional sub-diffusion equations on an unbounded domain using artificial boundary conditions, J. Comput. Phys., 231 (2012), 2865–2879. doi: 10.1016/j.jcp.2011.12.028. doi: 10.1016/j.jcp.2011.12.028
![]() |
[13] |
M. Tamsir, N. Dhiman, D. Nigam, A. Chauhan, Approximation of Caputo time-fractional diffusion equation using redefined cubic exponential B-spline collocation technique, AIMS Mathematics, 6 (2021), 3805–3820. doi: 10.3934/math.2021226. doi: 10.3934/math.2021226
![]() |
[14] | J. Shen, X. M. Gu, Two finite difference methods based on an H2N2 interpolation for two-dimensional time fractional mixed diffusion and diffusion-wave equations, Discrete Cont. Dyn. B, 2021. doi: 10.3934/dcdsb.2021086. |
[15] |
X. M. Gu, Y. L. Zhao, X. L. Zhao, B. Carpentieri, Y. Y. Huang, A note on parallel preconditioning for the all-at-once solution of Riesz fractional diffusion equations, Numer. Math. Theor. Meth. Appl., 14 (2021), 893–919. doi: 10.4208/nmtma.OA-2020-0020. doi: 10.4208/nmtma.OA-2020-0020
![]() |
[16] |
X. M. Gu, H. W. Sun, Y. L. Zhao, X. Zheng, An implicit difference scheme for time-fractional diffusion equations with a time-invariant type variable order, Appl. Math. Lett., 120 (2021), 107270. doi: 10.1016/j.aml.2021.107270. doi: 10.1016/j.aml.2021.107270
![]() |
[17] |
Y. Xu, Y. Zhang, J. Zhao, Backward difference formulae and spectral galerkin methods for the riesz space fractional diffusion equation, Math. Comput. Simulat., 166 (2019), 494–507. doi: 10.1016/j.matcom.2019.07.007. doi: 10.1016/j.matcom.2019.07.007
![]() |
[18] |
X. Gao, B. Yin, H. Li, Y. Liu, Tt-m FE method for a 2D nonlinear time distributed-order and space fractional diffusion equation, Math. Comput. Simulat., 181 (2021), 117–137. doi: 10.1016/j.matcom.2020.09.021. doi: 10.1016/j.matcom.2020.09.021
![]() |
[19] |
X. Zheng, H. Wang, An error estimate of a numerical approximation to a hidden-memory variable-order space-time fractional diffusion equation, SIAM J. Numer. Anal., 58 (2020), 2492–2514. doi: 10.1137/20M132420X. doi: 10.1137/20M132420X
![]() |
[20] |
X. Zheng, H. Wang, A Hidden-memory variable-order time-fractional optimal control model: Analysis and approximation, SIAM J. Control Optim., 59 (2021), 1851–1880. doi: 10.1137/20M1344962. doi: 10.1137/20M1344962
![]() |
[21] |
X. Zheng, H. Wang, Optimal-order error estimates of finite element approximations to variable-order time-fractional diffusion equations without regularity assumptions of the true solutions, IMA J. Numer. Anal., 41 (2021), 1522–1545. doi: 10.1093/imanum/draa013. doi: 10.1093/imanum/draa013
![]() |
[22] |
H. Wang, X. Zheng, Wellposedness and regularity of the variable-order time-fractional diffusion equations, J. Math. Anal. Appl., 475 (2019), 1778–1802. doi: 10.1016/j.jmaa.2019.03.052. doi: 10.1016/j.jmaa.2019.03.052
![]() |
[23] |
Y. Luchko, Initial-boundary-value problems for the one-dimensional time-fractional diffusion equation, Fract. Calc. Appl. Anal., 15 (2012), 141–160. doi: 10.2478/s13540-012-0010-7. doi: 10.2478/s13540-012-0010-7
![]() |
[24] |
B. Jin, B. Li, Z. Zhou, Numerical analysis of nonlinear subdiffusion equations, SIAM J. Numer. Anal., 56 (2018), 1–23. doi: 10.1137/16M1089320. doi: 10.1137/16M1089320
![]() |
[25] |
X. Zheng, H. Wang, Wellposedness and regularity of a variable-order space-time fractional diffusion equation, Anal. Appl., 18 (2020), 615–638. doi: 10.1142/S0219530520500013. doi: 10.1142/S0219530520500013
![]() |
[26] |
H. Fu, M. K. Ng, H. Wang, A divide-and-conquer fast finite difference method for space-time fractional partial differential equation, Comput. Math. Appl., 73 (2017), 1233–1242. doi: 10.1016/j.camwa.2016.11.023. doi: 10.1016/j.camwa.2016.11.023
![]() |
[27] | F. M. Salama, N. H. M. Ali, Fast O(N) hybrid method for the solution of two dimensional time fractional cable equation, Compusoft, 8 (2019), 3453–3461. |
[28] |
F. M. Salama, N. H. M. Ali, Computationally efficient hybrid method for the numerical solution of the 2D time fractional advection-diffusion equation, Int. J. Math. Eng. Manag., 5 (2020), 432–446. doi: 10.33889/IJMEMS.2020.5.3.036. doi: 10.33889/IJMEMS.2020.5.3.036
![]() |
[29] |
X. M. Gu, S. L. Wu, A parallel-in-time iterative algorithm for Volterra partial integro-differential problems with weakly singular kernel, J. Comput. Phys., 417 (2020), 109576. doi: 10.1016/j.jcp.2020.109576. doi: 10.1016/j.jcp.2020.109576
![]() |
[30] |
X. L. Lin, M. K. Ng, H. W. Sun, A multigrid method for linear systems arising from time-dependent two-dimensional space-fractional diffusion equations, J. Comput. Phys., 336 (2017), 69–86. doi: 10.1016/j.jcp.2017.02.008. doi: 10.1016/j.jcp.2017.02.008
![]() |
[31] |
J. Ren, Z. Z. Sun, W. Dai, New approximations for solving the caputo-type fractional partial differential equations, Appl. Math. Model., 40 (2016), 2625–2636. doi: 10.1016/j.apm.2015.10.011. doi: 10.1016/j.apm.2015.10.011
![]() |
[32] |
N. A. Khan, S. Ahmed, Finite difference method with metaheuristic orientation for exploration of time fractional partial differential equations, Int. J. Appl. Comput. Math., 7 (2021), 1–22. doi: 10.1007/s40819-021-01061-y. doi: 10.1007/s40819-021-01061-y
![]() |
[33] |
A. Ahmadian, S. Salahshour, M. Ali-Akbari, F. Ismail, D. Baleanu, A novel approach to approximate fractional derivative with uncertain conditions, Chaos Soliton. Fract., 104 (2017), 68–76. doi: 10.1016/j.chaos.2017.07.026. doi: 10.1016/j.chaos.2017.07.026
![]() |
[34] |
F. M. Salama, N. H. M. Ali, N. N. Abd Hamid, Fast O(N) hybrid Laplace transform-finite difference method in solving 2D time fractional diffusion equation, J. Math. Comput. Sci., 23 (2021), 110–123. doi: 10.22436/jmcs.023.02.04. doi: 10.22436/jmcs.023.02.04
![]() |
[35] |
N. H. M. Ali, L. M. Kew, New explicit group iterative methods in the solution of two dimensional hyperbolic equations, J. Comput. Phys., 231 (2012), 6953–6968. doi: 10.1016/j.jcp.2012.06.025. doi: 10.1016/j.jcp.2012.06.025
![]() |
[36] |
N. H. Mohd Ali, A. Mohammed Saeed, Convergence analysis of the preconditioned group splitting methods in boundary value problems, Abstr. Appl. Anal., 2012 (2012), 867598. doi: 10.1155/2012/867598. doi: 10.1155/2012/867598
![]() |
[37] |
L. M. Kew, N. H. M. Ali, New explicit group iterative methods in the solution of three dimensional hyperbolic telegraph equations, J. Comput. Phys., 294 (2015), 382–404. doi: 10.1016/j.jcp.2015.03.052. doi: 10.1016/j.jcp.2015.03.052
![]() |
[38] |
A. Saudi, J. Sulaiman, Robot path planning using four point-explicit group via nine-point laplacian (4EG9L) iterative method, Procedia Engineering, 41 (2012), 182–188. doi: 10.1016/j.proeng.2012.07.160. doi: 10.1016/j.proeng.2012.07.160
![]() |
[39] |
N. H. M. Ali, A. M. Saeed, Preconditioned modified explicit decoupled group for the solution of steady state navier-stokes equation, Appl. Math. Inform. Sci., 7 (2013), 1837–1844. doi: 10.12785/amis/070522. doi: 10.12785/amis/070522
![]() |
[40] |
M. A. Khan, N. H. M. Ali, N. N. Abd Hamid, A new fourth-order explicit 270 group method in the solution of two-dimensional fractional rayleigh–stokes problem for a heated generalized second-grade fluid, Adv. Diff. Equ., 2020 (2020), 598. doi: 10.1186/s13662-020-03061-6. doi: 10.1186/s13662-020-03061-6
![]() |
[41] |
N. Abdi, H. Aminikhah, A. H. Sheikhani, J. Alavi, M. Taghipour, An efficient explicit decoupled group method for solving two–dimensional fractional Burgers' equation and its convergence analysis, Adv. Math. Phys., 2021 (2021), 6669287. doi: 10.1155/2021/6669287. doi: 10.1155/2021/6669287
![]() |
[42] |
N. Abdi, H. Aminikhah, A. R. Sheikhani, High-order rotated grid point iterative method for solving 2D time fractional telegraph equation and its convergence analysis, Comp. Appl. Math., 40 (2021), 54. doi: 10.1007/s40314-021-01451-4. doi: 10.1007/s40314-021-01451-4
![]() |
[43] |
F. M. Salama, N. H. M. Ali, N. N. Abd Hamid, Efficient hybrid group iterative methods in the solution of two-dimensional time fractional cable equation, Adv. Differ. Equ., 2020 (2020), 257. doi: 10.1186/s13662-020-02717-7. doi: 10.1186/s13662-020-02717-7
![]() |
[44] |
N. Moraca, Bounds for norms of the matrix inverse and the smallest singular value, Linear Algebra Appl., 429 (2008), 2589–2601. doi: 10.1016/j.laa.2007.12.026. doi: 10.1016/j.laa.2007.12.026
![]() |
[45] |
A. Ali, N. H. M. Ali, On skewed grid point iterative method for solving 2d hyperbolic telegraph fractional differential equation, Adv. Differ. Equ., 2019 (2019), 303. doi: 10.1186/s13662-019-2238-6. doi: 10.1186/s13662-019-2238-6
![]() |
1. | Uday Singh, Numerical investigation of unconditionally stable spline function for three-dimensional time-fractional telegraph equations, 2022, 9, 26667207, 100180, 10.1016/j.rico.2022.100180 | |
2. | Abdul Hamid Ganie, Abdulkafi Mohammed Saeed, Sadia Saeed, Umair Ali, Muhammad Irfan, The Rayleigh–Stokes Problem for a Heated Generalized Second-Grade Fluid with Fractional Derivative: An Implicit Scheme via Riemann–Liouville Integral, 2022, 2022, 1563-5147, 1, 10.1155/2022/6948461 | |
3. | UMAIR ALI, MUHAMMAD NAEEM, FARAH AINI ABDULLAH, MIAO-KUN WANG, FOUAD MOHAMMAD SALAMA, ANALYSIS AND IMPLEMENTATION OF NUMERICAL SCHEME FOR THE VARIABLE-ORDER FRACTIONAL MODIFIED SUB-DIFFUSION EQUATION, 2022, 30, 0218-348X, 10.1142/S0218348X22402538 | |
4. | Nurathirah Sulaiman, Jumat Sulaiman, Mohammad Khatim Hasan, Samsul Ariffin Abdul Karim, 2022, Chapter 22, 978-3-031-04027-6, 341, 10.1007/978-3-031-04028-3_22 | |
5. | Fouad Mohammad Salama, Nur Nadiah Abd Hamid, Umair Ali, Norhashidah Hj. Mohd Ali, Fast hybrid explicit group methods for solving 2D fractional advection-diffusion equation, 2022, 7, 2473-6988, 15854, 10.3934/math.2022868 | |
6. | Fouad Mohammad Salama, Umair Ali, Ajmal Ali, Numerical Solution of Two-Dimensional Time Fractional Mobile/Immobile Equation Using Explicit Group Methods, 2022, 8, 2349-5103, 10.1007/s40819-022-01408-z | |
7. | Anam Naz, Umair Ali, Ashraf Elfasakhany, Khadiga Ahmed Ismail, Abdullah G. Al-Sehemi, Ahmed A. Al-Ghamdi, An Implicit Numerical Approach for 2D Rayleigh Stokes Problem for a Heated Generalized Second Grade Fluid with Fractional Derivative, 2021, 5, 2504-3110, 283, 10.3390/fractalfract5040283 | |
8. | Wajaree Weera, Chantapish Zamart, Zulqurnain Sabir, Muhammad Asif Zahoor Raja, Afaf S. Alwabli, S. R. Mahmoud, Supreecha Wongaree, Thongchai Botmart, Fractional Order Environmental and Economic Model Investigations Using Artificial Neural Network, 2023, 74, 1546-2226, 1735, 10.32604/cmc.2023.032950 | |
9. | Renu Choudhary, Satpal Singh, Devendra Kumar, A higher order unconditionally stable numerical technique for multi-term time-fractional diffusion and advection–diffusion equations, 2024, 43, 2238-3603, 10.1007/s40314-024-02688-5 | |
10. | Fouad Mohammad Salama, Faisal Fairag, On numerical solution of two-dimensional variable-order fractional diffusion equation arising in transport phenomena, 2024, 9, 2473-6988, 340, 10.3934/math.2024020 | |
11. | Shi-Ping Tang, Yu-Mei Huang, A fast preconditioning iterative method for solving the discretized second-order space-fractional advection–diffusion equations, 2024, 438, 03770427, 115513, 10.1016/j.cam.2023.115513 | |
12. | Mahdi Ahmadinia, Mokhtar Abbasi, Parisa Hadi, An intelligent non-uniform mesh to improve errors of a stable numerical method for time-tempered fractional advection–diffusion equation with weakly singular solution, 2024, 80, 0920-8542, 26280, 10.1007/s11227-024-06442-w | |
13. | Jianbing Hu, Analyzing the Approximate Error and Applicable Condition of the Fractional Reduced Differential Transform Method, 2024, 16, 2073-8994, 912, 10.3390/sym16070912 | |
14. | Fouad Mohammad Salama, Alla Tareq Balasim, Umair Ali, Muhammad Asim Khan, Efficient numerical simulations based on an explicit group approach for the time fractional advection–diffusion reaction equation, 2023, 42, 2238-3603, 10.1007/s40314-023-02278-x | |
15. | Bi-Yun Zhu, Ai-Guo Xiao, Xue-Yang Li, An efficient numerical method on modified space-time sparse grid for time-fractional diffusion equation with nonsmooth data, 2023, 94, 1017-1398, 1561, 10.1007/s11075-023-01547-4 | |
16. | Fouad Mohammad Salama, On Numerical Simulations of Variable-Order Fractional Cable Equation Arising in Neuronal Dynamics, 2024, 8, 2504-3110, 282, 10.3390/fractalfract8050282 | |
17. | Mohamed Abdelsabour Fahmy, Roqia Abdullah A. Jeli, A New Fractional Boundary Element Model for Anomalous Thermal Stress Effects on Cement-Based Materials, 2024, 8, 2504-3110, 753, 10.3390/fractalfract8120753 |
Algorithm 1: MHEG method utilizing the Gauss-Seidel iterative solver |
1: Branch the mesh points of the discretized solution domain into the following categories:
● Iterative points of type ⧫. ● Direct points of types ◻ and ○. 2: Arrange all the ⧫ points into groups of four mesh points as depicted in Figure 2. 3: Set an initial guess for the solution at the current time level. 4: For each four-point group at the current time level, iterate the intermediate solutions at ⧫ mesh points utilizing (ˆuk+1,n+1i,jˆuk+1,n+1i+2,jˆuk+1,n+1i+2,j+2ˆuk+1,n+1i,j+2)=1η(η1η2η3η4η2η1η4η3η3η4η1η2η4η3η2η1)(rhsi,jrhsi+2,jrhsi+2,j+2rhsi,j+2), where rhsi,j=q14uk+1i−2,j+q24uk+1i,j−2+uki,j+(r−1)τu0i,j+rτfk+1i,j,rhsi+2,j=q14uk+1i+4,j+q24uk+1i+2,j−2+uki+2,j+(r−1)τu0i+2,j+rτfk+1i+2,j,rhsi+2,j+2=q14uk+1i+4,j+2+q24uk+1i+2,j+4+uki+2,j+2+(r−1)τu0i+2,j+2+rτfk+1i+2,j+2,rhsi,j+2=q14uk+1i−2,j+2+q24uk+1i,j+4+uki,j+2+(r−1)τu0i,j+2+rτfk+1i,j+2, and perform the Gauss-Seidel solver (uk+1,n+1i,juk+1,n+1i+2,juk+1,n+1i+2,j+2uk+1,n+1i,j+2)=ω(ˆuk+1,n+1i,jˆuk+1,n+1i+2,jˆuk+1,n+1i+2,j+2ˆuk+1,n+1i,j+2)+(1−ω)(uk+1,ni,juk+1,ni+2,juk+1,ni+2,j+2uk+1,ni,j+2), where n is the iteration number and ω=1 is the relaxation parameter of the Gauss-Seidel solver. η, η1, η2, η3 and η4 are as defined before in this section. 5: Test the convergence. If convergence is achieved, go to step 6. Otherwise, step 4 is repeated until convergence is attained. 6: The solutions on the rest of mesh points of types ◻ and ○ are evaluated directly once in the following manner: ● For ◻ points, the approximating PDE (2.2) is discretized by using skewed finite difference operators. Such difference operators are derived by rotating the standard mesh an angle 45∘ clockwise and applying the Taylor series expansion afterwards [45]. This will lead to the following skewed difference scheme for (2.2): uk+1i,j=1(1+(r−1)τ+q1+q2)[q12(uk+1i+1,j−1+uk+1i−1,j+1)+q22(uk+1i+1,j+1+uk+1i−1,j−1)+uki,j+(r−1)τu0i,j+rτfk+1i,j]. ● For ○ points, the difference scheme formula defined by (2.6) is employed. 7: If the final time level is reached, go to step 8. Otherwise, adopt the solution values of the last time level as an initial guess for the next time level and go to step 4. 8: Display the numerical solutions. |
Method | Per iteration | After convergence | Total operations |
HSP | 13σ2∗Ite | - | 13σ2∗Ite |
MHEG | 4.5(σ−1)2∗Ite | 3.25(3σ2+2σ−1) | 4.5(σ−1)2∗Ite+3.25(3σ2+2σ−1) |
α | τ=h | Method | CPU time (sec) | Ite | Total operations | Emax |
0.55 | ||||||
1/6 | HSP | 0.062 | 27 | 8,775 | 1.6737E-03 | |
MHEG | 0.015 | 2 | 417 | 1.7983E-03 | ||
1/10 | HSP | 0.968 | 56 | 58,968 | 9.5405E-04 | |
MHEG | 0.140 | 11 | 4,013 | 1.0764E-03 | ||
1/14 | HSP | 3.937 | 86 | 188,942 | 6.0571E-04 | |
MHEG | 0.359 | 16 | 12,097 | 7.6144E-04 | ||
1/18 | HSP | 10.062 | 117 | 439,569 | 3.8681E-04 | |
MHEG | 0.734 | 22 | 28,269 | 5.7584E-04 | ||
0.75 | ||||||
1/6 | HSP | 0.062 | 26 | 8,450 | 2.3916E-03 | |
MHEG | 0.015 | 2 | 417 | 2.4750E-03 | ||
1/10 | HSP | 0.937 | 53 | 55,809 | 1.4477E-03 | |
MHEG | 0.109 | 10 | 3,725 | 1.5593E-03 | ||
1/14 | HSP | 3.109 | 80 | 175,760 | 9.9823E-04 | |
MHEG | 0.250 | 15 | 11,449 | 1.1488E-03 | ||
1/18 | HSP | 8.062 | 107 | 401,999 | 6.9546E-04 | |
MHEG | 0.640 | 21 | 27,117 | 9.0926E-04 | ||
0.95 | ||||||
1/6 | HSP | 0.062 | 25 | 8,125 | 2.9204E-03 | |
MHEG | 0.015 | 2 | 417 | 2.9789E-03 | ||
1/10 | HSP | 0.656 | 50 | 52,650 | 1.7502E-03 | |
MHEG | 0.062 | 10 | 3,725 | 1.8611E-03 | ||
1/14 | HSP | 2.562 | 75 | 164,775 | 1.1935E-03 | |
MHEG | 0.203 | 14 | 10,801 | 1.3410E-03 | ||
1/18 | HSP | 7.937 | 100 | 375,700 | 8.2842E-04 | |
MHEG | 0.515 | 19 | 24,813 | 1.0461E-03 |
α | τ=h | Method | CPU time (sec) | Ite | Total operations | Emax |
0.55 | ||||||
1/18 | HSP | 12.671 | 170 | 638,690 | 6.3352E-03 | |
MHEG | 0.859 | 29 | 36,333 | 6.8017E-03 | ||
1/22 | HSP | 31.906 | 221 | 1,266,993 | 5.0040E-03 | |
MHEG | 2.359 | 38 | 72,833 | 5.5331E-03 | ||
1/26 | HSP | 64.234 | 271 | 2,201,875 | 4.0351E-03 | |
MHEG | 3.390 | 46 | 125,485 | 4.6340E-03 | ||
1/30 | HSP | 118.625 | 320 | 3,498,560 | 3.2629E-03 | |
MHEG | 7.937 | 55 | 202,425 | 3.8992E-03 | ||
0.75 | ||||||
1/18 | HSP | 12.640 | 156 | 586,092 | 1.0239E-02 | |
MHEG | 0.625 | 27 | 34,029 | 1.0665E-02 | ||
1/22 | HSP | 28.828 | 199 | 1,140,867 | 8.4919E-03 | |
MHEG | 1.343 | 34 | 65,633 | 9.0001E-03 | ||
1/26 | HSP | 57.593 | 242 | 1,966,250 | 7.1948E-03 | |
MHEG | 3.359 | 42 | 115,117 | 7.8080E-03 | ||
1/30 | HSP | 102.688 | 284 | 3,104,972 | 6.1807E-03 | |
MHEG | 7.328 | 49 | 181,257 | 6.9414E-03 | ||
0.95 | ||||||
1/18 | HSP | 12.359 | 143 | 537,251 | 1.1950E-02 | |
MHEG | 0.562 | 25 | 31,725 | 1.2403E-02 | ||
1/22 | HSP | 25.968 | 182 | 1,043,406 | 9.7841E-03 | |
MHEG | 1.328 | 31 | 60,233 | 1.0285E-02 | ||
1/26 | HSP | 51.796 | 219 | 1,779,375 | 8.1651E-03 | |
MHEG | 3.203 | 38 | 104,749 | 8.8048E-03 | ||
1/30 | HSP | 94.796 | 255 | 2,787,915 | 6.9147E-03 | |
MHEG | 6.890 | 44 | 163,617 | 7.6882E-03 |
Test problem | τ | Emax | Computational order |
1 | 1/10 | 9.9154E-04 | 0.9835 |
1/20 | 5.0149E-04 | ||
2 | 1/10 | 1.1829E-02 | 0.9658 |
1/20 | 6.0570E-03 |
Algorithm 1: MHEG method utilizing the Gauss-Seidel iterative solver |
1: Branch the mesh points of the discretized solution domain into the following categories:
● Iterative points of type ⧫. ● Direct points of types ◻ and ○. 2: Arrange all the ⧫ points into groups of four mesh points as depicted in Figure 2. 3: Set an initial guess for the solution at the current time level. 4: For each four-point group at the current time level, iterate the intermediate solutions at ⧫ mesh points utilizing (ˆuk+1,n+1i,jˆuk+1,n+1i+2,jˆuk+1,n+1i+2,j+2ˆuk+1,n+1i,j+2)=1η(η1η2η3η4η2η1η4η3η3η4η1η2η4η3η2η1)(rhsi,jrhsi+2,jrhsi+2,j+2rhsi,j+2), where rhsi,j=q14uk+1i−2,j+q24uk+1i,j−2+uki,j+(r−1)τu0i,j+rτfk+1i,j,rhsi+2,j=q14uk+1i+4,j+q24uk+1i+2,j−2+uki+2,j+(r−1)τu0i+2,j+rτfk+1i+2,j,rhsi+2,j+2=q14uk+1i+4,j+2+q24uk+1i+2,j+4+uki+2,j+2+(r−1)τu0i+2,j+2+rτfk+1i+2,j+2,rhsi,j+2=q14uk+1i−2,j+2+q24uk+1i,j+4+uki,j+2+(r−1)τu0i,j+2+rτfk+1i,j+2, and perform the Gauss-Seidel solver (uk+1,n+1i,juk+1,n+1i+2,juk+1,n+1i+2,j+2uk+1,n+1i,j+2)=ω(ˆuk+1,n+1i,jˆuk+1,n+1i+2,jˆuk+1,n+1i+2,j+2ˆuk+1,n+1i,j+2)+(1−ω)(uk+1,ni,juk+1,ni+2,juk+1,ni+2,j+2uk+1,ni,j+2), where n is the iteration number and ω=1 is the relaxation parameter of the Gauss-Seidel solver. η, η1, η2, η3 and η4 are as defined before in this section. 5: Test the convergence. If convergence is achieved, go to step 6. Otherwise, step 4 is repeated until convergence is attained. 6: The solutions on the rest of mesh points of types ◻ and ○ are evaluated directly once in the following manner: ● For ◻ points, the approximating PDE (2.2) is discretized by using skewed finite difference operators. Such difference operators are derived by rotating the standard mesh an angle 45∘ clockwise and applying the Taylor series expansion afterwards [45]. This will lead to the following skewed difference scheme for (2.2): uk+1i,j=1(1+(r−1)τ+q1+q2)[q12(uk+1i+1,j−1+uk+1i−1,j+1)+q22(uk+1i+1,j+1+uk+1i−1,j−1)+uki,j+(r−1)τu0i,j+rτfk+1i,j]. ● For ○ points, the difference scheme formula defined by (2.6) is employed. 7: If the final time level is reached, go to step 8. Otherwise, adopt the solution values of the last time level as an initial guess for the next time level and go to step 4. 8: Display the numerical solutions. |
Method | Per iteration | After convergence | Total operations |
HSP | 13σ2∗Ite | - | 13σ2∗Ite |
MHEG | 4.5(σ−1)2∗Ite | 3.25(3σ2+2σ−1) | 4.5(σ−1)2∗Ite+3.25(3σ2+2σ−1) |
α | τ=h | Method | CPU time (sec) | Ite | Total operations | Emax |
0.55 | ||||||
1/6 | HSP | 0.062 | 27 | 8,775 | 1.6737E-03 | |
MHEG | 0.015 | 2 | 417 | 1.7983E-03 | ||
1/10 | HSP | 0.968 | 56 | 58,968 | 9.5405E-04 | |
MHEG | 0.140 | 11 | 4,013 | 1.0764E-03 | ||
1/14 | HSP | 3.937 | 86 | 188,942 | 6.0571E-04 | |
MHEG | 0.359 | 16 | 12,097 | 7.6144E-04 | ||
1/18 | HSP | 10.062 | 117 | 439,569 | 3.8681E-04 | |
MHEG | 0.734 | 22 | 28,269 | 5.7584E-04 | ||
0.75 | ||||||
1/6 | HSP | 0.062 | 26 | 8,450 | 2.3916E-03 | |
MHEG | 0.015 | 2 | 417 | 2.4750E-03 | ||
1/10 | HSP | 0.937 | 53 | 55,809 | 1.4477E-03 | |
MHEG | 0.109 | 10 | 3,725 | 1.5593E-03 | ||
1/14 | HSP | 3.109 | 80 | 175,760 | 9.9823E-04 | |
MHEG | 0.250 | 15 | 11,449 | 1.1488E-03 | ||
1/18 | HSP | 8.062 | 107 | 401,999 | 6.9546E-04 | |
MHEG | 0.640 | 21 | 27,117 | 9.0926E-04 | ||
0.95 | ||||||
1/6 | HSP | 0.062 | 25 | 8,125 | 2.9204E-03 | |
MHEG | 0.015 | 2 | 417 | 2.9789E-03 | ||
1/10 | HSP | 0.656 | 50 | 52,650 | 1.7502E-03 | |
MHEG | 0.062 | 10 | 3,725 | 1.8611E-03 | ||
1/14 | HSP | 2.562 | 75 | 164,775 | 1.1935E-03 | |
MHEG | 0.203 | 14 | 10,801 | 1.3410E-03 | ||
1/18 | HSP | 7.937 | 100 | 375,700 | 8.2842E-04 | |
MHEG | 0.515 | 19 | 24,813 | 1.0461E-03 |
α | τ=h | Method | CPU time (sec) | Ite | Total operations | Emax |
0.55 | ||||||
1/18 | HSP | 12.671 | 170 | 638,690 | 6.3352E-03 | |
MHEG | 0.859 | 29 | 36,333 | 6.8017E-03 | ||
1/22 | HSP | 31.906 | 221 | 1,266,993 | 5.0040E-03 | |
MHEG | 2.359 | 38 | 72,833 | 5.5331E-03 | ||
1/26 | HSP | 64.234 | 271 | 2,201,875 | 4.0351E-03 | |
MHEG | 3.390 | 46 | 125,485 | 4.6340E-03 | ||
1/30 | HSP | 118.625 | 320 | 3,498,560 | 3.2629E-03 | |
MHEG | 7.937 | 55 | 202,425 | 3.8992E-03 | ||
0.75 | ||||||
1/18 | HSP | 12.640 | 156 | 586,092 | 1.0239E-02 | |
MHEG | 0.625 | 27 | 34,029 | 1.0665E-02 | ||
1/22 | HSP | 28.828 | 199 | 1,140,867 | 8.4919E-03 | |
MHEG | 1.343 | 34 | 65,633 | 9.0001E-03 | ||
1/26 | HSP | 57.593 | 242 | 1,966,250 | 7.1948E-03 | |
MHEG | 3.359 | 42 | 115,117 | 7.8080E-03 | ||
1/30 | HSP | 102.688 | 284 | 3,104,972 | 6.1807E-03 | |
MHEG | 7.328 | 49 | 181,257 | 6.9414E-03 | ||
0.95 | ||||||
1/18 | HSP | 12.359 | 143 | 537,251 | 1.1950E-02 | |
MHEG | 0.562 | 25 | 31,725 | 1.2403E-02 | ||
1/22 | HSP | 25.968 | 182 | 1,043,406 | 9.7841E-03 | |
MHEG | 1.328 | 31 | 60,233 | 1.0285E-02 | ||
1/26 | HSP | 51.796 | 219 | 1,779,375 | 8.1651E-03 | |
MHEG | 3.203 | 38 | 104,749 | 8.8048E-03 | ||
1/30 | HSP | 94.796 | 255 | 2,787,915 | 6.9147E-03 | |
MHEG | 6.890 | 44 | 163,617 | 7.6882E-03 |
Test problem | τ | Emax | Computational order |
1 | 1/10 | 9.9154E-04 | 0.9835 |
1/20 | 5.0149E-04 | ||
2 | 1/10 | 1.1829E-02 | 0.9658 |
1/20 | 6.0570E-03 |