A partial integro-differential equation (PIDE) is an equation composed of a partial-differential term and an integral term. In the recent past, the solution of partial-integro differential equations has attracted attention and motivated research in the field in view of applications in mechanics, biology and finance     . We notice that in the past, research on integro-differential problems has focused on one-dimensional problems in the framework of ordinary differential equations. On the other hand, parabolic multi-dimensional problems with Volterra type integral terms have been considered  . Furthermore, independently of these topics, the problem of fast computation of Fredholm operators in multi-dimensions has been investigated. However, much less is known on the numerical analysis of multi-dimensional elliptic PIDEs with Fredholm integral terms. In  , a one dimensional PIDE with a convolution kernel is solved through conversion of the PIDE to an ordinary differential equation and the use of the inverse Laplace transform. The work in  develops a moving mesh finite-difference method for a PIDE that involves approximating the time dependent mapping of the coordinate transformation by a piecewise quadratic polynomial in space and piecewise linear functions in time. In   , compact finite-differences for one-dimensional PIDEs are studied. Additional results on high-order schemes for integro differential equations (IDE) can be found in  . The research in  , is devoted to an iterated Galerkin method for PIDE in one-dimesion; see  . Further, the work  considers the numerical solution of linear IDE using projection methods. The work in  investigates a Tau method with Chebychev and Legendre basis to find the numerical solutions of Fredholm integro-differential equations where the differential part is replaced by its operational Tau representation. We remark that the methodologies referred above are designed for one-dimensional problems and their complexity for multi-dimensional problems may become prohibitive.
The purpose of this work is to contribute to this field of research with the development and analysis of a methodology that is appropriate for multi-dimen- sional PIDE problems. We present a second-order accurate fast multigrid scheme to solve elliptic problems of the following form
where A represents an elliptic operator with given boundary conditions and . Our approach is to combine a multigrid scheme for elliptic problems with the multigrid kernel approximation strategy developed in  .
For this purpose, we discretize our PIDE problem by finite-differences and quadrature rules and analyse the stability and accuracy of the resulting scheme in the case of A being the minus Laplace operator that is combined with a Fredholm Hilbert-Schmidt integral operator.
It is well-known that a multigrid scheme solves elliptic problems with optimal computational complexity. However, this is in general not true if a straight- forward implementation of the integral term is considered. For this reason, with- in the multigrid framework, we investigate the multigrid kernel approximation stra- tegy proposed in  , where it is demonstrated that it is possible to approximate a Fredholm integral term with accuracy while reducing the complexity of its calculation from to , where n is the number of grid points.
Our work is organized as follows. In Section 2, we discuss the theory of an elliptic Fredholm partial integro-differential equation, proving existence and uni- queness of solutions. In Section 3, we discuss the finite-difference discretization of our PIDE problem and prove second-order accuracy of the numerical solution. In Section 4, we illustrate our multigrid solution process including the techniques in  for approximating the Fredholm integral operator. Our resulting PIDE multigrid solution procedure is analysed by local Fourier analysis in Section 5. In Section 6, results of numerical experiments are presented that successfully validate the theoretical estimates and the effectiveness of the proposed PIDE solution procedure. A section on conclusion completes this work.
2. An Elliptic Fredholm Partial Integro-Differential Equation
We consider the following PIDE problem
where is a two-dimensional, convex and bounded domain with a boundary or a rectangle. We denote and . We consider and a symmetric positive semi-definite Hilbert Schmidt kernel
, such that , and the following holds
We have the following theorem.
Theorem 2.1 Let be a Hilbert Schmidt kernel. The integral operator given by
defines a bounded mapping of into itself, with the Hilbert Schmidt norm .
Proof. From Tonelli’s theorem, is a measurable function of x and its L2-norm can be determined by the Cauchy-Schwarz inequality. Let , we have
Remark 2.1 From Schur’s test  , since the kernel k is a measurable function, it satisfies the following conditions
Then the integral operator defines a bounded mapping and .
With this preparation, we can prove the following.
Theorem 2.2 There exist a unique function that solves (2.1)-(2.2).
Proof. The proof is straightforward by using the Lax-Milgram theorem and the properties of the kernel.
3. Discretization of the Elliptic PIDE Problem
We discretize (2.1)-(2.2) using finite differences and the Simpson’s rule    . For simplicity, we assume that and such that we can evaluate these functions on grid points. Specifically, we consider and N is an integer with . We denote and
. Let be the mesh size. We denote the mesh points
#Math_37#, and . These grid points define the following grid
Later, we consider a sequence of nested uniform grids , where for .
For grid functions v and w defined on , we introduce the discrete L2-scalar product
with associated norm . The negative Laplacian with homogeneous Dirichlet boundary conditions is approximated by the five-point stencil and is denoted by . Given continuous functions in are approximated by grid functions defined through their values at the grid points. Thus the right-hand side of (2.1) in is represented by , if (otherwise by local average), and similarly for the kernel function.
Further, we introduce the following finite-difference operators. The forward finite-difference operator is given by
The backward finite-difference operator is as follows
With these operators, we can define the norm as . Notice that the bracket denotes summation up to N in the given direction , resp. ; see  . With this preparation, we have
The integral term of the elliptic PIDE in two-dimensions is written explicitly as follows
Using the Simpson’s rule, we have the following approximation of this integral operator
where represents the coefficients of the quadrature rule. In the case of the Simpson’s rule, we have
We refer to the Formula (6) as the full-kernel (FK) evaluation.
We need the following lemma.
Lemma 3.1 The positivity of the Hilbert Schmidt operator stated in (2.3) is preserved after discretization.
Proof. Consider the following function where is a suitable approximation of the Dirac delta function as , e.g., a narrow Gaussian. Inserting this function in (2.3) we have
Therefore by continuity, as the above integral tends to . Thus, we obtain .
The Simpson’s rule provides a fourth-order accurate approximation of the integral as follows
for any sufficiently smooth .
With the setting above, we write the finite-difference approximation of (2.1)- (2.2) as follows
where denotes the numerical approximation to u. Further, the integral function (6) evaluated at is given as follows
Notice that for functions that are zero on the boundary, summation can be restricted to the interior grid points. The double superscript in indi- cates that the pairs of indices refer to the mesh .
Next, we investigate the stability and accuracy of (8). For this purpose, we use the numerical analysis framework in  . We denote .
We need the following lemma, see also  .
Lemma 3.2 Suppose is a function defined on with on the boundary; then the following holds
Proof. Using the results of Lemma 3.3, we have
Lemma 3.3 Suppose U is a function defined on with on the boun- dary; then there exists a constant , which is independent of and , such that the following discrete Poincaré-Friedrichs inequality holds
for all such ; see  .
Remark 3.1 From (3.6) and (3.7), we obtain , where .
Theorem 3.4 The scheme (3.5) is stable in the sense that
Proof. We have
We conclude this section with the following theorem.
Theorem 3.5 Suppose and is a Hilbert Schmidt ker- nel, and assume that the weak solution to (2.1)-(2.2) belongs to ; then the solution to (3.5) approximates with second-order accuracy as follows
where c is a positive constant independent of h. In particular .
Proof. The proof uses Theorem 3.4 and the fact that the truncation error of (3.5) is of second order. This proof follows exactly the same reasoning as in Theorem 2.26 in  .
4. A Multigrid Scheme for Elliptic PIDE Problems
Our multigrid solution procedure for solving the discrete elliptic PIDE problem (3.5) is based on the full approximation storage (FAS) framework    and the multigrid fast integration technique presented in  . Notice that, although in this work we consider linear problems such that a linear multigrid scheme is well suited, our focus is on the nonlinear FAS framework in view of future applications (nonlinear problems, differential inequalities).
To illustrate our multigrid strategy, we first focus on the two-grid case, which involves the fine grid and the coarse grid , where .
In , consider the discretized PIDE Equation (3.5) as follows
where denotes the solution to this linear problem.
The main idea of any multigrid strategy for solving (4.1) is to combine a basic iterative method that is efficient in reducing short-wavelength errors of the approximate solution to (4.1), with a coarse-grid correction of the fine-grid long- wavelength solution’s errors that is obtained solving a coarse problem.
We denote the smoothing scheme with S. Specifically, when S is applied to (4.1), with a starting approximation , it results in . The smoothing property is such that the solution error has smaller higher-frequency modes than the error .
In the multigrid solution process, starting with an initial approximation and applying S to (4.1) -times, we obtain the approximate solution .
Now, the desired (smooth) correction to , to obtain the exact solution, is defined by . Equivalently, this correction can be defined as the solution to
where is the residual associated to .
Next, notice that the structure of and the smoothness of the error function allow to represent (4.2) on the coarse grid . On this grid, is repre- sented in terms of coarse variables as follows
where represents the restriction of to the coarse grid by means of the direct injection operator denoted with .
With this preparation, it appears natural to approximate (4.2) on the coarse grid as follows
Notice that this equation can be re-written as where . The term is the so-called fine-to-coarse defect correction.
Now, suppose to solve (4.4) to obtain . Then we can compute , which represents the coarse-grid approximation to . Notice that while and need not to be smooth, their difference is expected to be smooth by construction, and therefore it can be accurately interpolated on the fine grid to obtain an approximation to that is used to correct . This procedure defines the following coarse-grid correction step
In order to damp the high-frequency errors that may arise through the coarse- grid correction, a post-smoothing is applied.
Notice that in (4.4) a restriction operator is applied to the residual, while in (4.5) an interpolation operator is applied to the error function. We choose the following inter-grid transfer operators  
The transfer operators given by (4.6) are of 4th-order.The 4th-order interpo- lation operator is symmetric and accesses two grid points on either side of the interpolated point. To approximate at the grid point and , an asymmetric interpolation that accesses one point on one side and three points on the other side is used. The asymmetric fourth-order interpolation is given by
Notice that our choice of higher-order interpolation and restriction operators, as in  , appears advantageous for the fast integration technique that we discuss next.
The fast integration strategy   aims at performing integration mostly on coarser grids and to interpolate the resulting integral function to the original fine grid where this function is required. To illustrate this technique in the one- dimensional case, denote with the grid points on the grid with mesh size h, and with the grid points on the grid with mesh size . Notice that for .
Now, suppose that the kernel and are sufficiently smooth. (For the case of singular kernels see  .) In , the integral
is approximated by . On the
other hand, in the strategy of  , the kernel is approximated by , where the interpolation operator may be equal to . With this setting, we have
where is obtained by coarsening of . In particular, using straight injection for boundary values or the full-weighted restriction in (4.6), we have . Now, we go a step further and consider the coarse integral function
. This function is evaluated on the coarse grid and, from
the calculation above it is clear that it is equal to for all . Therefore we obtain the following approximation to the integral function on the fine grid
In one dimension, the summation complexity on the coarse grid is of order operations, which may still be large. However, assuming that the kernel is sufficiently smooth and using the fact that the coarse-grid summation has the same structure of the fine-grid summation, the coarsening-summation proce- dure just described can be applied recursively, until a grid is reached with grid points. On this grid the summation is then actually performed, requiring operations. Further, the computational effort of the restriction and of the interpolation is , where is the order of interpolation ( for linear interpolation). Therefore the order of total work required to obtain the (approximated) summation is operations. Notice that choosing , the quadrature error on using Simpson’s rule is . This gives an overall order of accuracy for the whole procedure.
Now, we can illustrate the multigrid procedure considering a sequence of nested grids (levels) of mesh size , indexed by the level number , where denotes the finest level. First, we summarize the fast inte- gration (FI) technique in Algorithm 1, where we perform full-kernel evaluation when a level , with given depth , is reached.
Next, we discuss our smoothing scheme. Our approach is to implement a Gauss-Seidel step for the Laplace operator, without updating the integral part of the equation operator. It can be appropriately called a Gauss-Seidel-Picard iteration, where the integral is evaluated using the FI scheme before the Gauss-Seidel step starts. In the one-dimensional case, our smoothing scheme is given by Algorithm 2.
Our multigrid scheme is given in Algorithm 3. Notice that this algorithm describes one cycle of the multigrid procedure that is repeated many times until a convergence criterion is satisfied. In Algorithm 3, the parameter is called the cycle index and it is the number of times the same multigrid procedure is applied to the coarse level. A V-cycle occurs when and a W-cycle results when .
5. Local Fourier Analysis
In this section, we investigate convergence of the two-grid version of our FAS-FI multigrid solution procedure using local Fourier analysis (LFA)     . In order to ease notation, we consider a one-dimensional case and use h and H indices to denote variables on the fine and coarse grid, respectively. For the LFA investigation, we assume that the kernel of the integral term is translational invariant in the sense that and require that de- cays rapidly to zero as becomes large. With these assumptions the stencil of our PIDE operator can be cast in the standard LFA framework. However, trea- ting the fast-kernel evaluation in this framework results too cumbersome. On the other hand, numerical experiments show that, apart of the different complexity, the convergence of our multigrid scheme with FI and with FK evaluation are very similar. Therefore we analyze our two-grid scheme with the latter procedure.
We apply the local Fourier analysis to the two-grid operator given by
where different pre- and post-smoothing steps are considered. The coarse grid operator is given by .
The local Fourier analysis considers infinite grids, , and there- fore the influence of boundary conditions is not taken into account. Nevertheless, LFA is able to provide sharp estimates of multigrid convergence factors. This analysis is based on the function basis
For any low frequency , we consider the high frequency mode given by
We have for and . We also have on for and .
The two components and are called harmonics. For a given , the two-dimensional space of harmonics is defined by
For each and a translational invariant kernel, we assume that the space is invariant under the action of . In fact in our case, the stencil of the discrete PIDE operator is defined by constant coefficients that do not depend on the choice of origin of the infinite grid. Now, we study the action of on the following function
Specifically, we determine how the coefficients , and , are transformed under the action of the two-grid operator. This requires to calculate the Fourier symbols of the components that enter in the construction of this operator.
First, we derive the Fourier symbol of our smoothing operator. For this purpose, we introduce the following Fourier representation of the solution error before and after one smoothing step. We drop the index as we assume invariance of under the action of the smoothing operator. We have
where and denote the error amplitude after and iterations of the smoother. Notice that , and thus for the coefficient of the Fourier mode, we have
Now, consider the following point-wise definition of our GSP iteration applied to our discretized PIDE on . We have
At this point, recall that and introduce the index . Since we assume that the element becomes very small as becomes large, we truncate the sum in (5.4) and consider the following equation
where we assume that the partial sum provides a sufficiently accurate approximation of the integral term on . Specifically, assuming that and requiring that the (double pre- cision machine epsilon) for , one should choose . However, in practice, a much smaller results in accurate LFA estimates.
Next, in (5.5), we insert (5.3) and obtain
which can be re-written as follows
Now, comparing the coefficients of equal frequency modes on both sides of (23), we obtain
Therefore an appropriate estimate of the smoothing factor of our GSP scheme is given by
With this definition, we obtain a smoothing factor of our GSP scheme given by for all mesh sizes. This is done by inspection of the function .
Next, in order to investigate the two-grid convergence factor, we construct the Fourier symbol of the two-grid operator. For this purpose, we derive the Fourier symbol for , applied to a generic vector with jth component given by . We have
Therefore we obtain
Now, recall that on the fine grid, we distinguish on the two harmonics. There- fore we have the following operator symbols acting on the vector of the two harmonics
On the coarse grid, we have the following
For the restriction operator, we have the following 
For the interpolation operator, we obtain
Now, we are able to compute the two-grid convergence factor as follows
where denotes the spectral radius of the matrix .
In Table 1, we report the values of the two-grid converge factor given by (5.12) for different numbers of pre- and post-smoothing steps, . These values are computed by inspection of the function , which is evaluated using MATLAB to compute the eigenvalues of the matrix on a fine grid of values, .
Further in the same table, we compare these values with the value of the observed convergence factor given by
Table 1. Estimated and observed multigrid convergence factors.
This numerical convergence factor represents the asymptotic ratio of reduction of the -norm of the residual between two multigrid cycles. These calculations refer to the choice . As shown in Table 1, the LFA estimates of the multigrid convergence factor are accurate. (The same values of are obtained with L ranging from 20 to 400.)
6. Numerical Experiments
In this section, we present results of numerical experiments to validate our FAS-FI multigrid strategy and the theoretical estimates. We demonstrate that our FAS-FI scheme has computational complexity (M denotes the total number of grid points on the finest grid) and provides second-order accurate solutions.
Our first purpose is to validate our accuracy estimates for the discretization scheme used. For this purpose, we consider an elliptic PIDE problem with a Gaus- sian convolution kernel in two dimensions as follows
where and the Gaussian kernel
To investigate the order of accuracy of the discretization scheme, we construct
an exact solution to (6.1) by choosing . With this
choice, the right-hand side of (6.1) is given by
The Dirichlet boundary is also given by the chosen u.
Using the exact solution above, we can validate the accuracy of our finite- differences and Simpson’s quadrature schemes. In Table 2, we report the values of the norm of the solution errors on different grids. We obtain second-order accuracy as predicted.
Next, we investigate the FI scheme. For this purpose, we consider the integral term in (3.3), and compute the norm . Notice that we can evaluate exactly, while is computed using the full-kernel (FK) evaluation formula (6) and the FI technique involving different depths, . For the FI scheme performs FK evaluation.
In Table 3, we report the values of corresponding to different working levels and different depths. Because we use a fourth-order quadrature formula, we can see an increase of accuracy of a factor 16 by halving the mesh
Table 2. L2−norm error using full kernel approximation.
Table 3. Errors for 2D integral evaluation using 4th−order interpolation and different depths.
size and using the FK scheme. On the other hand, increasing the depth of the FI scheme, this scaling factor deteriorates. However, since the truncation error corre- sponding to the Laplace operator is of second-order, the reduction of accuracy due to the use of the FI scheme with the fourth-order quadrature does not affect the overall solution accuracy of the PIDE problem as shown in Table 4.
Next, we validate our FAS-FI solution procedure. One main issue is how the accuracy of the solution obtained with the FAS-FI scheme is affected by the approximation of the integral due to the FI procedure. For this purpose, in Table 4, we compare the norm of the solution errors obtained with a FAS scheme with FK calculation and with our FAS scheme including the FI technique. We see a moderate degradation of the quality of the numerical solution while increasing the depth. On the other hand, we notice that a second-accurate solution is obtained by choosing d corresponding to the first before the coarsest grid.
For the same experiments as in Table 4, we show large speed up in com- putational time in Table 5. Further, in Figure 1, we demonstrate that the compu- tational complexity of our multigrid procedure is and is the total number of grid points. In Figure 2, we depict the convergence history of the norm of the residuals at a given working level using different numbers of pre- and post-smoothing steps, , , and 5 V-cycle iterations.
We complete this section, presenting results of experiments with a singular kernel, and consider an elliptic PIDE in one dimension of the following form
Table 4. Errors of FAS solution with FK and FI integral evaluation after 5 V-cycles.
Table 5. CPU time (secs.) of FAS solution with 5 V-cycles. In bold are the values of CPU time actually involved in the multigrid solution scheme.
Figure 1. Computational complexity of the FAS-FI method; .
where for . We further assume homogeneous Dirichlet boundary conditions.
We implement the FAS-FI scheme for this PIDE problem whose integral term has a singular kernel with one isolated singularity. On the singularity point, we cannot evaluate the kernel directly. However, we can estimate the integral using its values on neighbouring points. If the singularity is on one of the grid, we
use local averaging . In Figure 3, we depict the
observed multigrid computational complexity when solving the singular kernel problem and see that complexity appears to match or even improve on the typical estimate . Further, in Figure 4 the convergence history of the
Figure 2. Convergence history of the FAS-FI scheme with different , (green) to (red) along 5 V-cycles of FAS; .
Figure 3. Computational complexity of the FAS-FI scheme for the PIDE with a singular kernel.
Figure 4. Convergence history of the FAS-FI scheme for a PIDE with singular kernel.
multigrid scheme with different pre- and post-smoothing schemes applied to the singular kernel case is presented.
An efficient multigrid finite-differences scheme for solving elliptic Fredholm par- tial integro-differential equations (PIDE) was developed and investigated. This scheme combines a FAS multigrid scheme for elliptic problems with a multilevel fast integration technique. Theoretical estimates of second-order solution accuracy and LFA multigrid convergence estimates were presented. These estimates were confirmed by results of numerical experiments.
This work was supported in part by the German Academic Exchange Service (DAAD), the European Union under Grant Agreement Nr. 304617 “Multi-ITN STRIKE-Novel Methods in Computational Finance” and by the Würzburg-Wro- claw Center for Stochastic Computing. This publication was funded by the Ger- man Research Foundation (DFG) and the the University of Würzburg in the fun- ding programme Open Access Publishing.
 Amstrong, N.J., Painter, K.J. and Sheratt, J.A. (2006) A Continuum Approach to Modelling Cell-Cell Adhension. Journal of Theoretical Biology, 243, 98-113.
 Li, H., Di, L., Ware, A. and Yuan, G. (2014) The Applications of Partial Integro Differential Equations Related to Adaptive Wavelet Collocation Methods for Viscosity Solutions to Jump-Diffusion Models. Applied Mathematics and Computation, 246, 316-335.
 Itkin, A. (2016) Efficient Solution of Backward Jump-Diffusion Partial Integro-Differential Equations with Splitting and Matrix Exponentials. Journal of Computational Finance, 19, 29-70.
 Thorwe, J. and Bhalekar, S. (2012) Solving Partial Integro Differential Equation Using Laplace Transform Method. American Journal of Computational and Applied Mathematics, 2, 101-104.
 Ma, J., Jiang, Y. and Xiang, K. (2009) On a Moving Mesh Method for Solving Partial Integro-Differential Equations. Journal of Computational Mathematics, 27, 713-728.
 Soliman, A.F., El-Asyed, A.M.A. and El-Azab, M.S. (2012) On the Numerical Solution of Partial Integro-Differential Equations. Mathematical Sciences Letters, 1, 71-80.
 Liz, E. and Nieto, J.J. (1996) Boundary Value Problems for Second Order Integro-Differential Equations of Fredholm Type. Journal of Computational and Applied Mathematics, 72, 215-225.
 Hu, Q. (1998) Interpolation Correction for Collocation Solutions of Fredholm Integro-Differential Equations. Mathematics of Computation, 67, 987-999.
 Hosseini, S.M. and Shahmorad, S. (2003) Tau Numerical Solution of Fredholm Integro-Differential Equations with Arbitrary Polynimial Bases. Applied Mathematical Modelling, 27, 145-154.
 Brandt, A. and Lubrecht, A.A. (1990) Multilevel Matrix Multiplication and Fast Solution of Integral Equations. Journal of Computational Physics, 90, 348-370.
 Hemker, P.W. (1990) On the Order of Prolongations and Restrictions in Multigrid Procedures. Journal of Computational and Applied Mathematics, 32, 423-429.