Multi-Wavelets Galerkin Method for Solving the System of Volterra Integral Equations

Abstract

In this work, an efficient algorithm is proposed for solving the system of Volterra integral equations based on wavelet Galerkin method. This problem is reduced to a set of algebraic equations using the operational matrix of integration and wavelet transform matrix. For linear type, the computational effort decreases by thresholding. The convergence analysis of the proposed scheme has been investigated and it is shown that its convergence is of order $O\left(2^{-Jr}\right)$, where J is the refinement level and r is the multiplicity of multi-wavelets. Several numerical tests are provided to illustrate the ability and efficiency of the method.

1. Introduction

In this paper, we study and construct a novel numerical algorithm for the system of Volterra integral equations of the second kind

$$\mathbf{u}\left(x\right)=\mathbf{f}\left(x\right)+{\int }_0^x\mathbf{g}(x,t,\mathbf{u}\left(t\right))dt,x\in \mathsf{\Omega }:=[0,1],$$

(1)

where $\mathbf{f}:\mathsf{\Omega }\to {\mathbb{R}}^n$$(n\in \mathbb{N})$ is a given real-valued continuous function, $\mathbf{u}:\mathsf{\Omega }\to {\mathbb{R}}^n$ is the unknown function that will be determined and the function $\mathbf{g}:S\to {\mathbb{R}}^n$ with $S=\left\{\right(x,t):x,t\in \mathsf{\Omega }\}$ is a given linear or nonlinear function of $\mathbf{u}$ which satisfies the following Lipschitz condition with respect to the third variable: for all $x,t\in [0,1]$ and for all ${\mathbf{u}}_1,{\mathbf{u}}_2\in {\mathbb{R}}^n$,

$$|\mathbf{g}(x,t,{\mathbf{u}}_1\left(t\right))-\mathbf{g}(x,t,{\mathbf{u}}_2\left(t\right))|\le A|{\mathbf{u}}_1-{\mathbf{u}}_2|.$$

(2)

Therefore, the functions $\mathbf{f}$ and $\mathbf{g}$ are considered so that the Equation (1) has a unique solution.

Equation (1) is the general form of second-order Volterra integral equation and appears in scientific applications in chemistry, engineering, mathematics, and physics [1,2,3,4]. Numerical and analytical solutions of linear and nonlinear Volterra integral equations have been investigated in many papers. A useful method to solve such equations is the Adomian decomposition method. This method was used to investigate the existence and uniqueness of solutions of this type of equation [5,6]. One of the best paper which utilizes the multi-wavelets for solving integro-differential equations was presented by Saray [7]. In [7], an efficient algorithm was proposed for solving the Volterra integro-differential equation. This method outperforming former approaches. Golbabai et al. [8] developed a general method based on radial basis function networks to solve the system of Volterra integral equations. The modified homotopy perturbation method for solving this type of equation has been proposed by Aminikhah et al. [5,9]. Kılıçman et al. [10] used Simpson’s $3/8$ rule to solve this equation. Aguilar and Brunner used collocation techniques based on spline polynomials [11]. The umbral calculus and the Laplace transform methods were used as solution approaches as well [12].

Wavelets and specially multi-wavelets Galerkin method represent an efficient way to solve a variety of equations, including ordinary differential equations (ODEs), partial differential equations (PDEs), and integral equations [7,13,14]. Due to the discrete-time characterization of wavelet coefficient decay, the sparse form of the coefficients matrices arises. This property is very useful to reduce the computational cost. In this work, we aim to solve the system of the Volterra integral equation using Alpert’s multi-wavelets by exploiting the above-mentioned property. Some results are formally proved and supported by numerical experiments.

The paper is structured as follows. A brief introduction of the Alpert’s multi-wavelets is provided in Section 2. In Section 3, the wavelet Galerkin method is used to approximate the solution of the problem, and the convergence analysis is investigated. Some numerical experiments are performed to illustrate the efficiency and accuracy of the proposed method.

2. Alpert’s Multi-Wavelets and Multiresolution Analysis

Alpert et al. [13,15] introduced a class of multi-wavelets for $L^2$, which are indexed by a parameter $r\ge 0$ and built via Lagrange polynomials of degree less than r. These multi-wavelets are piecewise polynomials that are locally supported and orthonormal. The multiresolution analysis (MRA) framework, introduced and developed by Mallat [16] and Meyer [17], is useful to construct these bases.

According to MRA, a set of primal scaling functions $\{{\varphi }_{0,0}^0,\cdots ,{\varphi }_{0,0}^{r-1}\}$ is introduced for primal subspace $V_0^r\in L^2[0,1]$. By translation and dilation of primal scaling functions $\left\{{\varphi }^k\right\},k=0,\cdots ,r-1$, we determine a space $V_j^r$,

$$V_j^r=Span\{{\varphi }_{j,b}^k:={\mathcal{D}}_{2^j}{\mathcal{T}}_b{\varphi }^k,b\in {\mathcal{B}}_j,k=0,1,\cdots ,r-1\},$$

where ${\mathcal{D}}_a$ and ${\mathcal{T}}_b$ are the dilation and translation operators, respectively such that for a given function h, ${\mathcal{D}}_{a}h\left(x\right)=a^{\frac{1}{2}}h\left(ax\right)$ and ${\mathcal{T}}_{b}h\left(x\right)=h(x-b)$, also ${\mathcal{B}}_j:=\{0,1,\cdots ,2^j-1\}$ for $j\in {\mathbb{Z}}^{+}\cup \left\{0\right\}$. Therefore, ${\varphi }_{j,b}^k$ is a polynomial of degree less than k which is restricted to $I_{j,b}=[x_{j,b},x_{j,b+1}]$ where $x_{j,b}:=2^{-j}b$ and $\mathsf{\Omega }:=[0,1]={\bigcup }_{b\in {\mathcal{B}}_j}{I}_{j,b}$.

For a fixed integer $J\ge 0$, the orthogonal projection ${\mathcal{P}}_J^{r}h$ of $h\in L^2[0,1]$ onto $V_J^r$ is determined by

$$h\approx {\mathcal{P}}_J^r\left(h\right)=\sum _{b\in {\mathcal{B}}_J}\sum _{k=0}^{r-1}\langle h,{\varphi }_{J,b}^k\rangle {\varphi }_{J,b}^k.$$

(3)

The coefficients $h_{J,b}^k=\langle f,{\varphi }_{J,b}^k\rangle$ are determined by $h_{J,b}^k={\int }_{I_{J,b}}h\left(x\right){\varphi }_{J,b}^k\left(x\right)dx$ [18,19,20]. To avoid computing integrals, the Gauss–Legendre Quadrature are applied as follows

