# Fourier spectral methods for nonlocal models

Efficient and accurate spectral solvers for nonlocal models in any spatial dimension are presented. The approach we pursue is based on the Fourier multipliers of nonlocal Laplace operators introduced in a previous work. It is demonstrated that the Fourier multipliers, and the eigenvalues in particular, can be computed accurately and efficiently. This is achieved through utilizing the hypergeometric representation of the Fourier multipliers in which their computation in n dimensions reduces to the computation of a 1D smooth function given in terms of _2F_3. We use this representation to develop spectral techniques to solve periodic nonlocal time-dependent problems. For linear problems, such as the nonlocal diffusion and nonlocal wave equations, we use the diagonalizability of the nonlocal operators to produce a semi-analytic approach. For nonlinear problems, we present a pseudo-spectral method and apply it to solve a Brusselator model with nonlocal diffusion. Accuracy and efficiency of the spectral solvers are compared against a finite-difference solver.

## Authors

• 2 publications
• 3 publications
06/10/2022

### Linear peridynamics Fourier multipliers and eigenvalues

A characterization for the Fourier multipliers and eigenvalues of linear...
10/22/2020

### An efficient spectral-Galerkin method for fractional reaction-diffusion equations in unbounded domains

In this work, we apply a fast and accurate numerical method for solving ...
05/30/2021

### Computation of Eigenvalues for Nonlocal Models by Spectral Methods

The purpose of this work is to study spectral methods to approximate the...
11/02/2021

### FC-based shock-dynamics solver with neural-network localized artificial-viscosity assignment

This paper presents a spectral scheme for the numerical solution of nonl...
01/03/2020

### Tuning Multigrid Methods with Robust Optimization

Local Fourier analysis is a useful tool for predicting and analyzing the...
09/01/2021

### Analysis of pseudo-spectral methods used for numerical simulation of turbulence

Global spectral analysis (GSA) is used as a tool to test the accuracy of...
09/30/2020

### A frequency-dependent p-adaptive technique for spectral methods

When using spectral methods, a question arises as how to determine the e...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

This work concerns Fourier spectral methods for nonlocal equations involving peridynamic-type nonlocal Laplace operators. The focus is on nonlocal operators , parametrized by a spatial nonlocality parameter and an integral kernel exponent , of the form [3]

 Lδ,βu(x)=cδ,β∫Bδ(x)u(y)−u(x)∥y−x∥βdy, (1)

where is a ball in and the scaling constant is defined by

 cδ,β := (12n∫Bδ(0)1∥z∥β−2dz)−1, (2) = 2(n+2−β)Γ(n2+1)πn/2δn+2−β.

Peridynamic operators of the form (1) have been used in different applications including nonlocal diffusion, digital image correlation, and nonlocal wave phenomena, see for example [4, 6, 14, 18, 15, 16]. Nonlocal equations that involve Laplace-type operators have been addressed in several mathematical and numerical studies including [17, 8, 2, 1]. These nonlocal Lapalace operators are motivated by the peridynamic theory for continuum mechanics [19, 20, 5]

and were first introduced in nonlocal vector calculus

[9].

Spectral methods for nonlocal equations have been developed in [10, 11, 21]. The work in [10] focused on spectral approximations for a 1-D nonlocal Allen–Cahn (NAC) equation with periodic boundary conditions, and showed that these numerical methods are asymptotically compatible in the sense that they provide convergent approximations to both nonlocal and local NAC models. As pointed out in [10, 11, 21], a main challenge in developing spectral methods for nonlocal equations is to accurately and efficiently compute the eigenvalues of nonlocal operators. Examples of integral kernels for which the nonlocal eigenvalues are computed explicitly are presented in [10]. The work in [11] provides a more general method for the accurate computations of the eigenvalues of nonlocal operators in 1, 2, and 3 dimensions for a class of radially symmetric kernels. The method in [11] is given by a hybrid algorithm that is based on reformulating the integral representations of the eigenvalues as series expansions and as solutions to ODEs. The work in [21] presents spectral algorithms for solving nonlocal diffusion models that involve a nonlocal Laplace-Beltrami operator on the unit sphere. These algorithms are based on the diagonalizability of nonlocal diffusion operators on the sphere in the basis of spherical harmonics.

