# An Effective Local Radial Basis Function Method for Solving the Delay Volterra Integral Equation of Nonvanishing and Vanishing Types

Authors: Neda Khaksari; Mahmoud Paripour; Nasrin Karamikabir
Journal: Journal of Mathematics
 (2022)

---

## Abstract

This paper presents a numerical method for solving a class of the delay Volterra integral equation of nonvanishing and vanishing types by applying the local radial basis function method. This method converts these types of integral equations into an easily solvable system of algebraic equations. To prove the method, we use the discrete collocation method and the local radial basis function method to approximate the delay Volterra integral equation. Also, we use the nonuniform Gauss–Legendre integration method to calculate the integral part appearing in the method. In addition, the existence, uniqueness, and convergence of the solution are investigated in this paper. Finally, some numerical examples are shown to observe the accuracy and effectiveness of the numerical method. Some problems have been plotted and compared with other methods. Obtained numerical results and their comparison with other methods show the reliability and accuracy of this method.


---

## Body

## 1. Introduction

The solution of the linear system of Volterra–Fredholm integral equations (VFIEs) has been a topic of a significant interest. Studies in the field of mathematical modeling, boundary value problems, and diverse biological and physical models led to superior explanations in the VFIEs. A survey of the formulation of such problems is shown in [1, 2] and the references therein. One of the main classes of mentioned areas is called as the delay Volterra integral equation (DVIE). Such a class has been found to be the best tool to model physical and chemical processes. DVIEs have been widely used for analysis and predictions in different areas of science such as biology, control, and electrodynamics [3–5]. These problems are applied in the mathematical modeling of different physical phenomena as dynamical systems in control theory, dynamics of multispan uniform continuous beams, heat conduction, signal processing, and the generalized voltage divider. In this paper, we are going to study the application of a locally radial basis function-discrete collocation (RBF-DC), a method with the Gauss–Legendre integration rule for obtaining a numerical solution of the following DVIE:(1)ut=gt+∫0tK1t,susds+∫0ptK2t,susds,t∈0,T,p∈0,1,u0=uT=0,where ut, K1t,s, and K2t,s on their domain are continuous functions. The DVIE explains a model of the particle motion in liquid and the crystallization of polymers in [6]. This model is a particular kind of the pantograph model and often emerges in several scientific models such as electrodynamics, number theory, and nonlinear dynamical systems. In recent years, due to such vital and important applications, there have been numerous studies in formulating some impressive numerical methods and their stability and convergence analysis. The convergence of the numerical solution of the delay Volterra integral equations has been provided by Zarebnia and Shiri [7]. In [8], the authors proposed a Bernoulli wavelet method for solving Volterra integral equations of the nonlinear fuzzy Hammerstein type including constant delay. Some papers studied the numerical solution of the DVIEs. In some cases, if DVFIEs do not have exact solutions, then we are interested in obtaining their approximate solution. Therefore, creating efficient numerical methods to obtain approximate solutions for the DVIEs is important. For example, in [9], the authors studied a geometric mesh and considered the corresponding collocation technique to solve the DVIE. Dastjerdi et al. [10] applied a moving least squares collocation method for DVIE. Ming et al. [3] studied a collocation method to solve the DVIE. In [11], a differential transform method was studied to solve the numerical solution of the DVIEs. In [12], proposed a sinc-collocation method for the solution of the DVIE. Some numerical methods have been used to approximate the solution of the DVIEs with different types of conditions, such as meshless methods [13], a local meshless method based on a radial basis function-finite difference [4, 5], a numerical method based on a radial basis function finite-difference [14], an efficient local meshless collocation algorithm [15], the Legendre spectral collocation method [16], an extrapolation method [17], a well-conditioned Jacobi spectral Galerkin method [18], and the Chebyshev spectral method [19]. The subject of radial basis functions was first studied by Hardy [20] and applied for topographical mapping. However, radial basis functions started to attract the attention of scientists, mathematicians, and researchers after their usage of them by Kansa [21] for finding approximate solutions to PDEs, and they are applied extensively to solve many problems coming from various fields of engineering and science.A function whose independent variable is a distance, or with a more precise definition, a function that takes the distance from the center pointxc to computing node x, as an independent variable is called the radial basis function. In the last few years, radial basis functions as an impressive meshless scheme have been used for interpolating a function over a set of scattered points [22–24]. Since the classical radial basis functions are global functions, the resultant coefficient matrix with respect to them will be full and large when many points in the domain are investigated for obtaining high-order accurate results. However, this coefficient matrix is nonsingular, but usually, it is ill-conditioned. To overcome this problem, the advantages of radial basis functions with local supports can be used which are called local radial basis functions studied in the papers [25–28]. To make the approximation function by local radial basis functions, the only geometrical data demanded are applied fallen within the local influence domain. Therefore, local radial basis functions require much less computational work than the globally supported local radial basis functions, and eventually, the interpolation matrix of local radial basis functions will be well conditioned. The main purpose of this paper is to demonstrate the efficiency and implementation of a local radial basis function method based on the discrete collocation method with a nonuniform Gauss–Legendre integration method on numerical solutions of the DVIEs. This paper shows a computational method using the discrete collocation method together with the local radial basis function approach as a locally meshless method. To calculate integrals shown in the proposed method, we demonstrate a particular integration rule based on the application of the composite nonuniform Gauss–Legendre quadrature formula. The error analysis and convergence analysis for our scheme are studied. To the best of our knowledge, there are no reported studies based on the local radial basis function method for the DVIEs. Unlike other numerical methods such as the finite element method and the finite difference method, in the local radial basis function method there is no requirement to make a mesh or grid, only a set of scattered nodes is needed instead of the mesh or grid. This characteristic of the local radial basis function method makes them desirable, and particularly in complex geometries, getting rid of creating grids is a very good achievement that shows great flexibility to local radial basis function methods. This method provides sparse differentiation matrices and better-conditioned interpolation matrices. Owing to better-conditioned interpolation matrices, it is possible to obtain more stable results and sparse differentiation matrices that reduce both computational cost and storage requirements. According to the description given, the implementation of our numerical method on computers is simple and computationally attractive. Also, our numerical method is only independent of the pairwise distances between points and could be easily extended to higher dimension problems.The remainder of the paper is organized as follows. In Section2, we prove the existence and uniqueness of the solution for equation (1). In Section 3, we obtain a numerical method for solving DVIEs based on the local radial basis function approximation in the collocation method. The convergence analysis of the solution is provided in Section 4. In Section 5, numerical results using the proposed method and discussion for some numerical examples are provided with tables and figures. The conclusion of the paper is reported in Section 6.


