# On the Multiwavelets Galerkin Solution of the Volterra–Fredholm Integral Equations by an Efficient Algorithm

Authors: H. Bin Jebreen
Journal: Journal of Mathematics
 (2020)

---

## Abstract

We develop the multiwavelet Galerkin method to solve the Volterra–Fredholm integral equations. To this end, we represent the Volterra and Fredholm operators in multiwavelet bases. Then, we reduce the problem to a linear or nonlinear system of algebraic equations. The interesting results arise in the linear type where thresholding is employed to decrease the nonzero entries of the coefficient matrix, and thus, this leads to reduction in computational efforts. The convergence analysis is investigated, and numerical experiments guarantee it. To show the applicability of the method, we compare it with other methods and it can be shown that our results are better than others.


---

## Body

## 1. Introduction

A mathematical model of the spatiotemporal development of an epidemic yields the following Volterra–Fredholm integral equation (VFIE):(1)ux=fx+∫0xk1x,s,usds+∫01k2x,s,usds,x∈Ω≔0,1,where the given functions f:Ω⟶ℝ and k1:S×ℝ⟶ℝ with S≔x,s:x,s∈Ω are assumed to be continuous functions. Furthermore, we consider k2≔px,shus to be integrable, where hx is a nonlinear function and px,s:S⟶ℝ is a continuous function. Besides, the given functions are selected so that equation (1) has a unique solution.The parabolic boundary value problems lead to these types of integral equations and widely arise from various physical and biological models. The VFIE also appears in the literature in mixed form as(2)ux=fx+∫0x∫01kt,susdsdt,where f and k are analytic functions. Many authors studied the mixed form of VFIE numerically. Among these, we can mention collocation method [1], projection method [2], spline collocation method [3], wavelet collocation method [4], Adomian decomposition method [5], and so on [6–8]. Among these studies, we focus on a paper that uses the multiwavelet Galerkin method to solve linear mixed VFIE as mentioned in [9]. In [9], the wavelet transform matrix and the operational matrix of integration are utilized to reduce the problem of linear mixed VFIE to a sparse linear system of algebraic equations. By searching among the sources, we can find a small number of papers in the field of numerical and analytical solutions to problem (1). In [10], the Lagrange collocation method is employed to solve this problem. Wang and Wang [11] applied the Taylor collocation method to find the numerical solution of the equation. Also, the convergence analysis is investigated for the proposed method. The Taylor collocation method was applied by Karamete and Sezer [12] to solve this equation as well as the high-order linear Fredholm–Volterra integrodifferential equations [13]. The Bell polynomials have been employed for solving this equation [14].The motivation of our work is to develop the multiwavelet Galerkin method to solve (1). We split the problem into two configurations, linear and nonlinear. After using the multiwavelet Galerkin method, both types reduce to the system of the linear and nonlinear algebraic equations, respectively. The interesting results arise in the linear type where thresholding is employed to decrease the nonzero entries of the coefficient matrix. This gives the sparse system. This property is very useful to reduce the computational cost. We use Alpert’s multiwavelet bases constructed in [15] following [16], and these bases have been used to solved PDEs, ODEs, and integral equations [9, 17–19]. These bases allow us to have high vanishing moments, compact support, and properties such as orthogonality and interpolation [15]. These characteristics of Alpert’s multiwavelets lead to a sparse representation of differential and integral operators [16, 20].This paper is organized as follows: In Section2, we briefly introduce Alpert’s multiwavelet bases. In Section 3, the Multiwavelets Galerkin method is used to solve this proble, and the conditions for convergence of the proposed method are discussed. Section 5 contains some numerical results to confirm the validity and efficiency of the proposed method, and Section 6 contains a few concluding remarks.


## 2. Alpert’s Multiwavelet Bases

