Spectral Element Method for Fractional Klein – Gordon Equations Using Interpolating Scaling Functions

,


Introduction
Our aim in this study is to implement and develop the Galerkin method to approximate the solution of the fractional Klein-Gordon equation: subjected to the following conditions: where γ 1 and γ 2 are constants, qðx; tÞ is a continuous function, and N is a known function that fulfills the following Lipschitz condition: where L N is referred to the Lipschitz constant.The fractional derivative is of the Caputo fractional derivative (CFD) type, such that we shall introduce it in the sequel.
For many years, positive integer derivatives have been used in partial differential equations (PDEs).But, in recent years, we have noticed a trend toward using fractional derivatives to model physical phenomena.It seems like there are many advantages to this approach, including greater accuracy and more flexibility in the types of systems that can be modeled.Meanwhile, some applications of these equations can be mentioned such as colored noise [1], bioengineering [2][3][4], solid mechanics [5], continuum and statistical mechanics [6], earthquakes [7], anomalous transport [8], economics [9], fluid-dynamic traffic model [10], and engineering and natural sciences [11][12][13].There are some analytical methods to solve these types of equations [14][15][16].However, when the equations become more complicated, these methods no longer work.So, numerical approaches can address this deficiency.Here, we mention some of these methods, including finite difference method [17], collocation method [18,19], Galerkin method [20,21], finite element method [22], kernelbased pseudo-spectral method [23], nonuniform difference schemes [24], and fractional differential transform method [25].
It is common knowledge that the Klein-Gordon equations, both linear and nonlinear, are widely used when it comes to modeling various physical phenomena.They have been used to model everything from solitons to condensed matter physics, as well as classical and quantum mechanics.In 1926, two physicists named Klein [26] and Gordon [27] introduced this equation to describe relativistic electrons.According to the importance of these equations, studying and finding their solution can be very useful.To this end, through the utilization of different numerical techniques, we have able to solve these equations, such as the He's method [28,29], radial basis functions [30], decomposition method [31], meshless method [32], Galerkin method [33], Tau method [34], differential transform method [35], b-spline scaling function [36], and Adomians decomposition scheme [37].
The fractional Klein-Gordon equation has been garnering significant attention from scientists as of late.In fact, fractional Klein-Gordon equation is the generalization of fractional Klein-Gordon equation of integer order.On the other hand, this equation describes nonlocal relationships in space and time using power-law memory kernels.According to the definition of fractional derivative, it is worth noting that the fractional derivative at a point depends on all values of the function.Thus, it is expected that the fractional derivative operation involves some sort of boundary conditions, involving information on the function further out [38].As we know, when the fractional order approaches an integer, the fractional derivative converges to the ordinary derivative.So, we can conclude that the fractional differential equations are close to the real modeling of phenomena and the differential equations of integer order are the limits of them.Several numerical schemes have been developed by mathematicians to aid in solving this equation.It is remarkable to witness the rapid advancements being made in this field.The Sinc-Chebyshev based collocation method described in Nagy's [39] study is an effective means of solving the nonlinear version of this equation.I came across a fascinating article by Golmankhaneh et al. [40] on using the homotopy perturbation method to solve the problem.In Singh et al.'s [41] study, the authors employed Chebyshev polynomials of the third kind to implement the collocation method for finding the solution to this problem.Cubic B-spline functions are used to solve this problem by a numerical method [42].One of the approach that has been found to be efficient for finding numerical solutions for equations of this type is the meshless method suggested by Gharian et al. [43].The problem at hand may be resolved through the utilization of generalized polynomials, as presented in Hassani et al.'s [44] study.
Interpolating scaling functions (ISFs) are a family of powerful bases that emerge from a multiresolution analysis (MRA).These bases are orthonormal and by dilation and translation generate a family of nested subspaces of L 2 ½0; 1.These bases are also essential for constructing Alpert's multiwavelets.These bases were first introduced by Alpert et al. [45] where PDEs were solved by an adaptive method using these bases.ISFs have shown their adaptability and usefulness in a wide range of mathematical problems, such as solving ordinary differential equations [46,47], PDEs [48], integral equations [49], and more.
In this paper, we use the spectral element method based on ISFs to solve the fractional Klein-Gordon equation.The fundamental idea behind this method is to consider the unknown solution as a linear combination of basis functions.Then, using the operational matrix of CFD and spectral element method, the desired equation reduces to a system of algebraic equations.Considering the flexibility of the ISFs in the selection of m and s parameters, the presented method solves this type of equation with appropriate accuracy.As you know, the main challenges in equations such as the Klein-Gordon equation are the nonlinearity and existence of the fractional derivative.When we use other types of bases, these two factors cause some limitations in implementing the spectral element method.However, since there is no need to calculate the integral in using the ISFs to approximate the functions, the problems of nonlinearity can be easily solved.About the fractional derivatives, we know that the equation may have a nonsmooth solution near the boundaries, which many numerical methods fail to overcome.The scaling functions, due to the properties that they inherit from multiresolution analysis, by increasing the scale parameter s can overcome this challenge.
In this article, we have organized the information into different sections.Section 2 provides an introduction to ISFs and their properties, as well as a matrix representation of CFD.In Section 3, we explain and execute the proposed method, which utilizes the spectral element method.Additionally, we present a theorem that demonstrates the convergence of this method.Section 4 contains results from various numerical experiments.Finally, we complete this work through the inclusion of a conclusion in Section 5.