## 2. Existence and Uniqueness

In this section, we show and state the existence and uniqueness of the solution for equation (1). Let C0,T,ℝ be the Banach space. Then, we consider an operator Π:C0,T,ℝ⟶C0,T,ℝ by(2)Πut=∫0tK1t,susds+∫0ptK2t,susds,t∈0,T,p∈0,1.Ifut=C0,T,ℝ is a fixed point of the operator Π, then u is a solution of equation (1). Assume ut,u¯t∈C0,T,ℝ; then, we have(3)Πut−Πu¯t=supt∈0,T∫0tK1t,sus−u¯sds+∫0ptK2t,sus−u¯sds≤K1∗+K2∗u−u¯=Ku−u¯,where K1∗=supt∈0,T∫0tK1t,sds and K2∗=supt∈0,T∫0ptK2t,sds. If K<1, therefore the operator Π is a contraction and equation (1) has a unique fixed point which shows that there exists a unique solution of equation (1).


## 3. Explanation of the Proposed Method

In this section, we show the numerical method based on the local radial basis function collocation method for solving equation (1). Also, in this section, the locally supported RBFs so-called LRBFs are used to approximate a function in any dimension. A function Φ:ℝd⟶ℝ is called to be radial if there exists a univariate function ϕ:0,∞⟶ℝ such that(4)Φt=ϕr,in which r=t. Assume Θ⊆∪i=1NΘi, where Θ1,…,ΘN are an open bounded cover for Θ and Θi⊆Θ centered at ti with the radius ri. Then, a function ut can be developed as the following series:(5)ut≃Ωiut=∑j∈Iicjiϕit−tj,t∈Θi,where Iii=1N is the set of indexes corresponding to points fallen within the influence domain Θi with the cardinal number Ii=ni which is shown in Figure 1. The most favorite radial functions [29] are reported in Table 1. The Gaussian and inverse multiquadrics are strictly positive definite radial basis functions on any ℝd [30]. It guarantees that the associated interpolation matrix in the expansion (5) is nonsingular.Figure 1
One-dimensional case of the considered stencil.Table 1
The most popular radial functions.
Gaussian kernel  (GA)Multiquadric kernel (MQ)Inverse multiquadric kernel (IMQ)Φt=e−ct2Φt=c2+t2Φt=1/c2+t2In relation (5), the unknown coefficients cji for j∈Ii are calculated by enforcing the interpolation conditions.(6)Ωiutj=uj,j∈Ii.So, we can rewrite equation (6) as the following matrix form:(7)Biℂi=Ui,in which(8)Ui=u1,…,uniT,ℂi=c1i,…,cniiT,and Bi is the real symmetric coefficient matrix of size ni×ni and given by(9)Bi=ϕ0ϕt1−t2…ϕt1−tniϕt1−t2ϕ0…ϕt2−tni⋮⋮⋱⋮ϕt1−tniϕt2−tni…ϕ0.It should be noted that if a radial basis functionΦ is strictly positive definite, then the associated interpolation matrix Bi is positive definite and so nonsingular. Let us consider(10)Yit=ϕt−t1,…,ϕt−tniT.Then, equation (5) and the matrix Bi can be displayed by(11)Ωiut=YitTℂi,Bi=Yit1,…,Yitni.LetBi be an inverse matrix, then from equation (7), we have(12)ℂi=Bi−1Ui.By replacing equation (12) into equation (11), we obtain(13)Ωiut=YitTBi−1Ui=ΨitTUi.By substituting equation (13) into equation (5), we obtain(14)Ωiut=∑j∈IiujΨjit,where(15)Ψjit=∑k∈IiBikj−1ϕit−tk,j∈Ii.To discretize the integral in equation (1) by using the proposed method and applying a numerical integration method, we are required to estimate the solution of DVIEs for any point in the solution domain. Thus, we illustrate the local radial basis function collocation method to approximate a function ut at an arbitrary point t∈Θ by the following linear combination:(16)ut≃ℙNut=∑i=1NuiΨjiy,t∈Θ,where Ψjit is called shape functions for the local radial basis function interpolation. In general, the shape functions can be represented as follows:(17)Ψjitj=δij,tj∈Θi,0,tj∈Θi.Then, we can assume that(18)ℙNut=∑i∈ItuiΨiit,t∈Θ, in which It=i:t∈Θi and the influence domain of the point t are considered by Θt=∪i∈ItΘi. The choice of the parameter c, named as the shape parameter, influences heavily the accuracy of the method [29]. For example in Gaussians, if we select a bigger shape parameter with a fixed number N, the approximation will be more accurate, while the condition number of associated coefficient matrices will be large. On the other hand for fixed values of the shape parameter c, the condition number grows with N. However, the selection of the optimal value for c is still under check, and some values are suggested by many authors. Hardy has suggested in his original work on IMQ and MQ interpolations ℝ2 using c=0.815d, where d=∑i=1Ndi/N and di is the distance from ti to its nearest neighbor [31]. We consider the following DVI operator Γ:C0,T⟶C0,T:(19)Γut=∫0tK1t,susds+∫0ptK2t,susds.Then, we can rewrite equation (1) as the following operator identity:(20)I−Γu=g.To solve equation (1) by using the numerical method, we need N nodal points tii=1N considered in the interval 0,T. Thus, for any arbitrary function ut, the function ut can be estimated by applying the local radial basis functions as(21)ut≃∑j=1Nuj¯Ψjjt,t∈0,T.To calculate the coefficientsuj¯j=1N by using the collocation method, we substitute equation (21) and put the nodal points tii=1N into equation (1) and we obtain(22)∑j=1Nuj¯Ψjjti−∑j=1Nuj¯∫0tiK1ti,sΨjjsds+∫0ptiK2ti,sΨjjsds=gti,i=1,…,N.Because the support ofΨjjt is tj−rj,tj+rj for t∈0,T and the function Ψjjt satisfy the Kronecker delta condition, then we have(23)ui¯−∑j=1Nuj¯∫σjtiK1ti,sΨjjsds+∫σjptiK2ti,sΨjjsds=gti,i=1,…,N,σj=max0,tj−rj.AssumeℙN:C0,T⟶spanΨ11t,Ψ22t,…,ΨNNt⊂C0,T be a local collocation operator which is defined by(24)ℙNut=∑j=1NcjΨjjt,t∈0,T,where the coefficients cjj=1N are calculated by(25)ℙNutj=utj,j=1,…,N.Taking the operatorℙN of equation (1) yields(26)I−ℙNΓu=ℙNg.To calculate the integral in equation (23), we can use the μN-point quadrature rule which is given by(27)∫efytdt=Δt2∑i=1μNνi∑l=1Myξil+o1M2μN,M>0,yt∈C2μNe,f,where Δt=f−e/M, ξil=ϱi+2l−1/2Δt and ϱi∈−1,1. Therefore, we approximate the integrals given in equation (23) by using the μN-point quadrature rule as(28)∫σjtIK1ti,sΨjjsds≃Δsjti2∑l=1M∑i=1μNϱiK1ti,ξijlψjjξijlti,where ξijlti=ϱi+2l−1/2Δsjti+σj and Δsjti=ti−σj/M. Also, we obtain(29)∫σjptiK2ti,sΨjjsds=p∫σjtiK2ti,pzΨjjpzdz≃pΔsjti2∑l=1M∑i=1μNϱiK2ti,pξijlΨjjpξijlti,where ξijlti=ϱi+2l−1/2Δsjti+σj and Δsjti=ti−σj/M. Therefore, by substituting equations (28) and (29) into equation (23), we obtain a system of algebraic equations as follows:(30)ui¯−∑j=1Nuj¯Δsjti2∑l=1M∑i=1μNϱiK1ti,ξijlΨjjξijlti+pΔsjti0∑i=1M∑i=1μNϱiK2ti,ξijlΨjjpξijlti=gti,i=1,…,N.By solving these algebraic equations, the unknown constantsui¯ can be determined and the numerical solution given in equation (21) can be computed.Definition 1.
LetA=t1,…,tN be a set of points, such that A⊂D⊆0,T. Then, the fill distance of A is given by(31)hA,D=supmint∈Dt−tj:j=1,…,N.Theorem 1 (see [32]).
LetΦ be a positive definite radial basis function with infinite smoothness and u∈ℵΦa,b that, ℵΦa,b is defined by(32)ℵΦa,b=u∈Ca,b∩L2a,b:u^Φ^∈L2a,b,where u^ and Φ^ are the Fourier transform of u and Φ, respectively. Then, for any l∈ℕ and hA¯l≤h0¯l, there exist constants h0¯l and cl>0 such that(33)|ut−ℙNut≤clhA¯lmaxt∈a,buℵΦDt,where hA¯l=maxt∈a,bhA,Dt.With a similar process in [32], we can study the complexity of the proposed algorithm.