$$h_{J,b}^k\approx 2^{-J/2}\sqrt{\frac{{\omega }_k}{2}}h\left(2^{-J}(\frac{{\tau }_k+1}{2}+b)\right),k=0,\dots ,r-1,b\in {\mathcal{B}}_J,$$

(4)

where ${\omega }_k$ and ${\tau }_k$ are the Gauss-Legendre Quadrature weights and the roots of Legendre polynomial of order r, respectively which are introduced in [21,22,23]. The projection ${\mathcal{P}}_J^{r}h$ converges to h if the function $h\in C^r\left(\mathsf{\Omega }\right)$ (r times continuously differentiable) [15]. ${\mathcal{P}}_J^{r}h$ approximates h with mean error bounded as follows

$$\parallel {\mathcal{P}}_J^r\left(h\right)-h\parallel \le 2^{-Jr}\frac{2}{4^{r}r!}\underset{x\in \mathsf{\Omega }}{sup}\left|h^{\left(r\right)}\left(x\right)\right|.$$

(5)

Assume that ${\mathsf{\Phi }}_J^r:={[{\mathsf{\Phi }}_{r,J,0},\cdots ,{\mathsf{\Phi }}_{r,J,2^J-1}]}^T$, where ${\mathsf{\Phi }}_{r,J,b}:=[{\varphi }_{J,b}^0,\cdots ,{\varphi }_{J,b}^{r-1}]$. ${\mathsf{\Phi }}_J^r$ refers to the function vector called a multi-scaling function. In fact, ${\mathsf{\Phi }}_J^r$ is a function vector which includes the scaling function in the space $V_J^r$. By this definition, one can rewrite (3) viz, ${\mathcal{P}}_J^r\left(h\right)=H^T{\mathsf{\Psi }}_J^r$ where H is a vector with entries $H_{br+k+1}:=h_{J,b}^k$ and has dimension $N:=r{2}^J$.

This construction of multi-wavelets for $L^2\left(\mathsf{\Omega }\right)$ can be extended to the two–dimensional space including $L^2{\left(\mathsf{\Omega }\right)}^2$. Let us introduce the space $V_J^{r,2}:=V_J^r\times V_J^r$ which is spanned by orthonormal bases $\{{\varphi }_{J,b}^k{\varphi }_{J,b^{\prime }}^{k^{\prime }}:b,b^{\prime }\in {\mathcal{B}}_J,k,k^{\prime }=0,1,\cdots ,r-1\}.$ Therefore, any function $h\in L^2{\left(\mathsf{\Omega }\right)}^2$ can be projected onto the $V_J^{r,2}$ by the projection ${\mathcal{P}}_J^r$ via,

$$h\approx {\mathcal{P}}_J^{r}h=\sum _{b\in {\mathcal{B}}_j}\sum _{k^{\prime }=0}^{r-1}\sum _{k=0}^{r-1}\sum _{b^{\prime }\in {\mathcal{B}}_j}{H}_{rb+(k+1),r{b}^{\prime }+(k^{\prime }+1)}{\varphi }_{J,b}^k\left(x\right){\varphi }_{J,b^{\prime }}^{k^{\prime }}\left(y\right)={{\mathsf{\Phi }}_J^r}^T\left(x\right)H{\mathsf{\Phi }}_J^r\left(y\right),$$

(6)

where H is an $(N\times N)$ matrix whose elements are obtained by

$$H_{rb+(k+1),r{b}^{\prime }+(k^{\prime }+1)}\approx 2^{-J}\sqrt{\frac{{\omega }_k}{2}}\sqrt{\frac{{\omega }_{k^{\prime }}}{2}}h\left(2^{-J}({\widehat{\tau }}_k+b),2^{-J}({\widehat{\tau }}_{k^{\prime }}+b^{\prime })\right),$$

(7)

where ${\widehat{\tau }}_k=({\tau }_k+1)/2$. By the following theorem, it is possible to bound the error for such projection, if the function h is sufficiently smooth.

Theorem 1

([15]). Suppose that the function $h:{[0,1]}^2\to {\mathbb{R}}^2$ has continuous partial derivatives of order r and mixed partial derivative of order $2r$. Then

$$\parallel {\mathcal{P}}_J^{r}h-h\parallel \le {\mathcal{M}}_{max}\frac{2^{1-rJ}}{4^{r}r!}\left(2+\frac{2^{1-Jr}}{4^{r}r!}\right),$$

(8)

where

$${\mathcal{M}}_{max}=max\left\{\underset{\xi \in [0,1)}{sup}|\frac{{\partial }^r}{\partial x^r}h(\xi ,y)|,\underset{\eta \in [0,1)}{sup}|\frac{{\partial }^r}{\partial y^r}h(x,\eta )|,\underset{{\xi }^{\prime },{\eta }^{\prime }\in [0,1)}{sup}|\frac{{\partial }^{2r}}{\partial x^r\partial y^r}h({\xi }^{\prime },{\eta }^{\prime })|\right\}.$$

As the subspaces $V_j^r$ are nested, there exist complementary orthogonal subspaces $W_j^r$ such that

$$V_{j+1}^r=V_j^r⨁W_j^r,j\in \mathbb{Z}\cup \left\{0\right\},$$

(9)

here and in the following ⨁ denotes orthogonal sums. There is a family of other bases such that the dilations and translations of these bases span the complementary subspaces $W_j^r$, namely,

$$W_j^r=Span\{{\psi }_{j,b}^k:={\mathcal{D}}_{2^j}{\mathcal{T}}_b{\psi }^k,b\in {\mathcal{B}}_j,k=0,1,\cdots ,r-1\}.$$

The functions ${\psi }_{j,b}^k$ are called multi-wavelets. Due to (9), the multi-scale decomposition can be inductively found, $V_J^r=V_0^r\oplus \left({\oplus }_{j=0}^{J-1}{W}_j^r\right)$. This decomposition gives rise to the multi-scale projection operator ${\mathcal{M}}_J^r$ that maps $L^2\left(\mathsf{\Omega }\right)$ onto $V_J^r$ via

$$h\approx {\mathcal{M}}_J^r\left(h\right)=({\mathcal{P}}_0^r+\sum _{j=0}^{J-1}{\mathcal{Q}}_j^r)\left(h\right),$$

(10)

where ${\mathcal{Q}}_j^r$ is the orthonormal projection operator that maps $L^2\left(\mathsf{\Omega }\right)$ onto $W_j^r$. In fact, by using multi-scale projection operator, any function $h\in L^2\left(\mathsf{\Omega }\right)$ can be approximated by multi-wavelets of higher levels $W_j^r$, $j=0,1,\cdots ,J-1$ and the multi-scaling functions of the coarse space $V_0^r$ viz,

