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

Abstract

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.

1. 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.

The systems of integral and/or integro-differential equations are usually difficult to solve analytically, in particular systems of Volterra integral non-linear equations or with variable coefficients; so a numerical method is often needed. In such cases, it is required to approximate the solutions; and many different numerical techniques have been developed and presented during decades of research, with appropriate combinations of numerical integration and interpolation procedures (see the references [3,6], among others). In order to approximate numerically the solution of general integral equations, the predominant technique have been the use of some kind of piecewise constant basis functions (PCBFs) (see for example [7], among many others); Chebyshev polynomial ([8] and others). However, after a long period of time many other techniques have attracted much attention recently; like wavelets theory, started with the introduction of Haar function in 1910, and from 1990’s (see [9]) also many wavelet type methods have been applied for solving integral equations. Haar wavelets, despite its relative simplicity, have many valued properties: as its compact support and orthogonality properties. So they can be used for the solution of differential and integro-differential equations related with signal and image processing, for example. They have been also used to solve linear and nonlinear integral equations by Aziz et. al. [10], Babolian et. al. [11], Lepik [12], Maleknejad et al. [13], Farshid Mirzaae [14], among others. More recently, several numerical methods based on different triangular type and delta orthogonal functions were designed for approximating the solution of integral and/or integro-differential Volterra equations (see for example [15,16,17], and the references therein). All these publications have demonstrated and revealed that these techniques based on PCBF and wavelets are effective to obtain the solution of such integral equations.

Particularly, systems of linear integral equations, and their exact or approximate solutions, are of great importance in science and engineering. There are several numerical methods for solving systems of linear Volterra integral equations of the second kind, and they have been often solved by classical numerical and analytical methods: such as Galerkin and Finite Element methods, collocation and spectral methods, Taylor or Power series and expansion methods, transforming the equations into a linear or nonlinear system of algebraic equations, and so on. However, new methods also have been applied to solve them, like the homotopy perturbation method [18], Adomian decomposition method (and many others) [19], use of Legendre wavelets [20] or hybrid Legendre and block pulse functions [21], Chebyshev polynomials [22,23], etc. Berenguer et al. [24,25] have solved them with the aid of a combination of analytical methods and bi-orthogonal systems in Banach spaces, Sahn et al. [26] have used Bessel polynomials method, Malnekad et al. [15] have employed delta basis functions (DBFs), Balakumar et al. [27] have applied the block-pulse functions method, Li-Hong et al. [28] have applied reproducing kernel method. Furthermore, there are also expansion methods for integral equations such as El-gendi’s and Wolfe’s methods (see for example [29]). Additionally, the approximate solutions of systems of integral equations that usually appear in problems of physics, biology and engineering are based on numerical integration methods: such as Euler–Chebyshev or Runge–Kutta methods (see for example [30]).

Concerning many other possible techniques to solve these types of integral equations, Draidi and Qatanani [31] implemented a product Nystrom and sinc-collocation methods to solve Volterra integral equations with Carleman kernel; also Issa, Qatanani and Daraghmeh [32] used a Taylor expansion and the variational iteration methods to give an approximate solution of Volterra integral equations of the second kind. Aggarwal et al. [33] and Chauhan [34] used different integral transformations for obtaining the solutions of VIEs of second kind. Mahgoub [35] solved constant coefficient linear differential equations by defining the called Sawi transformation, but many other authors exploited this idea, or other appropriate transforms, to deal with these types of integral or integro-differential equations.

Next, we are going to cite the most recent references from the last 3 or 4 years. In [36] the authors present an approximation solution of system of Volterra integral equations of second kind in an analytical way, using an Adomian decomposition method in Mathematica. In [37] the authors propose a numerical algorithm based on Monte Carlo method for 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]):

$$\mathit{x}\left(t\right)=\mathit{f}\left(t\right)+{\int }_0^t\mathit{k}(t,s)\mathit{x}\left(s\right)ds,0\le s\le t\le 1,$$

(1)

where

$$\begin{array}{c}\mathit{x}\left(t\right)={(x_1\left(t\right),\dots ,x_n\left(t\right))}^{\top },\hfill \\ \mathit{f}\left(t\right)={(f_1\left(t\right),\dots ,f_n\left(t\right))}^{\top },\hfill \\ \mathit{k}(t,s)={\left(k_{ij}(t,s)\right)}_{1\le i,j\le n}.\hfill \end{array}$$

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

$$\mathit{f}=(I-K)\mathit{x},$$

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.

2. Notations and Preliminaries

Let ${\mathbb{R}}_0^{+}=\{x\in \mathbb{R}:x\ge 0\}$; and for $n\ge 1$, we denote by ${\langle ·\rangle }_n$ and ${\langle ·,·\rangle }_n$ the Euclidean norm and the inner product in ${\mathbb{R}}^n$.

On the other hand, for $m\ge 1$, we designate by $H^m((0,1);{\mathbb{R}}^n)$ the Sobolev space of order m of (classes of) functions $\mathit{u}\in L^2((0,1);{\mathbb{R}}^n)$ together with all j-th derivative functions ${\mathit{u}}^{\left(j\right)}$ of order $j\le m$, in the sense of distributions. This space is equipped with

