 # Numerical solution of a class of third-kind Volterra integral equations using Jacobi wavelets

We propose a spectral collocation method, based on the generalized Jacobi wavelets along with the Gauss-Jacobi quadrature formula, for solving a class of third-kind Volterra integral equations. To do this, the interval of integration is first transformed into the interval [-1,1], by considering a suitable change of variable. Then, by introducing special Jacobi parameters, the integral part is approximated using the Gauss-Jacobi quadrature rule. An approximation of the unknown function is considered in terms of Jacobi wavelets functions with unknown coefficients, which must be determined. By substituting this approximation into the equation, and collocating the resulting equation at a set of collocation points, a system of linear algebraic equations is obtained. Then, we suggest a method to determine the number of basis functions necessary to attain a certain precision. Finally, some examples are included to illustrate the applicability, efficiency, and accuracy of the new scheme.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

In this paper, we consider the following Volterra integral equation (VIE)

 tβu(t)=f(t)+∫t0(t−x)−ακ(t,x)u(x)dx,t∈[0,T], (1)

where , , , with a continuous function , and is continuous on the domain and is of the form

 κ(t,x)=xα+β−1κ1(t,x),

where is continuous on . The existence of the term in the left-hand side of equation (1) gives it special properties, which are distinct from those of VIEs of the second kind (where the left-hand side is always different from zero), and also distinct from those of the first kind (where the left-hand side is constant and equal to zero). This is why in the literature they are often mentioned as VIEs of the third kind. This class of equations has attracted the attention of researchers in the last years. The existence, uniqueness, and regularity of solutions to (1) were discussed in Sonia1 . In that paper, the authors have derived necessary conditions to convert the equation into a cordial VIE, a class of VIEs which was studied in detail in Vainikko1 ; Vainikko2 . This made possible to apply to equation (1) some results known for cordial equations. In particular, the case is of special interest, because in this case, if , the integral operator associated to (1) is not compact and it is not possible to assure the solvability of the equation by classical numerical methods. In Sonia1 , the authors have introduced a modified graded mesh and proved that with such mesh the collocation method is applicable and has the same convergence order as for regular equations.

In Nemati , two of the authors of the present paper have applied a different approach, which consisted in expanding the solution as a series of adjusted hat functions and approximating the integral operator by an operational matrix. The advantage of that approach is that it reduces the problem to a system of linear equations with a simple structure, which splits into subsystems of three equations, making the method efficient and easy to implement Nemati . A limitation of this technique is that the optimal convergence order () can be attained only under the condition that the solution satisfies , which is not the case in many applications.

It is worth remarking here that there is a close connection between equations of class (1) and fractional differential equations MyID:410 . Actually, the kernel of (1) has the same form as the one of a fractional differential equation, and if we consider the case , then the integral operator is the Riemann–Liouville operator of order . Therefore, it makes sense to apply to this class of equations numerical approaches that have recently been applied with success to fractional differential equations and related problems MyID:410 .

One of these techniques were wavelets, a set of functions built by dilation and translation of a single function , which is called the mother wavelet. These functions are known as a very powerful computational tool. The term wavelet was introduced by Jean Morlet about 1975 and the theory of the wavelet transform was developed by him and Alex Grossmann in the 80s Grossmann ; Grossmann1 . Some developments exist concerning the multiresolution analysis algorithm based on wavelets Daubechies2 and the construction of compactly supported orthonormal wavelet basis Mallat . Wavelets form an unconditional (Riesz) basis for , in the sense that any function in can be decomposed and reconstructed in terms of wavelets Daubechies3 . Many authors have constructed and used different types of wavelets, such as B-spline XLi , Haar Chen , Chebyshev YLi , Legendre Jafari , and Bernoulli Rahimkani wavelets. The advantage of employing wavelets, in comparison with other basis functions, is that when using them one can improve the accuracy of the approximation in two different ways: (i) by increasing the degree of the mother function (assuming that it is a polynomial); (ii) by increasing the level of resolution, that is, reducing the support of each basis function.

We underline that the application of wavelets has special advantages in the case of equations with non-smooth solutions, as it is the case of equation (1). In such cases, increasing the degree of the polynomial approximation is not a way to improve the accuracy of the approximation; however, such improvement can be obtained by increasing the level of resolution.