In this work, we present efficient and accurate spectral solvers for nonlocal models in dimensions. We pursue the Fourier multipliers’ approach introduced in [3]. The Fourier multipliers of are given by the integral representation

 m(ν)=cδ,β∫Bδ(0)cos(ν⋅z)−1∥z∥βdz, (3)

for . A key observation in our approach is that this -dimensional integral in (3) is realized through a 1D formula for which the multipliers satisfy

 m(r)=−r22F3(1,n+2−β2;2,n+22,n+4−β2;−14r2δ2), (4)

where . Noting that the hypergeometric function is smooth in its argument, it follows from (4) that evaluating the multipliers in reduces to computing this 1D smooth function. Numerical methods for computing hypergeometric functions such as are available in the literature, see for example [7]. We adopt the implementation of that is provided by Python’s mpmath library [12]. Efficiency and accuracy study of using mpmath to evaluate the multipliers is presented in Section 3.1. In addition, we provide a method for the accelerated evaluation of the eigenvalues of in Section 3.2.

Utilizing the computations of the eigenvalues of , we demonstrate spectral techniques to solve periodic nonlocal time-dependent problems. For linear problems, we diagonalize the nonlocal operators in the basis of trigonometric polynomials and provide semi-analytic solutions. We apply this in Section 4.1 to solve nonlocal diffusion and nonlocal wave equations in 2D. For nonlinear problems, we present a pseudo-spectral method in Section 4.2 and apply it to solve nonlocal diffusion in the Brusselator. In Section 4.3, we compare the spectral solver with a finite-difference solver for the nonlocal wave equation.

## 2 Fourier multipliers

For completeness of the presentation, we include a summary of the Fourier multipliers’ results, developed in [3], relevant to the computational spectral methods presented in Sections 3 and 4. The multipliers of

are defined through the Fourier transform by

 Lδ,βu(x)=1(2π)n∫Rnm(ν)^u(ν)eiν⋅xdν, (5)

where is given by the integral representation (3). The hypergeometric representation of the multipliers is provided by the following.

###### Theorem 1.

Let , and . The Fourier multipliers can be written in the form

 m(ν)=−∥ν∥22F3(1,n+2−β2;2,n+22,n+4−β2;−14∥ν∥2δ2). (6)

The smoothness of the function in its parameters, allows us to extend formula  (6) for all . Therefore, we utilize  (6) to extend the definition of the multipliers and, consequently, the definition of the operator acting on the class of Schwartz functions through the inverse Fourier transform in (5) to the case when (with ).

A characterization of the asymptotic behavior of the multipliers, and in particular that of the eigenvalues of , is given by the following result.

###### Theorem 2.

Let , and . Then, as ,

 m(ν)∼⎧⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪⎩−2n(n+2−β)δ2(n−β)+2(2δ)n+2−βΓ(n+4−β2)Γ(n+22)(n−β)Γ(β2)∥ν∥β−nif β≠n,−2nδ2(2log∥ν∥+log(δ24)+γ−ψ(n2))if β=n,

where is Euler’s constant and is the digamma function.

Near , moreover, (6) can be used to approximate the multipliers as . Figure 1 provides a comparison, in the case , of the multipliers for various choices of and . Observe that when or , the transition to the asymptotic behavior is more gradual than when both coefficients take values away from these special cases.

## 3 Eigenvalues on periodic domains

The eigenvalues of can be used to develop pseudo-spectral solvers for partial integro-differential equations involving on an -dimensional (flat) torus. In particular, suppose we consider as an operator on the periodic torus

 T=n∏i=1[0,ℓi],with ℓi>0,i=1,2,…,n.

For any , define

 να=(2πα1/ℓ1,2πα2/ℓ2,…,2παn/ℓn)Tandϕα(x)=eiνα⋅x.

The functions form a complete set in . Moreover,

 Lδ,βϕα(x)=(cδ∫Bδ(0)eiνα⋅z−1∥z∥βdz)ϕα(x)=m(να)ϕα(x),