• the semi–inner products, for any $\mathit{u},\mathit{v}\in H^m((0,1);{\mathbb{R}}^n)$,

  $${(\mathit{u},\mathit{v})}_j={\int }_0^1{\langle {\mathit{u}}^{\left(j\right)}\left(t\right),{\mathit{v}}^{\left(j\right)}\left(t\right)\rangle }_{n}dt,0\le j\le m,$$
• the corresponding semi–norms ${\displaystyle {\left|\mathit{u}\right|}_j={(\mathit{u},\mathit{u})}_j^{\frac{1}{2}}}$, for $0\le j\le m$,
• the inner product ${\displaystyle {\left((\mathit{u},\mathit{v})\right)}_m=\sum _0^m{(\mathit{u},\mathit{v})}_j}$,
• and the corresponding norm ${\displaystyle {\parallel \mathit{u}\parallel }_m={\left((\mathit{u},\mathit{u})\right)}_m^{\frac{1}{2}}}$.

For any $1\le i\le n$, let ${\mathit{k}}_i$ be a given function of the Sobolev vectorial functions space $H^m((0,1)\times (0,1);{\mathbb{R}}^n)$ and consider the matrix valued function

$$\mathit{k}(t,s)={\left({\mathit{k}}_i(t,s)\right)}_{1\le i\le n}\in H^m((0,1)\times (0,1);{\mathbb{R}}^{n,n}),$$

together with the associated integral operator

$$K\mathit{u}\left(t\right)={\left({\int }_0^t{\langle {\mathit{k}}_i(t,s),\mathit{u}\left(s\right)\rangle }_{n}ds\right)}_{1\le i\le n},t\in (0,1),\forall \mathit{u}\in H^m((0,1);{\mathbb{R}}^n).$$

Let ${\mathbb{R}}^{n,p}$ be the space of real matrices of n lines and p columns, equipped with the inner product

$${\langle A,B\rangle }_{n,p}=\sum _{i=1}^n\sum _{j=1}^p{a}_{ij}{b}_{ij},\forall A={\left(a_{ij}\right)}_{\genfrac{}{}{0pt}{}{1\le i\le n}{1\le j\le p}},B={\left(b_{ij}\right)}_{\genfrac{}{}{0pt}{}{1\le i\le n}{1\le j\le p}}\in {\mathbb{R}}^{n,p},$$

and the corresponding norm ${\displaystyle {\langle A\rangle }_{n,p}={\langle A,A\rangle }_{n,p}^{\frac{1}{2}}}$.

3. 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 $\varphi :{\mathbb{R}}_0^{+}\to \mathbb{R}$, a subset $\Omega \subset {\mathbb{R}}^d$, $d\ge 1$, and a point $\mathit{\xi }\in \Omega$, the radial function defined on Ω from the function ϕ with center $\mathit{\xi }$ is the continuous function ${\Phi }_{\mathit{\xi }}:\Omega \to R$ given by

$${\Phi }_{\mathit{\xi }}\left(\mathit{x}\right)=\varphi \left({\langle \mathit{x}-\mathit{\xi }\rangle }_d\right).$$

Then ${\Phi }_{\mathit{\xi }}$ only depends of the distance to $\mathit{\xi }$.

Definition 2.

Given a centers set $\Xi =\{{\mathit{\xi }}_1,\dots ,{\mathit{\xi }}_N\}$ the linear space generated by the functions

