# Pseudospectral Method Based on Müntz–Legendre Wavelets for Solving the Abel Integral Equation

Authors: Ioannis Dassios; Fairouz Tchier; F. M. O. Tawfiq
Journal: Journal of Mathematics
 (2022)

---

## Abstract

This paper deals with the numerical solution of the Abel integral equation based on Müntz–Legendre wavelets. To this end, the Abel integral operator is represented by Müntz–Legendre wavelets as an operational matrix. To find this matrix, we use the similarity between the Abel integral operator and the fractional integral operator. The proposed method can be easily used to solve weakly singular Volterra integral equations. We have proved the convergence of the proposed method. To demonstrate the ability and accuracy of the method, some numerical examples are presented.


---

## Body

## 1. Introduction

In this paper, we focus our attention on constructing and applying Müntz–Legendre (M-L) wavelets that will be used as the basis in the pseudospectral method to solve the famous Abel integral equation(1)ux−Aαkx=fx,in which Abel’s integral operator Aαkx of order 0<α<1 is defined in [1] as follows:(2)Aαkx≔∫0xkx,s,usx−s−αdx,x∈0,1.Here, givenΩ=0,1, fx, and kx,s,us are assumed to be continuous functions on Ω and S with S=x,s,u:x,s∈0,1,u∈ℝ. Further, we suppose that the kernel function kx,s,us is equal to the form gx,sus. In other words, the desired equation is assumed to be linear.The Abel equation is a special case of the integral equations with the weakly singular kernel that was first introduced by Abel. In investigating the generalization of the tautochrone problem, he introduced this equation [2]. This equation appears widely in modeling many physical problems, such as nuclear physics, X-ray radiography, fluid flow [3], scattering theory, plasma diagnostics, semiconductors, physical electronics, and nonlinear diffusion[1, 4]. Given this equation’s wide application, solving this equation is very important. But one cannot always solve the equation analytically, and we need to use numerical methods for it.Among the many papers that have considered the numerical solution of this equation, we can mention some of them. Saadatmandi and Dehghan [5] utilized the collocation method based on shifted Legendre polynomials. Piessens and Verbaeten [6] introduced a numerical method based on Chebyshev polynomials, and after approximating the unknown solution based on these bases, they obtained the solution as a sum of hypergeometric functions. Using the Bernstein operational matrix, Singh et al. [7] introduced a stable numerical method to solve this problem. In [8], we can find the integrable solution of the Abel integral equation under certain conditions, and also the sufficient and necessary conditions for the existence of this solution are presented. In [9], the authors proposed the Laplace transform method to solve the problem, where they assumed that the solution would be differentiable and continuous. Saray [10] introduced a novel and efficient method based on Alpert’s multiwavelets. In this work, after introducing the sparse representation of the Abel integral operator, the Abel integral equation is reduced to a sparse system of linear algebraic equations in the linear form, and this causes a reduction in time and computational costs. In [11], the unbounded solutions of the nonlinear Abel integral equations are investigated. Li and Zhao [12] used Mikusinski’s operator of fractional order to solve the Abel integral equation.The outline of this article is as follows: the M-L wavelets are constructed in Section2, and then the Abel integral operator is represented based on these bases. In Section 3, the Abel integral equation is solved by using the pseudospectral method based on M-L wavelets. This section contains the error analysis and the conditions for convergence are investigated. To demonstrate the efficiency and accuracy of the presented method, some numerical examples are given in Section 4.


## 2. Müntz–Legendre Wavelets

