Multigrid Solution of an Elliptic Fredholm Partial Integro-Differential Equation with a Hilbert-Schmidt Integral Operator

An efficient multigrid finite-differences scheme for solving elliptic Fredholm partial integro-differential equations (PIDE) is discussed. This scheme combines a second-order accurate finite difference discretization of the PIDE problem with a multigrid scheme that includes a fast multilevel integration of the Fredholm operator allowing the fast solution of the PIDE problem. Theoretical estimates of second-order accuracy and results of local Fourier analysis of convergence of the proposed multigrid scheme are presented. Results of numerical experiments validate these estimates and demonstrate optimal computational complexity of the proposed framework.


Introduction
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 [1] [2] [3] [4].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 [5].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 [6], 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 [7] 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 [8] [9], compact finite-differences for one-dimensional PIDEs are studied.Additional results on high-order schemes for integro differential equations (IDE) can be found in [10].The research in [11], is devoted to an iterated Galerkin method for PIDE in one-dimesion; see [12].Further, the work [13] considers the numerical solution of linear IDE using projection methods.The work in [14] 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-dimensional 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 d Ω ⊂  .Our approach is to combine a multigrid scheme for elliptic problems with the multigrid kernel approximation strategy developed in [15].
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 straightforward implementation of the integral term is considered.For this reason, within the multigrid framework, we investigate the multigrid kernel approximation strategy proposed in [15], where it is demonstrated that it is possible to approximate a Fredholm integral term with ( ) accuracy while reducing the complexity of its calculation from ( ) , 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 uniqueness 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 [15] 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.

An Elliptic Fredholm Partial Integro-Differential Equation
We consider the following PIDE problem where 2 , x y ∈ Ω ⊂  is a two-dimensional, convex and bounded domain with a 2 C boundary or a rectangle.We denote Γ = ∂Ω and Ω = Ω ∪ Γ .We consider ( ) Ω and a symmetric positive semi-definite Hilbert Schmidt kernel ( ) , and the following holds We have the following theorem.
Proof.From Tonelli's theorem, ( )( ) ( ) ( )  is a measurable function of x and its L 2 -norm can be determined by the Cauchy-Schwarz inequality.

Let
( ) Remark 2.1 From Schur's test [16], since the kernel k is a measurable function, it satisfies the following conditions ( ) ( ) ess sup , d , ess sup , d .
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 ( ) ( ) Proof.The proof is straightforward by using the Lax-Milgram theorem and the properties of the kernel.

Discretization of the Elliptic PIDE Problem
We discretize (2.1)-(2.2) using finite differences and the Simpson's rule [17] [18] [19].For simplicity, we assume that ( ) ∈ Ω such that we can evaluate these functions on grid points.Specifically, we consider and N is an integer with 2 N ≥ .We denote ( ) ( ) be the mesh size.We denote the mesh points ( ) . These grid points define the following grid Later, we consider a sequence of nested uniform grids { } 0 , where for ∈   .For grid functions v and w defined on h Ω , we introduce the discrete L 2 -scalar product with associated norm ( ) . The negative Laplacian with homogeneous Dirichlet boundary conditions is approximated by the five-point stencil and is denoted by h −∆ .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 h Ω is represented by , ∈ Ω (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 1 h H norm as ( ) Notice that the bracket ] denotes summation up to N in the given direction 1 x , resp. 2 x ; see [19].With this preparation, we have .
The integral term of the elliptic PIDE in two-dimensions is written explicitly as follows ( )( ) ( ) ( ) where ( ) ( ) ( ) , r l m r l r m =   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 ( ) ( ) L suitable approximation of the Dirac delta function as ( ) 0 x →  , e.g., a narrow Gaussian.Inserting this function in (2.3) we have N Therefore by continuity, as ( ) 0 x →  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 v .With the setting above, we write the finite-difference approximation of (2.1)-(2.2) as follows in , where denotes the numerical approximation to u.Further, the integral function ( 6) evaluated at ( ) Notice that for functions that are zero on the boundary, summation can be restricted to the interior grid points.The double superscript hh in , cates that the pairs of indices ( ) , lm ij refer to the mesh h Ω .Next, we investigate the stability and accuracy of (8).For this purpose, we use the numerical analysis framework in [19].We denote h h h A = −∆ +  .We need the following lemma, see also [19].

Lemma 3.2 Suppose U is a function defined on h
Ω with 0 U = on the boundary; then the following holds ( ) Using the results of Lemma 3.3, we have

Lemma 3.3 Suppose U is a function defined on h
Ω with 0 U = on the boundary; then there exists a constant * ρ , which is independent of U and h , such that the following discrete Poincaré-Friedrichs inequality holds ( ) for all such U ; see [19].
Remark 3.1 From (3.6) and (3.7), we obtain ( ) We conclude this section with the following theorem.
Theorem 3.5 Suppose ( ) Ω and ( ) , 0 k x y ≥ is a Hilbert Schmidt kernel, and assume that the weak solution u to (2.1)-(2.2) belongs to ( ) ; then the solution U to (3.5) approximates u 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 [19].

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 [20] [21] [22] and the multigrid fast integration technique presented in [15].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 h Ω and the coarse grid H Ω , where 2 H h = .In h Ω , consider the discretized PIDE Equation (3.5) as follows , where h U 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 longwavelength 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 . Equivalently, this correction can be defined as the solution to ( ) where Notice that this equation can be re-written as β is the so-called fine-to-coarse defect correction.Now, suppose to solve (4.4) to obtain H U .Then we can compute Ĥ , which represents the coarse-grid approximation to h e .Notice that while H U and ˆH h h U   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 h e that is used to correct h U  .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 H h  is applied to the residual, while in (4.5) an interpolation operator h H  is applied to the error function.We choose the following inter-grid transfer operators [23] [24] [ ] [ ] 1 9 9 1 and : 1 0 9 16 9 0 1 .16 32 The transfer operators given by (4.6) are of 4th-order.The 4th-order interpolation operator is symmetric and accesses two grid points on either side of the interpolated point.To approximate at the grid point