LetJ∈ℤ+∪0. We consider the uniform finite discretizations Ω≔0,1=∪b∈ℬXJ,b where the subintervals XJ,b≔xb,xb+1 are determined by the point xb≔b/2J with ℬ≔0,…,2J−1. For k∈ℛ≔0,1,…,r−1, we introduce the subspace VJr as a space of piecewise polynomial bases of degree less than multiplicity parameter r that is spanned by(3)VJr≔spanϕj,bk≔D2jTbϕk,b∈ℬj,k∈ℛ⊂L2Ω,r≥0,where D and T are the dilation and translation operators, respectively, and ϕkk∈ℛ are the primal interpolating scaling bases introduced by Alpert et al. [16]. Given nodes τkk∈ℛ, which are the roots of Legendre polynomial of degree r, the interpolating scaling bases are defined as(4)ϕkt=2ωkLk2t−1,t∈0,1,0,otherwise,where Lktk∈ℛ are the Lagrange interpolating polynomials at the point τkk∈ℛ and ωkk∈ℛ are the Gauss–Legendre quadrature weights [16, 18]. These bases form an orthonormal bases on Ω with respect to the L2-inner product. Due to the definition of the space VJr, the spaces Vjrj∈ℤ+∪0 have dimension N:=2Jr and obviously are nested:(5)V0r⊂V1r⊂⋯⊂VJr⊂⋯⊂L2Ω.Hence, we consider the complement subspaceWJr of VJr in VJ+1r such that(6)VJ+1r=VJr⊕WJr,WJr⊥VJr,where ⊕ denotes orthogonal sums. According to (6), the space VJ may be inductively decomposed to(7)VJr=V0r⊕⊕j=0J−1Wjr.The complement subspaceWJr has dimension 2Jr and is spanned by multiwavelet bases ψkk∈ℛ, as(8)WJr=spanψJ,bk≅W2JWbψk:b∈ℬJ,k∈ℛ.Because Alpert’s multiwavelets are completely introduced in [16], we avoid this and refer the readers to [15, 16, 19].Every functionp∈L2Ω can be represented in the form(9)p≈VJrp=∑b∈ℬJ∑k∈ℛpJ,bkϕJ,bk,where PJr is the orthogonal projection that maps L2Ω onto the subspace VJr. To find the coefficients pJ,bk that are determined by p,ϕJ,bk=∫XJ,bfxϕJ,bkxdx, we shall compute the integrals. We apply the r-point Gauss–Legendre quadrature by a suitable choice of the weights ωk and nodes τk for k∈ℛ to avoid these integrals [16, 19], via(10)pJ,bk≈2−J/2ωk2p2−Jτk+12+b,k∈ℛ,b∈ℬJ.Convergence analysis of the projectionPJrp is investigated for the r-times continuously on differentiable function p∈ℂrΩ:(11)PJrp−p≤2−Jr24rr!supx∈0,1prx.For the full proof of this approximation and further details, we refer the readers to [15]. Thus, we can conclude that PJrp converges to p with rate of convergence O2−Jr.Assume that the vector functionΦJr:=Φr,J,0,…,Φr,J,2J−1T with Φr,J,b:=ϕJ,b0,…,ϕJ,br−1 includes the scaling functions and is called multiscaling function. Approximation (9) may be rewritten using the vector P that includes entries Pbr+k+1:=pJ,bk as follows:(12)ΦJrp=PTΦJr,where P is an N-dimensional vector. The building blocks of these bases construction can be applied to approximate a higher-dimensional function. To this end, one can introduce the two-dimensional subspace VJr,2:=VJr×VJr⊂L2Ω2 that is spanned by(13)ϕJ,bkϕJ,b′k′:b,b′∈ℬJ,k,k′∈ℛ.Thus, by this assumption, to derive an approximation of the functionp∈L2Ω2 by the projection operator PJr, we have(14)px,y≈ΦJrpx,y=ΦJrTxPΦJry,where the components of the square matrix P of order N are obtained by(15)Prb+k+1,rb′+k′+1≈2−Jωk2ωk′2p2−Jτ^k+b,2−Jτ^k′+b′,where τ^k=τk+1/2. Consider the 2r-th partial derivatives of p:Ω2⟶ℝ to be continuous. Utilizing this assumption, the error of this approximation can be bounded as follows:(16)PJrp−p≤ℳmax21−rJ4rr!2+21−Jr4rr!,where ℳmax is a constant.Given orthogonal projection operatorQjr=Pj+1r−Pjr that maps L2Ω onto Wjr, the multiscale projection operator ℳJr can be represented as(17)ℳJr=Q0r+∑j=0J−1Qjr,and consequently, any function p∈L2Ω can be approximated as a linear combination of multiwavelets and single-scale interpolating scaling functions as(18)p≈ℳJrp=∑k=0r−1p0,0kϕ0,0k+∑j=0J−1∑b∈ℬj∑k∈ℛp˜j,bkψj,bk,where(19)p0,0k≅p,ϕ0,0k,p˜j,bk≅p,ψj,bk.Note that we can compute the coefficientsp0,0k by using (10). But in many cases, multiwavelet coefficients from zero up to higher-level J−1 must be evaluated numerically. To avoid this, we use multiwavelet transform matrix TJ, introduced in [18, 19]. This matrix connects multiwavelet bases and multiscaling functions, via(20)ΨJr=TJΦJr,where ΨJr:=Φr,0,b,Ψr,0,b,Ψr,1,b,…,Ψr,J−1,bT is a vector with the same dimension ΦJr (here Ψr,j,b:=ψj,b0,…,ψj,br−1). This representation helps to rewrite equation (18) as to form(21)p≈ℳJrp=P˜JTΨJr,where we have the N-dimensional vector P˜J whose entries are p0,0k and p˜j,bk and is given by employing the multiwavelet transform matrix TJ as P˜J=TJPJ.The multiwavelet coefficients (details) become small when the underlying function is smooth (locally) with increasing refinement levels. If the multiwavelet bases haveNψr vanishing moments [19, 21], then details decay at the rate of 2−JNψr [17]. Because vanishing moments of Alpert’s multiwavelets are equal to r, consequently, one can obtain p˜J,bk≈O2−Jr. This allows us to truncate the full wavelet transforms while preserving most of the necessary data. Thus, we can set to zero all details that satisfy a certain constraint ε using thresholding operator Cε:(22)CεP˜J=P¯J,and the elements of P¯J are determined by(23)p¯j,bk≅p˜j,b,kj,b,k∈Dε,0,else,b∈ℬj,j=0,…,J−1,k∈ℛ,where Dε:=j,b,k:p˜j,bk>ε. Now, we can bound the approximation error after thresholding [17] via(24)PJrp−PJ,DεrpL2Ω≤Cthrε,where PJ,Dεrp is the projection operator after thresholding with the threshold ε and Cthr>0 is a constant independent of Jandε.


