RBF Based Localized Method for Solving Nonlinear Partial Integro-Differential Equations

In this work, a numerical scheme is constructed for solving nonlinear parabolictype partial-integro differential equations. The proposed numerical scheme is based on radial basis functions which are local in nature like finite difference numerical schemes. The radial basis functions are used to approximate the derivatives involved and the integral is approximated by equal width integration rule. The resultant differentiation matrices are sparse in nature. After spatial approximation using RBF the partial integro-differential equations reduce to the system of ODEs. Then ODEs system can be solved by various types of ODE solvers. The proposed numerical scheme is tested and compared with other methods available in literature for different test problems. The stability and convergence of the present numerical scheme are discussed.


Introduction
In the fields of sciences and engineering, physical systems which depends on time and space can be formulated by partial differential equations. Some time the physical system may not be accurately model by such formulations at a given time by ignoring effect of the system in some past times. Specifically in the fields of nuclear reactors, thermoplastic and heat flow there is need to include the memory effect to the system. This addition always occurs as an integral term in basic differential equation resulting a partial integro-differential equation. Various types of partial integro-differential equations available in the literature using for different physical systems. In the present work, we study the parabolic-type Volterra partial integro-differential equations [Avazzadeh, Rizi, Ghaini et al. (2012)]. These type of equations have application in reaction-diffusion problems [Engler (1983)], compression of viscoelastic media [Visser (1997)] and nuclear reactor dynamics [Pachpatte (1983)]. The analytical solution can be available in few cases, so researchers have introduced different techniques for the approximating solutions of PIDEs. Some of these efficient numerical methods include the finite difference method [Tang (1993)], the Quintic B-spline collocation method [Zhang, Han and Yang (2013)], the Variational iteration method [Nawaz (2011)], the Quasi-wavelet with numerical method [Yang, Xu and Zhang (2011)], the spectral method [Fakhar-Izadi and Dehghan (2011)], the Method of lines [Kauthen (1992)], the Finite Element Methods [Yanik and Fairweather (1988)], the orthogonal spline collocation method [Yan and Fairweather (1992)], the Galerkin methods [Zhang, Lin, Lin et al. (2001)] and the θ-weighted with Radial Basis Functions [Avazzadeh, Rizi, Ghaini et al. (2012)]. Hardy [Hardy (1971)] introduced a numerical method using multiquadrics (MQ) radial basis functions. In year 1982 [Franke (1982)] a comparison had been made among different numerical methods and the Hardy's MQ radial basis functions method outclass all the other methods regarding accuracy, stability and ease of implementations. Tarwater [Tarwater (1985)] demonstrated the effect of shape parameter in MQ on the solution accuracy. In the work [Carlson and Foley (1991)] it has suggested that a small value of the shape parameter be used if the function varies rapidly, but a large value be used if the function has large curvature. In 1990 Kansa extended the method of Hardy's for solving numerically various types differential equations vary efficiently [Kansa and Multiquadrics (1990)]. The convergence analysis have been discussed by many authors (see for example [Micchelli (1984); Madych and Nelson (1990); Franke and Schaback (1998)]). The most advantage of utilizing the RBF for the approximation of PDEs is its effortlessness, pertinence to various PDEs, and viability in managing with multi-dimensional issues and complicated domains. In most of the cases the global method returned the differentiation matrices asymmetric and dense, which needed large amount of data and computations time. To overcome these issues the authors in Tolstykh et al. [Tolstykh (2000); Shu, Ding and Yeo (2003); Vertnik and Šarler (2006)] developed a local meshless technique where the kernel based interpolant in small sub-domains centered around each center is used for system matrices. This idea has been extended to develop different sorts of very efficient numerical methods and has been used effectively to a wide range of unsolved issues (see for example [Sarra (2012), Yao, Šarler and Chen (2011)]). Various other modification of local RBF methods and their application can be found in the literature for example local RBF method for Darcy flow by Kosec et al. [Kosec and Šarler (2008)], local RBF-based differential quadrature method for incompressible Navier-Stokes equations by Shu et al. [Shu, Ding and Yeo (2005)], H-adaptive local radial basis function collocation meshless method by Kosec et al. [Kosec and Šarler (2011)], the meshless local Petrov Galerkin (MLPG) method by Atluri et al. [Atluri and Shen (2005)], Stable calculation of Gaussian-based RBF-FD stencils by Fornberg et al. [Fornberg, Lehto and Powell (2013)], Scattered node compact finite difference-type formulas generated from radial basis functions by Wright et al. [Wright and Fornberg (2006)]. In this work, we use such type approach to solve the nonlinear parabolic-type Volterra partial integro-differential equations [Avazzadeh, Rizi, Ghaini et al. (2012)].

