Numerical Methods for a Class of Nonlocal Diffusion Problems with the Use of Backward SDEs This material is based upon work supported in part by the U.S. Air Force of Scientific Research under grant numbers FA9550-11-1-0149 and 1854-V521-12; by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research, Applied Mathematics program under contract numbers ERKJ259, and ERKJE45; by the National Natural Science Foundations of China under grant numbers 91130003 and 11171189; by Natural Science Foundation of Shandong Province under grant number ZR2011AZ002; and by the Laboratory Directed Research and Development program at the Oak Ridge National Laboratory, which is operated by UT-Battelle, LLC, for the U.S. Department of Energy under Contract DE-AC05-00OR22725.

# Numerical Methods for a Class of Nonlocal Diffusion Problems with the Use of Backward SDEs ††thanks: This material is based upon work supported in part by the U.S. Air Force of Scientific Research under grant numbers FA9550-11-1-0149 and 1854-V521-12; by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research, Applied Mathematics program under contract numbers ERKJ259, and ERKJE45; by the National Natural Science Foundations of China under grant numbers 91130003 and 11171189; by Natural Science Foundation of Shandong Province under grant number ZR2011AZ002; and by the Laboratory Directed Research and Development program at the Oak Ridge National Laboratory, which is operated by UT-Battelle, LLC, for the U.S. Department of Energy under Contract DE-AC05-00OR22725.

G. Zhang Department of Computational and Applied Mathematics, Oak Ridge National Laboratory, Oak Ridge, TN 37831(zhangg@ornl.gov, webstercg@ornl.gov).    W. Zhao School of Mathematics, Shandong University, Jinan 250100, China (wdzhao@sdu.edu.cn).    C. G. Webster22footnotemark: 2    M. Gunzburger Department of Scientific Computing, Florida State University (gunzburg@fsu.edu).
###### Abstract

We propose a novel numerical approach for nonlocal diffusion equations [8] with integrable kernels, based on the relationship between the backward Kolmogorov equation and backward stochastic differential equations (BSDEs) driven by Lèvy processes with jumps. The nonlocal diffusion problem under consideration is converted to a BSDE, for which numerical schemes are developed and applied directly. As a stochastic approach, the proposed method does not require the solution of linear systems, which allows for embarrassingly parallel implementations and also enables adaptive approximation techniques to be incorporated in a straightforward fashion. Moreover, our method is more accurate than classic stochastic approaches due to the use of high-order temporal and spatial discretization schemes. In addition, our approach can handle a broad class of problems with general nonlinear forcing terms as long as they are globally Lipchitz continuous. Rigorous error analysis of the new method is provided as several numerical examples that illustrate the effectiveness and efficiency of the proposed approach.

Key words. backward stochastic differential equation with jumps, nonlocal diffusion equations, superdiffusion, compound Poisson process, -scheme, adaptive approximation

## 1 Introduction

A diffusion process is deemed nonlocal when the associated underlying Lèvy process does not only consist of Brownian motions. The features of nonlocal diffusion has been verified experimentally to be present in a wide variety of applications, including, e.g., contaminant flow in groundwater, plasma physics, and the dynamics of financial markets. A comprehensive survey of nonlocal diffusion problems is given in [15]. In this work, we consider a partial-integral diffusion equation (PIDE) representation [6, 11, 8] of linear nonlocal diffusion problems. Since it is typically difficult to obtain analytical solutions of such problems, numerical solutions are highly desired in practical applications. There are mainly two types of numerical methods for nonlocal diffusion equations under consideration. The first is the family of deterministic approaches, e.g., finite-element-type methods [6, 11] and collocation methods, whereas the second can be classified as stochastic approaches, e.g., continuous-time random walk (CTRW) methods [5, 10, 15].

