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

Duncan Kioi Gathungu, Alfio Borzì

Institut für Mathematik, Universität Würzburg, Würzburg, Germany.

DOI: 10.4236/am.2017.87076                

Abstract

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.

Keywords

Elliptic Problems, Fredholm Operator, Multigrid Schemes, Finite Differences, Numerical Analysis

1. 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-dimen- sional PIDE problems. We present a second-order accurate fast multigrid scheme to solve elliptic problems of the following form

$Au\left(x\right)+\underset{\Omega }{{\displaystyle \int }}k\left(x\mathrm{,}y\right)u\left(y\right)\text{d}y=f\left(x\right)\text{in}\Omega \mathrm{,}$

where A represents an elliptic operator with given boundary conditions and $\Omega \subset {ℝ}^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 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 [15] , where it is demonstrated that it is possible to approximate a Fredholm integral term with $\mathcal{O}\left(h^{2s}\right)$ accuracy while reducing the complexity of its calculation from $\mathcal{O}\left(n^2\right)$ to $\mathcal{O}\left(sn\right)$ , 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 [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.

2. An Elliptic Fredholm Partial Integro-Differential Equation

We consider the following PIDE problem

$-\Delta u\left(x\right)+\underset{\Omega }{{\displaystyle \int }}k\left(x\mathrm{,}y\right)u\left(y\right)\text{d}y=f\left(x\right)\text{in}\Omega \mathrm{,}$ (2.1)

$u\left(x\right)=0\text{on}\Gamma \mathrm{,}$ (2.2)

where $x\mathrm{,}y\in \Omega \subset {ℝ}^2$ is a two-dimensional, convex and bounded domain with a $C^2$ boundary or a rectangle. We denote $\Gamma =\partial \Omega$ and $\bar{\Omega }=\Omega \cup \Gamma$ . We consider $f\in L^2\left(\Omega \right)$ and a symmetric positive semi-definite Hilbert Schmidt kernel

$k\in L^2\left(\Omega \times \Omega \right)$ , such that $\underset{\Omega }{{\displaystyle \int }}\underset{\Omega }{{\displaystyle \int }}{\left|k\left(x,y\right)\right|}^2\text{d}x\text{d}y<\infty$ , and the following holds

$\underset{\Omega }{{\displaystyle \int }}\underset{\Omega }{{\displaystyle \int }}k\left(x\mathrm{,}y\right)v\left(x\right)v\left(y\right)\text{d}x\text{d}y\ge \mathrm{0,}\text{for}\text{all}v\in L^2\left(\Omega \right)\mathrm{.}$ (2.3)

We have the following theorem.

Theorem 2.1 Let $k\in L^2\left(\Omega \times \Omega \right)$ be a Hilbert Schmidt kernel. The integral operator $\mathcal{I}$ given by

$\left(\mathcal{I}u\right)\left(x\right)=\underset{\Omega }{{\displaystyle \int }}k\left(x\mathrm{,}y\right)u\left(y\right)\text{d}y\mathrm{,}x\in \Omega ,$

defines a bounded mapping of $L^2\left(\Omega \right)$ into itself, with the Hilbert Schmidt norm $\Vert \mathcal{I}\Vert \le {\Vert k\Vert }_2$ .

Proof. From Tonelli’s theorem, $\left(\mathcal{I}u\right)\left(x\right)=\underset{\Omega }{{\displaystyle \int }}k\left(x\mathrm{,}y\right)u\left(y\right)\text{d}y$ is a measurable function of x and its L²-norm can be determined by the Cauchy-Schwarz inequality. Let $u\in L^2\left(\Omega \right)$ , we have

$\begin{array}{c}{\Vert \mathcal{I}u\Vert }_2^2=\underset{\Omega }{{\displaystyle \int }}{\left|\mathcal{I}u\left(x\right)\right|}^2\text{d}x=\underset{\Omega }{{\displaystyle \int }}{\left|\underset{\Omega }{{\displaystyle \int }}k\left(x\mathrm{,}y\right)u\left(y\right)\text{d}y\right|}^2\text{d}x\\ \le \underset{\Omega }{{\displaystyle \int }}\left(\underset{\Omega }{{\displaystyle \int }}{\left|k\left(x\mathrm{,}y\right)\right|}^2\text{d}y\right)\left(\underset{\Omega }{{\displaystyle \int }}{\left|u\left(y\right)\right|}^2\text{d}y\right)\text{d}x\\ =\underset{\Omega }{{\displaystyle \int }}\underset{\Omega }{{\displaystyle \int }}{\left|k\left(x\mathrm{,}y\right)\right|}^2{\Vert u\Vert }_2^2\text{d}y\text{d}x={\Vert k\Vert }_2^2{\Vert u\Vert }_2^2<\infty .\end{array}$

Hence $\mathcal{I}u\in L^2\left(ℝ\right)$ .

Remark 2.1 From Schur’s test [16] , since the kernel k is a measurable function, it satisfies the following conditions

${\xi }_1=\underset{x\in ℝ}{\text{ess}\text{sup}}\underset{\Omega }{{\displaystyle \int }}\left|k\left(x\mathrm{,}y\right)\right|\text{d}y<\infty \mathrm{,}{\xi }_2=\underset{y\in ℝ}{\text{ess}\text{sup}}\underset{\Omega }{{\displaystyle \int }}\left|k\left(x\mathrm{,}y\right)\right|\text{d}x<\infty \mathrm{.}$

Then the integral operator $\mathcal{I}$ defines a bounded mapping and ${\Vert \mathcal{I}\Vert }_2\le {\left({\xi }_1{\xi }_2\right)}^{\frac{1}{2}}$ .

With this preparation, we can prove the following.

Theorem 2.2 There exist a unique function $u\in H_0^1\left(\Omega \right)\cap H^2\left(\Omega \right)$ 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 [17] [18] [19] . For simplicity, we assume that $k\in C\left(\bar{\Omega }\times \bar{\Omega }\right)$ and $f\in C\left(\bar{\Omega }\right)$ such that we can evaluate these functions on grid points. Specifically, we consider $\Omega =\left(a,b\right)\times \left(a,b\right)$ and N is an integer with $N\ge 2$ . We denote $x=\left(x_1,x_2\right)$ and

$y=\left(y_1,y_2\right)$ . Let $h=\frac{b-a}{N-1}$ be the mesh size. We denote the mesh points

#Math_37#, and $x_{2j}=a+\left(j-1\right)h,j=1,\cdots ,N$ . These grid points define the following grid

${\Omega }_h=\left\{x_{ij}=\left(x_{1i},x_{2j}\right)\in {ℝ}^2:i,j=2,\cdots ,N-1\right\}\cap \Omega .$

Later, we consider a sequence of nested uniform grids ${\left\{{\Omega }_{h_{\mathcal{l}}}\right\}}_{h_{\mathcal{l}}>0}$ , where $N=N_{\mathcal{l}}=2^{\mathcal{l}}+1$ for $\mathcal{l}\in \mathbb{N}$ .

For grid functions v and w defined on ${\Omega }_h$ , we introduce the discrete L²-scalar product

${\left(v,w\right)}_h=h^2{\displaystyle \underset{x\in {\Omega }_h}{\sum }}v\left(x\right)w\left(x\right),$

with associated norm ${\Vert v\Vert }_h=\sqrt{{\left(v,v\right)}_h}$ . The negative Laplacian with homogeneous Dirichlet boundary conditions is approximated by the five-point stencil and is denoted by $-{\Delta }_h$ . Given continuous functions in $\Omega$ are approximated by grid functions defined through their values at the grid points. Thus the right-hand side of (2.1) in ${\Omega }_h$ is represented by $f_{ij}^h=f\left(x_{1i},x_{2j}\right)$ , if $f\in C\left(\Omega \right)$ (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

$D_{x_1}^{+}u\left(x_{1i},x_{2j}\right)\equiv \frac{u\left(x_{1i+1},x_{2j}\right)-u\left(x_{1i},x_{2j}\right)}{h}.$ (3.1)

The backward finite-difference operator is as follows

$D_{x_1}^{+}u\left(x_{1i}\mathrm{,}{x}_{2j}\right)\equiv \frac{u\left(x_{1i}\mathrm{,}{x}_{2j}\right)-u\left(x_{1i-1}\mathrm{,}{x}_{2j}\right)}{h}\mathrm{.}$ (3.2)

With these operators, we can define the $H_h^1$ norm as ${\Vert u\Vert }_{1,h}={\left({\Vert u\Vert }_h^2+\left|\right|D_{x_1}^{-}u]{|}_{x_1}^2+||D_{x_2}^{-}u]{|}_{x_2}^2\right)}^{\frac{1}{2}}$ . Notice that the bracket $\mathrm{<mo\; stretchy="false">\; ]\; </mo>}$ denotes summation up to N in the given direction $x_1$ , resp. $x_2$ ; see [19] . With this preparation, we have