As we know, multiresolution analysis (MRA) is a significant procedure for constructing wavelets. According to MRA, a family of nested subspaces exists such that they satisfy certain circumstances [13](3)0⊂⋯⊂V−1⊂V0⊂V1⊂⋯⊂L2Λ,where Λ is equal to ℝ or any bounded interval.Recently, the M-L wavelets have been used to solve some equations, such as fractional optimal control problems [14], fractional pantograph differential equations [15], fractional differential equations [16], and multiorder fractional differential equations [17]. To solve the Abel equation, we first briefly introduce the M-L wavelets as follows.GivenJ∈ℕ0, assume that the subspace VJ∈L2Λ is spanned by the scaled and translated version of a set of bases, which are called multiscaling functions, viz.,(4)VJ=spanϕJ,bn:b∈ℬ,n∈ℛ,where ℬ≔0,1,…,2J−1 and ℛ≔0,1,…,r−1 with r∈ℕ. The parameter J is called refinement level and r is the multiplicity parameter. In the sequel, we intend to introduce the functions ϕJ,bn.Motivated by [17], we denote the M-L polynomials Lnx as(5)Lnx=∑k=0nlk,nxλk,x∈Ω,where λk≔kμ:μ∈ℝ,k=0,…,n and the coefficient lk,n is obtained by(6)lk,n≔∏i=0n−1λk+λi+1∏i=0,i≠knλk−λi.It can be easily shown that these polynomials satisfy the orthogonality requirements and form an orthogonal system, via(7)Lnx,Ln′x=∫01LnxLn′xdx=δn′,n2λn+1,n≥n′,where δm′,m denotes the Kronecker symbol and is given by(8)δn′,n≔1,n′=n,0,n≥n′.Considering the definition ofLnx, one can introduce the M-L wavelets [17], via(9)ϕJ,bn=2J/22λn+1Ln2Jx−b,b2J≤x≤b+12J,0,otherwise.Due to the definition of M-L wavelets, one can introduce the projection operatorP that maps any function u∈L2Ω onto VJ as follows:(10)ux≈Pux=∑b=02J−1∑n=0r−1ub,nϕJ,bnx=UTΦx,where and throughout the paper, the superscript T is used for the matrix transpose. Here Φx is a vector function of dimension N=2Jr whose br+n+1-th element is ϕJ,bnx, and the br+n+1-th element of the vector U is evaluated by(11)ub,n=u,ϕJ,bn=∫01uxϕJ,bnxdx.It follows from [15] that one may be able to bound the projection error P in the sense of Sobolev norms.Lemma 1 (see [15]).
Givenn≥0, assume that r>m. If u∈HmΩ, then(12)u−PuL2Ω≤cr−1−m2J−1−mumL2Ω,and for s≥1, we have(13)u−PuHsΩ≤cr−12s−1/2−m2J−1s−mumL2Ω,in which HmΩ is the Sobolev space and the related norm is determined by(14)uHmΩ=∑j=0mujL2Ω21/2.


### 2.1. Representation of Abel Integral Operator in Müntz–Legendre Wavelets

