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

: In this work, an efﬁcient 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 ( 2 − Jr ) , where J is the reﬁnement level and r is the multiplicity of multi-wavelets. Several numerical tests are provided to illustrate the ability and efﬁciency of the method.


Introduction
In this paper, we study and construct a novel numerical algorithm for the system of Volterra integral equations of the second kind u(x) = f(x) + x 0 g(x, t, u(t))dt, x ∈ Ω := [0, 1], (1) where f : Ω → R n (n ∈ N) is a given real-valued continuous function, u : Ω → R n is the unknown function that will be determined and the function g : S → R n with S = {(x, t) : x, t ∈ Ω} is a given linear or nonlinear function of u which satisfies the following Lipschitz condition with respect to the third variable: for all x, t ∈ [0, 1] and for all u 1 , u 2 ∈ R n , |g(x, t, u 1 (t)) − g(x, t, u 2 (t))| ≤ A|u 1 − u 2 |. (2) 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.

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 ≥ 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 {φ 0 0,0 , . . . , φ r−1 0,0 } is introduced for primal subspace V r 0 ∈ L 2 [0, 1]. By translation and dilation of primal scaling functions {φ k }, k = 0, . . . , r − 1, we determine a space V r j , where D a and T b are the dilation and translation operators, respectively such that for a given function h, D a h(x) = a 1 2 h(ax) and T b h(x) = h(x − b), also B j := {0, 1, . . . , 2 j − 1} for j ∈ Z + ∪ {0}. Therefore, φ k j,b is a polynomial of degree less than k which is restricted to For a fixed integer J ≥ 0, the orthogonal projection P r The coefficients h k [18][19][20]. To avoid computing integrals, the Gauss-Legendre Quadrature are applied as follows where ω k and τ 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 P r J h converges to h if the function h ∈ C r (Ω) (r times continuously differentiable) [15]. P r J h approximates h with mean error bounded as follows whereτ k = (τ 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 → R 2 has continuous partial derivatives of order r and mixed partial derivative of order 2r. Then where As the subspaces V r j are nested, there exist complementary orthogonal subspaces W r j such that 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 r j , namely, The functions ψ k j,b are called multi-wavelets. Due to (9), the multi-scale decomposition can be inductively found, V r J = V r 0 ⊕ (⊕ J−1 j=0 W r j ). This decomposition gives rise to the multi-scale projection operator M r J that maps where Q r j is the orthonormal projection operator that maps L 2 (Ω) onto W r j . In fact, by using multi-scale projection operator, any function h ∈ L 2 (Ω) can be approximated by multi-wavelets of higher levels W r j , j = 0, 1, . . . , J − 1 and the multi-scaling functions of the coarse space V r 0 viz, where Note that the single-scale coefficients h k 0,0 can be determined by (4). However, for evaluating the multi-wavelets coefficientsh k j,b 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 where Using the vector function Ψ r J and (11), we can write

Thresholding
Alpert's multi-wavelets provide vanishing moments of order N k Furthermore, Alpert's multi-wavelets are uniformly bounded concerning to L ∞ and L 1 , i.e., The vanishing moments and normalization (16) imply that the detail coefficientsh k J,b become small when the underlying function is locally smooth. Therefore it is possible to obtain [25] So the detail coefficients decay at the rate of 2 −JN 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 T D ε is introduced by where D ε := {(J, b, k) : |h k J,b | > ε} and the elements of H J are defined bȳ 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 A D ε be the approximation operator A D ε := M r J −1 T D ε M r J . 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 ε j =ā j−J ε withā > 1. Then the approximation error concerning to the set of significant details D ε is uniformly bounded concerning to L q (Ω), q ∈ [1, ∞], i.e., for some constant C thr > 0 independent of J, ε. Here P r J h and P r J,D ε h are the projections according to (11) corresponding to the coefficientsH J and A D εH J .

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 u(x) of Equation (1) can be expanded using multi-scale projection operator M r J based on multi-wavelets as follows where ⊗ is the Kronecker product andŨ = (Ũ T 1 , . . . ,Ũ T n ) T is a (1 × nr2 J ) vector whose elements are n unknown sub-vectorsŨ i with a dimension of (r2 J × 1) such that One can imagine two types of equations, linear and nonlinear. For the linear type, the i-th component of the vector function g(x, t, u(t)) has the form and it can be approximated by multi-scale operator, i.e., where A ji (j, i = 1 : n) are (r2 J × r2 J ) matrices. Integrating from 0 to x, we get where A ji (j, i = 1 : n) are (r2 J × r2 J ) matrices and the rest are (r2 J × 1) vectors. But if the j-th component of vector function g(x, t, u(t)) is nonlinear, one can consider the following expansion where G j is a (r2 J × 1) vector whose elements are nonlinear algebraic equations. In view of the Equations (23) and (24), and using operational matrix I ψ of integration for multi-wavelets introduced in [7,14,15], one can write = 0, linear, nonlinear.
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 e J = u − M r J (u) where u and M r J (u) are the exact and approximate solutions of nonlinear system (1), respectively. Let X be an open set in R and let g : Ω × X → R be a function such that g(x, t, u(x)) ∈ C r (Ω) for any u ∈ X and the condition (2) is satisfied. Furthermore, presume that f ∈ C r (Ω). Furthermore, assume that the residual e :=r − (r r J ) where the residualr is specified as and r r J is introduced in (27).
If u ∈ C r (Ω), and the method used to solve system (28) is convergent then one has where κ is a positive constant. Consequently, e → 0 when J → ∞.
where κ = A I ψ 2 ∞ . Now Equation (5) leads to the desired result

Theorem 3.
Let the condition of Theorem 2 be valid. Suppose that 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 where η is a small positive number that desire to zero.
Since Alpert multi-wavelets are orthonormal, one can write Therefore one can bound the error of u − M r J (u) via According to the theorem hypotheses, the methods used to solve the obtained system are convergent. So η := M r J (u) − u J will be very small. Inequality (31) is obtained using (5) and (32), i.e.,

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 where u and M r J (u) 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 ε are set to zero. Finally, one can find the rate of sparsity S ε which is defined by [28] where N 0 is the total number of elements and N ε the number of elements remaining after thresholding.

Example 1.
Let us run the proposed method on the following linear Volterra integral equation [5,29] The exact solution is u(x) = (e −x , 2 sin(x)). The effect of thresholding on L 2 -errors and sparsity percentage is reported in Table 1 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 ε = 10 −4 , 10 −2 . Table 1. Effects of parameters r, J and ε on sparsity and L 2 -error for Example 1.

Example 2.
Let us consider the following system of Volterra integral equation as a further example The solution is reported in [5,30] and is 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 The exact solutions of this system are u 1 = sin(x) and u 2 = cos(x) [8]. Table 2 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 2 -errors decrease. In Tables 3 and 4, 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 Tables 3 and 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.

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.