implying that

is an eigenfunction of

with eigenvalue . Thus, (6) provides an alternative to the methods established in [11] for computing the eigenvalues of on ; any method for accurately evaluating the general hypergeometric function can be used to evaluate the eigenvalues in any spatial dimension and for any choices of and .

### 3.1 Efficiency and accuracy using Python’s mpmath library

One implementation of the hypergeometric function is provided by mpmath [12], a Python library for arbitrary precision arithmetic. Table 1 summarizes an efficiency and accuracy study using this library. For this study, we computed, for certain choices of , and , the values of with at equispaced values in the interval . (The interval was chosen in order to avoid relying too heavily on special values of at rational points.) Table 1 reports the results of this study. In each row, we report the average time (in seconds) required to compute each of the multipliers using (6). We also report the relative error at as compared against the value of the integral in (3), computed to 30 digits of precision using Mathematica [23]. This study was performed on a 2.7GHz laptop computer, and mpmath’s default precision (dps=15) was used.

### 3.2 Accelerated eigenvalue approximations

As demonstrated in Table 1, eigenvalues of

can be computed at a rate of a few thousand per second, or around 200,000 per minute. Considering the rapid convergence rate of Fourier-based solvers for smooth problems, this is probably sufficient for many problems. However, for large 2D or 3D problems, or non-smooth problems that require many Fourier modes, it may be desirable to accelerate the computation. In this section, we present one way to do this using high-accuracy Fourier interpolation combined with cubic spline interpolation.

This technique can be used effectively because is radially symmetric, only depending on the magnitude as described in (4). Suppose we wish to evaluate for many values of in an interval . Through the change of variables

 r=K2(1+cosθ), (7)

we can transform into a periodic function of . Moreover, since is an analytic function, so is , so has a Fourier series

 m(θ)=∞∑k=−∞^mkeikθ

with rapidly decaying Fourier coefficients . Figures 2 and 3 show two examples.

The technique for accelerating the evaluation of proceeds as follows. First, we sample at the points for . The parameter should be chosen sufficiently large that the coefficients are very small for . Now let be an even integer. can be resampled with high accuracy on the points for

using zero padding. (One simply uses an FFT on the

sample points, pads the coefficients with zeros in the highest frequency modes, and transforms back with an inverse FFT.) Because of the efficiency of the FFT, this can be done very quickly. Finally, by inverting the transformation (7), we obtain a highly accurate approximation of at the points for , the Chebyshev nodes of the second kind in the interval . If is chosen sufficiently large, the cubic spline interpolation of these points will provide an accurate approximation of throughout the interval.

Table 2 demonstrates the effectiveness of this approach for the two cases shown in Figures 2 and 3. Along with the relevant coefficients, the “prep time” row indicates the time required to build the interpolation data (including sampling the function on the -point grid, the FFT-based interpolation from the -point grid to the -point grid, and the construction of the cubic spline interpolant on the -point grid), the “avg interp time” row reports the average time required to evaluate using the cubic interpolation after the interpolant is constructed, and the “avg direct time” row reports the average time required to evaluate using (6) directly. The row titled “max error” reports the maximum value of where is the interpolated value and is the value from the formula in (6) at 10,000 uniformly spaced points in . Again, this study was performed on a 2.7GHz laptop computer using mpmath to compute the values of . The FFT and cubic interpolation were provided by scipy [13]. As can be seen from this example, by spending a few seconds generating interpolation data, it is possible to increase the speed of computing the multipliers by several orders of magnitude.

## 4 Spectral methods for time-dependent periodic problems

In this section, we demonstrate techniques for using the eigenvalues of to solve time-dependent problems with nonlocal spatial operators. We also compare these solvers with one of the finite-difference-like solvers described in [22].

### 4.1 Semi-analytic solutions to linear equations