In this subsection, we consider the Abel operator as a fractional integral operator, and after representing the fractional integral operator in M-L wavelets as an operational matrix, we find a representation of the Abel integral operator in M-L wavelets. To this end, it is necessary to define some concepts about fraction calculation.Definition 1 (see [18]).
Letu∈L1a,b. The Riemann–Liouville fractional integral operator ℐaα of order α∈ℝ+ is determined by(15)ℐaαux≔1Γα∫axx−sα−1usds,x∈a,b,where Γα represents the gamma function.Remark 1.
It can be verified that the fractional integral of the functionsxκ is given by(16)ℐaαxκ=Γκ+1Γκ+α+1xκ+α.Lemma 2 (cf Lemma 2.1(a). see [19]).
The fractional integration operatorsℐaα are bound in Lpa,b for 1≤p≤∞ as follows:(17)ℐaαup≤Kup,K≔b−aαΓα+1.It follows from [19] that if u∈L1a,b, then the function ℐaαu itself belongs to L1a,b. We recall that the Abel integral operator of order α∈0,1 is determined by(18)Aαux=∫0xx−s−αusds.There is a similarity between the Abel integral operatorAα and the Riemann–Liouville fractional integral operator ℐ0ββ≔1−α, viz,(19)Aα=Γαℐa1−α.Thus, we can use the Riemann–Liouville fractional integral operatorℐ0β instead of the Abel integral operator Aα.According to the definition of M-L wavelets, the fractional integral operatorℐ0β acting on the vector function Φx may be written as an expansion of M-L wavelets Φx, i.e.,(20)Pℐ0βΦx=IβxΦx,β∈0,1,where Iβx is an N×N matrix and is famous as the operational matrix of fractional integral for the M-L wavelets.Before we look at how to calculate the elements of the aforementioned matrix, it is intransitive to introduce the piecewise fractional-order Taylor functions. LetJ∈ℤ+∪0 be a fixed number, the piecewise fractional-order Taylor functions can be defined in the following form:(21)ψJ,bn=tλn,b2J≤x≤b+12J,0,otherwise.b∈ℬ,n∈ℳ.By introducing the square matrixT of dimension N×N, whose i,j-th element is computed by(22)Ti,j=Φix,Ψjx=∫01Φix,Ψjxdx,i,j=1,…,N.We can expand any elements of vector functionΦx (M-L wavelets) by the piecewise fractional-order Taylor functions Ψx, viz,(23)Φx=T−1Ψx.Here the matrixT is a N×N matrix and it is called the transformation matrix. In the sequel, we assume that Q is a vector of dimension r whose i-th element is xλi. Thus, it is easy to show that(24)Ψx=Q,…,QT.It follows from (16) that one can find the i-th element of ℐ0βΨx, via(25)ℐ0βΨix=Γλi+1Γλi+β+1xλi+β.This gives rise to introduce a diagonal matrixIΦ,βx, such that(26)ℐ0βΨx=IΨ,βxΨx.It is worth noting that this matrix expresses the fractional integral of functionsΨ as a combination of themselves and has the following form:(27)IΨ,βx=diagPβx,…,Pβx,where Pβx≔xβHℐ0βQx=PβxQx and H is a diagonal matrix of the form(28)Hi,i=Γλi+1Γλi+β+1.Now, we are able to introduce the operational matrix of fractional integral for the M-L wavelets via(29)PJℐ0βΦx=PJℐ0βT−1Ψx=T−1IΨ,βxΨx=T−1IΨ,βxTΦx.Thus, we get(30)Iβx≔T−1IΨ,βxT.


## 2.1. Representation of Abel Integral Operator in Müntz–Legendre Wavelets