Interpolating Scaling Bases
Before implementing the algorithm of the proposed scheme, it is worthwhile to give the bases' construction briefly.To this end, according to our knowledge of these bases, these bases are piecewise polynomials of the degree less than a constant number, called multiplicity parameter m [45,50,51].The main configuration of these bases is made by Lagrange polynomials L l ðxÞ, which are used in their construction from the roots of Legendre polynomials P m ðxÞ.Suppose that these roots are indicated by fτ l ; l ¼ 0; 1; ⋯; m − 1g.Furthermore, assume that the Gauss-Legendre quadrature weights ω l , l ¼ 0; 1; ⋯; m − 1, determined by ω l : ¼ 2=ðrP 0 m ðτ l ÞP m−1 ðτ l ÞÞ, is given.The ISFs are specified by the following equation: Motivated by L 2 -inner product, the ISFs consist of orthonormal bases for the finite-dimensional space: By dividing the interval Ω : ¼ ½0; 1 into subintervals Advances in Mathematical Physics which s 2 N ∪ f0g is given and B s : ¼ f0; 1; ⋯; 2 s − 1g, we are now ready to introduce the subspaces that generate a MRA.Recall that MRA consists of nested subspaces that satisfy some conditions cf [52].Here, we designate these subspaces by the following equation: where T k and D 2 s denote 9the translation and dilation operators, respectively, and are determined by the following equation: To approximate any function f , using an orthonormal projection P m s ð f Þ: L 2 ðΩÞ→A m s , it can be mapped onto A m s via a finite sum: where s ≥ 0 is a fixed integer number and the coefficients h f ; ϕ l s; k i can be computed by the following equation: As you see, computing these integrals lead to high computational cost, so to avoid this, the interpolating property of the bases could be a problem solver.This property gives rise to calculating these coefficients without integration, i.e.: , consisting of all scaling functions in A m s .According to this presentation, Equation ( 9) may be rewritten as follows: where the elements of N-dimensional vector F s are determined by ½ F s kmþlþ1 ¼ f l s; k , in which N ¼ m; 2 s .The superscript "T" in this text indicates the meaning of "transpose." There is an error in this approximation.In order to determine an estimate for the maximum error, we can refer to the Lemma 1, [50].
Lemma 1.Given m 2 N, assume that f : ½0; 1→R be of class C m , m times continuously differentiable function.The error of approximation ( 9) is estimated by the following equation: To extend the approximation to the 2D space A m; 2 s ⊂ L 2 ðΩÞ 2 , which is generated by fϕ l s; k ðxÞϕ l 0 s; k 0 ðtÞ : l; l 0 ¼ 0; …m − 1; k; k 0 2 B s g, the projection P m s can be utilized to approximate the 2D function f ðx; tÞ, viz.: in which the coefficients F mkþðlþ1Þ; mk 0 þðl 0 þ1Þ and l; l 0 ¼ 0; …m − 1, k; k 0 2 B s are computed by the following equation: But, similar to the 1D case, to avoid calculating these integrals, they can be found using the following equation: Lemma 2 (cf Theorem 1, Saray et al. [53]).Given m 2 N, assume that f : Ω × Ω → R be a sufficiently smooth function.
The approximation error ( 14) can be limited by the following equation: where: Advances in Mathematical Physics 2.1.Matrix Representation of CFD.Before introducing a matrix representation of CFD ( C D η 0 ), it is necessary to determine some preliminaries about fractional calculus.Additionally, since the matrix representation of the CFD is derived from the matrix representation of the fractional integral (FI), it is necessary to incorporate this matrix.Definition 1.Given η 2 R þ , the FI operator I η 0 of order η is specified by the following equation: where ΓðηÞ indicates the Gamma function.
Given the power function, its fractional integration is also a power function, that is: The norm of the FI operator is bounded, as motivated by Kilbas et al. [54].Lemma 3.There is an estimation of the bound of the fractional integral operator I η 0 in L q ð½0; 1Þ, viz.: The matrix representation of the fractional integral operator is denoted by I η .In order to determine I η , it is necessary to once again estimate the fractional integral of the vector function Φ m s ðxÞ using ISFs, viz.: where I η is a N × N square matrix.The objective is to identify the elements present in this matrix.To proceed, we have to take help from the following lemma.
Lemma 4 (cf, Asadzadeh and Saray [55]).The Lagrange polynomials L l for nodes fτ l g l¼0; 1; …; m−1 may be written as follows: Based on Equation (11), it is possible to compute the elements of matrix I η using the following equation: The elements of matrix I η is previously obtained by Asadzadeh and Saray [55].They proved this matrix is an upper triangular matrix in which entries can be computed by considering the three following cases.
(1) Given k 0 <k, we can easily prove that: 26) leads to the following equation: According to the closed form ISFs and Lemma 4, we can obtaine the following equation: Now, the elements situated on the diagonal I η are specified by the following equation: Furthermore, the beta function is denoted by B, and λ : So, it can be show that: Using the same argument as in Case (2), we get the following equation: Applying the binomial expansion of ð2λy − 1Þ m−1−n as follows: the entries of I η above the main diagonal can be specified the following equation: Assume that AC η ð½0; 1Þ is a space of functions such that: As we know, if wðxÞ 2 AC η ½0; 1, then the CFD: exists for almost every x 2 ½0; 1.Now, our objective is to find an operational matrix D η for the operator c D η 0 such that it satisfies the following equation: To find this matrix, motivated by Equation ( 37), the operational matrix I η can be utilized as follows: where the matrix D is used to represent the derivative operation, as explained in Alpert et al.'s [45] study.Consequently, the matrix D η can be determined by the following equation:

Methods Description
Recall that our objective is to introduce and implement an effective numerical algorithm for Equation (1).The approach used in this numerical algorithm involves implementing the Galerkin method.The method begins by expanding the unknown solution wðx; tÞ based on ISFs, viz.: or equivalently: in which we have a square matrix, W, of order N, with coefficients that are currently unknown.Substituting Equation (42) into Equation (1) leads to the following equation: Using Equation (42) and applying the operational matrices D and D η , we can approximate all terms of Equation ( 43), as follows.
(i) Taking CFD from Equation (41), and using the matrix D η , one can write the following equation: (ii) For the second one, namely, ∂ 2 ∂x 2 w m s ðx; tÞ, we use the matrix D, as follows: (iii) Given N m : ¼ Nðx; t; w m s ðx; tÞÞ, let us imagine a matrix F that has N × N dimensions and meets the following condition: Note that the entries of the matrix F consist of the linear or nonlinear equations of the unknowns W i; j .
(iv) Using the operator P m s , the given function qðx; tÞ can also map into the space A m s , viz.: Putting Equations ( 44)-( 47) back into Equation ( 43) gives rise to introducing the residual: Applying the Galerkin method requires fulfilling the following equation: or equivalent to the following equation: Because the ISFs are orthonormal bases, the aforementioned equation leads to a linear or nonlinear system: Now, the functions f 0 , f 1 , f 2 , and f 3 must be mapped to space A m s to apply the boundary and initial conditions.So, we have the following equation: By rewriting Equation ( 51) using the following equation: To gain the unknown matrix W, the generalized minimal residual method (GMRES method) [56] and Newton's method are used to determine the unknown vector W for the linear and nonlinear forms, respectively.

Convergence Verification
Theorem 1. Assume that the functions w and f are sufficiently smooth.Furthermore, if f satisfies the Lipschitz condition (4), thus the presented method for Equation ( 1) is convergent when the parameters m or s tend to infinity.
Proof.Subtracting Equation ( 1) from the following equation: the residual is determined by the following equation: Taking the norm from both sides and using the triangle inequality, we get the following equation: Given e ¼ w − w m s .By utilizing Lemmas 2, 3 and Equation (37), it can be verified that all expressions (57) have bounds, viz.: s Þðx; tÞjj 2 , we have the following equation: 6Advances in Mathematical Physics where C∂ κ w ∂t κ is a constant and can be obtained by replacing f with ∂ κ w ∂t κ in Lemma 2.
(ii) Employing Lemma 2 and the Lipschitz condition (4) gives rise to the following equation: where C w is constant.
(iii) For other expressions, using Lemma 2, it is straightforward to establish the following equation: To proceed, by substituting Equations ( 58)-( 60) in Equation (57), we obtain the following equation: ; C N m ; C q g.Therefore:

Illustrative Numerical Examples
This section includes some examples to showcase the effectiveness of the proposed method.To illustrate the results and make a global view of the present method and its efficiency, sometimes, the absolute errors: and L 2 error: are reported in tables or plotted in figures.
Example 1.We dedicate the first example to the following equation: subjected to the following conditions: in which: wðx; tÞ ¼ t 2 sinðxÞ provides the precise solution for this equation [57].
The approximate solution and corresponding absolute error for different choices of m, taking s ¼ 2 and η ¼ 1:5, are plotted in Figure 1. Figure 2 illustrates the accuracy of the presented method with m ¼ 8, s ¼ 2, and η ¼ 1:5.The L 2error is also tabulated in Table 1 with different choices of m, taking η ¼ 1:5 and s ¼ 2. From our observation, the results clearly showcase the effectiveness and precision of the method that was presented.To compare the present method with the method presented in Dehghan et al.'s [57]   Advances in Mathematical Physics show the efficiency and accuracy of the method, Table 2 is reported.
Example 2. We devote this example to the following equation: with the following boundary and initial conditions: The precise solution for this equation is assigned by wðx; tÞ ¼ t 2 ðe − e x ÞsinðxÞ [58].
To provide evidence of the method's accuracy, Table 3 is tabulated.According to our analysis, when the parameter m increases, the error must be reduced.Table 4 is reported to show the accuracy at different points.To this end, we can find the absolute error in this table.Compared to other existing methods, it is worth mentioning that the proposed scheme provides better accuracy than the high-order difference method [58].This superiority is reported in Table 5 with different eta.The approximate solution and absolute error are illustrated in Figure 3 for different m.
Example 3. The third example is dedicated to the nonlinear fractional Klein-Gordon equation: The precise solution for this equation is assigned by wðx; tÞ ¼ x 3 t 3 [59].
By setting η ¼ 1:5, m ¼ 4, and s ¼ 2, the L 2 -errors and maximum absolute error are reported in Table 6 at different times.To compare the presented method with others, Table 7 is reported.Compared to other methods including clique polynomials [59], tension spline approach [60] and radial basis functions [30], the presented method gives better accuracy.Figure 4 is provided to demonstrate the approximate solution and corresponding absolute error using the presented method with m ¼ 4, s ¼ 2, and η ¼ 1:8.

Conclusion
Due to the use of fractional derivatives and nonlinearity, solving the fractional Klein-Gordon equation is a challenging task.This work uses the spectral element method to solve the problem by mapping it to an approximation subspace of L 2 ½0; 1 through an orthonormal projection.The ISFs that possess specific properties and satisfy the MRA are made use in this method.To this end, a matrix representation of the CFD of ISFs is presented.We prove the method is convergent, and several examples have been prepared to verify the findings of this investigation.The results speak for themselves, and it is evident that this method presents the effective scheme that solve the problem in the accuracy manner.Additionally, when compared to other methods, it outperforms them in terms of accuracy and speed.These findings are truly remarkable and show the immense potential of this approach.The method proposed here has the potential to solve both fractional and nonfractional equations with ease.Its simplicity in implementation, combined with its high efficiency and significant accuracy, make it a strong candidate for solving the same equations.
To summarize, the following key points can be mentioned: (i) The method has a simple and practical structure, making it easy to extend for solving various fractional equations.(ii) The presented method can solve any kind of equation due to its flexibility in choosing scale and multiplicity parameters.(iii) The speed of calculations increases because there is no need for integration due to the properties of the bases used.(iv) Increasing the value of the scale parameter can lead to more appropriate accuracy for problems with nonsmooth solutions at the beginning and end points of the solution interval.
In the future, we plan to extend our numerical approaches for solving generalized fractional models, including the timefractional diffusion equation [61], the time-fractional mobileimmobile advection-dispersion equation [62], the timefractional mobile/immobile transport model [63], and multiterm boundary value problem of variable order [64], etc.

TABLE 1 :
In Example 1, the L 2 errors were calculated at various times using different values of m, with η set to 1:5 and s set to 2.

TABLE 2 :
The proposed method with the implicit RBF meshless method are compared in Example 1.

TABLE 3 :
The L 2 errors at various times with different choices of m, taking η ¼ 1:75 and s ¼ 2 for Example 2.

TABLE 4 :
The results obtained by present method compared with the high-order difference method for Example 2.