## 4. Convergence Analysis

In this section, we will check and study the convergence analysis of the DVIE based on the proposed method.Theorem 2.
Suppose thatut on 0,T be a differentiable function, such that u′t≤L. Then, over every subinterval In=n−1T/m,nT/m we have(34)u−un2=O1m3,where ut and un are the exact solution of equation (1) and its local radial basis series expansion, respectively.Proof.
By using the mean value theorem for integral, we obtain(35)u−un2=∫n−1T/mnT/mu−un2dt=Tmun−uη2,η∈In.
Also, for the functionunt, by applying the mean value theorem, we have(36)un=m∫n−1T/mnT/mutdt=mTmuη1,η1∈In.
By substituting equation (36) into equation (35) and applying the mean value theorem, we get(37)u−un2=Tmunη1−uη2≤TL2mη−η12≤T3L2m3.
If we consideret=∑n=1mu−un, we obtain(38)et22=∫0T∑n=1mu−un2dt=∫0T∑n=1mu−un2dt+2∑n<j∫0Tu−unu−ujdt.
Since forn≠j we have In∩Ij=∅, then,(39)et22=∑n=1m∫0Tu−un2dt=∑n=1mu−un2.
By substituting equation (37) into equation (39), we obtain(40)et22≤T3L2m2.
To prove the convergence analysis of equation (1), we consider the error function et=u−um, that u and um are the exact solution and approximate solution of equation (1). Then, we have(41)ut=gt+∫0tK1t,susds+∫0ptK2t,susds,(42)umt=gt+∫0tK1t,sumsds+∫0ptK2t,sumsds+rt,where rt is the remainder term. By subtracting equation (42) from equation (41), we obtain the following relation:(43)rt≤∫0tK1t,sus−umsds+∫0ptK2t,sus−umsds≤u−um∫0tK1t,sds+∫0ptK2t,sds.
Whenm⟶∞, from equation (40), we have u−um⟶0, then rt⟶0. Therefore, we show the convergence analysis of equation (1).


