Numerical Solution of Linear Volterra Integral Equation Systems of Second Kind by Radial Basis Functions

In this paper we propose an approximation method for solving second kind Volterra integral equation systems by radial basis functions. It is based on the minimization of a suitable functional in a discrete space generated by compactly supported radial basis functions of Wendland type. We prove two convergence results, and we highlight this because most recent published papers in the literature do not include any. We present some numerical examples in order to show and justify the validity of the proposed method. Our proposed technique gives an acceptable accuracy with small use of the data, resulting also in a low computational cost.


Introduction
A considerable large amount of research literature and books on the theory and applications of Volterra's integral equations have emerged over many decades since the apparition of Volterra's book "Leçons sur les équations intégrales et intégro-différentielles" [1] in 1913.
The applications include elasticity, plasticity, semi-conductors, scattering theory, seismology, heat and mass conduction or transfer, metallurgy, fluid flow dynamics, chemical reactions, population dynamics, and oscillation theory, among many others (see for example [2]). Other important references more related with the numerics of this type of equation are [3,4].
In fact, Volterra integral equations (VIEs) appear naturally when we try to transform an initial value problem into integral form, so that the solution of this integral equation is usually much easier to obtain than the original initial value problem. In the same way, some nonlinear Volterra integral equations are equivalent to an initial-value problem for a system of ordinary differential equations (ODEs). So, some authors (like for example [5]) have sought to exploit this connection for the numerical solution of the integral equations as well, since very effective ODE codes are widely available.
Volterra integral equations arise in many usual applications of technology, engineering and science in general: as in population dynamics, the spread of epidemics, some Dirichlet problems in potential theory, electrostatics, mathematical modeling of radioactive equilibrium, the particle transport problems of astrophysics and reactor theory, radiative energy and/or heat transfer problems, other general heat transfer problems, oscillation of strings and membranes, the problem of momentum representation in quantum mechanics, etc. However, many other complex problems of mathematics, chemistry, biology, astrophysics and mechanics, can be expressed in the terms of Volterra integral equations. Moreover, some practical problems, where impulses arise naturally (like in population dynamics or many biological applications) or are caused by some control system (like electric circuit problems and simulations of semiconductor devices) can be modeled by a differential equation, an integral equation, an integro-differential equation, or a system of these equations all combined. approximating solutions of the system of Volterra integral equations. In [38] the authors develop a numerical technique for the solution of 2D Volterra integral equations based on a discretization method by using two-dimensional Bernstein's approximations. In [39] the authors discussed the solution of linear Volterra integral equations of second kind using Mohand transform. In [40] the authors propose Bernstein polynomials to present effective solution for the second kind linear Volterra integral equations with delay. In [41] the author presents a method to solve numerically Volterra integral equations of the first kind with separable kernels.
In this work, we will present some specific variational methods adopted to study and approximate systems of linear Volterra integral equations with the aid of Radial Basis Functions (RBFs) of Wendland type. Wendland functions are compactly supported radial basis functions, which makes calculations with them quite simple. However, the general Wendland family of functions are defined recursively, and to determine the actual functions to use in any software implementation many calculations had to be done by hand or with the aid of some symbolic software (see for example [42]). There are for the moment just a few articles dealing with this type of techniques, like for example [43][44][45][46][47]; so we think there is still a lot to investigate in this regard.
Our goal in this work is to devise an appropriate approach procedure that is capable of solving this type of problem in a precise and efficient way. We consider then the linear Volterra equations system of the second kind as follows (see for example [1]): where . . , f n (t)) , k(t, s) = (k ij (t, s)) 1≤i,j≤n .
We assume that (1) has a unique continuous solution for appropriate functions f . In any case, the equations system (1) can be re-written in operator form as an equation of second kind where K is an integral operator and I denotes the identity operator. It is usual to impose certain assumptions on compactness on the operator K (see [48], Section 2.8.1) in order to establish the existence and uniqueness of the solution of (1), that we will assume throughout the work. Moreover, in [49] the authors proposed another method to solve second kind Fredholm integral equation systems, but the discrete functional space chosen in that article has been the space of spline functions. While at first glance it might seem that both works are similar, especially in the way they are presented, the two methods are totally different, not only be the fact that the discretization spaces are different (so we have adapted the notations accordingly), while the proofs (except the very preliminary ones, that can be also adapted), above all the proofs of the convergence results, are completely different, due to their greater complexity.
The outline of the paper is as follows. In Section 2 we briefly recall some notations and preliminaries. Section 3 is devoted to establish the discretization space as a radial basis functions space. The formulation of the minimization problem is realized in Section 4 and two equivalent variational problems are given. Section 5 is devoted to prove two convergence results. Section 6 deals with the description of the computation algorithm of the discrete problem solution. In Section 7 we present some numerical experiments and finally, in Section 8 we establish the conclusions of the work.