In a recent work Nemati2 , Legendre wavelets were applied to the numerical solution of fractional delay-type integro-differential equations. In the present paper we will apply a close technique to approximate the solution of (1).

The paper is organized as follows. Section 2 is devoted to the required preliminaries for presenting the numerical technique. In Section 3, we give some error bounds for the best approximation of a given function by a generalized Jacobi wavelet. Section 4 is concerned with the presentation of a new numerical method for solving equations of type (1). In Section 5, we suggest a criterion to determine the number of basis functions. Numerical examples are considered in Section 6 to confirm the high accuracy and efficiency of this new numerical technique. Finally, concluding remarks are given in Section 7.

## 2 Preliminaries

In this section, we present some definitions and basic concepts that will be used in the sequel.

### 2.1 Jacobi wavelets

The Jacobi polynomials , , , are a set of orthogonal functions with respect to the weight function

 w(ν,γ)(t)=(1−t)ν(1+t)γ,

with the following orthogonality property:

 ∫1−1w(ν,γ)(t)P(ν,γ)i(t)P(ν,γ)j(t)dt=h(ν,γ)iδij,

where is the Kronecker delta and

 h(ν,γ)i=2ν+γ+1Γ(ν+i+1)Γ(γ+i+1)i!(ν+γ+2i+1)Γ(ν+γ+i+1).

The Jacobi polynomials include a variety of orthogonal polynomials by considering different admissible values for the Jacobi parameters and . The most popular cases are the Legendre polynomials, which correspond to , Chebyshev polynomials of the first-kind, which correspond to , and Chebyshev polynomials of the second-kind, which correspond to .

We define the generalized Jacobi wavelets functions on the interval as follows:

 ψ(ν,γ)n,m(t)=⎧⎪⎨⎪⎩2k2√1h(ν,γ)mTP(ν,γ)m(2kTt−2n+1),n−12k−1T≤t

where is the level of resolution, , , is the degree of the Jacobi polynomial, and is the normalized time. The interested reader can refer to Boggess ; Walnut for more details on wavelets. Jacobi wavelet functions are orthonormal with respect to the weight function

 w(ν,γ)k(t)=⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩w(ν,γ)1,k(t),0≤t<12k−1T,w(ν,γ)2,k(t),12k−1T≤t<22k−1T,⋮w(ν,γ)2k−1,k(t),2k−1−12k−1T≤t

where

 w(ν,γ)n,k(t)=w(ν,γ)(2kTt−2n+1),n=1,2,…,2k−1.

An arbitrary function may be approximated using Jacobi wavelet functions as

 u(t)≃Ψ(ν,γ)k,M(t)=2k−1∑n=1M∑m=0un,mψ(ν,γ)n,m(t),

where

 un,m=⟨u(t),ψ(ν,γ)n,m(t)⟩w(ν,γ)k=∫T0w(ν,γ)k(t)u(t)ψ(ν,γ)n,m(t)dt=∫n2k−1Tn−12k−1Tw(ν,γ)n,k(t)u(t)ψ(ν,γ)n,m(t)dt.

For a given function , the Gauss–Jacobi quadrature formula is given by

 ∫1−1(1−t)ν(1+t)γu(t)dt=N∑l=1ωlu(tl)+RN(u),

where , , are the roots of , , , are the corresponding weights given by (see Shen ):

 ωl=2ν+γ+1Γ(ν+N+1)Γ(γ+N+1)N!Γ(ν+γ+N+1)(ddtP(ν,γ)N(tl))2(1−t2l), (2)

and is the remainder term which is as follows:

 RN(u)=2ν+γ+2N+1N!Γ(ν+N+1)Γ(γ+N+1)Γ(ν+γ+N+1)(ν+γ+2N+1)(Γ(ν+γ+2N+1))2×u(2N)(η)(2N)!,η∈(−1,1). (3)

According to the remainder term (3), the Gauss–Jacobi quadrature rule is exact for all polynomials of degree less than or equal to . This rule is valid if possesses no singularity in . It should be noted that the roots and weights of the Gauss–Jacobi quadrature rule can be obtained using numerical algorithms (see, e.g., Pang ).

## 3 Best approximation errors

The aim of this section is to give some estimates for the error of the Jacobi wavelets approximation of a function

in terms of Sobolev norms and seminorms. With this purpose, we extend to the case of Jacobi wavelets some results which were obtained in Canuto for the best approximation error by Jacobi polynomials in Sobolev spaces. The main result of this section is Theorem 9, which establishes a relationship between the regularity of a given function and the convergence rate of its approximation by Jacobi wavelets.