## 5. Illustrative Examples

This part shows some numerical examples for demonstrating the efficiency and accuracy of the proposed method. Also, a comparison of the proposed method with previously known methods is shown. For numerical examples in this work, we use the Gaussian (GA), inverse multiquadric (IMQ), and multiquadric (MQ) in which the applicability and efficiency of them heavily depend on the shape parametersc>0. Since local radial basis functions have much more freedom in selecting the shape parameter, we choose middle values c=0.1×N for GAs and c=5/N for IMQs and MQs. Also, in numerical examples, we use the linear l=1 and quadratic l=2 basis functions and the Gaussian weight functions for approximating integrals in the method. For assessing the performance and efficiency of the proposed method, we define the maximum of the absolute error and the mean error on a domain D⊂ℝd as follows:(44)eN22=∫Dut−uNt2dt,eN∞=maxt∈Dut−uNt,where ut is the exact solution and uNt is the numerical solution of ut obtained by using the proposed method. Also, the convergence rate of the presented method has been defined by(45)Ratio=lneN∞−lneN′∞lnhA′¯l−lnhA¯l.Example 1.
Consider the following DVIE:(46)ut=sint+cost+e−tcospt+1−1+∫0tusds+∫0pte−tusds,t∈0,π,p∈0,1,u0=uπ=0.
The exact solution of equation (45) is ut=sint. Numerical solutions by using the previous numerical methods have difficulties in solving these kinds of DVIEs, but we can easily calculate the numerical solution for this example by applying GAs and MQs displayed in this paper based on some random nodes over the interval 0,π. The distribution of nodes selected randomly on the interval 0,π is depicted for some N in Figure 2. We first use the proposed method for solving this example by choosing N=10,25,50, and p=0.5. Numerical solution and exact solution graphs by using the technique described in Section 3 for various values of N are shown in Figure 3. Figure 4 shows the graphs of the absolute error function with different values of N. The absolute error function with different values of N for this problem using GA (left) and MQ (right) is shown in Figure 4. Figure 5 demonstrates the obtained errors by using the proposed method for various values of N. Table 2 compares the absolute error function of the numerical solutions calculated by applying the proposed method when N=10,25,50, and p=0.5.Figure 2
Node distribution of numerical for Example1 with different values of N.Figure 3
The graphs of the numerical solutions by using the proposed method for Example1 with different values of N and the exact solution.Figure 4
The absolute error function graph by using the proposed method for Example1 with different values of N.Figure 5
The obtained errors by using the proposed method for Example1 with different values of N.Table 2
The absolute errors of the numerical solutions by using the proposed method for Example1.
NtOur methodGAMQceN2RatioceN2Ratio100.40.444.30527e−10-1.114.30527e−09—250.60.632.7551e−112.230.792.7551e−101.48500.80.703.2300e−1310.260.703.2300e−127.25Example 2.
Consider the DVFIE:(47)ut=1+e−pt2−∫0tusds+12∫0ptusds,t∈0,10,p∈0,1,u0≠0,u10≠0.
For this example, the exact solution isut=e−t. We solved this example by applying our suggested technique demonstrated in Section 3 by selecting different values of N and p=0.2. For this example, applying GAs and IMQs displayed in this paper based on some random nodes over the interval 0,10. Figure 6 represents the numerical solutions from the proposed method and the exact solution for different values of N. Figure 7 shows the absolute error in 0,10 by selecting different values of N. The absolute error function with different values of N for this problem using GA (left) and IMQ (right) is shown in Figure 7. Table 3 represents the maximum absolute error results of the studied method with various values of N. Figure 8 demonstrates the obtained errors by using the proposed method for various values of N. A comparison of our suggested technique and the numerical results shown in [3, 12] for different values of N is given in Table 4. The obtained numerical results by using the method presented in the current paper are better than the numerical results shown in [3, 12]. We have also compared the CPU times by applying our suggested technique for different numbers of N in Figure 9. The convergence rates of the proposed method are much higher than the convergence rates shown in [3, 12]. Also, the CPU times of the proposed method are lower than the numerical results shown in [3, 12].Figure 6
The graphs of the numerical solutions by using the proposed method for Example2 with different values of N and the exact solution.Figure 7
The absolute error function graph by using the proposed method for Example2 with different values of N.Table 3
The absolute errors of the numerical solutions by using the proposed method for Example2.
NtOur methodGAIMQceN2RatioceN2Ratio102.90.441.00e−16—1.111.00e−15—253.90.6307.720.7905.76504.80.70012.340.70010.74Figure 8
The obtained errors by using the proposed method for Example2 with different values of N.Table 4
Comparison of the obtained max erroreN∞ value for Example 2 using our method and those presented in the articles [3, 12].
Our methodMethod in  [12]Method in  [3]N=10N=25N=50N=10N=20N=100N=2002.4499e−088.9768e−105.9271e−131.18e−64.48e−81.67e−61.06e−7Figure 9
CPU times for Example2 with different values of N.


## 6. Conclusion

This paper used the local radial basis function method to solve a class of the delay Volterra integral equation of nonvanishing and vanishing types with given initial and boundary conditions. The method used locally supported radial basis functions made on scattered points as a basis in the collocation method. By using the nonuniform Gauss–Legendre quadrature, we gave a numerical formula for obtaining integrals in the method. This method has reduced the solution of the delay Volterra integral equation to the solution of a linear system of algebraic equations. This method does not need any background approximation cells. Some numerical examples are shown to demonstrate that the proposed method is very accurate and effective. The results obtained are compared with the errors presented by other methods published in the scientific literature in this regard. The numerical results of our proposed method are more accurate than those presented in the articles indicated in Table4 and the References.



---

*Source: 1527399-2022-06-22.xml*