The deterministic approaches are extensions of existing methods for local partial differential equations (PDEs) that incorporate schemes to discretize the underlying integral operators. The recently developed nonlocal vector calculus [9] provides helpful tools that allow one to analyze nonlocal diffusion problems in a similar manner to analyzing local PDEs; see [8, 11] for details. However, in the context of implicit time-stepping schemes, the nonlocal operator may result in severe computational difficulties coming from the dramatic deterioration of the sparsity of the stiffness matrices required by the underlying linear systems. On the other hand, stochastic approaches, e.g., CTRW methods, are based on the relation between nonlocal diffusion and a class of Lèvy jump processes. Although stochastic methods do not require the solution of linear systems, and the simulations of all paths of the random walk can be easily parallelized, they also have some drawbacks. Typically, most stochastic methods are sampling-based Monte Carlo approaches, suffering from slow convergence, thus requiring a very large number of samples to achieve small errors. Even worse, in the case of a general nonlinear forcing term, the nonlocal diffusion equation is no longer the master equation of the underlying jump processes.

In this paper, we propose a novel stochastic numerical scheme for the linear nonlocal diffusion problems studied in [5, 6, 10, 9, 11, 8] based on the relationship between the PIDEs and a certain class of backward stochastic differential equations (BSDEs) with jumps. The existence and uniqueness of solutions for nonlinear BSDEs with and without jumps have been proved in [17] and [1], respectively. Since then, BSDEs have become important tools in probability theory, stochastic optimal control and mathematical finance [7]. Unlike BSDEs driven by Brownian motions, there are very few numerical schemes proposed for BSDEs with jumps, and most of the schemes are focussed on time discretizations only. For example, a numerical scheme based on Picard’s method was proposed in [14], and a forward Euler scheme was proposed in  [2, 3] where the convergence rate was proved to be .

Our focus will be on the BSDEs corresponding to nonlocal diffusion equations with integrable jump kernels in unbounded domains. Since the PIDEs of interest do not possess local convection and diffusion terms, the corresponding BSDEs can be simplified, such that, the underlying Lèvy processes are merely compound Poisson processes. Also, motivated by various engineering applications discussed in [8], our goal is to develop and analyze computationally efficient numerical schemes that can achieve high-order accuracy in both temporal and spatial discretization, rather than focus on the scalability of our schemes in solving extremely high-dimensional problems (e.g., in option pricing). From this perspective, although the BSDEs under consideration is a special case of the equations studied in [2], high-order fully-discrete schemes and the corresponding error analysis are still missing in the literature. In fact, based on our previous work, it is particularly challenging to recover the classic finite element convergence rates in the spatial discretization, even for the simplified BSDEs. Thus, the main contributions of this paper are summarized as follows:

• Development of stable high-order matrix-free numerical schemes for the BSDEs and the corresponding nonlocal diffusion problems;

• Analysis of the approximation error of the proposed semi-discrete and fully-discrete schemes, and;

• Numerical demonstration of the capabilities of handling multiple types of integrable kernels, e.g., symmetric, non-symmetric, or singular kernels.

Specifically, the BSDEs will be discretized, in the temporal domain, using a -scheme extended from our previous works [21, 20] for BSDEs without jumps. In particular, the cases where , and correspond to forward Euler, Crank-Nicolson and backward Euler schemes, respectively. As discussed in [21], a quadrature rule adapted to the underlying stochastic processes is critical to approximate all the conditional mathematical expectations in our numerical scheme. Thus, we propose a general formulation of the quadrature rule for estimating the conditional expectations with respect to the compound Poisson processes, and a specific form of a certain PIDE can be determined based on regularities of the kernel and forcing term. In §LABEL:sec:ex, Gauss-Legendre, Gauss-Jacobi and Newton-Cotes rules are substituted into the proposed quadrature rule to approximate different kernels and forcing terms. In addition, since the total number of quadrature points grows exponentially with the number of time steps, we construct a piecewise Lagrange interpolation scheme based on a pre-determined spatial mesh, in order to evaluate the integrand of the expectations at all quadrature points.