Notations and Preliminaries
Let R + 0 = {x ∈ R : x ≥ 0}; and for n ≥ 1, we denote by · n and · , · n the Euclidean norm and the inner product in R n .
On the other hand, for m ≥ 1, we designate by H m ((0, 1); R n ) the Sobolev space of order m of (classes of) functions u ∈ L 2 ((0, 1); R n ) together with all j-th derivative functions u (j) of order j ≤ m, in the sense of distributions. This space is equipped with • the semi-inner products, for any u, v ∈ H m ((0, 1); R n ), • the corresponding semi-norms |u| j = (u, u) • and the corresponding norm u m = ((u, u)) 1 2 m . For any 1 ≤ i ≤ n, let k i be a given function of the Sobolev vectorial functions space H m ((0, 1) × (0, 1); R n ) and consider the matrix valued function together with the associated integral operator Let R n,p be the space of real matrices of n lines and p columns, equipped with the inner product 1≤j≤p ∈ R n,p , and the corresponding norm A n,p = A, A 1 2 n,p .

Discretization Space
For the remainder of the work, we are going to consider a space of finite dimension, where we will formulate and solve a discrete approximation problem. The discrete functional space we have chosen is the radial basis functions space with compact support, namely the radial basis function space generated by the Wendland functions (see [50]).

Definition 1. Given a continuous function
, and a point ξ ∈ Ω, the radial function defined on Ω from the function φ with center ξ is the continuous function Then Φ ξ only depends of the distance to ξ.

Definition 2.
Given a centers set Ξ = {ξ 1 , . . . , ξ N } the linear space generated by the functions is called a radial basis functions space.

Definition 3.
For a function u ∈ C([0, 1]; R n ), the radial basis function interpolating u on a set of distinct centers T N = {t 1 , . . . , t N } ⊂ [0, 1] is given by where φ : R + 0 → R is a continuous function and the coefficients α 1 , . . . , α N ∈ R n are determined by the interpolation conditions In [50] H. Wendland introduced a family of compactly supported radial basis functions in the following way: let the operator I and its inverse D for r ≥ 0 be given by where x denotes the largest integer less than or equal to x.
with a univariate polynomial p d,k of degree d 2 + 3k + 1. They possess continuous derivatives up to order 2k, and they are of minimal degree for a given constant factor, uniquely determined by this setting.
Thus, these functions are the natural candidates for interpolation by compactly supported radial basis functions, and they are called the Wendland's functions.
For the remainder of the work we suppose 0 ≤ k ≤ N − 1, and we take φ = φ 1,k in Definition 3. Table 1 shows the Wendland functions φ 1,k for k = 0, 1, 2, and its continuity order.
Let S N be the space of the restrictions of functions on [0, 1] of the functional space generated by the radial basis functions

Formulation of the Problem
We can define the operator ρ : H k+1 ((0, 1); R n ) → R N,n given by Let assume that f ∈ H k+1 ((0, 1); R n ) and consider the affine variety H N = {u ∈ S N : ρu = ( f (t i )) 1≤i≤N } and the linear subspace H 0 N = {u ∈ S N : ρu = 0 ∈ R N,n }. Proof. By adapting the notations, as in the proof of Proposition 4.1 of [49].
is an inner product on H k+1 ((0, 1); R n ) and its associated norm, given by [[u]] =< < u, u > > Definition 4. We say that u N ∈ H N is an approximating radial basis function relative to T N , ρ and f if u N is a solution of the following minimization problem: where J : H k+1 ((0, 1); R n ) → R is given by Theorem 2. Problem (4) has a unique solution u N ∈ H N which is the unique solution of the variational problem and, taking into account that H 0 N is a vector space, we obtain that ∀v ∈ H 0 N , < < u N , v > >= 0.
Therefore (5) holds. Finally, u N is the unique solution of (4) since Theorem 3. There exists a unique λ ∈ R N,n such that where u N is the unique solution of (5).
Proof. For i = 1, . . . , N, let us consider ϕ i ∈ S N the unique radial basis function determined by the interpolation conditions Let take v ∈ S N , and we consider the function that is ρw = 0 ∈ R N,n , and in fact w ∈ H 0 N . Thus, from Theorem 2, we have We notice Π : R n → R, for = 1, . . . , n, the projection application given by Π (x 1 , . . . , x n ) = x .