$\begin{array}{l}{\Delta }_{h}u\left(x_{1i},x_{2j}\right)=D_{x1}^{-}{D}_{x1}^{+}u\left(x_{1i},x_{2j}\right)+D_{x2}^{-}{D}_{x2}^{+}u\left(x_{1i},x_{2j}\right)\\ =\frac{u\left(x_{1i+1},x_{2j}\right)-2u\left(x_{1i},x_{2j}\right)+u\left(x_{1i-1},x_{2j}\right)}{h^2}+\frac{u\left(x_{1i},x_{2j+1}\right)-2u\left(x_{1i},x_{2j}\right)+u\left(x_{1i},x_{2j-1}\right)}{h^2}.\end{array}$

The integral term of the elliptic PIDE in two-dimensions is written explicitly as follows

$\left(\mathcal{I}u\right)\left(x\right)=\underset{\Omega }{{\displaystyle \int }}k\left(x\mathrm{,}y\right)u\left(y\right)\text{d}y\mathrm{,}x\in \Omega \mathrm{.}$

Using the Simpson’s rule, we have the following approximation of this integral operator

$\left({\mathcal{I}}^{h}u\right)\left(x\right)=h^2{\displaystyle \underset{l=1}{\overset{N}{\sum }}}{\displaystyle \underset{m=1}{\overset{N}{\sum }}}r\left(l,m\right)k\left(x,x_{lm}\right)u\left(x_{lm}\right),x\in {\Omega }_h,$ (3.3)

where $r\left(l\mathrm{,}m\right)=\tilde{r}\left(l\right)\tilde{r}\left(m\right)$ represents the coefficients of the quadrature rule. In the case of the Simpson’s rule, we have