Both theoretical analysis and numerical experiments show that our approach is advantageous for linear nonlocal diffusion equations with integrable kernels on unbounded domains. Compared to deterministic approaches, e.g., finite elements and collocation approaches, the proposed methods do not require the solution of dense linear systems. Instead, the PIDE can be solved independently at different spatial grid points on each time level, making it straightforward to incorporate parallel implementation and adaptive spatial approximation. Compared to stochastic approaches, our scheme is more accurate than the CTRW method due to the use of the -scheme for time discretization, the high-order quadrature rule, and piecewise Lagrange polynomial interpolation for spatial discretization. In addition, our method can handle a broad class of problems with general nonlinear forcing terms as long as they are globally Lipchitz continuous.

The outline of the paper is as follows. In §LABEL:sec:nlocBSDE, we introduce the mathematical description of the linear nonlocal diffusion equations and the corresponding class of BSDEs under consideration. In §LABEL:sec:scheme, we propose our numerical schemes for the BSDEs of interest, where the semi-discrete and the fully-discrete approximations are presented in §LABEL:sec:semi and §LABEL:sec:fully, respectively. Error estimates for the proposed scheme are proved in §LABEL:sec:err. Extensive numerical examples and comparisons to existing techniques are given in §LABEL:sec:ex, which are shown to be consistent with the theoretical results and reveal the effectiveness of our approach. Finally, several concluding remarks are given in §LABEL:sec:con.

## 2 Nonlocal diffusion models and BSDEs with jumps

This section is dedicated to describing the nonlocal diffusion problem that is the focus of this paper. In particular, in §LABEL:sec:nloc1 we introduce the definitions of a general nonlocal operator and diffusion equation. The relationship between the backward Kolmogorov equation, which is a generalization of the nonlocal diffusion equation of interest, and the BSDE with jumps, is described in §LABEL:sec:BSDE. Based on this equivalence, a simplified BSDE corresponding to the nonlocal diffusion problem of §LABEL:sec:nloc1 is also given at the end of §LABEL:sec:BSDE.

### 2.1 Nonlocal diffusion equations

Let us recall a nonlocal operator introduced in [8]. For a function , with and , we define the action of the linear operator on as

 Lu=∫Rd(u(t,x+e)−u(t,x))γ(e)de,∀(t,x)∈[0,T]×Rd, \hb@xt@.01(2.1)

where the properties of depend crucially on the kernel function . In this work, we focus on a particular class of kernel functions which are nonnegative and integrable, i.e.,

 γ(e)≥0∀e∈Rd and ∫Rdγ(e)de<∞. \hb@xt@.01(2.2)

Note that may be symmetric, i.e., for any , or non-symmetric, i.e., there exists such that . The operator exhibits nonlocal behavior because, for a fixed , the value of at a point requires information about at points ; this is contrasted with local operators, e.g., the value of at a point requires information about only in an infinitesimal neighborhood of . Our interest in the operator is due to its participation in the nonlocal diffusion equation

 ⎧⎪⎨⎪⎩∂u∂t(t,x)−Lu(t,x)=g(t,x,u) for (t,x)∈(0,T]×Rdu(0,x)=u0(x), for x∈Rd, \hb@xt@.01(2.3)

where is the initial condition and the forcing term may be a nonlinear function of and . We remark that (LABEL:nloc_diff) is a nonlocal Cauchy problem defined on the unbounded spatial domain . Also, in the context of kernel functions that are compactly supported in , the initial-boundary value nonlocal problem with volume constraints, i.e., constraints applied on a regions of nonzero measure in , has been well studied [6, 8, 11]. However, the well-posedness of the corresponding BSDEs with volume constraints is still an open challenge. As such, we do not impose volume constraints to (LABEL:nloc_diff) in this effort.

### 2.2 BSDEs and backward Kolmogorov equations