Convergence Result
Assume that f ∈ H k+1 ((0, 1); R n ) and k ∈ H k ((0, 1) × (0, 1); R n,n ), then there exists a unique solution x ∈ H k+1 ((0, 1); R n ) of (1). Moreover, the following convergence result is verified. Theorem 4. Suppose given f ∈ H k+1 ((0, 1); R n ) and k ∈ H k ((0, 1) × (0, 1); R n,n ). Let denote x ∈ H k+1 ((0, 1); R n ) the unique solution of (1) and u N ∈ H N the unique solution of (4). Suppose that the hypothesis (2) holds, where h is mentioned. Then, one has Proof. Let s x,T N be the interpolating radial basis function of x on T N from the Wendland function φ 1,k , then s x,T N ∈ S N . Thus J(u N ) ≤ J(s x,T N ), that also implies that In this case, we have From this, and that the operator (I − K) is linear and compact in the finite-dimensional space S N , and thus bijective, we can deduce that there exists C 1 > 0 verifying Taking into account (3), it is verified that there exists C 2 > 0 such and, from here and (8) we obtain that there exists C > 0 such that Thus, the family (u N ) N∈N is bounded in H k+1 ((0, 1); R n ), and consequently there exists a sequence (u N ) ∈N extracted from this family, and an element x * ∈ H k+1 ((0, 1); R n ) such that x * = lim →+∞ u N weakly in H k+1 ((0, 1); R n ). (9) Suppose that x * = x; then, from the continuous injection of H k+1 ((0, 1); R n ) into C([0, 1]; R n ), there exists γ > 0 and a nonempty interval ω ⊂ [0, 1] such that As this injection is compact, from (9) Thus, for any ≥ 0 and t ∈ ω it is verified On the other hand, as we are taking h → 0 along the whole process, using the density condition (2) we can assure that there exists ∈ N and t * ∈ ω such that t * ∈ T N ∩ ω and thus (I − K)u N (t * ) = (I − K)x(t * ).
The operator I − K, considering the hypotheses taken from the beginning, it is also a bijection in C((0, 1); R n ), and thus u N (t * ) = x(t * ), which is a contradiction with (10). Thus For any ∈ N it is verified Then, from (9) and the compact inclusion of H k+1 ((0, 1); R n ) into H k ((0, 1); R n ) (see for example [48]), one has lim Suppose now that u N − x k does not tend to 0 as h tends to 0; in this case, it would exist α > 0, and a sequence (u N ) ∈N such that However, the sequence (u N ) ∈N is bounded in H k+1 ((0, 1); R n ) and then, by reasoning as above, we deduce that from this sequence we can extract a subsequence convergent to x in H k ((0, 1); R n ), what contradicts (12). Thus Then, from here the result is obtained.

Numerical Examples
To check the validity of the described method for approximating the solution of Problem (1) we present some numerical experiments.
In order to show the accuracy of the method, we have computed two relative error estimations, given by the expressions which estimates how close u N is to the solution of (1) and which is an approximation of the relative error of u N with respect to x in L 2 ((0, 1); R n ) being {a 1 , . . . , a 1000 } ⊂ [0, 1] thousand distinct random points.
From Theorem 4 and Corollary 1, these relative error estimations E 1 and E 2 tend to 0 as h tends to 0.

Example 1.
We consider the following Volterra equation system of order 2 The exact solution is x 1 (t) = t, x 2 (t) = t 2 . Table 2 shows the relative error estimations for distinct values of N.
The exact solution is x 1 (t) = e t , x 2 (t) = e −t . Table 3 shows the relative error estimations for distinct values of N.
The exact solution is Table 4 shows the relative error estimations for distinct values of N.

Conclusions
We conclude that the above presented experiments (see Tables 1-4)) confirm the validity of the method and justify the convergence results given in Theorem 4 and Corollary 1. In fact, in all our experiments (see the Examples 1-3), by using small values of N, we obtain a significant good order of approximation using the relative errors E 1 and E 2 considered. So, our original goal to devise an appropriate variational procedure that is capable of solving this type of problems in a precise and efficient way has been completely accomplished.
As compared with the other recently published works, for example [36][37][38]40,41], they do not study convergence results. Likewise, our technique gives an acceptable accuracy with a small use of data, resulting also a low computational cost.
In ([37], Tables 1 and 2) the mean of the error is of the order 10 −5 . We have obtained the same order of error with only 50 points.
In [40] the authors use Bernstein polynomials and the degree of its approximation is of order 10 −4 in most of the tables. The same happens in ( [41], Table 4), it uses the simple block-by-block method and its degree of approximation is about 10 −3 .
In order to do more research on this topic in the future, among some of the open problems that we consider are: -a numerical comparison between our method and many others in the literature, - the theoretical study of the order of convergence of the presented method, -the adaptation of this procedure to find the numerical solution of the linear systems of 2D Volterra integral equations of the second kind.