We first introduce some notations that will be used in this paper. Suppose that is the space of measurable functions whose square is Lebesgue integrable in relative to the weight function . The inner product and norm of are, respectively, defined by

 ⟨u,v⟩w∗=∫baw∗(t)u(t)v(t)dt,∀ u,v∈L2w∗(a,b),

and

 ∥u∥L2w∗(a,b)=√⟨u,u⟩w∗.

The Sobolev norm of integer order in the interval , is given by

 ∥u∥Hrw∗(a,b)=(r∑j=0∥u(j)∥2L2w∗(a,b))12, (4)

where denotes the th derivative of and is a weighted Sobolev space relative to the weight function .

For ease of use, for some fixed values of , we set

 ψi,j(t):=ψ(ν,γ)i,j(t), w(t):=w(ν,γ)(t), wk(t):=w(ν,γ)k(t), wn,k(t):=w(ν,γ)n,k(t).

For starting the error discussion, first, we recall the following lemma from Canuto .

[See Canuto ] Assume that with and denotes the truncated Jacobi series of , where is the space of all polynomials of degree less than or equal to . Then,

 ∥u−LM(u)∥L2w(−1,1)≤CM−μ|u|Hμ;Mw(−1,1), (5)

where

 ∣u∣Hμ;Mw(−1,1)=⎛⎝μ∑j=min{μ,M+1}∥u(j)∥2L2w(−1,1)⎞⎠12 (6)

and is a positive constant independent of the function and integer . Also, for , one has

 ∥u−LM(u)∥Hrw(−1,1)≤CM2r−12−μ|u|Hμ;Mw(−1,1). (7)

Suppose that denotes the set of all functions whose restriction on each subinterval , , are polynomials of degree at most . Then, the following lemma holds.

Let , , be functions in with . Consider the function defined by for all . Then, for , we have

 ∥(¯un)(j)∥L2w(−1,1)=(2kT)12−j∥u(j)n∥L2wn,k(Ik,n).

Using the definition of the -norm and the change of variable , we have

 ∥¯un(j)∥2L2w(−1,1)=∫1−1w(t)|¯un(j)(t)|2dt=∫1−1w(t)∣∣∣un(j)(T2k(t+2n−1))∣∣∣2dt=∫n2k−1Tn−12k−1Twn,k(t′)(2kT)−2j∣∣u(j)n(t′)∣∣2(2kT)dt′=(2kT)1−2j∥u(j)n∥2L2wn,k(Ik,n),

which proves the lemma.

In order to continue the discussion, for convenience, we introduce the following seminorm for , , and , which replaces the seminorm (6) in the case of a wavelet approximation:

 ∣u∣Hr;μ;M;kwk(0,T)=⎛⎝μ∑j=min{μ,M+1}(2k)2r−2j∥u(j)∥2L2wk(0,T)⎞⎠12. (8)

If we choose such that , it can be easily seen that

 ∣u∣Hr;μ;M;kwk(0,T)=(2k)r−μ∥u(μ)∥L2wk(0,T). (9)

The next theorem provides an estimate of the best approximation error, when Jacobi wavelets are used, in terms of the seminorm defined by (8). Suppose that with and

 Ψk,M(u)=2k−1∑n=1M∑m=0un,mψn,m(t),

is the best approximation of based on the Jacobi wavelets. Then,

 ∥u−Ψk,M(u)∥L2wk(0,T)≤CM−μ|u|H0;μ;M;kwk(0,T) (10)

and, for ,

 ∥u−Ψk,M(u)∥Hrwk(0,T)≤CM2r−12−μ|u|Hr;μ;M;kwk(0,T), (11)

where in (10) and (11) the constant denotes a positive constant that is independent of and but depends on the length . Consider the function such that for all . Then, from (4) and Lemma 3 for , we have

 (12)

By setting in (12), we obtain

 ∥u−Ψk,M(u)∥2L2wk(0,T)≤C12k−1∑n=1(2k)−1∥¯un−(LM(¯un))∥2L2w(−1,1)≤C2M−2μ(2k)−12k−1∑n=1μ∑j=min{μ,M+1}∥∥¯u(j)n∥∥2L2w(−1,1)≤CM−2μμ∑j=min{μ,M+1}(2k)−2j2k−1∑n=1∥∥u(j)n∥∥2L2wn,k(In,k)=CM−2μμ∑j=min{μ,M+1}(2k)−2j∥∥u(j)∥∥2L2wk(0,T),