$\tilde{r}\left(l\right)=(\begin{array}{ll}\frac{1}{3}\hfill & \text{if}l=1,l=N\hfill \\ \frac{4}{3}\hfill & lmod2=0,l=2,\cdots ,N-1\hfill \\ \frac{2}{3}\hfill & \text{else}\text{.}\hfill \end{array}$

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 $v\left(x\right)={\displaystyle {\sum }_{m=1}^l}{\tilde{\delta }}_{ϵ}\left(x-x_{lm}\right)v_{lm}$ where ${\tilde{\delta }}_{\epsilon }\left(x\right)$ is a $L^2$ suitable approximation of the Dirac delta function as $ϵ\left(x\right)\to 0$ , e.g., a narrow Gaussian. Inserting this function in (2.3) we have

${\displaystyle \underset{l,m=1}{\overset{N}{\sum }}}{\displaystyle \underset{i,j=1}{\overset{N}{\sum }}}{v}_{lm}{v}_{ij}\underset{\Omega }{{\displaystyle \int }}\underset{\Omega }{{\displaystyle \int }}k\left(x,y\right){\tilde{\delta }}_{ϵ}\left(x\right)\left(x-x_{lm}\right){\tilde{\delta }}_{ϵ}\left(y\right)\left(y-y_{ij}\right)\text{d}x\text{d}y\ge 0.$ (3.4)

Therefore by continuity, as $ϵ\left(x\right)\to 0$ the above integral tends to $k\left(x_{lm}\mathrm{,}{y}_{ij}\right)$ . Thus, we obtain $\left({\mathcal{I}}^{h}u\mathrm{,}u\right)\ge 0$ .

The Simpson’s rule provides a fourth-order accurate approximation of the integral as follows

${\Vert {\mathcal{I}}^{h}v-\mathcal{I}v\Vert }_h=\mathcal{O}\left(h^4\right)\mathrm{,}$

for any sufficiently smooth $v$ .

With the setting above, we write the finite-difference approximation of (2.1)- (2.2) as follows

$-{\Delta }_{h}U+{\mathcal{I}}^{h}U=f^h\text{in}{\Omega }_h\mathrm{,}$ (3.5)

where $U=\left(U_{ij}\right)$ denotes the numerical approximation to u. Further, the integral function (6) evaluated at $\left(x_{1i}\mathrm{,}{x}_{2j}\right)$ is given as follows

${\left({\mathcal{I}}^{h}U\right)}_{ij}=h^2{\displaystyle \underset{lm}{\sum }}{k}_{ij\mathrm{,}lm}^{hh}{U}_{lm}\mathrm{.}$

Notice that for functions that are zero on the boundary, summation can be restricted to the interior grid points. The double superscript $hh$ in $k_{lm\mathrm{,}ij}^{hh}$ indi- cates that the pairs of indices $\left(lm\mathrm{,}ij\right)$ refer to the mesh ${\Omega }_h$ .

Next, we investigate the stability and accuracy of (8). For this purpose, we use the numerical analysis framework in [19] . We denote $A_h=-{\Delta }_h+{\mathcal{I}}^h$ .

We need the following lemma, see also [19] .

Lemma 3.2 Suppose $U$ is a function defined on ${\bar{\Omega }}_h$ with $U=0$ on the boundary; then the following holds

${\left(A_{h}U,U\right)}_h\ge \stackrel{N}{\underset{i=1}{{\displaystyle \sum }}}\stackrel{N-1}{\underset{j=1}{{\displaystyle \sum }}}{h}^2{\left|D_{x1}^{-}{U}_{ij}\right|}^2+\stackrel{N-1}{\underset{i=1}{{\displaystyle \sum }}}\stackrel{N}{\underset{j=1}{{\displaystyle \sum }}}{h}^2{\left|D_{x2}^{-}{U}_{ij}\right|}^2$

Proof. Using the results of Lemma 3.3, we have

$\begin{array}{c}{\left(AU\mathrm{,}U\right)}_h={\left(-D_{x1}^{+}{D}_{x1}^{-}U-D_{x2}^{+}{D}_{x2}^{-}U+{\mathcal{I}}^{h}U\mathrm{,}U\right)}_h\\ ={\left(-D_{x1}^{+}{D}_{x1}^{-}U\mathrm{,}U\right)}_h+{\left(-D_{x2}^{+}{D}_{x2}^{-}U\mathrm{,}U\right)}_h+{\left({\mathcal{I}}^{h}U\mathrm{,}U\right)}_h\\ \ge \stackrel{N}{\underset{i\mathrm{<mo>\; =\; </mo>1}}{{\displaystyle \sum }}}\stackrel{N-1}{\underset{j\mathrm{<mo>\; =\; </mo>1}}{{\displaystyle \sum }}}{h}^2{\left|D_{x1}{U}_{ij}\right|}^2+\stackrel{N-1}{\underset{i\mathrm{<mo>\; =\; </mo>1}}{{\displaystyle \sum }}}\stackrel{N}{\underset{j\mathrm{<mo>\; =\; </mo>1}}{{\displaystyle \sum }}}{h}^2{\left|D_{x2}{U}_{ij}\right|}^2\\ \ge \left|\right|D_{x1}^{-}U]{|}_{x1}^2+||D_{x2}^{-}U]{|}_{x2}^2\end{array}$ (3.6)

Lemma 3.3 Suppose U is a function defined on ${\bar{\Omega }}_h$ with $U=0$ on the boun- dary; then there exists a constant ${\rho }_{\mathrm{<mo>\; *\; </mo>}}$ , which is independent of $U$ and $h$ , such that the following discrete Poincaré-Friedrichs inequality holds

${\Vert U\Vert }_h^2\le {\rho }_{\mathrm{<mo>\; *\; </mo>}}\left(\left|\right|D_{x1}^{-}U]{|}_{x1}^2+||D_{x2}^{-}U]{|}_{x2}^2\right)$ (3.7)

for all such $U$ ; see [19] .

Remark 3.1 From (3.6) and (3.7), we obtain ${\left(-{\Delta }_{h}U+{\mathcal{I}}^{h}U\mathrm{,}U\right)}_h\ge {\rho }_0{\Vert U\Vert }_{\mathrm{1,}h}^2$ , where ${\rho }_0=1/\left(1+{\rho }^{*}\right)$ .

Theorem 3.4 The scheme (3.5) is stable in the sense that ${\Vert U^{\mathrm{<mo>\; *\; </mo>}}\Vert }_{\mathrm{1,}h}\le \frac{1}{{\rho }_0}{\Vert f^h\Vert }_h$

Proof. We have

$\begin{array}{c}{\rho }_0{\Vert U\Vert }_{\mathrm{1,}h}^2\le {\left(-{\Delta }_{h}U+{\mathcal{I}}^{h}U\mathrm{,}U\right)}_h={\left(f^h\mathrm{,}U\right)}_h\le \left|{\left(f^h\mathrm{,}U\right)}_h\right|\\ \le {\Vert f^h\Vert }_h{\Vert U\Vert }_h\le {\Vert f^h\Vert }_h{\Vert U\Vert }_{\mathrm{1,}h}\mathrm{,}\end{array}$

hence ${\Vert U^{\mathrm{<mo>\; *\; </mo>}}\Vert }_{\mathrm{1,}h}\le \frac{1}{{\rho }_0}{\Vert f^h\Vert }_h$ .

We conclude this section with the following theorem.

Theorem 3.5 Suppose $f\in L^2\left(\Omega \right)$ and $k\left(x\mathrm{,}y\right)\ge 0$ is a Hilbert Schmidt ker- nel, and assume that the weak solution $u$ to (2.1)-(2.2) belongs to $H_0^1\cap H^3\left(\Omega \right)$ ; then the solution $U$ to (3.5) approximates $u$ with second-order accuracy as follows

${\Vert u-U\Vert }_{\mathrm{1,}h}\le c{h}^2\mathrm{,}$