In this subsection, we consider the Abel operator as a fractional integral operator, and after representing the fractional integral operator in M-L wavelets as an operational matrix, we find a representation of the Abel integral operator in M-L wavelets. To this end, it is necessary to define some concepts about fraction calculation.Definition 1 (see [18]).
Letu∈L1a,b. The Riemann–Liouville fractional integral operator ℐaα of order α∈ℝ+ is determined by(15)ℐaαux≔1Γα∫axx−sα−1usds,x∈a,b,where Γα represents the gamma function.Remark 1.
It can be verified that the fractional integral of the functionsxκ is given by(16)ℐaαxκ=Γκ+1Γκ+α+1xκ+α.Lemma 2 (cf Lemma 2.1(a). see [19]).
The fractional integration operatorsℐaα are bound in Lpa,b for 1≤p≤∞ as follows:(17)ℐaαup≤Kup,K≔b−aαΓα+1.It follows from [19] that if u∈L1a,b, then the function ℐaαu itself belongs to L1a,b. We recall that the Abel integral operator of order α∈0,1 is determined by(18)Aαux=∫0xx−s−αusds.There is a similarity between the Abel integral operatorAα and the Riemann–Liouville fractional integral operator ℐ0ββ≔1−α, viz,(19)Aα=Γαℐa1−α.Thus, we can use the Riemann–Liouville fractional integral operatorℐ0β instead of the Abel integral operator Aα.According to the definition of M-L wavelets, the fractional integral operatorℐ0β acting on the vector function Φx may be written as an expansion of M-L wavelets Φx, i.e.,(20)Pℐ0βΦx=IβxΦx,β∈0,1,where Iβx is an N×N matrix and is famous as the operational matrix of fractional integral for the M-L wavelets.Before we look at how to calculate the elements of the aforementioned matrix, it is intransitive to introduce the piecewise fractional-order Taylor functions. LetJ∈ℤ+∪0 be a fixed number, the piecewise fractional-order Taylor functions can be defined in the following form:(21)ψJ,bn=tλn,b2J≤x≤b+12J,0,otherwise.b∈ℬ,n∈ℳ.By introducing the square matrixT of dimension N×N, whose i,j-th element is computed by(22)Ti,j=Φix,Ψjx=∫01Φix,Ψjxdx,i,j=1,…,N.We can expand any elements of vector functionΦx (M-L wavelets) by the piecewise fractional-order Taylor functions Ψx, viz,(23)Φx=T−1Ψx.Here the matrixT is a N×N matrix and it is called the transformation matrix. In the sequel, we assume that Q is a vector of dimension r whose i-th element is xλi. Thus, it is easy to show that(24)Ψx=Q,…,QT.It follows from (16) that one can find the i-th element of ℐ0βΨx, via(25)ℐ0βΨix=Γλi+1Γλi+β+1xλi+β.This gives rise to introduce a diagonal matrixIΦ,βx, such that(26)ℐ0βΨx=IΨ,βxΨx.It is worth noting that this matrix expresses the fractional integral of functionsΨ as a combination of themselves and has the following form:(27)IΨ,βx=diagPβx,…,Pβx,where Pβx≔xβHℐ0βQx=PβxQx and H is a diagonal matrix of the form(28)Hi,i=Γλi+1Γλi+β+1.Now, we are able to introduce the operational matrix of fractional integral for the M-L wavelets via(29)PJℐ0βΦx=PJℐ0βT−1Ψx=T−1IΨ,βxΨx=T−1IΨ,βxTΦx.Thus, we get(30)Iβx≔T−1IΨ,βxT.


## 3. Pseudospectral Method

To derive the numerical solution of the second kind of Abel integral equation based on the pseudospectral method, we can approximate the unknown solution with the projection operatorPJ, as follows:(31)u≈PJu=UTΨ,where U is a vector of dimension N, whose elements should be found. Note that the function fx and the kernel function kx,s,PJus can be approximated in the same manner, i.e.,(32)fx≈PJu=FTΨ,kx,s,PJus≈PJkx,s,PJus=ΨTxKΨs,where F∈ℝN and K∈ℝN×N. By substituting (32) into the integral part of the Abel integral equation (1), we obtain(33)PJΨxTK˜∫0xx−s−αΨsds=PJrΨxTK˜IαΨx=K^Ψx.We now substitute equations (31)–(33) into Abel integral equation (1) and simplify to get(34)rx≔U−F+K^TΨx=0,where rx is the residual function that our goal is to reduce to zero. Let xi be a number of points in Ω, we select the solution that satisfies the collocation condition rxi=0, where xi are called the collocation points. In this paper, we use the shifted Chebyshev and Legendre polynomials zeros as collocation points. The collocation method gives rise to a system of linear or nonlinear algebraic equations. One can derive the unknown coefficients U after solving this system.


### 3.1. Error Analysis