where we have used (5) and Lemma 3. This completes the proof of (10). Furthermore, using (12) for and , we get:

 ∥u−Ψk,M(u)∥2Hrwk(0,T) ≤C1(2k)2r−12k−1∑n=1r∑j=0∥∥¯u(j)n−(LM(¯un))(j)∥∥2L2w(−1,1) =C1(2k)2r−12k−1∑n=1∥¯un−LM(¯un)∥2Hrw(−1,1) ≤C2M4r−1−2μ(2k)2r−12k−1∑n=1μ∑j=min{μ,M+1}∥∥¯u(j)n∥∥2L2w(−1,1) =C2M4r−1−2μ(2k)2r−1μ∑j=min{μ,M+1}2k−1∑n=1∥∥¯u(j)n∥∥2L2w(−1,1) ≤CM4r−1−2μμ∑j=min{μ,M+1}(2k)2r−2j2k−1∑n=1∥∥u(j)n∥∥2L2wn,k(Ik,n) =CM4r−1−2μμ∑j=min{μ,M+1}(2k)2r−2j∥∥u(j)∥∥2L2wk(0,T),

where we have used (4), (7), and Lemma 3. Therefore, we have proved (11).

We can also obtain estimates for the Jacobi wavelets approximation in terms of the -norm. With , if we combine (9) with (10), we get

 ∥u−Ψk,M(u)∥L2wk(0,T)≤CM−μ2−μk∥u(μ)∥L2wk(0,T);

from (9) and(10), we obtain

 ∥u−Ψk,M(u)∥Hrwk(0,T)≤CM2r−12−μ(2k)r−μ∥u(μ)∥L2wk(0,T),r≥1.

## 4 Method of solution

In this section, we propose a method for solving the VIE (1). To this end, by using a suitable change of variable, we transform the interval of the integral to . Suppose that

 s=2(xt)−1,ds=2tdx.

Therefore, the equation (1) is transformed into the following integral equation:

 tβu(t)=f(t)+(t2)1−α∫1−1(1−s)−ακ(t,t2(s+1))u(t2(s+1))ds. (13)

In order to compute the integral part of (13), we set and as the Jacobi parameters and use the Gauss–Jacobi quadrature rule. Then, we have

 tβu(t)=f(t)+(t2)1−αN∑l=1ωlκ(t,t2(sl+1))u(t2(sl+1)), (14)

where are the zeros of and are given using (2) as

 ωl=21−α(ddxP(−α,0)N(sl))2(1−s2l),l=1,2,…,N.

We consider an approximation of the solution of (1) in terms of the Jacobi wavelets functions as follows:

 u(t)≃2k−1∑i=1M∑j=0ui,jψ(ν,γ)i,j(t), (15)

where the Jacobi wavelets coefficients are unknown. In order to determine these unknown coefficients, we substitute (15) into (14) and get

 (16)

In this step, we define the following collocation points

 tn,m=T2k(τm+2n−1),n=1,2,…,2k−1, m=0,1,…,M,

where , , are the zeros of . Therefore, are the shifted Gauss–Jacobi points in the interval , corresponding to the Jacobi parameters and . By collocating (16) at the points , we obtain

 2k−1∑i=1M∑j=0[tβn,mψ(ν,γ)i,j(tn,m)−(tn,m2)1−αN∑l=1ωlκ(t,tn,m2(sl+1))×ψ(ν,γ)i,j(tn,m2(sl+1))]ui,j=f(tn,m). (17)

By considering and , in the above equation, we have a system of linear algebraic equations that can be rewritten as the following matrix form:

 AU=F, (18)

where

 U=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣u1,0u1,1⋮u1,M⋮u2k−1,0u2k−1,1⋮u2k−1,M⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦,F=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣f(t1,0)f(t1,1)⋮f(t1,M)⋮f(t2k−1,0)f(t2k−1,1)⋮f(t2k−1,M)⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦,

and the entries of the rows of the matrix are the expressions in the bracket in (17), which vary corresponding to the values of and , i.e., the coefficients of , , , for , that are all nonzero and positive. Since the functions are orthonormal and the nodes are pairwise distinct, the matrix is nonsingular. Therefore, (18) is unique solvable. By solving this system using a direct method, the unknown coefficients are obtained. Finally, an approximation of the solution of (1) is given by (15).

