A Meshfree Approach for Solving Fractional Galilei Invariant Advection–Diffusion Equation through Weighted–Shifted Grünwald Operator

: Fractional Galilei invariant advection–diffusion (GIADE) equation, along with its more general version that is the GIADE equation with nonlinear source term, is discretized by coupling weighted and shifted Grünwald difference approximation formulae and Crank–Nicolson technique. The new version of the backward substitution method, a well-established class of meshfree methods, is proposed for a numerical approximation of the consequent equation. In the present approach, the ﬁnal approximation is given by the summation of the radial basis functions, the primary approximation, and the related correcting functions. Then, the approximation is substituted back to the governing equations where the unknown parameters can be determined. The polynomials, trigonometric functions, multiquadric, or the Gaussian radial basis functions are used in the approximation of the GIADE. Moreover, a quasilinearization technique is employed to transform a nonlinear source term into a linear source term. Finally, three numerical experiments in one and two dimensions are presented to support the method.


Introduction
The problem of diffusion has drawn significant interest and wide recognition in the scientific community in the past decade [1][2][3][4][5].Fractional calculus has played an important role in obtaining advection-diffusion equations.The method of separation of variables, as well as the theory of eigenvalues and eigenfunctions when constructing a solution to an advection-diffusion equation, is studied in [6].The analytical and numerical solution of a space-time fractional advection-diffusion equation with first-order accuracy in the temporal direction and second-order accuracy in the spatial direction is discussed in [7].Moreover, the enhanced unconditionally positive finite difference method investigates the solutions to the linear and nonlinear advection-diffusion-reaction (ADR) to reduce the degree of freedom of the ADR equations [8].A possible realization of the fractional Galilei invariant advection-diffusion (GIADE) equation is the transport of material in an external velocity or force field [9,10].More precisely, we study the following fractional GIADE equation with the nonlinear source term      ∂u(x,t) ∂t + ν∇u(x, t) − ∇ 2 u(x, t) = 0 D 1−γ t ∇ 2 u(x, t) + f (x, t) + G(x, t, u(x, t)), x ∈ Ω, u(x, 0) = ϕ(x), x ∈ Ω, u(x, t) = ψ(x), x ∈ ∂Ω, (1) where 0 < t ≤ T, ∇ = ∂ ∂x + ∂ ∂y , the Laplacian operator is defined as ∇ 2 = ∂ 2 ∂x 2 + ∂ 2 ∂y 2 , ϕ(x), ψ(x) and f (x, t) are specific functions, and ν is a real constant with the nonlinear source term G(x, t, u(x, t)).Ω is a domain in the plane R 2 or R with the boundary ∂Ω.
The symbol 0 D 1−γ t is the Riemann-Liouville fractional derivative operator of order (0 < γ < 1) with respect to t with the starting point at t = 0 and is defined as ( where Γ(.)is the gamma function.Chen et al. [11] have presented an implicit difference approximation with second-order spatial accuracy and first-order temporal accuracy for a variable-order GIADE equation in one dimension.Zaky et al. [12] have studied the existence and uniqueness solutions to nonlinear tempered fractional boundary value problems.Moreover, they have proved that the convergence rate for nonsmooth solutions can be enhanced by using a suitable smoothing transformation [13].In [9], the fractional GIADE equation has been investigated by interpolating a stabilized element-free Galerkin method with the order of convergence O(τ 2 + r m+1 ).Moreover, in the numerical results, the effect of the advection coefficient and the smoothness of the numerical solution has been studied.There is also a wealth of literature regarding numerical solutions to the GIADE equation [14,15].Therefore, it is more challenging to consider numerical schemes and related numerical analyses for the problem in a more general form.This article aims to describe, for the first time, a backward substitution method (BSM) for numerical approximation of a fractional problem with a coupled weighted and shifted Grünwald difference (WSGD) approximation formula and the Crank-Nicolson (C-N) technique.First, the standard WSGD formula has been discussed in detail and successfully applied to solve fractional equations [16].In [17], stability and consistency results for higher-order Grünwald-type formulae have been used to space fractional partial differential equations (PDEs).Finite difference schemes and the quasicompact difference schemes by shifted WSGD operators for discretizing the Riemann-Liouville fractional derivative have been established to solve the fractional diffusion equations numerically [18,19].However, the 3-WSGD operator that was applied in [19] failed to solve the time-dependent fractional equations with unconditional stability.Wang [20] took advantage of the WSGD operator and analyzed a finite difference scheme for solving the nonlinear complex fractional Ginzburg-Landau equation.Meerschaert et al. [21,22] explored consistency and stability for numerical schemes for fractional PDEs using a shifted Grünwald formula to approximate the fractional derivative.Hence, we adopt the shifted formula to discretize the time-fractional derivative.In real applications, analytical methods cannot work well on most PDEs, but there have been many techniques to solve PDEs covered by Laplace transform method [23], the Green function method [24,25], factorization method [26,27] and so on [28,29].This paper makes an important contribution to the literature by expressing the solution by BSM at each time step utilizing the WSGD operator.Reutskiy [30] applied the BSM based on the Fourier series expansion for solving initial and boundary value problems for fractional PDEs in one dimension.In [31], the combination of the C-N method with the BSM for solving space-fractional advection-diffusion reaction equations in one dimension has been considered.Safari and Azarsa [32] took advantage of the Müntz polynomials and numerically solved nonlinear space-fractional PDEs with meshless BSM.
Because of the nonlocal operator of the fractional derivative, it is challenging to harness using the BSM as a meshfree method for the fractional problem.So, the algorithm of the method contains two stages.First, the WSGD operator to discretize the involved fractional derivative is used.Then, the expression of the approximation is obtained by the sum of two parts, which includes the free parameters and already satisfies all the boundary conditions of the original problem but does not satisfy the equation.The parameters are determined in the second stage by the backward substitution into the original equation.The first part of the approximate solution, which is called primary approximation, is obtained from the boundary conditions, and the second part are the corrected basis functions, which consist of trigonometric functions and radial basis functions (RBFs).
It should be noted here that, for the first time, the method for the numerical approximation of fractional problems in two dimensions and three dimensions in the near future is applied.The approach consists of the advantages of RBFs in reducing the condition number of the matrix [33,34].However, unlike the singular boundary method, the approach can be applied to problems in which equations contain a nonlinear term with variable coefficients [1,35].Furthermore, the method of fundamental solutions is applied to solve the mentioned problem [36,37].However, the truth of the real-world problems is that the determination of the optimal sources in the method has remained an open issue.
To solve the fractional GIADE equation, first, a mesh and discretization of the general equation in time is introduced, then the C-N scheme and the WSGD operator to discretize the involved fractional derivative (Section 2) are coupled.Next, in Section 3, the numerical scheme's BSM to solve the consequent equation is proposed.Finally, in Section 4, the results of some numerical experiments, including a two-dimensional GIADE equation with nonlinear source term highlighting the efficiency of the BSM is given.

The Time Discretization Approximation
In this section, we introduce the concepts that are applied for discretization of the time variable.For a fixed positive number τ = t n+1 − t n , the grid points in the time interval [0, T] are labeled t n = nτ.It is well known that the approximation of the first-order derivative is as follows: Second-order approximations of the Riemann-Liouville fractional derivative (the WSGD operators), which is important for the numerical approximation of fractional differential equations, are given as follows [18,38,39] 0 where ω γ j are the normalized Grünwald weights and the coefficients of the series of the function (1 − z) γ and can be evaluated as and Denoting the approximations of u(x, t n ), G(x, t n , u(x, t n )), and f (x, t n ) by u n , G n , and f n , respectively, omitting the truncation error and using the C-N technique for time discretization of Equation ( 1), we obtain the following equation for n ≥ 0 Using the Grünwald formula in Equation ( 4) and rearranging Equation ( 7), we obtain the following equations: Now, we apply the quasilinearization technique [40] which makes G n+1 to a linear function of u n+1 .Suppose that u n+1 0 are the given functions of x, which are the initial approximations of the corresponding exact values at the n + 1 steps.Then, we have the relations By substitution Equation ( 9) in Equation ( 8), we have We denote the left-hand side of the Equation ( 8) as L[u n+1 ], then the above equation can be written as with the boundary equation where Notice that for n = 0, from Equation ( 8) and initial condition u 0 = ϕ, we have We are now ready to introduce the BSM method for solving Equation (11) in the next section.We extend the BSM approach for governing Equation (11) in each step time.This extension is based on RBFs and a particular solution method.We later discuss the method for one-dimensional and two-dimensional cases.

Description of the BSM Method
Now, the meshless BSM method is described as a numerical tool for the solution of Equation (11).The approach for solving the inhomogeneous problem is to express the solution u n+1 by the sum of two parts of a particular solution u n+1 p and w n+1 as where u n+1 p (x) approximates the boundary condition Equation ( 11) and w n+1 (x) is the correcting solution of the problem with the homogeneous boundary condition as follows: To obtain the approximate solution, we consider the series where φ m (x) are the multiquadric or the Gaussian RBFs on Ω, ω m (x) is a defined correcting function that can be approximated by using basis functions, and {q m } M m=1 are unknown coefficients.Moreover, ω m (x) is taken in such a way that φ m (x) satisfy the boundary condition In order to approximate u n+1 p (x) and ω m (x), we chose a system of trigonometric functions in a one-dimensional problem ρ (x) = sin ρπ x+ζ 2ζ and a two-dimensional problem where {p m,ρ } K ρ=1 and {p M+1,ρ } K ρ=1 are the unknown coefficients which can be determined by enforcing the Equation ( 19) on the boundary nodes as where the number of boundary collocation nodes K is larger than the number of basis function K 1 .

L[w
Substituting a set of collocation points into Equation ( 21) , {q m } M m=1 can be achieved.

Numerical Results in Two Space Dimensions
In this section, we present numerical results for the fractional GIADE (1) with G(x, t, u(x, t)) = u 2 (x, t) and G(x, t, u(x, t)) = 0 on the domains, as shown in Figure 1. Figure 1a is the star domain Ω 1 which is discretized with uniform distribution.The boundary of the irregular domain Ω 1 is determined by r = 1 2m 2 m 2 + 2m + 2 − 2(m + 1)cos(mθ) , where 0 ≤ θ ≤ 2π and m are assumed to be 8. Figures 1b,c are shown to be irregular geometries where r = 0.5 sin 2 (5θ) + sin 2 (2θ) + cos 2 θ and r = 0.29(e sin θ sin 2 (2θ) + e cos θ cos 2 (2θ)) define the irregular domains Ω 2 and Ω 3 with an irregular distribution of nodes, respectively.In these examples, we report the following error indicators: where u and u M are analytical and numerical solutions, respectively.The computational order of the new technique is computed by in which E 1 and E 2 are errors corresponding to grids with mesh sizes h 1 and h 2 , respectively.
Remark 1.We have performed our calculations by applying MATLAB software with a Core i7-8565 PC with a 1.80 GHz CPU and 8-GB RAM.
Example 2. We consider the fractional GIADE (1) with G(x, t, u(x, t)) = u 2 (x, t).The initial value is given by u(x, 0) = 0, x ∈ Ω, u(x, t) = t 1+γ sin(πx) sin(πy), (x, y) ∈ ∂Ω, 0 < t ≤ T, We apply the new approach to obtain the solution on the domains Ω 1 , Ω 2 , and Ω 3 .To investigate the effect of the chosen basis function and domain on the accuracy of the proposed method, the results obtained using two different RBFs are shown in Table 1.When running our algorithms using the same parameters as in Table 1, we obtain the results depicted in Figure 2 for different shape parameters c on the domains Ω 1 .It can be observed that the errors for Multiquadric and Gaussian basis functions decrease with the increase in c, even in the case of the multiquadric RBF for advection coefficients ν = 2 and ν = 8, and α = 5 and α = 10.It should be remarked that the optimal shape parameter is very difficult to determine.However, for the best accuracy of the method that is sensitive to the shape parameter c, we chose c = 5 for the multiquadrics and the Gaussian RBF.The numerical and analytical solutions and E 2 error function for the Gaussian basis function are shown in Figures 3 and 4. The first row of Figure 3 is the numerical solution, and the second one is the analytical solution.It is easy to find that the E 2 error in the middle of the domain tends to decrease.We can even observe the effect of the distributed nodes inside the domains and the behavior of the E 2 error from this figure.Different factors affect the convergence rate of the proposed methods.Obviously, increasing the number of K and K 1 improves the accuracy of the method.Now, we test the problems with multiquadric basis function for two types of geometries to show the result of our method.The E 2 errors are depicted in Figure 5, which show a little better results as we increase the number of K or K 1 .Figure 6a shows E max and E 2 errors and the effect of the fractional order at (0.5, 0.5) in Ω 2 .We can see that the error drops with the increase in fractional order γ.Moreover, from Figure 6b, it can be observed that the E max error value has an increasing trend for different time t.
The values of E 2 errors yielded from the proposed method via different values of α and domains are reported for multiquadric basis functions in Figure 7.These figures confirm that the accuracy of the method improves by a small enough α.Example 3. The following fractional GIADE with an analytical solution is considered: We chose a suitable right-hand-side function f such that the exact solution to (26) is u(x, y, t) = t 2 (1 − x) σ (1 − y) σ [9].First, we investigated the effect of the time, time step, and domain on the accuracy of the proposed method.In Table 2, we list the evolution of the solution at time T = 1 and T = 1.5 with τ = 0.04.Next, we investigated the effect of the time steps and domain on the maximum errors.A comparison of numerical solutions using the BSM method and ISEFG method that is introduced in [9,41,42] is listed in Table 3.The stability of the present method and ISEFG method can be seen from the obtained results.However, the error reduction with decreasing time steps is negligible in our results.The values of E 2 errors yielded from the proposed method via different values of α and domains are reported for multiquadric basis functions in Figure 8.Using the same parameters as in Figure 8 and the best choice of α in every domain, graphs of E 2 errors for different values of K and K 1 for multiquadric functions are shown in Figure 9. Hence, we observe that the convergence rate will be better when K and K 1 are big enough.A contour plot for the behavior of numerical and analytical solutions and the E 2 error function for the multiquadric basis function are shown in Figures 10 and  11.We can even observe the effect of distributed nodes inside the domains and the behavior of the E 2 error from these figures.When running our algorithms using the same parameters as in Figure 9, different shape parameter c and the order of the Riemann-Liouville fractional derivatives γ and M and N do not play any significant role in the solution quality in this example.To verify the theoretical results and show the efficiency and applicability of the proposed schemes, as well as the extension of the proposed schemes to the corresponding one-dimensional problems and even three-dimensional problems, a fractional GIADE with a nonhomogeneous source term in one-space dimension is presented in the next example.
Example 4. We consider the fractional GIADE with a nonhomogeneous source term on the domain The exact solution of the problem ( 27) is u(x, t) = e x t 1+γ .Similar to the previous example, the boundary and initial condition are obtained from the exact solution.First, we solve the problem by using the different parameters.Figure 12 illustrates how the shape parameters c and α play a significant role in the solution quality.It can be observed that the errors for multiquadric and Gaussian basis functions increase with the increase in c, and the good shape parameter range is between 1 and 7.Moreover, from the second row of the figure, it can be observed that the error value has an increasing trend when γ is arising from γ = 0.5.Now, we investigate convergence in time and space.In Figure 13, we illustrate the effect of the number of collocation points for this problem with τ, γ, and M fixed.It should be remarked that error reduction is negligible in our results.To evaluate the approximate solution of the problem, we illustrate Figure 14 to show the effect of the fractional order in both time and space.So, the different shape parameter c and the order of Riemann-Liouville fractional derivative γ and K played a significant role in the solution quality in this example.Finally, a comparison of our method by means of multiquadric basis functions is listed.Table 4 shows all the methods are numerically stable in terms of E max error when the first-order finite difference scheme [11], ISEFG method [9], and BSM method are applied to the problem.In addition, in Figure 15, the evolution of the solution at space x = 0.2, 0.5, 0.8 with different time step is illustrated.Obviously, the present schemes are numerically stable.Example 5. We consider the fractional GIADE with a nonhomogeneous source term on the domain The initial and boundary conditions are as follows: The numerical solutions are listed in Table 5, with τ = 0.1, c = 0.1, α = 10, M = 69, N = 37, K = 2, and K 1 = 20.We define E n = max u n+1 M − u n M to measure the stability of the time discretization scheme where u n M = u M (x, t n ), n ≥ 1.The stability of the time discretization scheme at T = 1 is presented in Figure 16.According to the lack of known reference solutions of the example, we validate the convergence rates by producing a numerical reference solution at τ = 0.001.The convergence of time discretization scheme at T = 1 is presented in Figure 17.

Conclusions
We coupled the WSGD approximation formula and the C-N technique that succeed in discretization of the involved fractional derivative of the GIADE equation along with its more general version, which is the GIADE equation with nonlinear source term.Next, we proposed the new version of the BSM, a well-established class of meshfree methods, to solve the consequent equation.In the present approach, the final approximation is given by the summation of the RBFs; the primary approximation satisfies only the boundary conditions.Moreover, a quasilinearization technique was employed to transform a nonlinear source term into a linear source term.From the verification of the efficiency and applicability of the proposed schemes, as well as an extension of the proposed schemes to the corresponding one-and two-dimensional problems, we observed acceptable accuracy depending on the shape parameter, the grid size, the domain of the problem, and the other potential factors.

Figure 2 .
Figure 2. Comparison E 2 with respect to the shape parameter c on the domains Ω 1 at time T = 1 for Example 2 where τ = 0.1.

Table 1 .
E max and E 2 for the BSM method and the effect of the chosen basis function and domain at times T = 1 for Example 2 where τ

Table 2 .
E max and E 2 for the BSM method and the effect of time and domain forExample 3, where