= − 
Notice that our choice of higher-order interpolation and restriction operators, as in [15], appears advantageous for the fast integration technique that we discuss next.
The fast integration strategy [15] [25] 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 onedimensional case, denote with ( ) the grid points on the grid with mesh size h, and with ( ) the grid points on the grid with mesh Now, suppose that the kernel ( ) , k x y and u are sufficiently smooth.(For the case of singular kernels see [15].)In h Ω , the integral ( )( ) ( ) ( )  .On the other hand, in the strategy of [15], the kernel is approximated by , ,.
where the interpolation operator h H  may be equal to h H  .With this setting, we have where , where H U is obtained by coarsening of h U .In particular, using straight injection for boundary values or the full-weighted restriction in (4.6), we have 2 . 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 ( ) 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 procedure just described can be applied recursively, until a grid is reached with ( ) N  grid points.On this grid the summation is then actually performed, requiring ( ) operations.Further, the computational effort of the restriction  ) , where p is the order of interpolation ( 2 p = for linear interpolation).Therefore the order of total work required to obtain the (approximated) summation is ( ) Notice that choosing H h = , the quadrature error on using Simpson's rule is 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 1 γ = and a W-cycle results when 2 γ = .

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) [20] [21] [26] [27].
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 ( ) cays rapidly to zero as x y − becomes large.With these assumptions the stencil of our PIDE operator can be cast in the standard LFA framework.However, treating 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,  We have

{ }
For each θ and a translational invariant kernel, we assume that the space h E θ is invariant under the action of TG .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 H h TG on the following function θ ∈ − and 0,1 α = , 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 h E θ under the action of the smoothing operator.We have , 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 h G .We have At this point, recall that ( ) ( ) Since we assume that the element jl becomes very small as n becomes large, we truncate the sum in (5.4) and consider the following equation 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 estimate the smoothing factor of our GSP scheme is given by ( With this definition, we obtain a smoothing factor of our GSP scheme given by GSP 0.447 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 h h h A = −∆ +  , applied to a generic vector v with jth component given by e i j Now, recall that on the fine grid, we distinguish on the two harmonics.Therefore 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 [27] ( )   where ρ denotes the spectral radius of the 2 2 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,

Numerical Experiments
In this 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 Gaussian convolution kernel in two dimensions as follows , ,  , d d  , ,  u x y k x y t s u t s t s f x y where ( ) ( ) and the Gaussian kernel , , , exp 2 To investigate the order of accuracy of the discretization scheme, we construct an exact solution to (6.1) by choosing ( ) , exp 2 x y u x y . 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 finitedifferences 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 In Table 3  Table 3. Errors for 2D integral evaluation using 4 th −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 corresponding 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.

, we report the values of
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 computational time in Table 5.Further, in Figure 1, we demonstrate that the computational complexity of our multigrid procedure is   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 i x of the grid, we use local averaging ( ) . In Figure 3  multigrid scheme with different pre-and post-smoothing schemes applied to the singular kernel case is presented.

Conclusion
An efficient multigrid finite-differences scheme for solving elliptic Fredholm partial 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.


Using the Simpson's rule, we have the following approximation of this integral operator 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

(
can illustrate the multigrid procedure considering a sequence of nested grids (levels) k k h Ω = Ω of mesh size k h , indexed by the level number 1, , k l =  , where l 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 k l d = − , with given depth d , is reached.
and therefore 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 that the partial sum provides a sufficiently accurate approximation of the integral term on h G .Specifically, assuming that cision machine epsilon) for n L > , one should choose 1 L h ∝ .However, in practice, a much smaller L results in accurate LFA estimates.
are able to compute the two-grid convergence factor as follows

3
This numerical convergence factor represents the asymptotic ratio of reduction of the2  L -norm of the residual between two multigrid cycles.These calculations refer to the choice ( ) in Table1, the LFA estimates of the multigrid convergence factor are accurate.(The same values of TG µ are obtained with L ranging from 20 to 400.) the full-kernel (FK) evaluation formula (6) and the FI technique involving different depths, d l k = − .For 0 d = the FI scheme performs FK evaluation.
working levels l 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 of grid points.In Figure2, we depict the convergence history of the norm of the residuals at a given working level using different numbers of preand post-smoothing steps, 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

Figure 1 .
Figure 1.Computational complexity of the FAS-FI method;

Figure 4 .
Figure 4. Convergence history of the FAS-FI scheme for a PIDE with singular kernel.
smoothing property is such that the solution error

Table 1 .
Estimated and observed multigrid convergence factors.

Table 2 .
L 2 −norm error using full kernel approximation.

Table 4 .
Errors of FAS solution with FK and FI integral evaluation after 5 V-cycles.

Table 5 .
CPU (secs.) of solution with 5 V-cycles.In bold are the values of CPU time actually involved in the multigrid solution scheme.