We write the Abel integral equation (1) in the form(35)I−Ku=f,where K is a compact operator that maps any continuous function onto C0,1. As we said, our goal is to reduce the residual function rx to zero. Symbolically, we have(36)rx=I−KuJ−f,where uJ≔PJu. We note that PJrx=0 if and only if rxi=0 or equivalently,(37)PJI−KuJ=PJf.LetuJ is a solution of (37), then by applying PJuJ=uJ, equation (37) can be written as follows:(38)I−PJKuJ=PJf.Since both the original equations (1) and (38) are defined on C0,1, then for the error analysis, we compare them.Theorem 1.
Let us assume thatI−K:C0,1⟶C0,1 is a bijections operator. Further, let us assume(39)K−PJK⟶0,asr⟶∞.Then for all sufficiently largerr≤N, the operator I−PJK−1 exists as a bounded operator from C0,1 to C0,1. Moreover, it is uniformly bounded(40)supr≤NI−PJK−1<∞.For the solutions of (35) and (38)(41)u−uJ=I−PJK−1u−PJu,1I−PJKu−PJu≤u−uJ≤I−PJK−1u−PJu.Proof.
For details, refer to [20].
SinceK:C0,1⟶C0,1 is a compact operator and PJ is a bounded projection, such that PJu⟶0 as r⟶∞. Then motivated by [20] (Lemma 3.1.2), we have(42)K−PJK⟶0,asr⟶∞.
Therefore, the condition of Theorem1 is held.


## 3.1. Error Analysis

We write the Abel integral equation (1) in the form(35)I−Ku=f,where K is a compact operator that maps any continuous function onto C0,1. As we said, our goal is to reduce the residual function rx to zero. Symbolically, we have(36)rx=I−KuJ−f,where uJ≔PJu. We note that PJrx=0 if and only if rxi=0 or equivalently,(37)PJI−KuJ=PJf.LetuJ is a solution of (37), then by applying PJuJ=uJ, equation (37) can be written as follows:(38)I−PJKuJ=PJf.Since both the original equations (1) and (38) are defined on C0,1, then for the error analysis, we compare them.Theorem 1.
Let us assume thatI−K:C0,1⟶C0,1 is a bijections operator. Further, let us assume(39)K−PJK⟶0,asr⟶∞.Then for all sufficiently largerr≤N, the operator I−PJK−1 exists as a bounded operator from C0,1 to C0,1. Moreover, it is uniformly bounded(40)supr≤NI−PJK−1<∞.For the solutions of (35) and (38)(41)u−uJ=I−PJK−1u−PJu,1I−PJKu−PJu≤u−uJ≤I−PJK−1u−PJu.Proof.
For details, refer to [20].
SinceK:C0,1⟶C0,1 is a compact operator and PJ is a bounded projection, such that PJu⟶0 as r⟶∞. Then motivated by [20] (Lemma 3.1.2), we have(42)K−PJK⟶0,asr⟶∞.
Therefore, the condition of Theorem1 is held.


## 4. Numerical Examples