## 3. Multiwavelet Galerkin Method

In order to obtain multiwavelet Galerkin solution of (1), assume that solution can be approximated as an expansion of the Alpert’s multiwavelets, i.e.,(25)ux≈PJrux=UTΨJrx,where the N dimension vector U of unknowns must be specified. This solution is selected such that it satisfies (1) approximately. Also, it is obtained from the solution of the minimization problem(26)u−PJru=minz∈L2Ωu−z.SinceL2Ω is an inner product space with finite dimension, it can be shown that this minimization problem has a unique solution [22].Let us rewrite (1) in the operator form(27)I−ℱu=g,where gx:=fx+Vux with Vux:=∫0xk1x,s,usds. Furthermore, the operator ℱu:=∫01px,shusds is assumed to be compact on L2Ω to L2Ω and k1 is a given continuous function. Due to these assumptions, Vux is a continuous function and, consequently, g is also continuous function. Due to (11), the function gx can be approximated at a rate of at least 2−Jr:(28)gx≈PJrgx=FTΦJrx+K1TΦJrx=FTTJ−1ΨJrx+K1TTJ−1ΨJrx=F˜TΨJrx+UTA1ΨJrx,linear,C1ΨJrx,nonlinear,=G˜TΨJrx=ℳJrxx,where A1 and C1 are N×N matrices and the rest are N-dimensional vectors. Now, we introduce the residual in the approximation of (1):(29)rJx=uJx−ℳJrgx−ℳJrℱJrx,where uJ:=ℳJru and(30)ℳJrℱJruJ=ℳJr∫01ℳJrk2x,s,uJsds=ℳJr∫01ℳJrpx,shuJsds=∫01HTΨJrsΨJrTsPΨJrxds=HTPΨJrx=UTΨ2ΨJrx,linear,Ψ2ΨJrx,nonlinear.Symbolically,(31)rJ=ℳJrI−ℱJruJ−g.To finduJ, it requires that the approximate solution uJ satisfies(32)rJ,ΨJri=0,i=1,…,N.This is multiwavelet Galerkin’s method and yields a linear or nonlinear system that must be solved to obtain the approximate solution. Note thatℳJrz=0 if and only if z,ΨJri=0. Thus, we can rewrite (31) as(33)ℳJrrJ=0,or equivalently,(34)ℳJrI−ℱJruJ=ℳJrg.Note thatℳJruJ=uJ whenever uJ∈VJr. Due to this, equation (34) can be rewritten as(35)I−ℳJrℱJruJ=ℳJrg.According to (35), we obtain the system of linear or nonlinear equations. Due to the higher vanishing moments r and increasing refinement level J, for the linear type of this equation, we can discard coefficients by hard thresholding introduced in the previous section. We can reduce the computational efforts using proper methods for this type of system such as the GMRES method. The GMRES method is introduced by Saad and Shultz [23] for solving sparse and large linear systems. The GMRES generates an approximate solution whose residual norm is minimum by using a Krylov subspace. In this paper, we use restarted GMRES Algorithm 2 [24]. To use this method, we must first define Arnoldi’s algorithm. Arnoldi’ s procedure is an algorithm for building an orthogonal basis of the Krylov subspace κm. The N-th Krylov subspace is defined as follows:(36)κmΛ,w=spanw1,Λw1,…,Λm−1w1.Here, we assume that system (35) of the linear type to be of form ΛU=D and Wm is a N×m matrix with column vectors w1,…,wm. Also, H¯m is a m+1×m Hessenberg matrix whose nonzero entries hi,j are defined by Algorithm 1.Algorithm 1: Arnoldi”s algorithm.
(1)
Choose a vectorw1, such that w12=1(2)
Forj=1,2,…,m do(3)
Computehi,j=Λwj,wi for i=1,2,…,j(4)
wj=Λwj,wi for i=1,2,…,j(5)
hj+1,j=wj2(6)
Ifh2j+1,j=0 then stop(7)
wj+1=wj/hj+1,j(8)
End doAlgorithm 2: Restarted (GMRES).
(1)
Computer0=D−ΛU0, β=r02andw1=r0/β(2)
Generate the Arnoldi basis and the matrixH¯m using the Arnoldi algorithm starting with w1(3)
Computeym the minimizer of βe1−H˜my2 and Um=U0+Wmym(4)
If satisfied then stop, else setU0=Um and go to 1(1)
Computer0=D−ΛU0, β=r02andw1=r0/β(2)
Generate the Arnoldi basis and the matrixH¯m using the Arnoldi algorithm starting with w1(3)
Computeym the minimizer of βe1−H˜my2 and Um=U0+Wmym(4)
If satisfied then stop, else setU0=Um and go to 1To investigate the convergence analysis, one can prove thatℱJr−ℳJrℱJr=O2−Jr. Thus, ℱJr−ℳJrℱJr⟶0 when J⟶∞ because ℱ is a compact operator. Now, we can raise the convergence theorem.Theorem 1.
Letℱ be a compact operator and I−ℱ be injective. Assume that the sequence ℱJr:L2Ω⟶L2Ω is collectively compact and pointwise convergent to ℱ.Then,I−ℱJr−1 exists and is uniformly bounded. Also, the solution of (27) and (31) satisfy the error estimate(37)u−uJ≤I−ℱJr−11−ɛJ0I−ℱJr−1u−ℳJru+ℳJrℱ−ℳJrℱJru.Proof.
Because the sequenceℱJr converges pointwise to ℱ in L2Ω and is collectively compact, we conclude that(38)ℱJr−ℱℱJr⟶0,asJ⟶∞.
For all sufficiently largeJ, we have(39)I−ℱ−1ℱJr−ℱℱJr<1,and as a consequence of this, I−ℱJr is reversible. Note that the inverse operator I−ℱ exists due to the Riesz theorem. The investigation is based on the approximation of I−ℱJr by I−ℳJrℱJr,(40)I−ℳJrℱJr=I−ℱJr+ℱJr−ℳJrℱJr=I−ℱJrI+I−ℱJr−1ℱJr−ℳJrℱJr.
To prove the existence ofI−ℳJrℱJr−1, assume that(41)ɛJ0≡supJ≥J0ℱJr−ℳJrℱJr<1I−ℱJr−1.
Thus,I+I−ℱJr−1ℱJr−ℳJrℱJr−1 exists and is uniformly bounded due to geometric series theorem, i.e.,(42)I+I−ℱJr−1ℱJr−ℳJrℱJr−1≤11−ɛJ0I−ℱJr−1.
Using (42), I−ℳJrℱJr exists. Taking norm and using (42), we obtain(43)I−ℳJrℱJr−1≤I−ℱJr−11−ɛJ0I−ℱJr−1.
Applying the operatorℳJr to both sides of (27) and then rearranging, we can obtain(44)I−ℳJrℱJru=ℳJrg+u−ℳJru+ℳJrℱ−ℳJrℱJru.
Subtracting (35) from (44), we obtain(45)I−ℳJrℱJru−uJ=u−ℳJru+ℳJrℱ−ℳJrℱJru.
Taking norm and employing (43),(46)u−uJ≤I−ℱJr−11−ɛJ0I−ℱJr−1u−ℳJru+ℳJrℱ−ℳJrℱJru.
This is equivalent to (37). It is straightforward to show that ℳJrℱ−ℳJrℱJru⟶0 as J⟶∞. Assume that uJ be a sequence of continuous functions so that uJ⟶u as J⟶∞. Since the orthonormal projection ℳJr satisfies ℳJr=1, we can obtain(47)u−ℳJru≤u−uJ+uJ−ℳJruJr+ℳJru−uJ≤2u−uJr+uJr−ℳJruJr.
Thus, for each real numberɛ>0, there exists a number J0 such that, for every number J0≤J, one can write u−uJ≤ɛ/4. This then implies that(48)u−ℳJru≤ɛ2+uJr−ℳJruJr.
This implies thatu−ℳJru≤ɛ, for sufficiently large value of J and because ɛ is arbitrary, consequently, ℳJru⟶u as J⟶∞.