Development of the proposed method
The nonlinear parabolic-type Volterra partial integro-differential equation is of the type u t ðx; tÞ ¼ Luðx; tÞ þ Z t 0 jðx; t; s; uðx; sÞÞds þ f ðx; tÞ; (1) where x 2 & R d ; d ! 1, and t∈[0, T], and subject to the following initial and and the boundary conditions respectively Buðx; tÞ ¼ gðx; tÞ; x 2 @; where L and B are spatial operators. It is assume that the functions κ and f are continuous on {(x,t): x∈Ω, 0≤t≤T}. Let N be the number of points in the set {x i , i=1, …, N}∈Ω and Ω k be a sub-domain in the domain Ω contains m points for each center x k in Ω. The unknown function u is approximated by the linear combination of radial basis functions in each sub-domain Ω k , k=1, …, N (see for example ; Uddin, Minullah, Ali et al. (2015)]) is given by where r j =‖x − x j ‖, and x, x j ∈Ω k is the Euclidean norm between two centers x and x j and φ is a radial basis function defined for r≥0 and λ j are the expansion coefficients. From Eq. (4), N number of m×m systems of linear equations in the matrix form are given by and we can also approximate Luðx; tÞ by which can be represented by the matrix-vector form By eliminating the expansion coefficients λ k from Eqs. (6)-(8), we get where W k ¼ A k L ðA k Þ À1 , denote the corresponding weight for the kth center, hence for all centers k=1, …, N, we have from Eq. (9) where W is a N×N sparse differentiation matrix. Using these kernel based approximation of function u and Lu in the model Eq.
(1), we get the system of ODEs of Eq. (1), Advancing the solution in time, we write the solution of Eq. (11) in the form where U n , F n and K n denote the values at (x i ,t n =n δt, i=1, …, N), with step size δt and the ith element of K n can be computed by using the numerical trapezoidal rule as jðx i ; t n ; s q ; uðx i ; s q ÞÞ þ jðx i ; t n ; s n ; uðx i ; s n ÞÞg: (13) This is the required scheme for obtaining the numerical solution at any time level n. Initially we take U 0 =h 0 (x) from the given initial condition u(x,0)=h 0 (x). In the next section we will discuss the stability and convergence of the this scheme.
3 Stability and convergence of the scheme The scheme in Eq. (12) is a recurrence relation that allows us to advance the solution in time from t n =nδt to t n+1 =(n+1)δt. The value of the amplification matrix B=1+δtW depend on the ratio δt/h r , here r is the order of largest space derivative and h denote the distance between two nodes. Suppose u n be the exact solution of Eq.
(1) at time t n =n δt, then it follows that | D N u(x)−D N U(x)|≤Ch N |u| see [Fasshauer (2007); Uddin and Haq (2011)]. Further assume the numerical scheme in Eq. (12) is of order p in space, so we get We define the error at nth time level by ε n =u n −U n , By subtracting Eq. (12) from Eq. (14) we get By Lax-Richtmyer definition of stability the scheme Eq. (15) is stable if The result B≤ρ(B) is always true and B=ρ(B) if B is normal. Let us assume that the initail condition and the solution of the given integro-differential equation must be sufficiently smoothand h to be enough small. For keeping the values of δt/h r to be constant we must have δt → 0. so there exist a constant C such that ke nþ1 k kBke n k þ CððdtÞ 2 þ h p Þ; n¼0; 1; 2; …; T : Since initailly at n=0, U 0 =u 0 therefore ε 0 =0 and so by mathematical induction, we have ke nþ1 k ð1 þ kBk 2 þ … þ kBk nÀ1 ÞCððdtÞ 2 þ h p Þ; n¼0; 1; 2; …; T: By using Eq. (16), we have ke nþ1 k nCððdtÞ 2 þ h p Þ; n¼0; 1; 2; …; T: This shows the scheme is convergent.

Application of the method to problems
In this section, we apply the present numerical scheme to various problems 1-D and 2-D of the type given in Eq.
Problem 4.5 We consider an other set of function for approximate the problem in Eq. (20) by the current numerical scheme with the initial condition and boundary conditions given by uðx; 0Þ ¼ sin x; x 2 ½0; 1; uð0; tÞ ¼ 0; uð1; tÞ ¼ e Àt sinð1Þ; 0 t T ; where α(x, t)=x 5 +4x 2 and κ(x,s,t,u)=(x 3 t 3 +1) s 2 u 2 and f(x,t) can be found from the exact solution u(x, t)=e −t sin x, given in Avazzadeh et al. [Avazzadeh, Rizi, Ghaini et al. (2012)]. The present localized kernel based method is used and the solution is advanced in time with the time step δt. The results are shown in form of error norms L 2 and L ∞ respectively and are shown in Tabs. 6-7 and Fig. 5. The proposed method performed more accurate results than the method in Avazzadeh et al. [Avazzadeh, Rizi, Ghaini et al. (2012)]. Problem 4.6 Here we applied the local kernel based method for solving Eq. (20) having the exact solution given by u(x, t)=x cos(t) with initial condition u(x,0)=x, x∈[0,1], and with the boundary conditions u(0,t)=0, u(1,t)=cos(t), 0≤t≤T. The functions α(x, t)=cos 2 x+sin(x) and κ(x, s, t, u)=e x+s+t u 2 are used and the function f(x,t) can be found from the exact solution.
The error norms L 2 and L ∞ are shown in Tab. 8 and Fig. 6.

Conclusion
In this work, a local kernel based numerical scheme is developed for solving the nonlinear Volterra type integro-differential equations. The numerical scheme is local in nature like finite differences numerical scheme and the system matrices are small and well conditioned. The numerical scheme is stable for a large range of RBF shape parameters. Convergence and stability of the proposed numerical scheme have been established. A number of 1D and 2D integral equations have been solved to validate the present numerical scheme. The efficiency and applicability of the proposed method is well demonstrated and compared with some available results. The developed numerical scheme is an alternative for solving such types of models effetely and accurately.