where c is a positive constant independent of h. In particular ${\Vert u-U\Vert }_h\le c{h}^2$ .

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

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 [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 ${\Omega }_h$ and the coarse grid ${\Omega }_H$ , where $H=2h$ .

In ${\Omega }_h$ , consider the discretized PIDE Equation (3.5) as follows

$A_h{U}^h=f^h,$ (4.1)

where $U^h$ 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 $U^{h\mathrm{,}m-1}$ , it results in $U^{h\mathrm{,}m}=S\left(U^{h\mathrm{,}m-1}\mathrm{,}{f}^h\right)$ . The smoothing property is such that the solution error $e^{h,m}=U^h-U^{h,m}$ has smaller higher-frequency modes than the error $e^{h,m-1}=U^h-U^{h,m-1}$ .

In the multigrid solution process, starting with an initial approximation $U^{h\mathrm{,0}}$ and applying S to (4.1) ${\nu }_1$ -times, we obtain the approximate solution ${\tilde{U}}^h=U^{h,{\nu }_1}$ .

Now, the desired (smooth) correction $e_h$ to ${\tilde{U}}^h$ , to obtain the exact solution, is defined by $A_h\left({\tilde{U}}^h+e^h\right)=f^h$ . Equivalently, this correction can be defined as the solution to

$A_h\left({\tilde{U}}^h+e^h\right)=r^h+A_h{\tilde{U}}^h,$ (4.2)

where $r^h=f^h-A_h{\tilde{U}}^h$ is the residual associated to ${\tilde{U}}^h$ .

Next, notice that the structure of $A_h$ and the smoothness of the error function allow to represent (4.2) on the coarse grid ${\Omega }_H$ . On this grid, ${\tilde{U}}^h+e^h$ is repre- sented in terms of coarse variables as follows

$U^H={\hat{\mathbb{I}}}_h^H{\tilde{U}}^h+e^H\mathrm{,}$ (4.3)

where ${\hat{\mathbb{I}}}_h^H{\tilde{U}}^h$ represents the restriction of ${\tilde{U}}_h$ to the coarse grid by means of the direct injection operator denoted with ${\hat{\mathbb{I}}}_h^H$ .

With this preparation, it appears natural to approximate (4.2) on the coarse grid as follows

$A_H{U}^H={\mathbb{I}}_h^H\left(f^h-A_h{\tilde{U}}^h\right)+A_H{\hat{\mathbb{I}}}_h^H{\tilde{U}}^h\mathrm{.}$ (4.4)

Notice that this equation can be re-written as $A_H{U}^H={\mathbb{I}}_h^H{f}^h+{\beta }_h^H$ where ${\beta }_h^H=A_H{\hat{\mathbb{I}}}_h^H{\tilde{U}}^h-{\mathbb{I}}_h^H{A}_h{\tilde{U}}^h$ . The term ${\beta }_h^H$ is the so-called fine-to-coarse defect correction.

Now, suppose to solve (4.4) to obtain $U^H$ . Then we can compute $e^H=U^H-{\hat{\mathbb{I}}}_h^H{\tilde{U}}^h$ , which represents the coarse-grid approximation to $e^h$ . Notice that while $U^H$ and ${\hat{\mathbb{I}}}_h^H{\tilde{U}}^h$ 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 $e^h$ that is used to correct ${\tilde{U}}^h$ . This procedure defines the following coarse-grid correction step

$U^h={\tilde{U}}^h+{\mathbb{I}}_H^h\left(U^H-{\hat{\mathbb{I}}}_h^H{\tilde{U}}^h\right)\mathrm{.}$ (4.5)

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 ${\mathbb{I}}_h^H$ is applied to the residual, while in (4.5) an interpolation operator ${\mathbb{I}}_H^h$ is applied to the error function. We choose the following inter-grid transfer operators [23] [24]

${\mathbb{I}}_H^h\mathrm{<mo>\; :\; </mo>}=\frac{1}{16}\left[\begin{array}{cccc}-1& 9& 9& -1\end{array}\right]\text{and}{\mathbb{I}}_h^H\mathrm{<mo>\; :\; </mo>}=\frac{1}{32}\left[\begin{array}{ccccccc}-1& 0& 9& 16& 9& 0& -1\end{array}\right]\mathrm{.}$ (4.6)

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 $i=2$ and $i=N-1$ , 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

${\mathbb{I}}_H^h\mathrm{<mo>\; :\; </mo>}=\frac{1}{16}\left[\begin{array}{cccc}5& 15& -5& 1\end{array}\right]\mathrm{.}$

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 one- dimensional case, denote with $x_j=a+\left(j-1\right)h$ the grid points on the grid with mesh size h, and with $x_J=a+\left(J-1\right)H$ the grid points on the grid with mesh size $H=2h$ . Notice that $x_j=x_J$ for $j=2J-1$ .

Now, suppose that the kernel $k\left(x\mathrm{,}y\right)$ and $u$ are sufficiently smooth. (For the case of singular kernels see [15] .) In ${\Omega }_h$ , the integral

$\left(\mathcal{I}u\right)\left(x\right)=\underset{\Omega }{{\displaystyle \int }}k\left(x\mathrm{,}y\right)u\left(y\right)\text{d}y$ is approximated by ${\left({\mathcal{I}}^{h}U\right)}_i=h\underset{j}{{\displaystyle \sum }}{k}_{i\mathrm{,}j}^{hh}{U}_j$ . On the

other hand, in the strategy of [15] , the kernel is approximated by ${\bar{k}}_{i\mathrm{,}j}^{hh}={\left[{\bar{\mathbb{I}}}_H^h{k}_{i\mathrm{,.}}^{hH}\right]}_j$ , where the interpolation operator ${\bar{\mathbb{I}}}_H^h$ may be equal to ${\mathbb{I}}_H^h$ . With this setting, we have

$\begin{array}{c}{\left({\mathcal{I}}^h{U}^h\right)}_i\approx {\left({\bar{\mathcal{I}}}^h{U}^h\right)}_i=h{\displaystyle \underset{j}{\sum }}{\bar{k}}_{i\mathrm{,}j}^{hh}{U}_j^h=h{\displaystyle \underset{j}{\sum }}{\left[{\mathbb{I}}_H^h{k}_{i\mathrm{,.}}^{hH}\right]}_j{U}_j^h\\ =h{\displaystyle \underset{j}{\sum }}{k}_{i\mathrm{,}J}^{hH}{\left[{\left({\mathbb{I}}_H^h\right)}^{\text{T}}{U}^h\right]}_J=H\underset{j}{{\displaystyle \sum }}{k}_{i\mathrm{,}J}^{hH}{U}_J^H\text{where}{U}^H={\mathbb{I}}_h^H{U}^h\mathrm{,}\end{array}$

where $U^H$ is obtained by coarsening of $U^h$ . In particular, using straight injection for boundary values or the full-weighted restriction in (4.6), we have $U_J^H=U_{2J}^h$ . Now, we go a step further and consider the coarse integral function

${\left({\mathcal{I}}^H{U}^H\right)}_I=H{\displaystyle \sum }{k}_{I\mathrm{,}J}^{HH}{U}_J^H$ . This function is evaluated on the coarse grid and, from

the calculation above it is clear that it is equal to ${\left({\bar{\mathcal{I}}}^h{U}^h\right)}_i$ for all $i=2I-1$ . Therefore we obtain the following approximation to the integral function on the fine grid

${\mathcal{I}}^h{U}^h\approx {\mathbb{I}}_H^h\left({\mathcal{I}}^H{U}^H\right)\text{where}{U}^H={\mathbb{I}}_h^H{U}^h\mathrm{.}$ (4.7)

In one dimension, the summation complexity on the coarse grid is of order $\mathcal{O}\left(N^2/2\right)$ 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 $\mathcal{O}\left(\sqrt{N}\right)$ grid points. On this grid the summation is then actually performed, requiring $\mathcal{O}\left(N\right)$ operations. Further, the computational effort of the restriction ${�}_h^H{U}^h$ and of the interpolation ${\mathbb{I}}_H^h\left({\mathcal{I}}^H{U}^H\right)$ is $\mathcal{O}\left(2pN\right)$ , where $p$ is the order of interpolation ( $p=2$ for linear interpolation). Therefore the order of total work required to obtain the (approximated) summation is $\mathcal{O}\left(N\right)$ operations. Notice that choosing $H=\sqrt{h}$ , the quadrature error on using Simpson’s rule is $\mathcal{O}\left(H^4\right)=\mathcal{O}\left(h^2\right)$ . This gives an overall $\mathcal{O}\left(h^2\right)$ order of accuracy for the whole procedure.

Now, we can illustrate the multigrid procedure considering a sequence of nested grids (levels) ${\Omega }_k={\Omega }_{h_k}$ of mesh size $h_k$ , indexed by the level number $k=1,\cdots ,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.

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 $\gamma$ 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 $\gamma =1$ and a W-cycle results when $\gamma =2$ .

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) [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 $k\left(x,y\right)=k\left(\left|x-y\right|\right)$ and require that $k\left(\left|x-y\right|\right)$ de- cays rapidly to zero as $\left|x-y\right|$ 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

$T{G}_h^H=S_h^{{\nu }_2}\left[{\mathbb{I}}_h-{\mathbb{I}}_H^h{\left(A_H\right)}^{-1}{\mathbb{I}}_h^H{A}_h\right]S_h^{{\nu }_1}\mathrm{,}$ (5.1)

where different pre- and post-smoothing steps are considered. The coarse grid operator is given by $C{G}_h^H=\left[{\mathbb{I}}_h-{\mathbb{I}}_H^h{\left(A_H\right)}^{-1}{\mathbb{I}}_h^H{A}_h\right]$ .

The local Fourier analysis considers infinite grids, $G_h=\left\{jh\mathrm{,}j\in \mathbb{Z}\right\}$ , 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

${\varphi }_h\left(\theta \mathrm{,}x\right)={\text{e}}^{i\theta x/h}\mathrm{,}\theta \in \left(-\text{π}\mathrm{,}\text{π}\right]\mathrm{.}$

For any low frequency ${\theta }^0\in \left[-\text{π}/2\mathrm{,}\text{π}/2\right)$ , we consider the high frequency mode given by

${\theta }^1={\theta }^0-\text{signum}\left({\theta }^0\right)\text{π}.$ (5.2)

We have ${\varphi }_h\left({\theta }^0\mathrm{,}\cdot \right)={\varphi }_h\left({\theta }^1\mathrm{,}\cdot \right)$ for ${\theta }^0\in \left[-\text{π}/2\mathrm{,}\text{π}/2\right)$ and $x\in G_H$ . We also have ${\varphi }_h\left(\theta \mathrm{,}x\right)={\varphi }_H\left(2{\theta }^0\mathrm{,}x\right)$ on $G_H$ for $\theta ={\theta }^0$ and $\theta ={\theta }^1$ .

The two components ${\varphi }_h\left({\theta }^0\mathrm{,}\cdot \right)$ and ${\varphi }_h\left({\theta }^1\mathrm{,}\cdot \right)$ are called harmonics. For a given ${\theta }^0\in \left[-\text{π}/2\mathrm{,}\text{π}/2\right)$ , the two-dimensional space of harmonics is defined by

$E_h^{\theta }=span\left[{\varphi }_h\left({\theta }^{\alpha }\mathrm{,}\cdot \right)\mathrm{<mo>\; :\; </mo>}\alpha \in \left\{\mathrm{0,1}\right\}\right]\mathrm{.}$

For each $\theta$ and a translational invariant kernel, we assume that the space $E_h^{\theta }$ is invariant under the action of $T{G}_h^H$ . 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 $T{G}_h^H$ on the following function

$\psi \left(x_j\right)={\displaystyle \underset{\alpha \mathrm{,}\theta }{\sum }}{A}_{\theta }^{\alpha }{\varphi }_h\left({\theta }^{\alpha }\mathrm{,}{x}_j\right)\mathrm{,}{x}_j\in G_h\mathrm{.}$

Specifically, we determine how the coefficients $A_{\theta }^{\alpha }$ , $\theta \in \left[-\text{π}/2\mathrm{,}\text{π}/2\right)$ and $\alpha =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 $\alpha$ as we assume invariance of $E_h^{\theta }$ under the action of the smoothing operator. We have

$e_j^m=\underset{\theta }{{\displaystyle \sum }}{\text{E}}_{\theta }^m{\text{e}}^{i\theta j}\text{and}{e}_j^{m+1}=\underset{\theta }{{\displaystyle \sum }}{\text{E}}_{\theta }^{m+1}{\text{e}}^{i\theta j}\mathrm{,}$ (5.3)

where ${\text{E}}_{\theta }^m$ and ${\text{E}}_{\theta }^{m+1}$ denote the error amplitude after $m$ and $m+1$ iterations of the smoother. Notice that $e^{m+1}=S_h{e}^m$ , and thus for the coefficient of the $\theta$ Fourier mode, we have ${\text{E}}_{\theta }^{m+1}={\hat{S}}_h\left(\theta \right){\text{E}}_{\theta }^m\mathrm{.}$

Now, consider the following point-wise definition of our GSP iteration applied to our discretized PIDE on $G_h$ . We have

$\frac{2}{h^2}{e}_j^{m+1}-\frac{1}{h^2}{e}_{j-1}^{m+1}=\frac{1}{h^2}{e}_{j+1}^m-h\stackrel{\infty }{\underset{l=-\infty }{{\displaystyle \sum }}}{k}_{jl}{e}_l^m\mathrm{.}$ (5.4)

At this point, recall that $k\left(x\mathrm{,}y\right)=k\left(\left|x-y\right|\right)$ and introduce the index $n=j-l$ . Since we assume that the element $k_{jl}=k_{\left|j-l\right|}=k_{\left|n\right|}$ becomes very small as $n$ becomes large, we truncate the sum in (5.4) and consider the following equation

$\frac{2}{h^2}{e}_j^{m+1}-\frac{1}{h^2}{e}_{j-1}^{m+1}=\frac{1}{h^2}{e}_{j+1}^m-h\stackrel{L}{\underset{n=-L}{{\displaystyle \sum }}}{k}_{\left|n\right|}{e}_{j-n}^m.$ (5.5)

where we assume that the partial sum provides a sufficiently accurate approximation of the integral term on $G_h$ . Specifically, assuming that $k\left(\left|x-y\right|\right)=exp\left(-{\left|x-y\right|}^2\right)$ and requiring that the $k_{\left|n\right|}=O\left({10}^{-16}\right)$ (double pre- cision machine epsilon) for $\left|n\right|>L$ , one should choose $L\propto 1/h$ . However, in practice, a much smaller $L$ results in accurate LFA estimates.

Next, in (5.5), we insert (5.3) and obtain

$\begin{array}{l}\frac{2}{h^2}\underset{\theta }{{\displaystyle \sum }}{\text{E}}_{\theta }^{m+1}{\text{e}}^{i\theta j}-\frac{1}{h^2}\underset{\theta }{{\displaystyle \sum }}{\text{E}}_{\theta }^{m+1}{\text{e}}^{i\theta \left(j-1\right)}\\ =\frac{1}{h^2}\underset{\theta }{{\displaystyle \sum }}{\text{E}}_{\theta }^m{\text{e}}^{i\theta \left(j+1\right)}-h\stackrel{L}{\underset{n=-L}{{\displaystyle \sum }}}{k}_{\left|n\right|}\underset{\theta }{{\displaystyle \sum }}{\text{E}}_{\theta }^m{\text{e}}^{i\theta \left(j-n\right)}\mathrm{,}\end{array}$ (5.6)

which can be re-written as follows

$\underset{\theta }{{\displaystyle \sum }}{\text{E}}_{\theta }^{m+1}\left(\frac{2}{h^2}-\frac{1}{h^2}{\text{e}}^{-i\theta }\right){\text{e}}^{i\theta j}=\underset{\theta }{{\displaystyle \sum }}{\text{E}}_{\theta }^m\left(\frac{1}{h^2}{\text{e}}^{i\theta }-h\stackrel{L}{\underset{n=-L}{{\displaystyle \sum }}}{k}_{\left|n\right|}{\text{e}}^{-i\theta n}\right){\text{e}}^{i\theta j}\mathrm{.}$ (5.7)

Now, comparing the coefficients of equal frequency modes on both sides of (23), we obtain

${\hat{S}}_h\left(\theta \right)\mathrm{<mo>\; :\; </mo>}=\frac{{\text{E}}_{\theta }^{m+1}}{{\text{E}}_{\theta }^m}=\frac{{\text{e}}^{i\theta }-h^3\stackrel{L}{\underset{n=-L}{{\displaystyle \sum }}}{k}_{\left|n\right|}{\text{e}}^{-i\theta n}}{2-{\text{e}}^{-i\theta }}\mathrm{.}$ (5.8)

Therefore an appropriate estimate of the smoothing factor of our GSP scheme is given by

${\mu }_{\text{GSP}}=\underset{\frac{\text{π}}{2}\le \left|\theta \right|\le \text{π}}{\text{max}}\left(\left|{\hat{S}}_h\left(\theta \right)\right|\right)\mathrm{.}$ (5.9)

With this definition, we obtain a smoothing factor of our GSP scheme given by ${\mu }_{\text{GSP}}=0.447$ for all mesh sizes. This is done by inspection of the function $\left|{\hat{S}}_h\left(\theta \right)\right|$ .

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 $A_h=-{\Delta }_h+{\mathcal{I}}^h$ , applied to a generic vector $v$ with jth component given by $v_j=\underset{\theta }{{\displaystyle \sum }}{V}_{\theta }{\text{e}}^{i\theta j}$ . We have

${A_{h}v|}_j=\underset{\theta }{{\displaystyle \sum }}{V}_{\theta }\left(-\frac{{\text{e}}^{i\theta \left(j+1\right)}-2{\text{e}}^{i\theta j}+{\text{e}}^{i\theta \left(j-1\right)}}{h^2}+h\stackrel{L}{\underset{n=-L}{{\displaystyle \sum }}}{k}_{\left|n\right|}{\text{e}}^{i\theta \left(j-n\right)}\right)\mathrm{.}$

Therefore we obtain

${\hat{A}}_h\left(\theta \right)=-\frac{2\left(1-cos\left(\theta \right)\right)}{h^2}+h\stackrel{L}{\underset{n=-L}{{\displaystyle \sum }}}{k}_{\left|n\right|}{\text{e}}^{-i\theta n}\mathrm{.}$ (5.10)

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

${\hat{A}}_h\left(\theta \right)=\left[\begin{array}{cc}{\hat{A}}_h\left({\theta }^0\right)& 0\\ 0& {\hat{A}}_h\left({\theta }^1\right)\end{array}\right]\text{and}{\hat{S}}_h\left(\theta \right)=\left[\begin{array}{cc}{\hat{S}}_h\left({\theta }^0\right)& 0\\ 0& {\hat{S}}_h\left({\theta }^1\right)\end{array}\right]\mathrm{.}$

On the coarse grid, we have the following

${\hat{A}}_H\left(2\theta \right)=-\frac{2\left(1-cos\left(2\theta \right)\right)}{H^2}+H\stackrel{\frac{L}{2}}{\underset{n=-\frac{L}{2}}{{\displaystyle \sum }}}{k}_{\left|n\right|}{\text{e}}^{-i\left(2\theta \right)n}\mathrm{.}$ (5.11)

For the restriction operator, we have the following [27]

$\begin{array}{c}{\mathbb{I}}_h^H{\text{e}}^{i\theta j}=\frac{-{\text{e}}^{i\theta \left(j-3\right)}+9{\text{e}}^{i\theta \left(j-1\right)}+16{\text{e}}^{i\theta j}+9{\text{e}}^{i\theta \left(j+1\right)}-{\text{e}}^{i\theta \left(j+3\right)}}{32}\mathrm{,}\\ =\left(\frac{-{\text{e}}^{-3i\theta }+9{\text{e}}^{-i\theta }+16+9{\text{e}}^{i\theta }-{\text{e}}^{3i\theta }}{32}\right){\text{e}}^{i\theta j}\mathrm{,}\\ =\left(\frac{18cos\left(\theta \right)+16-2cos\left(3\theta \right)}{32}\right){\text{e}}^{i\theta j}\mathrm{.}\end{array}$

Hence

${\hat{\mathbb{I}}}_h^H\left(\theta \right)=\left[\begin{array}{cc}\frac{18cos\left({\theta }^0\right)+16-2cos\left(3{\theta }^0\right)}{32}& \frac{18cos\left({\theta }^1\right)+16-2cos\left(3{\theta }^1\right)}{32}\end{array}\right]\mathrm{.}$

For the interpolation operator, we obtain

$\begin{array}{c}{\mathbb{I}}_H^h{\text{e}}^{i\theta j}=\frac{-{\text{e}}^{i\theta \left(j-3\right)}+9{\text{e}}^{i\theta \left(j-1\right)}+16{\text{e}}^{i\theta \left(j\right)}+9{\text{e}}^{i\theta \left(j+1\right)}-{\text{e}}^{i\theta \left(j+3\right)}}{16}\mathrm{,}\\ =\left(\frac{18cos\left(\theta \right)+16-2cos\left(3\theta \right)}{16}\right){\text{e}}^{i\theta j}\mathrm{.}\end{array}$

Hence

${\hat{\mathbb{I}}}_H^h\left(\theta \right)=\left[\begin{array}{c}\frac{18cos\left({\theta }^0\right)+16-2cos\left(3{\theta }^0\right)}{16}\\ \frac{18cos\left({\theta }^1\right)+16-2cos\left(3{\theta }^1\right)}{16}\end{array}\right]\mathrm{.}$

Now, we are able to compute the two-grid convergence factor as follows

${\mu }_{TG}=\underset{\left|\theta \right|\le \frac{\text{π}}{2}}{\text{max}}\left\{\rho \left(\widehat{TG}{}_h^H\left(\theta \right)\right)\right\}\mathrm{,}$ (5.12)

where $\rho$ denotes the spectral radius of the $2\times 2$ matrix $\widehat{TG}{}_h^H\left(\theta \right)$ .

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, ${\nu }_1\mathrm{,}{\nu }_2$ . These values are computed by inspection of the function $\rho \left(\widehat{TG}{}_h^H\left(\theta \right)\right)$ , which is evaluated using MATLAB to compute the eigenvalues of the matrix $\widehat{TG}{}_h^H\left(\theta \right)$ on a fine grid of $\theta$ values, $-\text{π}/2\le \theta \le \text{π}/2$ .

Further in the same table, we compare these values with the value of the observed convergence factor given by

$\rho =\frac{{\Vert r_h^{m+1}\Vert }_{L_h^2}}{{\Vert r_h^m\Vert }_{L_h^2}}.$

Table 1. Estimated and observed multigrid convergence factors.

This numerical convergence factor represents the asymptotic ratio of reduction of the $L^2$ -norm of the residual between two multigrid cycles. These calculations refer to the choice $k\left(\left|x-y\right|\right)=exp\left(-{\left|x-y\right|}^2\right)$ . As shown in Table 1, the LFA estimates of the multigrid convergence factor are accurate. (The same values of ${\mu }_{TG}$ 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 $\mathcal{O}\left(MlogM\right)$ 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

$-\Delta u\left(x\mathrm{,}y\right)+\underset{\Omega }{{\displaystyle \int }}\underset{\Omega }{{\displaystyle \int }}k\left(x\mathrm{,}y\mathrm{,}t\mathrm{,}s\right)u\left(t\mathrm{,}s\right)\text{d}t\text{d}s=f\left(x\mathrm{,}y\right)\mathrm{,}$ (6.1)

where $\Omega =\left(-1,1\right)\times \left(-1,1\right)$ and the Gaussian kernel

$k\left(x\mathrm{,}y\mathrm{,}t\mathrm{,}s\right)=exp\left(-\frac{{\left(x-t\right)}^2+{\left(y-s\right)}^2}{2}\right)$ .

To investigate the order of accuracy of the discretization scheme, we construct

an exact solution to (6.1) by choosing $u\left(x,y\right)=\exp \left(-\frac{x^2+y^2}{2}\right)$ . With this

choice, the right-hand side of (6.1) is given by

$\begin{array}{c}f\left(x\mathrm{,}y\right)=2exp\left(-\frac{x^2+y^2}{2}\right)-\left(x^2+y^2\right)ex{p}^{\left(-\frac{x^2+y^2}{2}\right)}+\frac{\text{π}}{4}ex{p}^{\left(-\left(\frac{x^2}{4}+\frac{y^2}{4}\right)\right)}\\ \times \left[\text{erf}\left(1+\frac{y}{2}\right)+\text{erf}\left(\frac{2-y}{2}\right)\right]\left[\text{erf}\left(1-\frac{x}{2}\right)+\text{erf}\left(1+\frac{x}{2}\right)\right]\mathrm{.}\end{array}$

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 ${\Vert {\mathcal{I}}^{h}u-\mathcal{I}u\Vert }_h$ . Notice that we can evaluate $\mathcal{I}u$ exactly, while ${\mathcal{I}}^{h}u$ is computed using the full-kernel (FK) evaluation formula (6) and the FI technique involving different depths, $d=l-k$ . For $d=0$ the FI scheme performs FK evaluation.

In Table 3, we report the values of ${\Vert {\mathcal{I}}^{h}u-\mathcal{I}u\Vert }_h$ corresponding to different 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

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

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 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 $\mathcal{O}\left(MlogM\right)$ and $M=N^2$ 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, $\nu =1,\cdots ,10$ , $\nu ={\nu }_1+{\nu }_2$ , 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

$-\Delta u\left(x\right)+\underset{\Omega }{{\displaystyle \int }}log\left|x-y\right|u\left(y\right)\text{d}y=f\left(x\right)$

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; $M=N^2$ .

where $f\left(x\right)=1$ for $x\in \Omega :=\left(-1,1\right)$ . 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 $x_i$ of the grid, we

use local averaging $k\left(x_i\right)\approx \frac{1}{2}\left(k\left(x_{i-\frac{1}{2}}\right)+k\left(x_{i+\frac{1}{2}}\right)\right)$ . 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 $\mathcal{O}\left(MlogM\right)$ . Further, in Figure 4 the convergence history of the

Figure 2. Convergence history of the FAS-FI scheme with different $\nu ={\nu }_1+{\nu }_2$ , $\nu =1$ (green) to $\nu =10$ (red) along 5 V-cycles of FAS; $l=8,d=3$ .

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.

7. Conclusion

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.

Acknowledgements

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.

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	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. https://doi.org/10.1016/j.jtbi.2006.05.030
[2]	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. https://doi.org/10.1016/j.amc.2014.08.002
[3]	Hirsa, A. and Neftci, S.N. (2013) An Introduction to the Mathematics of Financial Derivatives. Academic Press, Cambridge.
[4]	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. https://doi.org/10.21314/JCF.2016.208
[5]	Yaser, R. and Khosrow, M. (2017) Numerical Solution of Partial Integro-Differential Equations by Using Projection Method Mediterranean. Journal of Mathematics, 14, 113.
[6]	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. https://doi.org/10.5923/j.ajcam.20120203.06
[7]	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. https://doi.org/10.4208/jcm.2009.09-m2852
[8]	Zhao, J. and Corless, R.M. (2006) Compact Finite Difference Method for Integro-Differential Equations. Applied Mathematics and Computation, 177, 271-288. https://doi.org/10.1016/j.amc.2005.11.007
[9]	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. https://doi.org/10.12785/msl/010109
[10]	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. https://doi.org/10.1016/0377-0427(95)00273-1
[11]	Volk, W. (1988) The Iterated Garlekin Method for Linear Integro-Differential Equations. Journal of Computational and Applied Mathematics, 21, 63-74. https://doi.org/10.1016/0377-0427(88)90388-3
[12]	Hu, Q. (1998) Interpolation Correction for Collocation Solutions of Fredholm Integro-Differential Equations. Mathematics of Computation, 67, 987-999. https://doi.org/10.1090/S0025-5718-98-00956-9
[13]	Volk, W. (1985) The Numerical Solutions of Linear Integro-Differential Equations by Projection Methods. Journal of Integral Equations, 9, 171-190.
[14]	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. https://doi.org/10.1016/S0307-904X(02)00099-9
[15]	Brandt, A. and Lubrecht, A.A. (1990) Multilevel Matrix Multiplication and Fast Solution of Integral Equations. Journal of Computational Physics, 90, 348-370. https://doi.org/10.1016/0021-9991(90)90171-V
[16]	Heli, C.E. (2008) Real Analysis II. Georgia Institute of Technology, Atlanta.
[17]	Hackbusch, W. (1993) Elliptic Differential Equations Theory and Numerical Treatment. Journal of Materials Science, 28, 5831-5835.
[18]	Stoer, J. and Bulirsch, R. (2002) Introduction to Numerical Analysis. Texts in Applied Mathematics, 79, 243. https://doi.org/10.1007/978-0-387-21738-3
[19]	Jovanovic, B.S. and Süli, E. (2014) Analysis of Finite Difference Schemes. Springer Series in Computational Mathematics, 46, 198. https://doi.org/10.1007/978-1-4471-5460-0
[20]	Brandt, A. (1977) Multi-Level Adaptive Solutions to Boundary-Value Problems. Mathematics of Computation, 31, 333-390. https://doi.org/10.1090/S0025-5718-1977-0431719-X
[21]	Trottenberg, U., Osterlee, C. and Schüller, A. (2001) Multigrid. Elsevier Academic Press, Amsterdam.
[22]	Wesseling, P. (1991) An Introduction to Multigrid Methods. John Wiley and Sons, New York.
[23]	Hemker, P.W. (1990) On the Order of Prolongations and Restrictions in Multigrid Procedures. Journal of Computational and Applied Mathematics, 32, 423-429. https://doi.org/10.1016/0377-0427(90)90047-4
[24]	Hackbusch, W. (1985) Multi-Grid Methods and Applications. Springer, Berlin. https://doi.org/10.1007/978-3-662-02427-0
[25]	Lubrecht, T. and Venner, C.H. (2000) Multi-Level Methods in Lubrication. Elsevier, Amsterdam.
[26]	Brandt, A. (1982) Guide to Multi Grid Development. The Weizmann Institute of Science, Rehovot.
[27]	Wienands, R. and Joppich, W. (1994) Practical Fourier Analysis for Multigrid methods. Chapman and Hall/CRC Press, Boca Raton.