$$h\approx {\mathcal{M}}_J^r\left(h\right)=\sum _{k=0}^{r-1}{h}_{0,0}^k{\varphi }_{0,0}^k+\sum _{j=0}^{J-1}\sum _{b\in {\mathcal{B}}_j}\sum _{k=0}^{r-1}{\tilde{h}}_{j,b}^k{\psi }_{j,b}^k,$$

(11)

where

$$h_{0,0}^k:=\langle h,{\varphi }_{0,0}^k\rangle ,{\tilde{h}}_{j,b}^k:=\langle h,{\psi }_{j,b}^k\rangle .$$

(12)

Note that the single-scale coefficients $h_{0,0}^k$ can be determined by (4). However, for evaluating the multi-wavelets coefficients ${\tilde{h}}_{j,b}^k$ of higher levels, in many cases, they have to be calculated numerically. To avoid such numerical computations, the wavelet transform matrix $T_J$ can be applied, as introduced in [14,24]. This matrix is useful to find the multi-wavelets by using the scaling functions

$${\mathsf{\Psi }}_J^r=T_J{\mathsf{\Phi }}_J^r,$$

(13)

where ${\mathsf{\Psi }}_J^r:={[{\mathsf{\Phi }}_{r,0,b},{\mathsf{\Psi }}_{r,0,b},{\mathsf{\Psi }}_{r,1,b},\cdots ,{\mathsf{\Psi }}_{r,J-1,b}]}^T$ and ${\mathsf{\Psi }}_{r,j,b}:=[{\psi }_{j,b}^0,\cdots ,{\psi }_{j,b}^{r-1}]$, $b\in {\mathcal{B}}_j$. Using the vector function ${\mathsf{\Psi }}_J^r$ and (11), we can write

$$h\approx {\mathcal{M}}_J^r\left(h\right)={{\tilde{H}}_J}^T{\mathsf{\Psi }}_J^r,$$

(14)

where ${\tilde{H}}_J$ is an N-dimensional vector with entries $h_{0,0}^k$ and ${\tilde{h}}_{j,b}^k$ for $b\in {\mathcal{B}}_j$, $j=0,\cdots ,J-1$ and $k=0,\cdots ,r-1$. Besides, it is obvious that ${\tilde{H}}_J={T_J^{-1}}^T{H}_J$.

Thresholding

Alpert’s multi-wavelets provide vanishing moments of order $N_{\psi }^k=k+r-1$ for $k=0,1,\cdots ,r-1$, i.e.,

$${\mathcal{N}}_p^k={\int }_{-\infty }^{\infty }{x}^p{\psi }_{0,0}^k\left(x\right)dx,0\le p<N_{\psi }^k,andk=0,1,\cdots ,r-1.$$

(15)

Furthermore, Alpert’s multi-wavelets are uniformly bounded concerning to $L_{\infty }$ and $L_1$, i.e.,

$$\parallel {\psi }_{J,b}^k{\parallel }_{L_{\infty }\left(\mathsf{\Omega }\right)}\lesssim 1,{\parallel {\psi }_{J,b}^k\parallel }_{L_1\left(\mathsf{\Omega }\right)}\lesssim 1.$$

(16)

The vanishing moments and normalization (16) imply that the detail coefficients ${\tilde{h}}_{J,b}^k$ become small when the underlying function is locally smooth. Therefore it is possible to obtain [25]

$${\tilde{h}}_{J,b}^k=|\langle h,{\psi }_{J,b}^k\rangle |\le \underset{P\in {\prod }_{N_{\psi }^k}}{inf}|\langle h-P,{\psi }_{J,b}^k\rangle |\lesssim 2^{-J{N}_{\psi }^k}{\parallel h\parallel }_{W^{1,N_{\psi }^k}\left(\mathsf{\Omega }\right)}.$$

(17)

So the detail coefficients decay at the rate of $2^{-J{N}_{\psi }^k}$ and in the regions where the function is smooth, most of the detail coefficients may be discarded when the refinement level J increases. This gives rise to thresholding. The thresholding operator ${\mathcal{T}}_{D_{\epsilon }}$ is introduced by

$${\mathcal{T}}_{D_{\epsilon }}({\tilde{H}}_J)={\mathcal{H}}_J,$$

(18)

where $D_{\epsilon }:=\{(J,b,k):|{\tilde{h}}_{J,b}^k|>\epsilon \}$ and the elements of ${\mathcal{H}}_J$ are defined by