## 4. Numerical Experiments

To verify the accuracy and efficiency of the proposed method, we consider a series of numerical examples. In the linear type of equation (1), we aim to generate a sparse matrix to reduce the computational costs. We illustrate the rate of sparsity Sε which is defined by [21,25](49)Sε=N0−NεN0×100%,where ε is the threshold (small positive number) and Nε and N0 are the number of elements remaining after thresholding and the total number of elements, respectively.Example 1.
Consider the linear VFIE given in [11] as(50)ux=e−x−exx−1+∫0xex+susds−∫01ex+susds.
The exact solution ise−x [11].
The effects of the multiplicity parameterr, the refinement level J, and thresholding with different thresholds are reported in Table 1 and Figure 1. The results confirm the theoretical claims and demonstrate the effectiveness of the method. Note that the L2 error decreases as parameters r and J increase due to the rate of convergence O2−Jr. We compare the error of the proposed method, the Lagrange collocation method [10], Taylor collocation method [11], and Taylor polynomial method [26] in Table 2. Due to Table 2, our method is flexible than other methods, and without changing the multiplicity parameter r, we can improve the results. Figure 2 illustrate the effect of thresholding with different threshold parameters ε on the coefficient matrix. It can be seen that the number of matrix elements decreases when the threshold parameter increases.Table 1
Effects of parametersr, J, and ε on sparsity and L2-error for Example 1.
rJWithout thresholdingε=10−5ε=10−3ε=10−1SεL2-errorSεL2-errorSεL2-errorSεL2-error5206.39e−924.81.20e−6542.74e−4912.96e−2302.00e−1059.84.87e−678.82.78e−496.52.96e−27205.97e−13431.72e−665.91.38e−497.57.98e−3304.71e−1572.42.65e−685.41.38e−493.87.98e−3Figure 1
Effects of the refinement levelJ and multiplicity parameter r for Example 1.
(a)(b)Table 2
Comparison of the error for Example1.
rPresented methodTaylor collocationTaylor polynomial methodLagrange collocation21.55e−37.87e−23.41e−27.87e−256.39e−96.23e−53.68e−46.23e−584.63e−151.89e−81.24e−51.77e−794.01e−162.35e−83.46e−77.21e−6Figure 2
The effect of thresholding with thresholdε=10−5 (left), ε=10−3 (middle), and ε=10−1 (right) on the coefficient matrix for Example 1.
(a)(b)(c)Example 2.
Let us consider the following VFIE:(51)ux=e2x−12cosx+ex+1−e22sinx+∫0xescosxusds−∫01essinxusds.
In Table3 and Figure 3, we show the effect of the parameters r, J, and ε on sparsity and L2-error. It is obvious that increasing the parameters r and J reduces the error. A comparison of the proposed method with other methods such as Taylor method [26] and the Lagrange collocation [10] is reported in Table 4. In Figure 4, we illustrate the effect of thresholding on the coefficient matrix by taking different threshold ε when r=8 and J=2.Table 3
Effects of parametersr, J, and ε on sparsity and L2-error for Example 2.
rJWithout thresholdingε=10−6ε=10−4ε=10−2SεL2-errorSεL2-errorSεL2-errorSεL2-error5201.73e−822.21.00e−651.51.74e−479.25.15e−3305.43e−1058.11.18e−676.21.74e−491.95.15e−36201.81e−1032.55.24e−1059.57.81e−582.11.37e−2302.84e−1265.31.02e−681.77.80e−593.71.37e−27201.62e−1239.46.78e−1265.19.95e−584.81.84e−2301.27e−1471.32.74e−784.79.95e−5951.84e−2Figure 3
Effects of the refinement levelJ and multiplicity parameter r for Example 2.
(a)(b)Table 4
Comparison of the error for Example2.
rPresented methodLagrange collocationTaylor method24.17e−31.06e−24.51e−251.73e−81.12e−62.98e−481.27e−149.35e−76.14e−7Figure 4
The effect of thresholding with thresholdε=10−5 (left), ε=10−4 (middle), and ε=10−3 (right) on the coefficient matrix for Example 2.
(a)(b)(c)Example 3.
Let us consider nonlinear VFIE (1) with(52)ux=sinx+1−esin1x3−x2−sin2x4−∫0xx−su2sds−∫01x3cosseusds.
The exact solution of this equation isux=sinx.
In Figure5, we illustrate the effect of parameter r and J on the L2-error and the error of approximation is plotted in Figure 6 taking r=7 and J=2.Figure 5
Effects of the refinement levelJ and multiplicity parameter r for Example 3.
(a)(b)Figure 6
Error of Example3.
(a)(b)Example 4.
Consider the nonlinear VFIE as(53)ux=fx+∫0xcosx+su2sds+∫01sins−xu2sds,with(54)fx=−2cos2x−sin2xe2x5+2e2sin−1+x5+e2cos−1+x5+cosx5+ex−sinx5.
The exact solution isux=ex.
Table5 is reported to show the efficiency and accuracy of the proposed method. We observe when the refinement level J and multiplicity parameter r increase, the L2-errors decrease. Figure 7 shows the error of proposed method on taking r=6 and J=2.Table 5
Effect of the refinement levelJ and multiplicity parameter r on L2-error for Example 4.
rJ=2J=3J=4J=5J=6J=724.18e−31.04e−32.61e−46.53e−51.68e−54.15e−638.77e−51.10e−51.67e−61.71e−72.19e−82.68e−941.38e−68.92e−85.09e−93.31e−102.11e−111.32e−1251.73e−85.71e−101.79e−115.28e−131.61e−145.16e−1661.81e−102.23e−124.64e−146.95e−161.15e−171.69e−19Figure 7
Error of Example4.
(a)(b)


## 5. Conclusion

We have employed the multiwavelet Galerkin method to solve the Volterra–Fredholm integral equations. To this end, the Volterra and Fredholm operators are represented in multiwavelet bases. Applying this method leads to a linear or nonlinear system of algebraic equations. In the linear type, we obtain a new sparse system using thresholding due to the decay in the wavelet coefficients. The convergence analysis is investigated, and one can show that the rate of convergence isO2−Jr. The numerical examples illustrate the efficiency and accuracy of the method.



---

*Source: 2672683-2020-11-11.xml*