In this section, some numerical examples are solved to show the validity and efficiency of the method. To do this, we carry out the Maple and MATLAB software simultaneously.Example 1.
For the first example, let us consider the linear Abel integral equation of the second kind with the kernel functionkt,x,ux≔1/10Γ1−αux, and ft=1. The exact solution is given in [10] as follows:(43)ut=E1−α,1t1−α10,where Eσ,β is the Mittag-Leffler function(44)Eσ,βz=∑l=0∞zlΓσl+β,σ,β,z∈ℝ,σ>0.
Table1 shows the absolute value of errors at different times xi when the collocation points are chosen to be the Legendre polynomial nodes. As we expected, when the r increases (the degree of the bases as well as the number of collocation points increases) the error will decrease. To show the effect of the parameter μ in the L2-error, we plot Figure 1 and report Table 2. Figure 2 illustrates the effect of the choosing nodes on the L2-error and also absolute value of errors at the Chebyshev nodes using different multiplicity r. Also, we can see the effect of the multiplicity parameter r on the L2-error and absolute error in Figure 2.Table 1
The absolute errors at the Legendre nodes for Example1.
xr=5r=7r=90.19.72e−42.78e−41.19e−40.22.41e−42.24e−41.30e−50.33.15e−41.22e−64.64e−50.42.93e−41.09e−44.83e−50.51.23e−57.38e−71.93e−60.61.94e−47.95e−63.14e−50.71.53e−46.80e−62.21e−50.87.94e−56.32e−56.04e−60.91.89e−44.44e−51.60e−51.04.09e−41.60e−47.86e−5Figure 1
Plot of theL2-errors taking μ=1/2 (a) and μ=1/3 (b) at the Legendre and Chebyshev nodes for Example 1.
(a)(b)Table 2
TheL2-errors at the Legendre and Chebyshev nodes taking different μ for Example 1.
rμ=1/2μ=1/3Chebyshev nodesLegendre nodesChebyshev nodesLegendre nodes27.57e−47.87e−44.03e−34.48e−331.69e−52.06e−59.34e−59.79e−543.82e−75.32e−79.81e−51.27e−458.58e−91.33e−83.37e−54.29e−561.89e−103.16e−101.76e−52.26e−574.06e−127.19e−129.89e−61.27e−5Figure 2
Plot of theL2-errors (a) taking the Legendre and Chebyshev nodes and the absolute error at the Chebyshev nodes (b) for Example 1.
(a)(b)Example 2.
The second example is devoted to the Abel integral equation (1) with(45)ft≔2t,kt,x,ux≔ux.
The exact solution is given byut=1−eπterfcπt [4].
To show the effect of the multiplicity parameterr and choosing the collocation points, we report Table 3. This Table also illustrates the effect of parameter μ on the L2-error. Due to Table 3, it is obvious that these three parameters have a direct effect on the L2-error such that when r increases, the error decreases. Also, choosing the Chebyshev nodes gives us a better result than Legendre nodes. In Figure 3, we demonstrate the absolute error when multiplicity parameter r increases taking μ=1/2 and different collocation points.Table 3
TheL2-errors at the Legendre and Chebyshev nodes taking different μ for Example 2.
rμ=1/2μ=1Chebyshev nodesLegendre nodesChebyshev nodesLegendre nodes23.38e−23.52e−25.68e−25.62e−239.21e−21.08e−23.15e−23.15e−242.57e−33.36e−31.95e−22.00e−257.29e−41.04e−31.31e−21.38e−262.06e−43.17e−49.40e−31.01e−275.77e−59.40e−57.04e−37.73e−3Figure 3
Plot of the absolute error at the Chebyshev nodes (a) and the Legendre nodes (b) taking differentr for Example 2.
(a)(b)Example 3.
Let us consider the following Abel integral equation:(46)ft≔sinx−43x3/2F121;54,74;−x24,kt,x,ux≔ux,in which F12 is the hypergeometric function defined in [21]. Also, the exact solution is ux=sinx.
Table4 shows the absolute value of errors at different times xi when the collocation points are chosen to be the Chebyshev polynomials nodes. Figure 4 illustrates the effect of selecting the Chebyshev and Legendre nodes. We can also see the effect of increasing the parameter r. It is observed that by increasing the parameter r, the error decreases.Table 4
The absolute errors at the Chebyshev nodes for Example3.
xr=6r=8r=100.11.42e−41.20e−61.60e−70.26.05e−62.45e−61.53e−90.32.39e−63.44e−74.50e−80.44.67e−56.55e−76.19e−80.58.29e−51.61e−67.59e−80.61.03e−41.87e−61.08e−70.71.23e−42.31e−61.48e−70.81.66e−43.36e−62.01e−70.92.40e−44.64e−62.76e−71.03.21e−46.34e−63.78e−7Figure 4
Plot of theL2-errors taking μ=1/2 at the Legendre and Chebyshev nodes for Example 3.


## 5. Conclusion

In this paper, we utilize an efficient algorithm based on the wavelet pseudospectral method to solve the well-known Abel integral equation. This method can easily be used to solve weakly singular Volterra integral equations, and this shows the ability of the proposed method. We have compared the method with other methods and shown that this method offers better results. We have proved the convergence of the proposed method. Given the construction of these bases and the role of the parameterμ, which can be polynomials with fractional powers, compared to other bases, if the exact solution or the known functions in the equation are of the fractional type, the proposed method will provide better results.



---

*Source: 2251623-2022-03-04.xml*