## 5 A criterion for choosing the number of wavelets

Now we discuss the choice of adequate values of and (number of basis functions). To do this, we suppose that and . Using the error of the Gauss–Jacobi quadrature rule given by (3), and substituting and in it, the exact solution of (1) satisfies the equation

 tβu(t)=f(t)+(t2)1−α[N∑l=1ωlκ(t,t2(sl+1))u(t2(sl+1))+RN(κu)],

where

 RN(κu)=2−α+2N+1(N!)2(Γ(−α+N+1))2(2N)!(−α+2N+1)(Γ(−α+2N+1))2×(t2)2N∂2N(κ(t,x)u(x))∂x2N|x=η,

for . Therefore, we have

 tβu(t)=f(t)+(t2)1−α(N∑l=1ωlκ(t,t2(sl+1))u(t2(sl+1)))+t1−α+2Nξα,N∂2N(κ(t,x)u(x))∂x2N|x=η,

where

 ξα,N=(N!)2(Γ(−α+N+1))2(2N)!(−α+2N+1)(Γ(−α+2N+1))2.

Suppose that . Then we obtain that

 u(t)=g(t)+2α−1t1−α−β(N∑l=1ωlκ(t,t2(sl+1))u(t2(sl+1)))+t1−α−β+2Nξα,N∂2N(κ(t,x)u(x))∂x2N|x=η. (19)

Let be the numerical solution of (1) obtained by the proposed method in Section 4. From the definition of the Jacobi wavelets, for the collocation points , , , we have

 Uk,M(tn,m)=M∑j=0un,jψ(ν,γ)n,j(tn,m),n=1,…,2k−1.

By definition, the restriction of the functions , , on the subinterval , which we denote here by , is smooth. Therefore, we can define

 ζn,j=maxx,t∈Ik,n∣∣ ∣ ∣∣∂2N(κ(t,x)ρ(ν,γ)n,j(x))∂x2N∣∣ ∣ ∣∣.

For a given , since all the collocation points, , , , are positive, using (19), we can choose and such that for all , the following criterion holds:

 ∣∣ ∣∣Uk,M(tn,m)−g(tn,m)−2α−1t1−α−βn,m(N∑l=1ωlκ(tn,m,tn,m2(sl+1))×Uk,M(tn,m2(sl+1)))∣∣∣+t1−α−β+2Nn,mξα,N∣∣ ∣∣M∑j=0un,jζn,j∣∣ ∣∣<ε.

## 6 Numerical examples

In this section, we consider three examples of VIEs of the third-kind and apply the proposed method to them. The weighted -norm is used to show the accuracy of the method. In all the examples, we have used in the Gauss–Jacobi quadrature formula and the following notation is used to show the convergence of the method:

 Ratio=e(k−1)e(k),

where ) is the -error obtained with resolution .

As the first example, we consider the following third-kind VIE, which is an equation of Abel type Sonia1 ; Nemati :

 t2/3u(t)=f(t)+∫t0√33πx1/3(t−x)−2/3u(x)dx,t∈[0,1],

where

 f(t)=t4712⎛⎝1−Γ(13)Γ(5512)π√3Γ(5912)⎞⎠.

The exact solution of this equation is , which belongs to the space . We have employed the method for this example with different values of , , and , and reported the results in Tables 1, 2 and Figure 1. Table 1 displays the weighted -norm of the error for three different classes of the Jacobi parameters, which include (second-kind Chebyshev wavelets), (Legendre wavelets), and (first-kind Chebyshev wavelets) with different values of and . Moreover, the ratio of the error versus is given in this table. It can be seen from Table 1 that the method converges faster in the case of the second-kind Chebyshev wavelets. In Table 2, we compare the maximum absolute error at the collocation points obtained by our method with the error of the collocation method introduced in Sonia1 , and the operational matrix method based on the adjusted hat functions Nemati . From this table, it can be seen that our method gives more accurate results with less collocation points (we have used ) than the method of Sonia1 () and also, has higher accuracy with a smaller number of basis functions (we have used ) than the method of Nemati (). In Figure 1, we show the error function obtained by the method based on the second-kind Chebyshev wavelets with , (left) and , (right).