Now we discuss the probabilistic representation of the solution of the nonlocal diffusion equation in (LABEL:nloc_diff). The nonlinear Feynman-Kac formula studied in [16] shows that BSDEs driven by Brownian motion provides a probabilistic representation of the solutions of a class of second-order quasi-linear parabolic PDEs. This result was extended in [1] to the case of BSDEs driven by general Lévy processes with jumps. Due to the inclusion of Lèvy jumps, the counterpart of a BSDE is a partial integral differential equation (PIDE), i.e., the backward Kolmogorov equation [13]. It turns out that such a PIDE is a generalization of the nonlocal diffusion equation in (LABEL:nloc_diff). Here we first introduce a general form of a BSDE with jumps and the backward Kolmogorov equation, then give a simplified form of the BSDE corresponding to the nonlocal diffusion equation (LABEL:nloc_diff) under consideration.

Let be a probability space with a filtration for a finite terminal time . We assume the filtration satisfies the usual hypotheses of completeness, i.e., contains all sets of -measure zero and maintains right continuity, i.e., . Moreover, the filtration is assumed to be generated by two mutually independent processes, i.e., one -dimensional Brownian motion and one -dimensional Poisson random measure on , where is equipped with its Borel field . The compensator of and the resulting compensated Poisson random measure are denoted by and respectively, such that is a martingale for all . We also assume that is a -finite measure on satisfying

 ∫E(1∧|e|2)λ(de)<+∞. \hb@xt@.01(2.4)

Based on the stochastic basis , we introduce the backward stochastic differential equation with jumps

 ⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩Xt=X0+∫t0b(s,Xs)ds+∫t0σ(s,Xs)dWs+∫t0∫Eβ(s,Xs−,e)~μ(de,ds)Yt=ξ+∫Ttf(s,Xs,Ys,Zs,Γs)ds−∫TtZsdWs−∫Tt∫EUs(e)~μ(de,ds), \hb@xt@.01(2.5)

where the solution is the quadruplet with , , and . The map is referred to as the drift coefficient, is the local diffusion coefficient, is the jump coefficient, is the generator of the BSDE, and the processes is defined by for a given bounded function . The terminal condition , which is an -measurable random vector, is assumed to be a function of , denoted by . Note that the integrals in (LABEL:eq:BSDE) with respect to the -dimensional Brownian motion and the -dimensional compensated Poisson measure are Itô-type stochastic integrals. The first equation in (LABEL:eq:BSDE) is the forward SDE system, which is a general Lévy process, while the second equation is a system of BSDEs driven by . The quadruplet is called an -adapted solution if it is -adapted, square integrable process which satisfies the BSDE in (LABEL:eq:BSDE). Under standard assumptions on the given functions , , , and , the existence and uniqueness of the solution of the system in (LABEL:eq:BSDE) with a nonlinear generator have been proved in [1].

Based on the extension of the nonlinear Feynman-Kac formula for BSDEs proposed in [1], the adapted solution of (LABEL:eq:BSDE) can be associated to the unique viscosity solution of the backward Kolmogorov equation, i.e.,

 ⎧⎪⎨⎪⎩∂v∂t(t,x)+˜Lv(t,x)+f(t,x,v,σ∇v,Bv)=0, for (t,x)∈[0,T)×Rdv(T,x)=φ(x), for x∈Rd. \hb@xt@.01(2.6)

In (LABEL:eq:Kol), is the terminal condition at the time , and the second-order integral-differential operator is of the form

 ˜Lv(t,x)=d∑i=1bi(t,x)∂v∂xi(t,x)+12d∑i,j=1(σσ⊤)i,j(t,x)∂2v∂xi∂xj(t,x) \hb@xt@.01(2.7) +∫E(v(t,x+β(t,x,e))−v(t,x)−d∑i=1∂v∂xi(t,x)β(t,x,e))λ(de),

with is an integral operator defined as

 Bv(t,x)=∫E[v(t,x+β(t,x,e))−v(t,x)]η(e)λ(de).