$$\{\varphi {\left(\langle ·-{\mathit{\xi }}_1\rangle \right)}_d,\dots ,\left\{\varphi {\left(\langle ·-{\mathit{\xi }}_N\rangle \right)}_d\right\}$$

is called a radial basis functions space.

Definition 3.

For a function $\mathit{u}\in C([0,1];{\mathbb{R}}^n)$, the radial basis function interpolating $\mathit{u}$ on a set of distinct centers ${\mathit{T}}_N=\{t_1,\dots ,t_N\}\subset [0,1]$ is given by

$${\mathit{s}}_{\mathit{u},{\mathit{T}}_N}\left(t\right)=\sum _{i=1}^N{\mathit{\alpha }}_i\varphi \left(\right|t-t_i\left|\right),t\in [0,1],$$

where $\varphi :{\mathbb{R}}_0^{+}\to \mathbb{R}$ is a continuous function and the coefficients ${\mathit{\alpha }}_1,\dots ,{\mathit{\alpha }}_N\in {\mathbb{R}}^n$ are determined by the interpolation conditions

$${\mathit{s}}_{\mathit{u},{\mathit{T}}_N}\left(t_i\right)=\mathit{u}\left(t_i\right),1\le i\le N.$$

In [50] H. Wendland introduced a family of compactly supported radial basis functions in the following way: let the operator $\mathcal{I}$ and its inverse $\mathcal{D}$ for $r\ge 0$ be given by

$$\begin{array}{c}{\displaystyle \left(\mathcal{I}\varphi \right)\left(r\right)={\int }_r^{\infty }t\varphi \left(t\right)dt,}\hfill \\ {\displaystyle \left(\mathcal{D}\varphi \right)\left(r\right)=-\frac{1}{r}{\varphi }^{\prime }\left(r\right),}\hfill \end{array}$$

for any differentiable function $\varphi :{\mathbb{R}}_0^{+}\to \mathbb{R}$.

Given the truncated power function ${\varphi }_{\ell }\left(r\right)={(1-r)}_{+}^{\ell }$, we define

$${\varphi }_{d,k}={\mathcal{I}}^k{\varphi }_{\lfloor \frac{d}{2}\rfloor +k+1},$$

where $\lfloor x\rfloor$ denotes the largest integer less than or equal to x.

Theorem 1. 

([50], Theorem 1.2) The functions ${\varphi }_{d,k}$ induce positive definite functions on ${\mathbb{R}}^d$ of the form

$${\varphi }_{d,k}\left(r\right)=\left\{\begin{array}{cc}{p}_{d,k}\left(r\right),\hfill & 0\le r\le 1,\hfill \\ 0,\hfill & r>1,\hfill \end{array}\right.$$

with a univariate polynomial $p_{d,k}$ of degree $\lfloor \frac{d}{2}\rfloor +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\le k\le N-1$, and we take $\varphi ={\varphi }_{1,k}$ in Definition 3.

Table 1. The Wendland functions ${\varphi }_{1,k}$ for $k=0,1,2$ and its continuity order.

k	Wendland Function	Continuity Order
$k=0$	${\varphi }_{1,0}\left(r\right)={(1-r)}_{+}$	$C^0$
$k=1$	${\varphi }_{1,1}\left(r\right)\doteq {(1-r)}_{+}^3(3r+1)$	$C^2$
$k=2$	${\varphi }_{1,2}\left(r\right)\doteq {(1-r)}_{+}^5(8r^2+5r+1)$	$C^4$

shows the Wendland functions ${\varphi }_{1,k}$ for $k=0,1,2$, and its continuity order.

Let

$$h=\underset{t\in [0,1]}{sup}\underset{1\le i\le N}{min}|t-t_i|.$$

(2)

From ([50], Theorem 2.1) we can affirm that ${\varphi }_{1,k}\in C^{2k}\left([0,1]\right)$ and the corresponding native space is $H^{k+1}\left([0,1]\right)$. Finally, from ([50], Theorem 2.1) and ([51], Theorem 4.1) we conclude that there exists $C>0$ such that

$$\parallel \mathit{u}-{\mathit{s}}_{\mathit{u},{\mathit{T}}_N}{\parallel }_{L_{\infty }((0,1);{\mathbb{R}}^n)}\le C{\parallel \mathit{u}\parallel }_{k+1}{h}^{k+\frac{1}{2}},\forall \mathit{u}\in H^{k+1}([0,1];{\mathbb{R}}^n),$$

and

$$|\mathit{u}-{\mathit{s}}_{\mathit{u},{\mathit{T}}_N}{|}_j\le C{h}^{k+1-j}{\parallel \mathit{u}\parallel }_{k+1},0\le j\le k+1,\forall \mathit{u}\in H^{k+1}([0,1];{\mathbb{R}}^n).$$

(3)

Let $S_N$ be the space of the restrictions of functions on $[0,1]$ of the functional space generated by the radial basis functions $\{{\varphi }_{1,k}\left(\right|·-t_1\left|\right),\dots ,{\varphi }_{1,k}\left(\right|·-t_N\left|\right)\}$ and ${\mathit{S}}_N={\left(S_N\right)}^n$. Then ${\mathit{S}}_N\subset H^{k+1}((0,1);{\mathbb{R}}^n)\cap C^{2k}([0,1];{\mathbb{R}}^n)$.

4. Formulation of the Problem

We can define the operator $\rho :H^{k+1}((0,1);{\mathbb{R}}^n)\to {\mathbb{R}}^{N,n}$ given by

$$\rho \mathit{v}={\left((I-K)\mathit{v}\left(t_i\right)\right)}_{1\le i\le N}.$$

Let assume that $\mathit{f}\in H^{k+1}((0,1);{\mathbb{R}}^n)$ and consider the affine variety ${\mathit{H}}_N=\{\mathit{u}\in {\mathit{S}}_N:\rho \mathit{u}={\left(\mathit{f}\left(t_i\right)\right)}_{1\le i\le N}\}$ and the linear subspace ${\mathit{H}}_N^0=\{\mathit{u}\in {\mathit{S}}_N:\rho \mathit{u}=\mathit{0}\in {\mathbb{R}}^{N,n}\}$.

Proposition 1.

The set ${\mathit{H}}_N$ is a nonempty closed bounded convex subset of ${\mathit{S}}_N$. Moreover it is an affine variety associated with the linear subspace ${\mathit{H}}_N^0$.

Proof. 

By adapting the notations, as in the proof of Proposition 4.1 of [49]. □

Lemma 1.

The application $<<·,·>>:H^{k+1}((0,1);{\mathbb{R}}^n)\times H^{k+1}((0,1);{\mathbb{R}}^n)\to \mathbb{R}$ defined by

$$<<\mathit{u},\mathit{v}>>={\langle \rho \mathit{u},\rho \mathit{v}\rangle }_{N,n}+{((I-K)\mathit{u},(I-K)\mathit{v})}_{k+1}$$

is an inner product on $H^{k+1}((0,1);{\mathbb{R}}^n)$ and its associated norm, given by ${\displaystyle \left[\left[\mathit{u}\right]\right]=<<\mathit{u},\mathit{u}>{>}^{\frac{1}{2}}}$, is equivalent to the usual Sobolev norm ${\parallel ·\parallel }_{k+1}$.

Proof. 

By adapting the notations as in the proof of Lemma 4.2 of [49] and using ([48], Theorem 7.3.12) the proof can be obtained. □

Definition 4.

We say that ${\mathit{u}}_N\in {\mathit{H}}_N$ is an approximating radial basis function relative to ${\mathit{T}}_N$, ρ and $\mathit{f}$ if ${\mathit{u}}_N$ is a solution of the following minimization problem:

$$Find {\mathit{u}}_N\in {\mathit{H}}_N such that \forall \mathit{v}\in {\mathit{H}}_N,J\left({\mathit{u}}_N\right)\le J\left(\mathit{v}\right),$$

(4)

where $J:H^{k+1}((0,1);{\mathbb{R}}^n)\to \mathbb{R}$ is given by

$$J\left(\mathit{v}\right)={\left|(I-K)\mathit{v}\right|}_{k+1}^2.$$

Theorem 2.

Problem (4) has a unique solution ${\mathit{u}}_N\in {\mathit{H}}_N$ which is the unique solution of the variational problem

$$\forall \mathit{v}\in {\mathit{H}}_N^0,{((I-K){\mathit{u}}_N,(I-K)\mathit{v})}_{k+1}=0.$$

(5)

Proof. 

From Proposition 1 and ([48], Theorem 3.4.3) we can deduce that there exists a unique ${\mathit{u}}_N\in {\mathit{H}}_N$, which is the projection of $\mathit{0}$ on ${\mathit{H}}_N$ such that

$$\left[\left[{\mathit{u}}_N\right]\right]\le \left[\left[\mathit{v}\right]\right],\forall \mathit{v}\in {\mathit{H}}_N$$

and verifying

$$\forall \mathit{w}\in {\mathit{H}}_N,<<-{\mathit{u}}_N,\mathit{w}-{\mathit{u}}_N>>\le 0,$$

that is

$$\forall \mathit{v}\in {\mathit{H}}_N^0,<<-{\mathit{u}}_N,\mathit{v}>>\le 0$$

and, taking into account that ${\mathit{H}}_N^0$ is a vector space, we obtain that

$$\forall \mathit{v}\in {\mathit{H}}_N^0,<<{\mathit{u}}_N,\mathit{v}>>=0.$$

Therefore (5) holds. Finally, ${\mathit{u}}_N$ is the unique solution of (4) since $J\left(\mathit{v}\right)={\left[\left[\mathit{v}\right]\right]}^2-{\langle \rho \mathit{f}\rangle }_{N,n}^2$, for any $\mathit{v}\in {\mathit{H}}_N$. □

Theorem 3.

There exists a unique $\lambda \in {\mathbb{R}}^{N,n}$ such that

$$\forall \mathit{v}\in {\mathit{S}}_N,{((I-K){\mathit{u}}_N,(I-K)\mathit{v})}_{k+1}+{\langle \lambda ,\rho \mathit{v}\rangle }_{N,n}=0,$$

(6)

where ${\mathit{u}}_N$ is the unique solution of (5).

Proof. 

For $i=1,\dots ,N$, let us consider ${\phi }_i\in S_N$ the unique radial basis function determined by the interpolation conditions

$${\phi }_i\left(t_j\right)={\delta }_{ij},\forall j=1,\dots N.$$

Let take $\mathit{v}\in {\mathit{S}}_N$, and we consider the function

$$\mathit{w}=\mathit{v}-\sum _{i=1}^N(I-K)\mathit{v}\left(t_i\right){\phi }_i,$$

then

$$(I-K)\mathit{w}\left(t_j\right)=(I-k)\mathit{v}\left(t_j\right)-\sum _{i=1}^N(I-K)\mathit{v}\left(t_i\right){\phi }_i\left(t_j\right)=0,\forall j=1,\dots ,N,$$

that is $\rho \mathit{w}=\mathit{0}\in {\mathbb{R}}^{N,n}$, and in fact $\mathit{w}\in {\mathit{H}}_N^0$. Thus, from Theorem 2, we have

$${((I-K){\mathit{u}}_N,(I-K)\mathit{w})}_{k+1}=0.$$

(7)

We notice ${\mathsf{\Pi }}_{\ell }:{\mathbb{R}}^n\to \mathbb{R}$, for $\ell =1,\dots ,n$, the projection application given by ${\mathsf{\Pi }}_{\ell }(x_1,\dots ,x_n)=x_{\ell }$.

Then, for $i=1,\dots ,N$, it is verified that

$$\begin{array}{c}\hfill {\displaystyle {\left((I-K){\mathit{u}}_N,(I-K)\mathit{v}\left(t_i\right){\phi }_i\right)}_{k+1}=\sum _{\ell =1}^n{\left({\mathsf{\Pi }}_{\ell }((I-K){\mathit{u}}_N,{\mathsf{\Pi }}_{\ell }\left((I-K)\mathit{v}\left(t_i\right){\phi }_i\right)\right)}_{k+1}}\\ \hfill {\displaystyle =\sum _{\ell =1}^n{\mathsf{\Pi }}_{\ell }\left((I-K)\mathit{v}\left(t_i\right)\right){\left({\mathsf{\Pi }}_{\ell }((I-K){\mathit{u}}_N,{\phi }_i\right)}_{k+1}.}\end{array}$$

Let denote ${\lambda }_{i\ell }=-{\left({\mathsf{\Pi }}_{\ell }((I-K){\mathit{u}}_N,{\phi }_i\right)}_{k+1}\in \mathbb{R}$ and $\lambda ={\left({\lambda }_{i\ell }\right)}_{\genfrac{}{}{0pt}{}{1\le i\le N}{1\le \ell \le n}}\in {\mathbb{R}}^{N,n}$, then

$$\begin{array}{c}{\displaystyle {\left((I-K){\mathit{u}}_N,(I-K)\mathit{w}\right)}_{k+1}=}\hfill \\ {\displaystyle {\left((I-K){\mathit{u}}_N,(I-K)\mathit{v}\right)}_{k+1}-\sum _{i=1}^N{\left((I-K){\mathit{u}}_N,(I-K)\mathit{v}\left(t_i\right){\phi }_i\right)}_{k+1}=}\hfill \\ {\displaystyle {\left((I-K){\mathit{u}}_N,(I-K)\mathit{v}\right)}_{k+1}+\sum _{i=1}^N\sum _{\ell =1}^n{\mathsf{\Pi }}_{\ell }\left((I-K)\mathit{v}\left(t_i\right)\right){\lambda }_{i\ell }=}\hfill \\ {\displaystyle {\left((I-K){\mathit{u}}_N,(I-K)\mathit{v}\right)}_{k+1}+{\langle \lambda ,\rho \mathit{v}\rangle }_{N,n}.}\hfill \end{array}$$

From (7), we conclude that there exists $\lambda ={(-{({\mathsf{\Pi }}_{\ell }\left((I-K){\mathit{u}}_n\right),{\phi }_i)}_{k+1})}_{\genfrac{}{}{0pt}{}{1\le i\le N}{1\le \ell \le n}}\in {\mathbb{R}}^{N,n}$ such that

$${\left((I-K){\mathit{u}}_N,(I-K)\mathit{v}\right)}_{k+1}+{\langle \lambda ,\rho \mathit{v}\rangle }_{N,n}=0$$

and (6) holds.

The uniqueness of $\lambda$ is immediate. □

5. Convergence Result

Assume that $\mathit{f}\in H^{k+1}((0,1);{\mathbb{R}}^n)$ and $\mathit{k}\in H^k((0,1)\times (0,1);{\mathbb{R}}^{n,n})$, then there exists a unique solution $\mathit{x}\in H^{k+1}((0,1);{\mathbb{R}}^n)$ of (1). Moreover, the following convergence result is verified.

Theorem 4.

Suppose given $\mathit{f}\in H^{k+1}((0,1);{\mathbb{R}}^n)$ and $\mathit{k}\in H^k((0,1)\times (0,1);{\mathbb{R}}^{n,n})$. Let denote $\mathit{x}\in H^{k+1}((0,1);{\mathbb{R}}^n)$ the unique solution of (1) and ${\mathit{u}}_N\in {\mathit{H}}_N$ the unique solution of (4). Suppose that the hypothesis (2) holds, where h is mentioned. Then, one has

$$\underset{h\to 0}{lim}{\parallel {\mathit{u}}_N-\mathit{x}\parallel }_k=0.$$

Proof. 

Let ${\mathit{s}}_{\mathit{x},{\mathit{T}}_N}$ be the interpolating radial basis function of $\mathit{x}$ on ${\mathit{T}}_N$ from the Wendland function ${\varphi }_{1,k}$, then ${\mathit{s}}_{\mathit{x},{\mathit{T}}_N}\in {\mathit{S}}_N$. Thus $J\left({\mathit{u}}_N\right)\le J\left({\mathit{s}}_{\mathit{x},{\mathit{T}}_N}\right)$, that also implies that

$$|(I-K){\mathit{u}}_N{|}_{k+1}\le {\left|(I-K){\mathit{s}}_{\mathit{x},{\mathit{T}}_N}\right|}_{k+1}.$$

In this case, we have

$$\left[\left[(I-K){\mathit{u}}_N\right]\right]\le \left[\left[(I-K){\mathit{s}}_{\mathit{x},{\mathit{T}}_N}\right]\right].$$

From this, and that the operator $(I-K)$ is linear and compact in the finite-dimensional space ${\mathit{S}}_N$, and thus bijective, we can deduce that there exists $C_1>0$ verifying

$$\parallel {\mathit{u}}_N{\parallel }_{k+1}\le C_1{\parallel {\mathit{s}}_{\mathit{x},{\mathit{T}}_N}\parallel }_{k+1}.$$

(8)

Taking into account (3), it is verified that there exists $C_2>0$ such

$$\parallel {\mathit{s}}_{\mathit{x},{\mathit{T}}_N}{\parallel }_{k+1}\le C_2{\parallel \mathit{x}\parallel }_{k+1}.$$

and, from here and (8) we obtain that there exists $C>0$ such that

$$\parallel {\mathit{u}}_N{\parallel }_{k+1}\le C{\parallel \mathit{x}\parallel }_{k+1}.$$

Thus, the family ${\left({\mathit{u}}_N\right)}_{N\in \mathbb{N}}$ is bounded in $H^{k+1}((0,1);{\mathbb{R}}^n)$, and consequently there exists a sequence ${\left({\mathit{u}}_{N_{\ell }}\right)}_{\ell \in \mathbb{N}}$ extracted from this family, and an element ${\mathit{x}}^{*}\in H^{k+1}((0,1);{\mathbb{R}}^n)$ such that

$${\mathit{x}}^{*}=\underset{\ell \to +\infty }{lim}{\mathit{u}}_{N_{\ell }} \mathrm{weakly} \mathrm{in} H^{k+1}((0,1);{\mathbb{R}}^n).$$

(9)

Suppose that ${\mathit{x}}^{*}\ne \mathit{x}$; then, from the continuous injection of $H^{k+1}((0,1);{\mathbb{R}}^n)$ into $C([0,1];{\mathbb{R}}^n)$, there exists $\gamma >0$ and a nonempty interval $\omega \subset [0,1]$ such that

$$\forall t\in \omega ,{\langle {\mathit{x}}^{*}-\mathit{x}\rangle }_n>\gamma .$$

As this injection is compact, from (9)

$$\exists {\ell }_0\in \mathbb{N},\forall \ell \ge {\ell }_0,{\langle {\mathit{u}}_{N_{\ell }}\left(t\right)-{\mathit{x}}^{*}\left(t\right)\rangle }_n\le \frac{\gamma }{2}.$$

Thus, for any $\ell \ge {\ell }_0$ and $t\in \omega$ it is verified

$${\langle {\mathit{u}}_{N_{\ell }}\left(t\right)-\mathit{x}\left(t\right)\rangle }_n\ge {\langle {\mathit{x}}^{*}\left(t\right)-\mathit{x}\left(t\right)\rangle }_n-{\langle {\mathit{u}}_{N_{\ell }}\left(t\right)-{\mathit{x}}^{*}\left(t\right)\rangle }_n>\frac{\gamma }{2}.$$

(10)

On the other hand, as we are taking $h\to 0$ along the whole process, using the density condition (2) we can assure that there exists $\ell \in \mathbb{N}$ and $t_{\ell }^{*}\in \omega$ such that $t_{\ell }^{*}\in T_{N_{\ell }}\cap \omega$ and thus

$$(I-K){\mathit{u}}_{N_{\ell }}\left(t_{\ell }^{*}\right)=(I-K)\mathit{x}\left(t_{\ell }^{*}\right).$$

The operator $I-K$, considering the hypotheses taken from the beginning, it is also a bijection in $C((0,1);{\mathbb{R}}^n)$, and thus ${\mathit{u}}_{N_{\ell }}\left(t_{\ell }^{*}\right)=\mathit{x}\left(t_{\ell }^{*}\right)$, which is a contradiction with (10). Thus ${\mathit{x}}^{*}=\mathit{x}$.

For any $\ell \in \mathbb{N}$ it is verified

$$\parallel {\mathit{u}}_{N_{\ell }}{-\mathit{x}\parallel }_k^2=\parallel {\mathit{u}}_{N_{\ell }}{\parallel }_k^2+{\parallel \mathit{x}\parallel }_k^2-2{({\mathit{u}}_{N_{\ell }},\mathit{x})}_k.$$

Then, from (9) and the compact inclusion of $H^{k+1}((0,1);{\mathbb{R}}^n)$ into $H^k((0,1);{\mathbb{R}}^n)$ (see for example [48]), one has

$$\underset{\ell \to +\infty }{lim}{\parallel {\mathit{u}}_{N_{\ell }}-\mathit{x}\parallel }_k=0.$$

(11)

Suppose now that $\parallel {\mathit{u}}_N{-\mathit{x}\parallel }_k$ does not tend to 0 as h tends to 0; in this case, it would exist $\alpha >0$, and a sequence ${\left({\mathit{u}}_{N_{\ell }^{\prime }}\right)}_{\ell \in \mathbb{N}}$ such that

$$\forall \ell \in \mathbb{N},\parallel {\mathit{u}}_{N_{\ell }^{\prime }}{-\mathit{x}\parallel }_k>\alpha .$$

(12)

However, the sequence ${\left({\mathit{u}}_{N_{\ell }^{\prime }}\right)}_{\ell \in \mathbb{N}}$ is bounded in $H^{k+1}((0,1);{\mathbb{R}}^n)$ and then, by reasoning as above, we deduce that from this sequence we can extract a subsequence convergent to $\mathit{x}$ in $H^k((0,1);{\mathbb{R}}^n)$, what contradicts (12). Thus

$$\underset{h\to 0}{lim}{\parallel {\mathit{u}}_N-\mathit{x}\parallel }_k=0.$$

□

Corollary 1.

Under the conditions of Theorem 4 one has

$$\underset{h\to 0}{lim}{\parallel \mathit{f}-(I-K){\mathit{u}}_N\parallel }_k=0.$$

Proof. 

From Theorem 4 and the continuity of the operator $I-K$ we have

$$\underset{h\to 0}{lim}(I-K){\mathit{u}}_N=(I-K)\mathit{x}=\mathit{f} \mathrm{in} H^k((0,1);{\mathrm{R}}^n).$$

Then, from here the result is obtained. □

6. Computation

Let us compute the unique solution of (6). The solution of problem (5) can be expressed by

$${\mathit{u}}_N=\sum _{i=1}^N{\mathit{\alpha }}_i{\varphi }_{1,k}\left(\right|·-t_i\left|\right),$$

with ${\mathit{\alpha }}_1,\dots ,{\mathit{\alpha }}_N\in {\mathbb{R}}^n$.

Consider the basis $\{{\mathit{B}}_1,\dots ,{\mathit{B}}_{Nn}\}$ of the space ${\mathit{S}}_N$ given, for $\ell =1,\dots Nn$, by

$${\mathit{B}}_{\ell }\left(t\right)={\varphi }_{1,k}\left(\right|t-t_i\left|\right){\mathit{e}}_j,$$

being $i=quotient(\ell -1,n)+1$ and $j=\ell -(i-1)n$.

Then, the solution of (5) can be expressed by

$${\mathit{u}}_N=\sum _{\ell =1}^{Nn}{\alpha }_{\ell }{\mathit{B}}_{\ell },$$

with ${\alpha }_1,\dots ,{\alpha }_{Nn}\in \mathbb{R}$.

By replacing in (6), we have

$$\sum _{\ell =1}^{Nn}{\alpha }_{\ell }{\left((I-K){\mathit{B}}_{\ell },(I-K)\mathit{v}\right)}_{k+1}+{\langle \lambda ,\rho \mathit{v}\rangle }_{N,n}=0,\forall \mathit{v}\in {\mathit{S}}_N,$$

subject to the restrictions

$$\sum _{\ell =1}^{Nn}{\alpha }_{\ell }(I-K){\mathit{B}}_{\ell }\left(t_i\right)=\mathit{f}\left(t_i\right),i=1,\dots ,N.$$

Taking $\mathit{v}={\mathit{B}}_j$, for $j=1,\dots ,Nn$, we obtain a linear system of order $2Nn$ with unknowns ${\alpha }_1,\dots ,{\alpha }_{Nn},{\lambda }_1,\dots ,{\lambda }_{Nn}\in \mathbb{R}$, that can be expressed in matrix form as follows:

$$\left(\begin{array}{cc}\mathcal{C}\hfill & \mathcal{D}\hfill \\ {\mathcal{D}}^{\top }\hfill & 0\hfill \end{array}\right)\left(\begin{array}{c}\mathit{\alpha }\hfill \\ \mathit{\lambda }\hfill \end{array}\right)=\left(\begin{array}{c}\mathit{0}\hfill \\ \mathit{F}\hfill \end{array}\right),$$

with

$$\begin{array}{c}\mathcal{C}={\left({((I-K){\mathit{B}}_{\ell },(I-K){\mathit{B}}_j)}_{k+1}\right)}_{\genfrac{}{}{0pt}{}{1\le \ell \le Nn}{1\le j\le Nn}},\hfill \\ \mathcal{D}={\left(d_{ij}\right)}_{\genfrac{}{}{0pt}{}{1\le i\le Nn}{1\le j\le Nn}},\hfill \\ \mathit{\alpha }={({\alpha }_1,\dots ,{\alpha }_{Nn})}^{\top },\mathit{\lambda }={({\lambda }_1,\dots ,{\lambda }_{Nn})}^{\top },\hfill \\ \mathit{F}={\left(f_i\right)}_{1\le i\le Nn},\hfill \end{array}$$

being, for $i=1\dots ,Nn$ and $j=1,\dots ,Nn$,

$$d_{ij}={\mathsf{\Pi }}_{\ell }\left((I-K){\mathit{B}}_r\left(t_s\right)\right),$$

with $r=quotient(i-1,n)+1$, $s=quotient(j-1,n)+1$, $\ell =j-(s-1)n$ and for $i=1,\dots ,Nn$,

$$f_i={\mathsf{\Pi }}_{\ell }\left(\mathit{f}\left(t_s\right)\right),$$

with $s=quotient(i-1,n)$ and $\ell =i-(s-1)n$.

7. 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

$$E_1=\frac{1}{1000}\sum _{i=1}^{1000}{\langle \mathit{f}\left(a_i\right)-(I-K){\mathit{u}}_N\left(a_i\right)\rangle }_n,$$

which estimates how close ${\mathit{u}}_N$ is to the solution of (1) and

$$E_2=\sqrt{\frac{{\displaystyle \sum _{i=1}^{1000}{\langle {\mathit{u}}_N\left(a_i\right)-\mathit{x}\left(a_i\right)\rangle }_n^2}}{{\displaystyle \sum _{i=1}^{1000}{\langle \mathit{x}\left(a_i\right)\rangle }_n^2}}},$$

which is an approximation of the relative error of ${\mathit{u}}_N$ with respect to $\mathit{x}$ in $L^2((0,1);{\mathrm{R}}^n)$ being $\{a_1,\dots ,a_{1000}\}\subset [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.

Moreover, in all the examples, the discrete space that we use to calculate the approximated solution ${\mathit{u}}_N$ is the radial basis function space constructed from the Wendland function ${\varphi }_{1,1}$ and the centers set ${\displaystyle {\mathit{T}}_N=\{t_i=\frac{i}{N},i=0,\dots ,N\}}$.

In order to compute the numerical integrals, we have used the following quadrature formula (see [52])

$$\begin{array}{ccc}{\displaystyle {\int }_a^{b}g\left(t\right)dt}\hfill & \approx \hfill & {\displaystyle \sum _{i=6}^{n-3}g\left({\xi }_i\right)+h\left(\frac{206}{1575}(g\left({\xi }_1\right)+g\left({\xi }_{n+2}\right))+\frac{107}{128}(g\left({\xi }_2\right)+g\left({\xi }_{n+1}\right))+\right.}\hfill \\ {\displaystyle }\hfill & & {\displaystyle \frac{6019}{5760}(g\left({\xi }_3\right)+g\left({\xi }_n\right))+\frac{9467}{9600}(g\left({\xi }_4\right)+g\left({\xi }_{n-1}\right))+}\hfill \\ & & {\displaystyle \left.\frac{13,469}{13,440}(g\left({\xi }_5\right)+g\left({\xi }_{n-2}\right))\right),}\hfill \end{array}$$

where ${\displaystyle h=\frac{b-a}{n}}$ and

$${\xi }_1=a,{\xi }_{n+2}=b,{\xi }_i=a+\frac{2i-1}{2}h,i=2,\dots ,n+1.$$

This formula has an error order of $O\left(h^6\right)$ for $g\in C^6\left([a,b]\right)$.

Example 1.

We consider the following Volterra equation system of order 2

$$\left\{\begin{array}{c}{\displaystyle x_1\left(t\right)-{\int }_0^t({(t-s)}^3{x}_1\left(s\right)+{(t-s)}^2{x}_2\left(s\right))ds=t-\frac{t^5}{12},}\hfill \\ {\displaystyle x_2\left(t\right)-{\int }_0^t({(t-s)}^4{x}_1\left(s\right)+{(t-s)}^3{x}_2\left(s\right))ds=t^2-\frac{t^6}{20}.}\hfill \end{array}\right.$$

The exact solution is

$$x_1\left(t\right)=t,x_2\left(t\right)=t^2.$$

Table 2. Computed relative error estimations for Example 1 from some values of N.

N	${\mathit{E}}_{\mathbf{1}}$	${\mathit{E}}_{\mathbf{2}}$
5	$2.1868\times {10}^{-2}$	$3.1058\times {10}^{-2}$
10	$3.6034\times {10}^{-3}$	$4.8048\times {10}^{-3}$
20	$6.2683\times {10}^{-4}$	$8.2990\times {10}^{-4}$
30	$2.0727\times {10}^{-4}$	$3.0254\times {10}^{-4}$
40	$1.0215\times {10}^{-4}$	$1.2509\times {10}^{-4}$
50	$6.4520\times {10}^{-5}$	$9.2824\times {10}^{-5}$

shows the relative error estimations for distinct values of N.

Example 2.

We consider the following Volterra equation system of order 2

$$\left\{\begin{array}{c}{\displaystyle x_1\left(t\right)-{\int }_0^t(e^{t-s}{x}_1\left(s\right)+e^{t+s}{x}_2\left(s\right))ds=e^t(1-2t),}\hfill \\ {\displaystyle x_2\left(t\right)-{\int }_0^t(-e^{t-s}{x}_1\left(s\right)+e^{t+s}{x}_2\left(s\right))ds=e^{-t}.}\hfill \end{array}\right.$$

The exact solution is

$$x_1\left(t\right)=e^t,x_2\left(t\right)=e^{-t}.$$

Table 3. Computed relative error estimations for Example 2 from some values of N.

N	${\mathit{E}}_{\mathbf{1}}$	${\mathit{E}}_{\mathbf{2}}$
5	$3.0586\times {10}^{-2}$	$2.5854\times {10}^{-2}$
10	$6.4473\times {10}^{-3}$	$3.7229\times {10}^{-3}$
20	$1.1610\times {10}^{-3}$	$6.3689\times {10}^{-4}$
30	$4.4048\times {10}^{-4}$	$2.2905\times {10}^{-4}$
40	$1.5159\times {10}^{-4}$	$1.1068\times {10}^{-4}$
50	$9.9079\times {10}^{-5}$	$6.3629\times {10}^{-5}$

shows the relative error estimations for distinct values of N.

Example 3.

We consider the following Volterra equation system of order 3

$$\left\{\begin{array}{c}{\displaystyle x_1\left(t\right)-{\int }_0^t(x_1\left(s\right)+t{x}_3\left(s\right))ds=-t+t^2,}\hfill \\ {\displaystyle x_2\left(t\right)-{\int }_0^t((t+s)x_1\left(s\right)+x_2\left(s\right)+(t-s)x_3\left(s\right))ds=1-t-\frac{t^4}{2},}\hfill \\ {\displaystyle x_3\left(t\right)-{\int }_0^t((-t-s)x_1\left(s\right)-x_2\left(s\right)+(-t+s)x_3\left(s\right))ds=-t-t^2+\frac{t^4}{2}.}\hfill \end{array}\right.$$

The exact solution is

$$x_1\left(t\right)=t^2,x_2\left(t\right)=1,x_3\left(t\right)=-t^2.$$

Table 4. Computed relative error estimations for Example 3 from some values of N.

N	${\mathit{E}}_{\mathbf{1}}$	${\mathit{E}}_{\mathbf{2}}$
5	$2.0024\times {10}^{-2}$	$3.5705\times {10}^{-2}$
10	$2.4457\times {10}^{-3}$	$5.2296\times {10}^{-3}$
20	$2.9878\times {10}^{-4}$	$7.5222\times {10}^{-4}$
30	$7.5462\times {10}^{-5}$	$2.6518\times {10}^{-4}$
40	$2.6133\times {10}^{-5}$	$1.1834\times {10}^{-4}$
50	$1.0932\times {10}^{-5}$	$7.1453\times {10}^{-5}$

shows the relative error estimations for distinct values of N.

8. Conclusions

We conclude that the above presented experiments (see Table 1, Table 2, Table 3, Table 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.