On a periodic torus , the peridynamic heat and wave equations take the form,

 {ut=Lδ,βuin (0,∞)×T,u(0,x)=u0(x), (8)

and

 ⎧⎪⎨⎪⎩utt=Lδ,βuin (0,∞)×T,u(0,x)=u0(x),ut(0,x)=v0(x), (9)

respectively. The solution to (8) can be written as

 u(t,x)=∑α∈Zneλαt^u0,αeiνα⋅x, (10)

where and are the Fourier coefficients for . Similarly, the solution to (9) can be written as

 u(t,x)=∑α∈Zn^uα(t)eiνα⋅x,where^uα(t)=12(^u0,α+^v0,α√λα)e√λαt+12(^u0,α−^v0,α√λα)e−√λαt, (11)

where (as before), and and are the Fourier coefficients for and respectively. Linear nonlocal problems like these can be solved semi-analytically on using the Fast Fourier Transform (FFT) to approximate the Fourier coefficients and the initial data. To do so, we simply apply the FFT to the initial data, apply one of the solution formulas above, and then use the inverse FFT to find the solution at the specified time.

#### 4.1.1 Example: the nonlocal heat equation

Consider, for example, the nonlocal heat equation in (8) on the 2D torus with initial data . Although this is not truly a periodic function, it is indistinguishable from a smooth periodic function at machine precision. Figure 4 shows the solutions in -space (with the positive direction receding back into the page and toward the upper right) as computed by the methods described in this section. In all cases, the spatial domain was discretized using a grid with points, which is sufficient to capture the Fourier spectrum of to machine precision. The Fourier coefficients for were approximated by FFT. The temporal domain was discretized using 250 equally spaced points. For each choice of operator and each time level, (10) was used to compute the solution in the frequency domain and transformed back to the spatial domain by inverse FFT. To aid in visualization, values of are omitted from the 3D plots; the outer “shell” of each plot corresponds to the level surface . A wedge running parallel to the axis has been removed in order to show values of near .

Note that, as could be inferred from Figure 1, when , the close agreement between and the multipliers of the Laplacian ensures that the smoothed square bump decays similarly by rapidly transitioning into a circular bump and dissipating away. On the other hand, when is relatively large, the dependence on becomes more pronounced. For , the solution values spread more slowly than the classical heat equation. In fact, for , the solution does not appear to spread at all. On the other hand, when , the solution spreads and decays more rapidly than in the classical heat equation.

A similar phenomenon can be seen in Figure 5. This example is identical to the previous except in the initial condition, which can be written in polar form as

 u(t,r,θ)=exp⎡⎣−(r4(1+0.3sin(5θ)))8⎤⎦.

The figure shows how the level lines change as a function of time. Again, the results for are quite similar to the results from the classical heat equation. When is larger, though, one again sees how the choice of parameter affects the diffusion phenomenon as compared to classical diffusion: values of slow the diffusion while values of accelerate the diffusion.

#### 4.1.2 Example: the nonlocal wave equation

As a second example, consider the 2D nonlocal wave equation (9) on the torus with Gaussian initial conditions , . Once again, this function is periodic to machine precision. Figure 6 shows the solutions at time for several choices of nonlocal operator (as well as the classical wave equation). In all cases, the spatial domain was discretized using 400 points in each direction. (The Gaussian can be resolved using fewer points, but the final images are not as clear.) The solution at time was then produced using (11) and the FFT. As in Section 4.1.1, we observe that when , the behavior is quite similar to the classical wave equation while the cases with larger exhibit more sensitivity to the parameter . Again, appears as a threshold with displaying slower wave propagation than the classical case and displaying faster.

### 4.2 Pseudo-spectral methods

For nonlinear problems (or even problems with varying coefficients), the semi-analytic method from the previous section tend not to work well. The main difficulty arises from the fact that products of functions, under Fourier transform, give rise to convolutions that prevent the decoupling of the Fourier modes present in the linear case.

Nonetheless, as shown in [11], the fact that the trigonometric polynomials are eigenfunctions for in the periodic setting allows the development of fast and accurate pseudo-spectral (or Fourier collocation) numerical methods. Essentially, whenever we need to apply to a function, we first pass to the frequency domain using an FFT, then apply the operator there (which amounts to simply scaling each Fourier mode by the appropriate eigenvalue) and then invert the FFT to return to the spatial domain. Thus is applied efficiently in the frequency domain, while function multiplication can be applied efficiently in the spatial domain.

#### 4.2.1 Nonlocal diffusion in the Brusselator

As an example, consider the following nonlocal Brusselator reaction-diffusion model

 {ut=DuLδ,βu+a−(b+1)u+u2v,vt=DvLδ,βv+bu−u2v,

where we have simply replaced the Laplacian by its nonlocal version.

Figure 7 shows several solutions to the 1D periodic system approximated using the pseudospectral method described in Section 4.2. For this example, we chose parameter values

 Du=0.0625,Dv=0.12,a=3,andb=11.

These equations were solved on the periodic spatial interval for times in the interval . The intial conditions chosen were

 u(0,x) =a(1+12sin(π10x)) v(0,x) =ba+110cos(3π5x).

The numerical solution was advanced using a standard RK4 timestepping scheme combined with a spectral filter (zeroing the coefficients of the 1/3 largest wave numbers) to prevent aliasing errors. A CFL condition of the form was chosen empirically from numerical tests. The periodic interval was discretized using points, resulting in a temporal discretization using points. This discretization is sufficiently fine that the plots in the figure are not altered if the number of spatial nodes is halved.

The complexity of these figures illustrates again the importance of utilizing highly accurate methods as well as the nontrivial role played by in the nonlocal operator. Small differences in this parameter can yield significant differences in solutions.

### 4.3 Comparison with a finite difference method

In this section, we compare the pseudo-spectral method in 4.2 with a finite difference approximation in 1D. Although the pseudo-spectral method described above is very accurate for smooth, periodic problems, it does not show the same accuracy in non-periodic situations. For this reason, spectral methods have typically been avoided in numerical solvers, with finite difference and finite element methods being the preferred treatment. These solvers come with their own difficulties, however. When refining the spatial mesh size to improve accuracy, one finds that, in order for to remain fixed, the number of discretization points (or finite element cells) used in computing at a point must grow like , leading to increased computation complexity as goes to . To avoid this blow-up in computational cost, a common technique is to fix the ratio as . The problem with this technique can be seen in Figure 8.

In this figure, we consider the finite difference operator found in [22], using the scaling with . This operator approximates on an equispaced grid using a finite difference formula of the form

 Lδ,βu(x)≈AβD,0u(x)=a0u(x)+r∑j=1aj[u(x+jΔx)+u(x−jΔx)].

When treated as a finite difference operator on a periodic interval, the trigonometric polynomials are eigenfunctions of the discrete operator. To see this, let . Then

 AβD,0ϕk(x)=ϕk(x)(a0+r∑j=1aj[e2πikjΔx/L+e−2πikjΔx/L])=ϕk(x)(a0+2r∑j=1ajcos(2πkjΔx/L)).

The left and center plots in Figure 8 show the eigenvalues (upper dashed curve) of associated with for the cases and respectively. In these plots, and . Shown for comparison are the true eigenvalues of (solid curve) as well as the eigenvalues for the 1D Laplacian operator (lower dashed curve). Notice that the shape of the eigenvalue curves do not change, indicating a certain scale invariance. What this indicates is that the finite difference method with fixed ratio is never accurate for high frequency modes, exactly where the eigenvalues of differ significantly from those of the 1D Laplacian. Refining the spatial discretization in order to better resolve high-frequency modes simply has the effect of shifting those modes into the Laplacian-like part of the spectrum.

The way to remedy this problem, of course, is to fix and widen the finite difference stencil as . However, truly capturing the “interesting” part of the spectrum typically requires very large stencil sizes. Consider, for example, the right-most plot of Figure 8. The figure shows the first 30 eigenvalues of for , and along with the corresponding eigenvalues of the finite difference approximations for stencils with radii , and . Evidently, correctly capturing the spectrum of requires finite difference stencils spanning hundreds of points. Such stencils are impractical in numerical solvers.

As another way to see how the solutions when is fixed can be misleading, consider Figure 9. These plots show the results of a periodic nonlocal wave equation

 utt=Lδ,βu,(x,t)∈(0,20)×(0,40) u(0,x)=exp(−(x−10)8),ut(0,x)=0.

Note that at and , the initial conditions are zero in double precision. For this example, rather than choosing a time stepper with fixed step size, we have instead used the solve_ivp method of scipy [13], which implements and adaptive RK45 timestep algorithm.

The top row shows the results of a solver based on the finite difference operator with , . The bottom row shows the corresponding solutions for a pseudo-spectral grid points. The dashed curve in all plots shows the true solution evaluated using the FFT and (11) on a grid with nodes, which is sufficient to resolve the initial data with high accuracy.

From the figure, it is clear that the finite difference solver only approximates the true solution in the limit as and both go to zero; in other words, only when the solution to the nonlocal equation is sufficiently close to the classical wave equation. In particular, this means that the standard finite difference and finite element solvers have extreme difficulty resolving nonlocal behaviors when the solution varies on a length scale much smaller than , for example, the solutions shown in Figure 10, which demonstrates that, for large very large might be insufficient to correctly capture nonlocal wave phenomena.

This observation is quantified in Table 3, which summarizes the results displayed in Figures 9 and 10. For each test, the table provides the value of for the spectral method, and the value of and for the finite difference method. Also displayed are the associated maximum absolute errors at time (as compared to the exact solution evaluated using the FFT and (11)) along with the total CPU time required to obtain the solution (including the time to compute the multipliers in the spectral case and time to compute the finite difference stencils in the finite difference case). In all cases, the spectral method is able to attain an error of less than with only grid points in under 4 seconds. On the other hand, the finite differences approach with fixed reaches an error just under using points and 22 seconds. Worse, in the case of fixed , the finite differences method with points and obtains an error just under and requires nearly 17 minutes. This again shows the difficulty inherent in solving peridynamics problems accurately using finite difference or finite element approaches; one must either scale with the spatial discretization size (thus making the method accurate only when ) or increase the stencil size or element order when refining the discretization, leading to very expensive computations.

## 5 Discussion

The hypergeometric representation of the Fourier multipliers (6) has shown to be a useful tool for the analysis of nonlocal equations [3]. The usefulness of this representation for computations of nonlocal equations has been demonstrated in this paper. Specifically, formula (6) allows for efficient and accurate computations of the multipliers in dimensions, without the need for integration, as shown in Section 3. In addition, the fact that the nonlocality only appears as a parameter in the argument of the function in (6), leads to the uncoupling of and the grid size . This enables the development of efficient spectral methods in which does not scale with as demonstrated in Section 4. Moreover, extending the definition of the operator , through extending the definition of the multipliers, to all provides a unified approach for the computations of local and nonlocal models considered in this paper. This is demonstrated in Section 4; in particular, the solutions displayed in Figures 4, 5, (and  6), which correspond to local () or nonlocal () equations have been computed using the same diffusion spectral solver (wave spectral solver), respectively. Therefore, coupling local and nonlocal models is seamless using this approach.

The computational experiments in Section 4 indicates that, in addition to , the parameter is another locality/nonlocality parameter for which smaller values of correspond to more nonlocal behavior while larger values correspond to more local behavior. In the context of diffusion, in Section 4.1.1, the computational experiments suggest that has a critical value at corresponding to classical diffusion, with subdiffusion occurring when , and superdiffusion occurring when . This observation combined with the smoothness of the multipliers’ representation (6) in provide justification to extending to larger values of . Lastly, these observations provide a way to interpret and set apart the roles of these two nonlocality parameters with as the parameter that determines the type of phenomena in a given model while affects the scaling in the given model.

## References

• [1] Aksoylu, B., Beyer, H. R., and Celiker, F. Application and implementation of incorporating local boundary conditions into nonlocal problems. Numerical Functional Analysis and Optimization 38, 9 (2017), 1077–1114.
• [2] Aksoylu, B., and Parks, M. L. Variational theory and domain decomposition for nonlocal problems. Applied Mathematics and Computation 217, 14 (2011), 6498–6515.
• [3] Alali, B., and Albin, N. Fourier multipliers for peridynamic laplace operators. arXiv preprint arXiv:1810.11877v2 (2018). Submitted.
• [4] Bobaru, F., and Duangpanya, M. The peridynamic formulation for transient heat conduction. International Journal of Heat and Mass Transfer 53, 19-20 (2010), 4047–4059.
• [5] Bobaru, F., Foster, J. T., Geubelle, P. H., and Silling, S. A. Handbook of peridynamic modeling. CRC press, 2016.
• [6] Burch, N., and Lehoucq, R. Classical, nonlocal, and fractional diffusion equations on bounded domains. International Journal for Multiscale Computational Engineering 9, 6 (2011).
• [7] NIST Digital Library of Mathematical Functions. http://dlmf.nist.gov/, Release 1.0.17 of 2017-12-22. F. W. J. Olver, A. B. Olde Daalhuis, D. W. Lozier, B. I. Schneider, R. F. Boisvert, C. W. Clark, B. R. Miller and B. V. Saunders, eds.
• [8] Du, Q., Gunzburger, M., Lehoucq, R. B., and Zhou, K. Analysis and approximation of nonlocal diffusion problems with volume constraints. SIAM review 54, 4 (2012), 667–696.
• [9] Du, Q., Gunzburger, M., Lehoucq, R. B., and Zhou, K. A nonlocal vector calculus, nonlocal volume-constrained problems, and nonlocal balance laws. Mathematical Models and Methods in Applied Sciences 23, 03 (2013), 493–540.
• [10] Du, Q., and Yang, J. Asymptotically compatible fourier spectral approximations of nonlocal allen–cahn equations. SIAM Journal on Numerical Analysis 54, 3 (2016), 1899–1919.
• [11] Du, Q., and Yang, J. Fast and accurate implementation of fourier spectral approximations of nonlocal diffusion operators and its applications. Journal of Computational Physics 332 (2017), 118–134.
• [12] Johansson, F., et al.

mpmath: a Python library for arbitrary-precision floating-point arithmetic (version 0.18)

, December 2013.
• [13] Jones, E., Oliphant, T., Peterson, P., et al. SciPy: Open source scientific tools for Python, 2001–. [Online; accessed 3/20/2018].
• [14] Lehoucq, R. B., Reu, P. L., and Turner, D. Z. A novel class of strain measures for digital image correlation. Strain 51, 4 (2015), 265–275.
• [15] Madenci, E., and Oterkus, E. Peridynamic theory. In Peridynamic Theory and Its Applications. Springer, 2014, pp. 19–43.
• [16] Oterkus, S., Madenci, E., and Agwai, A. Peridynamic thermal diffusion. Journal of Computational Physics 265 (2014), 71–96.
• [17] Radu, P., Toundykov, D., and Trageser, J. A nonlocal biharmonic operator and its connection with the classical analogue. Archive for Rational Mechanics and Analysis 223, 2 (2017), 845–880.
• [18] Seleson, P., Gunzburger, M., and Parks, M. L. Interface problems in nonlocal diffusion and sharp transitions between local and nonlocal domains. Computer Methods in Applied Mechanics and Engineering 266 (2013), 185–204.
• [19] Silling, S. A. Reformulation of elasticity theory for discontinuities and long-range forces. Journal of the Mechanics and Physics of Solids 48 (2000), 175–209.
• [20] Silling, S. A., Epton, M., Weckner, O., Xu, J., and Askari, E. Peridynamic states and constitutive modeling. Journal of Elasticity 88, 2 (2007), 151–184.
• [21] Slevinsky, R. M., Montanelli, H., and Du, Q. A spectral method for nonlocal diffusion operators on the sphere. Journal of Computational Physics 372 (2018), 893–911.
• [22] Tian, X., and Du, Q. Analysis and comparison of different approximations to nonlocal diffusion and linear peridynamic equations. SIAM Journal on Numerical Analysis 51, 6 (2013), 3458–3482.
• [23] Wolfram Research, Inc. Mathematica, Version 11.3. Champaign, IL, 2018.