The functions , , , , in (LABEL:eq:L1) have the same definitions as in the BSDE in (LABEL:eq:BSDE). Under the condition that for a fixed , the solution of the BSDE for can be represented by

 ⎧⎪⎨⎪⎩Yt=v(t,Xt),Zt=σ(t,Xt)∇v(t,Xt),Ut=v(t,Xt−+β(t,Xt−,e))−v(t,Xt−), \hb@xt@.01(2.8)

where denotes the gradient of with respect to and . By comparing the operator in (LABEL:eq:L1) and the BSDE in (LABEL:eq:BSDE), we can see that the incremental consists of three components: the drift term , the Brownian diffusion , and the jump diffusion .

Now let us relate the BSDE in (LABEL:eq:BSDE) to the nonlocal diffusion equation in (LABEL:nloc_diff). Without loss of generality, we only consider the case when is a scalar function, i.e., by setting in (LABEL:eq:BSDE); the numerical schemes and analysis presented in this paper can be directly extended to the case of . Observing that the operator in (LABEL:eq:L1) is a generalization of the nonlocal operator in (LABEL:nloc_diff), we aim to find a simplified form of the BSDE which is the stochastic representation of the nonlocal diffusion equation in (LABEL:nloc_diff). To proceed, we define for and for , satisfying the condition in (LABEL:lambda). Due to the integrability assumption of in (LABEL:gamma), the jump diffusion term in (LABEL:eq:BSDE) is a compensated compound Poisson process. In this case, the compensated Poisson random measure can be represented by

 ~μ(de,dt)=μ(de,dt)−λρ(e)dedt, \hb@xt@.01(2.9)

where , with being the intensity of Poisson jumps and being the probability measure of the jump size satisfying . Then, the operator in (LABEL:eq:L1) can be simplified to the nonlocal operator in (LABEL:eq:L) by setting , i.e., removing the Brownian motion, and

 σ=0,β=eID(e),λ=∫Eγ(e)de, \hb@xt@.01(2.10) ρ(e)=1λγ(e),bi=∫Eeγ(e)de for i=1,…,d,

where is the characteristic function of the domain . In the context of , the drift coefficient is a constant vector which indeed helps cancel the first-order derivatives of in . On the other hand, substituting into the definition of , we have

 ∫t0bids=λ∫t0∫Eeρ(e)deds=λtE[e] for i=1,…,d, \hb@xt@.01(2.11)

which is the compensator of . Hence, is just a compound Poisson process under the Poisson measure without compensation. Then, by defining and , with and being the forcing term and initial condition in (LABEL:nloc_diff) respectively, we obtain the BSDE corresponding to the nonlocal diffusion equation in (LABEL:nloc_diff)

 \hb@xt@.01(2.12)

According to Theorem 2.1 in [1], the well-posedness of the BSDE in (LABEL:eq:BSDE2) requires the generator is globally Lipschitz continuous, i.e., there exists such that

 |f(t,x,y)−f(t,x′,y′)|≤K(|x−x′|+|y−y′|)

for all with and . In this case, there exists a unique solution where is the set of -adapted càdlàg processes such that

 ∥Yt∥2S2=E⎡⎣(sup0≤t≤T|Yt|)2⎤⎦<∞,

and is the set of mappings such that

 ∥Ut∥2L2(~μ)=E[∫T0∫EUt(e)2λ(de)dt]<∞.

To relate the solution to the viscosity solution of the nonlocal diffusion equation in (LABEL:nloc_diff), we assume and are continuous functions satisfying

 |f(t,x,0)|≤C(1+|x|p)and|φ(x)|≤C(1+|x|p),(t,x)∈[0,T]×Rd,

for some real . Then, under the condition that for a fixed , the solution for can be represented by

 Yt=u(T−t,Xt) and Ut=u(T−t,Xt−+e)−u(T−t,Xt−),

where is the viscosity solution of (LABEL:nloc_diff) satisfying

 |u(t,x)|≤C(1+|x|p),(t,x)∈[0,T]×Rd,