$${\overline{h}}_{j,b}^k:=\left\{\begin{array}{cc}{\tilde{h}}_{j,b}^k,\hfill & (j,b,k)\in D_{\epsilon },\hfill \\ 0,\hfill & \mathrm{else},\hfill \end{array}\right.b\in {\mathcal{B}}_j,j=0,\cdots ,J-1,k=0,\cdots ,r-1.$$

(19)

Note that the thresholding affects only the detail coefficients while the coarse scale coefficients remain unchanged.

The approximation error due to the thresholding can be estimated similarly to the classical wavelets. Let ${\mathcal{A}}_{D_{\epsilon }}$ be the approximation operator ${\mathcal{A}}_{D_{\epsilon }}:={{\mathcal{M}}_J^r}^{-1}{\mathcal{T}}_{D_{\epsilon }}{\mathcal{M}}_J^r$. The approximation error due to the thresholding can be bounded as stated by the following proposition [25].

Proposition 1.

(Approximation error). Let Ω be bounded and ${\epsilon }_j={\overline{a}}^{j-J}\epsilon$ with $\overline{a}>1$. Then the approximation error concerning to the set of significant details $D_{\epsilon }$ is uniformly bounded concerning to $L^q\left(\mathsf{\Omega }\right)$, $q\in [1,\infty ]$, i.e.,

$$\parallel {\mathcal{P}}_J^{r}h-{\mathcal{P}}_{J,D_{\epsilon }}^r{h\parallel }_{L^q\left(\mathsf{\Omega }\right)}\le C_{thr}\epsilon ,$$

(20)

for some constant $C_{thr}>0$ independent of $J,\epsilon$. Here ${\mathcal{P}}_J^{r}h$ and ${\mathcal{P}}_{J,D_{\epsilon }}^{r}h$ are the projections according to (11) corresponding to the coefficients $\widetilde{H_J}$ and ${\mathcal{A}}_{D_{\epsilon }}\widetilde{H_J}$.

3. Multi-Wavelets Galerkin Method

In this section, we use the wavelet Galerkin method to solve the system of the Volterra integral Equation (1). To this end, we will apply the interpolation property of scaling functions to reach an efficient algorithm. Assume that the solution $\mathbf{u}\left(x\right)$ of Equation (1) can be expanded using multi-scale projection operator ${\mathcal{M}}_J^r$ based on multi-wavelets as follows

$$\mathbf{u}\left(x\right)\approx {\mathcal{M}}_J^r\left(\mathbf{u}\right)\left(x\right)=({\mathcal{P}}_0^r+\sum _{j=0}^{J-1}{\mathcal{Q}}_j^r)\left(\mathbf{u}\right)\left(x\right)={\tilde{\mathbf{U}}}^T\otimes {\mathsf{\Psi }}_J^r\left(x\right),$$

(21)

where ⊗ is the Kronecker product and $\tilde{\mathbf{U}}={({\tilde{U}}_1^T,\cdots ,{\tilde{U}}_n^T)}^T$ is a $(1\times nr{2}^J)$ vector whose elements are n unknown sub-vectors ${\tilde{U}}_i$ with a dimension of $(r{2}^J\times 1)$ such that

$$u_i\left(x\right):={\tilde{U}}_i^T{\mathsf{\Psi }}_J^r\left(x\right),i=1,2,\cdots ,n.$$

One can imagine two types of equations, linear and nonlinear. For the linear type, the i-th component of the vector function $\mathbf{g}(x,t,\mathbf{u}(t\left)\right)$ has the form

$${\mathbf{g}}_j(x,t,\mathbf{u}\left(t\right)):=\sum _{i=1}^n{a}_{ji}(x,t)u_i\left(t\right),j=1,2,\cdots ,n,$$

(22)

and it can be approximated by multi-scale operator, i.e.,

$$\begin{array}{cc}\hfill \sum _{i=1}^n{a}_{ji}(x,t)u_i\left(t\right)& \approx \sum _{i=1}^n{\mathcal{P}}_J^r\left(a_{ji}{u}_i\right)(x,t)=\sum _{i=1}^n{{\mathsf{\Phi }}_J^r}^T\left(x\right)A_{ji}^T{\mathsf{\Phi }}_J^r\left(t\right)\hfill \\ \hfill & \sum _{i=1}^n{{\mathsf{\Psi }}_J^r}^T\left(x\right)T_J{A}_{ji}^T{T}_J^{-1}{\mathsf{\Psi }}_J^r\left(t\right),j=1,2,\cdots ,n,\hfill \end{array}$$

where $A_{ji}$$(j,i=1:n)$ are $(r{2}^J\times r{2}^J)$ matrices. Integrating from 0 to x, we get

$$\begin{array}{cc}\hfill {\int }_0^x{\mathbf{g}}_j(x,t,\mathbf{u}\left(t\right))dt& =\sum _{i=1}^n\underset{p_j\left(x\right)}{\underbrace{{{\mathsf{\Psi }}_J^r}^T\left(x\right)T_J{A}_{ji}^T{T}_J^{-1}{I}_{\psi }{\mathsf{\Psi }}_J^r\left(x\right)}}\hfill \\ \hfill & =\sum _{i=1}^n{P}_i^T{T}_J^{-1}{\mathsf{\Psi }}_J^r\left(x\right)=\sum _{i=1}^n{\tilde{U}}_i^T{\mathcal{A}}_{ji}{T}_J^{-1}{\mathsf{\Psi }}_J^r\left(t\right),\hfill \end{array}$$

(23)

where ${\mathcal{A}}_{ji}$$(j,i=1:n)$ are $(r{2}^J\times r{2}^J)$ matrices and the rest are $(r{2}^J\times 1)$ vectors. But if the j-th component of vector function $\mathbf{g}(x,t,\mathbf{u}(t\left)\right)$ is nonlinear, one can consider the following expansion

$$\begin{array}{cc}\hfill {\int }_0^x{g}_j(x,t,\mathbf{u}\left(t\right))dt& \approx {\mathcal{P}}_J^r\left({\int }_0^x{g}_j\left(x,t,{\mathcal{M}}_J^r\left(\mathbf{u}\right)\left(t\right)\right)dt\right)\hfill \\ \hfill & G_j^T{\mathsf{\Phi }}_J^r\left(x\right)=G_j^T{T}_J^{-1}{\mathsf{\Psi }}_J^r\left(x\right),j=1,2,\cdots ,n,\hfill \end{array}$$

(24)

where $G_j$ is a $(r{2}^J\times 1)$ vector whose elements are nonlinear algebraic equations. In view of the Equations (23) and (24), and using operational matrix $I_{\psi }$ of integration for multi-wavelets introduced in [7,14,15], one can write

$${\int }_0^x\mathbf{g}(x,t,\mathbf{u}\left(t\right))dt\approx \left\{\begin{array}{cc}{\tilde{\mathbf{U}}}^T{\Gamma }^T\otimes T_J^{-1}{\mathsf{\Psi }}_J^r\left(x\right),\hfill & \mathrm{linear},\hfill \\ {\mathbf{G}}^T\otimes T_J^{-1}{\mathsf{\Psi }}_J^r\left(x\right),\hfill & \mathrm{nonlinear},\hfill \end{array}\right.$$

(25)

with $\Gamma :={\left(\mathcal{A}\right)}_{ji}$, $(j,i=1:n)$ and $\mathbf{G}:={(G_1^T,G_2^T,\cdots ,G_n^T)}^T$.

Such an approximation can be considered for the j-th element of $\mathbf{f}$ viz,

$$f_j\left(x\right)\approx {\mathcal{P}}_J^r\left(f_j\right)\left(x\right)=F_j^T{\mathsf{\Phi }}_J^r\left(x\right)=F_j^T{T}_J^{-1}{\mathsf{\Psi }}_J^r\left(x\right),$$

and thus by putting $\mathbf{F}:={(F_1^T,F_2^T,\cdots ,F_n^T)}^T$ we have

$$\mathbf{f}\approx {\mathbf{F}}^T\otimes T_J^{-1}{\mathsf{\Psi }}_J^r\left(x\right).$$

(26)

Now, we introduce the residual as

$$r_J^r\left(x\right)={\tilde{\mathbf{U}}}^T\otimes {\mathsf{\Psi }}_J^r\left(x\right)-{\mathbf{F}}^T\otimes T_J^{-1}{\mathsf{\Psi }}_J^r\left(x\right)-\left\{\begin{array}{cc}{\tilde{\mathbf{U}}}^T{\Gamma }^T\otimes T_J^{-1}{\mathsf{\Psi }}_J^r\left(x\right),\hfill & \mathrm{linear},\hfill \\ {\mathbf{G}}^T\otimes T_J^{-1}{\mathsf{\Psi }}_J^r\left(x\right),\hfill & \mathrm{nonlinear}.\hfill \end{array}\right.$$

(27)

To apply the Galerkin method, it is necessary that $\langle r_J^r,{\mathsf{\Psi }}_J^r\rangle =0$. Thus we have

$${\tilde{\mathbf{U}}}^T-{\mathbf{F}}^T\otimes T_J^{-1}-\left\{\begin{array}{c}{\tilde{\mathbf{U}}}^T{\Gamma }^T\otimes T_J^{-1}\hfill \\ {\mathbf{G}}^T\otimes T_J^{-1}\hfill \end{array}\right.=0,\begin{array}{c}\mathrm{linear},\hfill \\ \mathrm{nonlinear}.\hfill \end{array}$$

(28)

By solving this system of linear and nonlinear equations using restarted generalized minimal residual method (GMRES) and Newton methods, respectively, we obtain the approximate solution of the Equation (1). Note that because we use the Galerkin method with orthogonal bases, such a system will have a unique solution [26].

Convergence Analysis

Theorem 2.

Suppose that${\mathbf{e}}_J=\mathbf{u}-{\mathcal{M}}_J^r\left(\mathbf{u}\right)$where$\mathbf{u}$and${\mathcal{M}}_J^r\left(\mathbf{u}\right)$are the exact and approximate solutions of nonlinear system (1), respectively. Let X be an open set in$\mathbb{R}$and let$g:\mathsf{\Omega }\times X\to \mathbb{R}$be a function such that$g(x,t,\mathbf{u}\left(x\right))\in C^r\left(\mathsf{\Omega }\right)$for any$\mathbf{u}\in X$and the condition (2) is satisfied. Furthermore, presume that$\mathbf{f}\in C^r\left(\mathsf{\Omega }\right)$. Furthermore, assume that the residual$\mathbf{e}:=\tilde{r}-\left(r_J^r\right)$where the residual$\tilde{r}$is specified as

$$\tilde{r}\left(x\right)=u\left(x\right)-f\left(x\right)-{\int }_0^{x}g(x,t,u\left(t\right))dt,$$

(29)

and$r_J^r$is introduced in (27).

If$\mathbf{u}\in C^r\left(\mathsf{\Omega }\right)$, and the method used to solve system (28) is convergent then one has

$${\parallel \mathbf{e}\parallel }_2^2\le \frac{2^{1-Jr}}{4^{r}r!}\left(|1+\sqrt{nr{2}^J}\kappa |\underset{x\in \mathsf{\Omega }}{sup}\left|\mathbf{u}\left(x\right)\right|+\underset{x\in \mathsf{\Omega }}{sup}\left|\mathbf{f}\left(x\right)\right|\right),$$

where κ is a positive constant. Consequently,$\mathbf{e}\to 0$when$J\to \infty$.

Proof. 

Using (27), (29) and the hypotheses of the theorem, we can write

$$\mathbf{e}\left(x\right)={\mathbf{e}}_J\left(x\right)-\left(\mathbf{f}\left(x\right)-{\mathcal{M}}_J^r\left(\mathbf{f}\right)\left(x\right)\right)-\left({\int }_0^x\mathbf{g}(x,t,\mathbf{u}\left(t\right))dt-{\mathcal{M}}_J^r({\int }_0^x\mathbf{g}(x,t,{\mathcal{M}}_J^r\left(\mathbf{u}\right))\left(t\right)dt)\right),$$

(30)

where $\mathbf{e}:=\tilde{r}-\left(r_J^r\right)$. Since the function $\mathbf{g}$ satisfies the Lipschitz condition (2), Equation (30) can be reduced to

$$\mathbf{e}\left(x\right)={\mathbf{e}}_J\left(x\right)-\left(\mathbf{f}\left(x\right)-{\mathcal{M}}_J^r\left(\mathbf{f}\right)\left(x\right)\right)-A{\int }_0^x{\mathbf{e}}_J\left(t\right)dt.$$

Now, suppose that

$$\mathbf{e}\approx \mathcal{E}\otimes {\mathsf{\Psi }}_J^r,{\mathbf{e}}_J\approx {\mathcal{E}}_J\otimes {\mathsf{\Psi }}_J^r,$$

where $\mathcal{E}$ and ${\mathcal{E}}_J$ are the $(1\times nr{2}^J)$ vectors and thus, one can write

$$\mathcal{E}\otimes {\mathsf{\Psi }}_J^r={\mathcal{E}}_J\otimes {\mathsf{\Psi }}_J^r-\left(\mathbf{f}\left(x\right)-{\mathcal{M}}_J^r\left(\mathbf{f}\right)\left(x\right)\right)-A{\mathcal{E}}_J\otimes I_{\psi }{\mathsf{\Psi }}_J^r,$$

Taking $L_2$-norm from both sides and using the triangle inequality yields

$$\begin{array}{cc}\hfill \parallel \mathcal{E}\otimes {\mathsf{\Psi }}_J^r{\parallel }_2^2& \le \parallel {\mathcal{E}}_J\otimes {\mathsf{\Psi }}_J^r{\parallel }_2^2+\parallel \mathbf{f}\left(x\right)-{\mathcal{M}}_J^r{\left(\mathbf{f}\right)\left(x\right)\parallel }_2^2+A{\parallel {\mathcal{E}}_J\otimes I_{\psi }{\mathsf{\Psi }}_J^r\parallel }_2^2\hfill \\ \hfill & =\parallel {\mathcal{E}}_J{\parallel }_2^2\parallel {\mathsf{\Psi }}_J^r{\parallel }_2^2+A\parallel {\mathcal{E}}_J{\parallel }_2^2\parallel I_{\psi }{\mathsf{\Psi }}_J^r{\parallel }_2^2+{\parallel \mathbf{f}\left(x\right)-{\mathcal{M}}_J^r\left(\mathbf{f}\right)\left(x\right)\parallel }_2^2\hfill \\ \hfill & \le \parallel {\mathcal{E}}_J{\parallel }_2^2\parallel {\mathsf{\Psi }}_J^r{\parallel }_2^2+A\parallel {\mathcal{E}}_J{\parallel }_2^2\parallel I_{\psi }{\parallel }_2^2\parallel {\mathsf{\Psi }}_J^r{\parallel }_2^2+{\parallel \mathbf{f}\left(x\right)-{\mathcal{M}}_J^r\left(\mathbf{f}\right)\left(x\right)\parallel }_2^2\hfill \end{array}$$

where the second row comes from theorem 8 in [27]. Since Alpert multi-wavelets are orthonormal, one can write

$$\begin{array}{cc}\hfill {\parallel \mathcal{E}\parallel }_2^2& \le \parallel {\mathcal{E}}_J{\parallel }_2^2+A\parallel {\mathcal{E}}_J{\parallel }_2^2\parallel I_{\psi }{\parallel }_2^2+{\parallel \mathbf{f}\left(x\right)-{\mathcal{M}}_J^r\left(\mathbf{f}\right)\left(x\right)\parallel }_2^2\hfill \\ \hfill & =\parallel {\mathcal{E}}_J{\parallel }_2^2\left(\parallel I_{nr{2}^J}{\parallel }_2^2+A{\parallel I_{\psi }\parallel }_2^2\right)+{\parallel \mathbf{f}\left(x\right)-{\mathcal{M}}_J^r\left(\mathbf{f}\right)\left(x\right)\parallel }_2^2.\hfill \end{array}$$

According to the previous section, when ${\mathsf{\Psi }}_J^r$ has high vanishing moments and the function h is smooth, $\langle h,{\mathsf{\Psi }}_J^r\rangle$ decays fast in $J\to \infty$. By means of vanishing moments of Alpert’s multi-wavelets and the matrix norms inequalities, we get

$${\parallel \mathcal{E}\parallel }_2^2\le \parallel {\mathcal{E}}_J{\parallel }_2^2|1+\sqrt{nr{2}^J}\kappa |+\parallel \mathbf{f}\left(x\right)-{\mathcal{M}}_J^r{\left(\mathbf{f}\right)\left(x\right)\parallel }_2^2,$$

where $\kappa =A\parallel I_{\psi }{\parallel }_{\infty }^2$. Now Equation (5) leads to the desired result

$${\parallel \mathcal{E}\parallel }_2^2\le \frac{2^{1-Jr}}{4^{r}r!}\left(|1+\sqrt{nr{2}^J}\kappa |\underset{x\in \mathsf{\Omega }}{sup}\left|\mathbf{u}\left(x\right)\right|+\underset{x\in \mathsf{\Omega }}{sup}\left|\mathbf{f}\left(x\right)\right|\right).$$

 ☐

Theorem 3.

Let the condition of Theorem 2 be valid. Suppose that${\mathbf{u}}_J$is the approximate solution obtained from solving (28) using restarted GMRES or Newton methods. If these methods solve (28) with proper accuracy, the error can be estimated from

$$\parallel \mathbf{u}-{\mathbf{u}}_J{\parallel }_2^2\le (1-A\parallel I_{\psi }{\parallel }_2^2{)}^{-1}\frac{2^{1-Jr}}{4^{r}r!}\underset{x\in \mathsf{\Omega }}{sup}\left|\mathbf{f}\left(x\right)\right|+\eta ,$$

(31)

where η is a small positive number that desire to zero.

Proof. 

Taking ${\mathcal{M}}_J^r\left(\mathbf{u}\right)$ as the approximate solution of (1) and ${\mathbf{u}}_J$ as the approximate solution obtained from solving (28) using restarted GMRES or Newton methods, the convergence can be concluded from

$$\parallel \mathbf{u}-{\mathbf{u}}_J\parallel \le \parallel \mathbf{u}-{\mathcal{M}}_J^r\left(\mathbf{u}\right)\parallel +\parallel {\mathcal{M}}_J^r\left(\mathbf{u}\right)-{\mathbf{u}}_J\parallel .$$

(32)

The approximate solution of (1) satisfies

$${\mathcal{M}}_J^r\left(\mathbf{u}\right)\left(x\right)={\mathcal{M}}_J^r\left(\mathbf{f}\right)\left(x\right)+{\mathcal{M}}_J^r\left({\int }_0^x\mathbf{g}(x,t,{\mathcal{M}}_J^r\left(\mathbf{u}\right)\left(t\right)dt\right).$$

(33)

Subtracting (33) from (1), and using the Lipschits condition (2), one can write

$${\mathbf{e}}_J\le \mathbf{f}-{\mathcal{M}}_J^r\left(\mathbf{f}\right)+A{\int }_0^x{\mathbf{e}}_{J}dt.$$

(34)

Let us consider ${\mathbf{e}}_J\approx {\mathcal{E}}_J\otimes {\mathsf{\Psi }}_J^r$ where ${\mathcal{E}}_J$ is the $(1\times nr{2}^J)$ vector and thus, we have

$${\mathcal{E}}_J\otimes {\mathsf{\Psi }}_J^r=\left(\mathbf{f}\left(x\right)-{\mathcal{M}}_J^r\left(\mathbf{f}\right)\left(x\right)\right)+A{\mathcal{E}}_J\otimes I_{\psi }{\mathsf{\Psi }}_J^r.$$

(35)

Taking $L_2$-norm from both sides of (35) and using Theorem 8 in [27] yields

$$\begin{array}{cc}\hfill \parallel {\mathcal{E}}_J\otimes {\mathsf{\Psi }}_J^r{\parallel }_2^2& \le \parallel \mathbf{f}\left(x\right)-{\mathcal{M}}_J^r{\left(\mathbf{f}\right)\left(x\right)\parallel }_2^2+A{\parallel {\mathcal{E}}_J\otimes I_{\psi }{\mathsf{\Psi }}_J^r\parallel }_2^2\hfill \\ \hfill & \le \parallel \mathbf{f}\left(x\right)-{\mathcal{M}}_J^r{\left(\mathbf{f}\right)\left(x\right)\parallel }_2^2+A\parallel {\mathcal{E}}_J{\parallel }_2^2{\parallel I_{\psi }{\mathsf{\Psi }}_J^r\parallel }_2^2\hfill \end{array}$$

Since Alpert multi-wavelets are orthonormal, one can write

$$\parallel {\mathcal{E}}_J{\parallel }_2^2\le \parallel \mathbf{f}\left(x\right)-{\mathcal{M}}_J^r{\left(\mathbf{f}\right)\left(x\right)\parallel }_2^2+A\parallel {\mathcal{E}}_J{\parallel }_2^2{\parallel I_{\psi }\parallel }_2^2.$$

(36)

Therefore one can bound the error of $\parallel \mathbf{u}-{\mathcal{M}}_J^r\left(\mathbf{u}\right)\parallel$ via

$$\parallel {\mathcal{E}}_J{\parallel }_2^2\le (1-A\parallel I_{\psi }{\parallel }_2^2{)}^{-1}{\parallel \mathbf{f}\left(x\right)-{\mathcal{M}}_J^r\left(\mathbf{f}\right)\left(x\right)\parallel }_2^2.$$

(37)

According to the theorem hypotheses, the methods used to solve the obtained system are convergent. So $\eta :=\parallel {\mathcal{M}}_J^r\left(\mathbf{u}\right)-{\mathbf{u}}_J\parallel$ will be very small. Inequality (31) is obtained using (5) and (32), i.e.,

$$\parallel \mathbf{u}-{\mathbf{u}}_J{\parallel }_2^2\le (1-A\parallel I_{\psi }{\parallel }_2^2{)}^{-1}\frac{2^{1-Jr}}{4^{r}r!}\underset{x\in \mathsf{\Omega }}{sup}\left|\mathbf{f}\left(x\right)\right|+\eta .$$

 ☐

4. Numerical Examples

In this section, some numerical experiments are considered to illustrate the convergence and efficiency of the proposed method. To this end, we report the $L^2$ errors of the solution which is defined by

$${\xi }_u:={\parallel u-{\mathcal{M}}_J^r\left(u\right)\parallel }_2={\left({\int }_{\mathsf{\Omega }}{|u\left(x\right)-{\mathcal{M}}_J^r\left(u\right)\left(x\right)|}^{2}dx\right)}^{1/2},$$

where u and ${\mathcal{M}}_J^r\left(u\right)$ are the exact and approximate solution of systems (1), respectively. In order to get the sparse coefficients matrix in the linear type, all the entries of this matrix that are less than a small positive number $\epsilon$ are set to zero. Finally, one can find the rate of sparsity $S_{\epsilon }$ which is defined by [28]

$$S_{\epsilon }=\frac{N_0-N_{\epsilon }}{N_0}\times 100\%.$$

where $N_0$ is the total number of elements and $N_{\epsilon }$ the number of elements remaining after thresholding.

Example 1.

Let us run the proposed method on the following linear Volterra integral equation [5,29]

$$\begin{array}{cc}\hfill u_1\left(x\right)& =-{\int }_0^x{e}^{-(s-x)}{u}_1\left(s\right)ds-{\int }_0^{x}cos(s-x)u_2\left(s\right)ds+cosh\left(x\right)+xsin\left(x\right),\hfill \\ \hfill u_2\left(x\right)& =-{\int }_0^x{e}^{s+x}{u}_1\left(s\right)ds-{\int }_0^{x}xcos\left(s\right)u_2\left(s\right)ds+2sin\left(x\right)+x({sin}^2\left(x\right)+e^x).\hfill \end{array}$$

The exact solution is$\mathbf{u}\left(x\right)=(e^{-x},2sin\left(x\right))$. The effect of thresholding on$L_2$-errors and sparsity percentage is reported in

Table 1. Effects of parameters r, J and $\epsilon$ on sparsity and $L^2$-error for Example 1.

		Without Thresholding		$\mathit{\epsilon }={10}^{-5}$		$\mathit{\epsilon }={10}^{-3}$
$\mathit{r}$	$\mathit{J}$	${\mathit{S}}_{\mathbf{\epsilon }}$	${\mathit{L}}^{\mathbf{2}}$-Error	${\mathit{S}}_{\mathbf{\epsilon }}$	${\mathit{L}}^{\mathbf{2}}$-Error	${\mathit{S}}_{\mathbf{\epsilon }}$	${\mathit{L}}^{\mathbf{2}}$-Error
3	2	6.25	${\xi }_{u_1}=3.24\times {10}^{-5}$	$24.48$	${\xi }_{u_1}=3.30\times {10}^{-5}$	$49.83$	${\xi }_{u_1}=2.42\times {10}^{-4}$
			${\xi }_{u_2}=8.41\times {10}^{-5}$		${\xi }_{u_2}=8.41\times {10}^{-5}$		${\xi }_{u_2}=3.04\times {10}^{-4}$
	3	$17.19$	${\xi }_{u_1}=4.05\times {10}^{-6}$	$53.30$	${\xi }_{u_1}=6.02\times {10}^{-6}$	$74.13$	${\xi }_{u_1}=2.41\times {10}^{-4}$
			${\xi }_{u_2}=1.05\times {10}^{-5}$		${\xi }_{u_2}=1.26\times {10}^{-5}$		${\xi }_{u_2}=5.15\times {10}^{-4}$
5	2	$6.25$	${\xi }_{u_1}=6.38\times {10}^{-9}$	$45.69$	${\xi }_{u_1}=4.62\times {10}^{-6}$	$70.69$	${\xi }_{u_1}=1.53\times {10}^{-4}$
			${\xi }_{u_2}=1.66\times {10}^{-8}$		${\xi }_{u_2}=1.28\times {10}^{-6}$		${\xi }_{u_2}=6.51\times {10}^{-4}$
	3	$17.19$	${\xi }_{u_1}=2.00\times {10}^{-10}$	$71.42$	${\xi }_{u_1}=1.24\times {10}^{-5}$	$86.75$	${\xi }_{u_1}=1.53\times {10}^{-4}$
			${\xi }_{u_2}=5.19\times {10}^{-10}$		${\xi }_{u_2}=8.54\times {10}^{-6}$		${\xi }_{u_2}=6.52\times {10}^{-4}$

for different values of r, J ad ε. To illustrate the effect of the refinement level J and the multiplicity parameter r on$L^2$error, Figure 1 is plotted. Figure 2 shows sparse matrix when$r=5$,$J=3$and$\epsilon ={10}^{-4},{10}^{-2}$.

Example 2.

Let us consider the following system of Volterra integral equation as a further example

$$\begin{array}{cc}\hfill u_1\left(x\right)& =-\frac{x^5}{3}-\frac{x^4}{4}+\frac{x^3}{3}+x+{\int }_0^x(x^2-s)(u_1\left(s\right)+u_2\left(s\right))ds\hfill \\ \hfill u_2\left(x\right)& =\frac{x^3}{2}-\frac{x^4}{3}+x^2-{\int }_0^{x}x(u_1\left(s\right)-u_2\left(s\right))ds.\hfill \end{array}$$

The solution is reported in [5,30] and is $\mathbf{u}=(x,x^2)$. To illustrate the effect of thresholding on the coefficients matrix obtained from proposed method, the graph in Figure 3 is provided. Figure 4 shows the effect of parameters r and J.

Example 3.

Let us consider the following nonlinear system

$$\begin{array}{cc}\hfill u_1\left(x\right)& =sin\left(x\right)-x+{\int }_0^x(u_1^2\left(s\right)+u_2^2\left(s\right))ds\hfill \\ \hfill u_2\left(x\right)& =cos\left(x\right)-\frac{1}{2}{sin}^2\left(x\right)+{\int }_0^x{u}_1\left(s\right)u_2\left(s\right)ds.\hfill \end{array}$$

The exact solutions of this system are$u_1=sin\left(x\right)$and$u_2=cos\left(x\right)$[8].

Table 2. Effect of the refinement level J and multiplicity parameter r on $L_2$-error for Example 3.

	r	$\mathit{J}=2$	$\mathit{J}=3$	$\mathit{J}=4$	$\mathit{J}=5$
$u_1$	3	$4.20\times {10}^{-5}$	$5.25\times {10}^{-6}$	$6.56\times {10}^{-7}$	$8.20\times {10}^{-8}$
	4	$4.04\times {10}^{-7}$	$2.53\times {10}^{-8}$	$1.58\times {10}^{-9}$	$9.86\times {10}^{-11}$
	5	$8.31\times {10}^{-9}$	$2.60\times {10}^{-10}$	$8.12\times {10}^{-12}$	$2.54\times {10}^{-13}$
$u_2$	3	$2.50\times {10}^{-5}$	$3.21\times {10}^{-6}$	$4.01\times {10}^{-7}$	$5.02\times {10}^{-8}$
	4	$6.62\times {10}^{-7}$	$4.14\times {10}^{-8}$	$2.59\times {10}^{-9}$	$1.62\times {10}^{-10}$
	5	$5.06\times {10}^{-9}$	$1.58\times {10}^{-10}$	$4.94\times {10}^{-12}$	$1.54\times {10}^{-13}$

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 L₂-errors decrease. In

Table 3. The comparison between absolute errors of Example 3 for $u_1$.

	HPM [31]	Method Based upon Discretization [32]	Simpson’s 3/8 Rule [10]	Bernstein Collocation Method [33]	Present Method $(\mathit{r}=7,\mathit{J}=2)$
x	$(\mathit{n}=\mathbf{5})$	$(\mathit{h}=\mathbf{200})$	$(\mathit{h}=\mathbf{0.025})$	$(\mathit{n}=\mathbf{10})$
0.1	$1.4\times {10}^{-07}$	$2.4\times {10}^{-05}$	$3.0\times {10}^{-10}$	$5.5\times {10}^{-10}$	$1.0\times {10}^{-12}$
0.2	$3.5\times {10}^{-06}$	$9.8\times {10}^{-05}$	$1.1\times {10}^{-09}$	$1.7\times {10}^{-10}$	$1.1\times {10}^{-12}$
0.3	$5.5\times {10}^{-05}$	$2.3\times {10}^{-04}$	$3.6\times {10}^{-09}$	$2.6\times {10}^{-10}$	$1.1\times {10}^{-12}$
0.4	$3.8\times {10}^{-04}$	$4.1\times {10}^{-04}$	$6.0\times {10}^{-09}$	$1.0\times {10}^{-10}$	$9.6\times {10}^{-13}$
0.5	$1.6\times {10}^{-03}$	$6.6\times {10}^{-04}$	$8.7\times {10}^{-09}$	$1.1\times {10}^{-10}$	$3.2\times {10}^{-12}$
0.6		$9.7\times {10}^{-04}$	$1.4\times {10}^{-08}$	$2.5\times {10}^{-10}$	$8.4\times {10}^{-13}$
0.7		$1.4\times {10}^{-03}$	$1.9\times {10}^{-08}$	$5.8\times {10}^{-10}$	$9.2\times {10}^{-13}$
0.8		$1.8\times {10}^{-03}$	$2.4\times {10}^{-08}$	$1.0\times {10}^{-08}$	$9.2\times {10}^{-13}$
0.9		$2.5\times {10}^{-03}$	$3.3\times {10}^{-08}$	$1.0\times {10}^{-07}$	$5.0\times {10}^{-12}$
1.0		$3.2\times {10}^{-03}$	$4.0\times {10}^{-08}$	$8.2\times {10}^{-07}$	$6.7\times {10}^{-11}$

and

Table 4. The comparison between absolute errors of Example 3 for $u_2$.

	HPM [31]	Method Based upon Discretization [32]	Simpson’s 3/8 Rule [10]	Bernstein Collocation Method [33]	Present Method $(\mathit{r}=7,\mathit{J}=2)$
x	$(\mathit{n}=\mathbf{5})$	$(\mathit{h}=\mathbf{200})$	$(\mathit{h}=\mathbf{0.025})$	$(\mathit{n}=\mathbf{10})$
0.1	$3.2\times {10}^{-07}$	$2.5\times {10}^{-04}$	$5.3\times {10}^{-10}$	$6.5\times {10}^{-11}$	$1.2\times {10}^{-13}$
0.2	$1.1\times {10}^{-05}$	$5.0\times {10}^{-04}$	$3.0\times {10}^{-10}$	$1.3\times {10}^{-10}$	$1.6\times {10}^{-13}$
0.3	$1.2\times {10}^{-04}$	$7.5\times {10}^{-04}$	$6.0\times {10}^{-10}$	$6.5\times {10}^{-12}$	$4.0\times {10}^{-13}$
0.4	$6.3\times {10}^{-04}$	$1.0\times {10}^{-04}$	$2.1\times {10}^{-09}$	$6.5\times {10}^{-11}$	$3.9\times {10}^{-13}$
0.5	$2.2\times {10}^{-03}$	$1.2\times {10}^{-03}$	$3.1\times {10}^{-09}$	$9.0\times {10}^{-11}$	$1.3\times {10}^{-12}$
0.6		$1.5\times {10}^{-03}$	$5.3\times {10}^{-09}$	$4.2\times {10}^{-11}$	$6.0\times {10}^{-13}$
0.7		$1.8\times {10}^{-03}$	$9.4\times {10}^{-09}$	$4.2\times {10}^{-11}$	$6.8\times {10}^{-13}$
0.8		$2.0\times {10}^{-03}$	$1.4\times {10}^{-08}$	$3.6\times {10}^{-09}$	$8.1\times {10}^{-13}$
0.9		$2.3\times {10}^{-03}$	$2.0\times {10}^{-08}$	$3.5\times {10}^{-08}$	$4.8\times {10}^{-13}$
1.0		$2.6\times {10}^{-03}$	$2.9\times {10}^{-08}$	$2.9\times {10}^{-07}$	$4.2\times {10}^{-11}$

, results are compared with other methods [10,31,32,33] in terms of the absolute errors. In this paper, J and r are the criteria for the discretization and the degree of polynomials used as a basis, respectively. Taking r = 7 and J = 2, the results of Table 3 and Table 4 indicate that the proposed method solves this equation better than others [10,31,32,33]. Furthermore, we reported the exact and numerical solution by Figure 5.

5. Conclusions

In this paper, we proposed the multi-wavelets Galerkin method to solve the linear and nonlinear Volterra integral equation. The convergence analysis and numerical simulations indicate that the proposed method gives a satisfactory approximation to the exact solution. Thresholding can be used to increase sparsity for a lower computational cost, without affecting the error in $L^2$. Using the interpolation property of Alpert’s multi-wavelets, the proposed method turns out to be fast and very competitive against state-of-the-art techniques. The main advantages of the proposed method are the lower computational cost and lower complexity.