again, for some real . Moreover, if and are bounded and uniformly continuous, then is also bounded and uniformly continuous.

Now we introduce the following notations which will be used throughout. Let be a -field generated by the stochastic process starting from the time-space point , and set . Denote by the mathematical expectation and the conditional mathematical expectation under the -field , i.e., . Particularly, for the sake of notational convenience, when , we use to denote , i.e., the mathematical expectation under the condition that .

## 3 Numerical schemes for BSDEs with jumps

In this section we propose numerical schemes for the BSDE in (LABEL:eq:BSDE2) in order to solve the nonlocal diffusion equation in (LABEL:nloc_diff). Specifically, a semi-discrete scheme for time discretization is studied in §LABEL:sec:semi, and a fully-discrete scheme is constructed in §LABEL:sec:fully by incorporating appropriate spatial discretization techniques.

### 3.1 The semi-discrete scheme

For the time interval , we introduce the partition

 Tτ={0=t0<⋯

with and . In the time interval for , under the condition that , the BSDE in (LABEL:eq:BSDE2) can be rewritten as

 Xs=x+∫stn∫Eeμ(de,dr) for tn≤s≤tn+1, \hb@xt@.01(3.2)
 Ytn=Ytn+1+∫tn+1tnf(s,Xs,Ys)ds−∫tn+1tn∫EUs(e)~μ(de,ds). \hb@xt@.01(3.3)

Taking the conditional mathematical expectation on both sides of (LABEL:s3:e1_2), due to the fact that for is a martingale, we obtain that

 Ytn=Extn[Ytn+1]+∫tn+1tnExtn[f(s,Xs,Ys)]ds, \hb@xt@.01(3.4)

where , and the integrand is a deterministic function of . To estimate the integral in (LABEL:s3:e2), several numerical methods for approximating integrals with deterministic integrands can be used. Here, we utilize the -scheme proposed in [20, 21], which yields

 Ytn=Extn[Ytn+1] +(1−θ)ΔtnExtn[f(tn+1,Xtn+1,Ytn+1)] \hb@xt@.01(3.5) +θΔtnf(tn,Xtn,Ytn)+Rn,

where and the residual is given by

 Rn=∫tn+1tn {Extn[f(s,Xs,Ys)]−θf(tn,Xtn,Ytn) \hb@xt@.01(3.6) −(1−θ)Extn[f(tn+1,Xtn+1,Ytn+1)]}ds.

Next we propose a semi-discrete scheme for the BSDE in (LABEL:eq:BSDE2) as follows: given the random variable as the terminal condition, for , under the condition that , the solution is approximated by satisfying

 ⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩Xn+1=Xn+∫Eeμ(de,Δt),Yn=Extn[Yn+1]+(1−θ)ΔtnExtn[f(tn+1,Xn+1,Yn+1)]+θΔtnf(tn,Xn,Yn), \hb@xt@.01(3.7)

where . By choosing different values for , we obtain different semi-discrete schemes, e.g., , and lead to forward Euler (FE), backward Euler (BE) and Crank-Nicolson (CN) schemes, respectively.

Note that since is the standard compound Poisson process, the probability distributions of or any incremental , with , are well known. Thus, one does not need to discretize so that, in (LABEL:eq:semi), denotes the exact solution evaluated at . In this case, in (LABEL:eq:Ry) is the local truncation error of the scheme (LABEL:eq:semi); a rigorous error analysis of is given in §LABEL:sec:err. Other time-stepping schemes, e.g., linear multi-step schemes [22], can be directly used in order to further improve the accuracy of time discretization. Note that in the general case where the coefficients and in (LABEL:eq:BSDE) and (LABEL:eq:L1) are functions of and , another time-stepping scheme [18] is also needed for to discretize , so that an extra local truncation error will be introduced into the semi-discrete scheme (LABEL:eq:semi). This is beyond the scope of this effort but will be the focus of future work.

### 3.2 The fully-discrete scheme

The semi-discrete scheme (LABEL:eq:semi) is not computationally feasible, because the involved conditional mathematical expectation is defined over the whole real space . Thus, an effective spatial discretization approach is also necessary in order to approximate . To proceed, we partition the -dimensional Euclidean space by

 Sdh=S1h1×S2h2×⋯×Sdhd,

where for is a partition of the one-dimensional real space , i.e.,

 Skhk={xki∣∣xki∈R,i∈Z,xki

such that and . For each multi-index , the corresponding grid point in is denoted by .

Now we turn to constructing a quadrature rule for approximating in (LABEL:eq:semi) for . Such expectation is defined with respect to the probability distribution of the incremental stochastic process , which is a compound Poisson process starting from the grid point . Here we take as an example and the proposed quadrature rule can be directly extended to approximate . It is well known that the number of jumps of within , denoted by , follows a Poisson distribution with intensity , and the size of each Poisson jump follows the distribution defined in (LABEL:eq:coeff). Thus, can be represented by

 Exitn[Yn+1] \hb@xt@.01(3.8) =∞∑m=0exp(−λΔtn)(λΔtn)mm!E[Yn+1(xi+m∑k=1ek)] =exp(−λΔtn)Yn+1(xi)+∞∑m=1exp(−λΔtn)(λΔtn)mm! ×∫E⋯∫EYn+1(xi+m∑k=1ek)(m∏k=1ρ(ek))de1⋯dem,

where for is the size of the -th jump. Observing that the probability of having jumps within is of order , the conditional mathematical expectation can be approximated by a truncation of (LABEL:eq:EY), i.e., the sum of a finite sequence, where the number of terms retained is determined according to the prescribed accuracy. For example, if we take in (LABEL:eq:semi), we expect a second-order convergence from the time discretization, i.e., the local truncation error in (LABEL:s3:e3) should be of order . In this case, assuming the intensity value is of order , the first three terms should be retained in (LABEL:eq:EY) in order to guarantee the error from the truncation of (LABEL:eq:EY) matches the order of . Hence, in the sequel, we denote by the approximation of by retaining the first terms, where indicates the number of jumps included in . An analogous notation is used to represent the approximation of by retaining jumps.

Next, we also need to approximate the -dimensional integral in (LABEL:eq:EY) for . This can be accomplished by selecting an appropriate quadrature rule base on the properties of and the smoothness of . A straightforward choice is to utilize Monte Carlo methods by drawing samples from , but this is inefficient because of slow convergence. To enhance the convergence rate, an alternative way is to use the tensor product of high-order one-dimensional quadrature rules, e.g., Newton-Cotes, Gaussian, etc.. In addition, for multi-dimensional problems, i.e., , by requiring second-order convergence , i.e., , sparse-grid quadrature rules [4, 12] can be applied to further improve the computational efficiency. Without loss of generality, for , we denote by and the set of weights and points, respectively, of the selected quadrature rule for estimating the -th integral in (LABEL:eq:EY). Then, the approximation of , denoted by is given by

 ˆExitn,My[Yn+1] =exp(−λΔt)Yn+1(xi) \hb@xt@.01(3.9)

Note that it is possible that the quadrature points in (LABEL:appro_EY) do not belong to the spatial grid . The simplest approach to manage this difficult is to add all quadrature points to the spatial grid at each time stage. However, this will result in an exponential growth of the total number of grid points as the number of time steps increases. Thus, we follow the same strategy as in [20, 21, 22], and construct piecewise Lagrange interpolating polynomials based on to approximate at the non-grid points. Specifically, at any given point , is approximated by

 Yn+1(x)≈Yn+1,p(x)=p+1∑j1=1⋯p+1∑jd=1(Yij1,…,ijdn+1,pd∏k=1∏1≤j≤p+1j≠jkxk−xkijxkijk−xkij),

where is a -th order tensor-product Lagrange interpolating polynomial. For , the interpolation points are the closest neighboring points of such that , for and , constitute a local tensor-product sub-grid around . The fully-discrete solution at the interpolation point is given by . Note that does not interpolate the semi-discrete solution due to the error . Therefore, the fully-discrete scheme of the BSDE in (LABEL:eq:BSDE2) is described as follows: given the random variable for , and for , solve the quantities for such that

 \hb@xt@.01(3.10)

where and the non-negative integers , indicate the Poisson jumps included in the approximations of and , respectively.

Note that, by using the quadrature rule, the approximation in (LABEL:appro_EY) is defined in a bounded domain for a fixed . As such, the approximate solution in any bounded subdomain of can be obtained by solving the fully-discrete scheme (LABEL:eq:full) in a bounded subset of . We also observe from (LABEL:eq:full) that at each point , the computation of only depends on and even though an implicit time-stepping scheme () is used. This means on each time stage can be computed independently, so that the difficulty of solving linear systems with possibly dense matrices is completely avoided. Instead, can be either computed explicitly (in case is linearly dependent on ), or obtained by solving a nonlinear equation (in case is nonlinear with respect to ). This feature makes it straightforward to develop massively parallel algorithms, which are, in terms of scalability, very similar to the CTRW method. Moreover, the scheme (LABEL:eq:full) is expected to outperform the CRTW method because it achieves high-order accuracy (discussed in §LABEL:sec:err) from the discretization and is also capable of solving problems that include a general nonlinear forcing term . In fact, the combination of accurate approximations and scalable computations are the key advantages of our approach, compared to the existing methods for the nonlocal diffusion problem in (LABEL:nloc_diff).

###### Remark 1

The nonlinear Feynman-Kac formula converted the effect of the nonlocal operator in (LABEL:eq:L) to the behavior of the corresponding compound Poisson process, which is accurately approximated by the conditional mathematical expectations . As such, there is no need to discretize the nonlocal operator in the fully-discrete scheme (LABEL:eq:full), so that stability does not require any CFL-type restriction on the time step; this is in contrast to classic numerical schemes for which explicit schemes do require a CFL condition. Moreover, according to the analysis in [20, 22], the approximation (LABEL:eq:full) is stable as long as the semi-discrete scheme is absolutely stable.

###### Remark 2

The total computational cost of the scheme (LABEL:eq:full) is dominated by the cost of approximating at each grid point , using the formula in (LABEL:appro_EY). For example, when solving a three-dimensional problem () and retaining two Lèvy jumps (,) our approach requires the approximation of multiple six-dimensional integrals. In this case, sparse-grid quadrature rules [19] can be used to alleviate the explosion in computational cost coming from the high-dimensional systems.

## 4 Error estimates

In this section, we analyze both the semi-discrete and fully-discrete schemes, given by (LABEL:eq:semi) and (LABEL:eq:full) respectively, for the BSDE in (LABEL:eq:BSDE2). Without loss of generality, we only consider the one-dimensional case () and set ; the following analysis can be directly extended to multi-dimensional cases. For simplicity, we also assume and are both uniform grids, i.e., for and for .

Note that the well-posedness of the BSDE in (LABEL:eq:BSDE2) only requires global Lipchitz continuity on with respect to and . However, in order to obtain error estimates, we need to impose stronger regularity on . To this end, we first introduce the following notation:

 Ckb(D1×⋯×DJ) = {ψ:J∏j=1Dj→R∣∣∣∂|α|ψ∂xα is % bounded and continuous for 0≤|α|≤k where α=(α1,…,αJ) with {αj}Jj=1⊂N and |α|=α1+⋯+αJ},

where . Next we present the following lemma that relates the smoothness of to the regularity of the solution of the nonlocal diffusion equation in (LABEL:nloc_diff).

###### Lemma 1

Under assumptions in (LABEL:gamma) and (LABEL:lambda), if , in (LABEL:eq:BSDE2) and