Fourier spectral methods for nonlocal models

TL;DR

This paper introduces efficient Fourier spectral methods for solving nonlocal models in multiple dimensions, leveraging hypergeometric representations to compute eigenvalues accurately, and applies these methods to both linear and nonlinear nonlocal PDEs.

Contribution

The paper develops a novel spectral solver framework for nonlocal models using hypergeometric functions, enabling accurate eigenvalue computation and efficient solutions in multiple dimensions.

Findings

Spectral methods outperform finite-difference solvers in accuracy and efficiency.

Eigenvalues of nonlocal operators can be computed accurately using hypergeometric functions.

The approach is successfully applied to linear and nonlinear nonlocal PDEs, including the Brusselator model.

Abstract

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.

Tables

Table 1. Table 1 : Efficiency and accuracy study, as described in Section  3.1 , using mpmath ’s implementation of F 3 2 subscript subscript 𝐹 3 2 {}_{2}F_{3} to compute m ​ ( ν ) 𝑚 𝜈 m(\nu) .

				rel err
$n$	$\beta$	$\delta$	avg time	$\Vert \nu \Vert =318\pi$
1	0.25	0.1	2.689e-04	1.216e-16
1	1.00	0.1	2.958e-04	2.401e-17
1	1.50	0.1	3.152e-04	8.259e-17
2	0.75	0.1	3.076e-04	7.624e-17
2	2.00	0.1	2.515e-04	8.525e-17
2	3.00	0.1	2.879e-04	3.949e-18
3	1.75	0.1	3.791e-04	2.442e-17
3	3.00	0.1	2.813e-04	2.960e-17
3	4.50	0.1	2.766e-04	3.990e-17

Table 2. Table 2 : Results of the interpolation test described in Section  3.2 .

$n$	2	2
$\beta$	0.5	2.3
$\delta$	1.2	0.4
$K$	1000	1000
$N$	1500	600
$M$	20000	10000
prep time	3.117 s	0.423 s
avg interp time	8.720e-8 s	1.354e-7 s
avg direct time	1.949e-3 s	8.099e-4 s
max error	3.758e-10	5.301e-11

Table 3. Table 3 : Comparison of accuracy and computing time for the spectral and finite difference methods described in Section  4.3 .

	spectral			finite differences
$\delta$	$N$	error	CPU (s)	$N$	$r$	error	CPU (s)
0.3000	2000	7.490e-06	3.46	1000	3	4.635e-01	1.06
0.1500	2000	6.522e-06	3.53	2000	3	3.293e-01	2.21
0.0750	2000	7.853e-06	3.53	4000	3	1.976e-01	6.77
0.0375	2000	9.464e-06	3.60	8000	3	9.423e-02	22.11
5.0000	2000	7.462e-06	3.77	2000	100	3.429e-01	17.33
				4000	200	1.934e-01	121.34
				8000	400	9.051e-02	1019.27

Equations

$$L^{\delta ,\beta }u(x)=c^{\delta ,\beta }{\int }_{B_{\delta }(x)}\frac{u(y)-u(x)}{\Vert y-x{\Vert }^{\beta }}dy,$$

$$c^{\delta ,\beta }$$

$$=$$

$$m(\nu )=c^{\delta ,\beta }{\int }_{B_{\delta }(0)}\frac{\cos (\nu \cdot z)-1}{\Vert z{\Vert }^{\beta }}dz,$$

$$m(r)=-r^2{}_2{F}_3\left(1,\frac{n+2-\beta }{2};2,\frac{n+2}{2},\frac{n+4-\beta }{2};-\frac{1}{4}{r}^2{\delta }^2\right),$$

$$L^{\delta ,\beta }u(x)=\frac{1}{(2\pi {)}^n}{\int }_{{\mathbb{R}}^n}m(\nu )\hat{u}(\nu )e^{i\nu \cdot x}d\nu ,$$

$$m(\nu )=-\Vert \nu {\Vert }^2{}_2{F}_3\left(1,\frac{n+2-\beta }{2};2,\frac{n+2}{2},\frac{n+4-\beta }{2};-\frac{1}{4}\Vert \nu {\Vert }^2{\delta }^2\right).$$

$$m(\nu )\sim \{\begin{array}{ll}-\frac{2n(n+2-\beta )}{{\delta }^2(n-\beta )}+2{\left(\frac{2}{\delta }\right)}^{n+2-\beta }\frac{\Gamma \left(\frac{n+4-\beta }{2}\right)\Gamma \left(\frac{n+2}{2}\right)}{(n-\beta )\Gamma \left(\frac{\beta }{2}\right)}\Vert \nu {\Vert }^{\beta -n}& \text{if β=n},\\ -\frac{2n}{{\delta }^2}\left(2\log \Vert \nu \Vert +\log \left(\frac{{\delta }^2}{4}\right)+\gamma -\psi (\frac{n}{2})\right)& \text{if β=n},\end{array}$$

$$\mathbf{T}=\prod _{i=1}^n[0,{\ell }_i],\text{with }{\ell }_i>0,i=1,2,\dots ,n.$$

$${\nu }_{\alpha }=(2\pi {\alpha }_1/{\ell }_1,2\pi {\alpha }_2/{\ell }_2,\dots ,2\pi {\alpha }_n/{\ell }_n{)}^T\text{and}{\varphi }_{\alpha }(x)=e^{i{\nu }_{\alpha }\cdot x}.$$

$$L^{\delta ,\beta }{\varphi }_{\alpha }(x)=\left(c_{\delta }{\int }_{B_{\delta }(0)}\frac{e^{i{\nu }_{\alpha }\cdot z}-1}{\Vert z{\Vert }^{\beta }}dz\right){\varphi }_{\alpha }(x)=m({\nu }_{\alpha }){\varphi }_{\alpha }(x),$$

$$r=\frac{K}{2}(1+\cos \theta ),$$

$$m(\theta )=\sum _{k=-\infty }^{\infty }{\hat{m}}_k{e}^{ik\theta }$$

$$\{\begin{array}{l}{u}_t=L^{\delta ,\beta }u\text{in }(0,\infty )\times \mathbf{T},\\ u(0,x)=u_0(x),\end{array}$$

$$\{\begin{array}{l}{u}_{tt}=L^{\delta ,\beta }u\text{in }(0,\infty )\times \mathbf{T},\\ u(0,x)=u_0(x),\\ u_t(0,x)=v_0(x),\end{array}$$

$$u(t,x)=\sum _{\alpha \in {\mathbb{Z}}^n}{e}^{{\lambda }_{\alpha }t}{\hat{u}}_{0,\alpha }{e}^{i{\nu }_{\alpha }\cdot x},$$

$$\begin{array}{rl}u(t,x)& =\sum _{\alpha \in {\mathbb{Z}}^n}{\hat{u}}_{\alpha }(t)e^{i{\nu }_{\alpha }\cdot x},\text{where}\\ {\hat{u}}_{\alpha }(t)& =\frac{1}{2}\left({\hat{u}}_{0,\alpha }+\frac{{\hat{v}}_{0,\alpha }}{\sqrt{{\lambda }_{\alpha }}}\right)e^{\sqrt{{\lambda }_{\alpha }}t}+\frac{1}{2}\left({\hat{u}}_{0,\alpha }-\frac{{\hat{v}}_{0,\alpha }}{\sqrt{{\lambda }_{\alpha }}}\right)e^{-\sqrt{{\lambda }_{\alpha }}t},\end{array}$$

$$u(t,r,\theta )=\exp \left[-{\left(\frac{r}{4(1+0.3\sin (5\theta ))}\right)}^8\right].$$

$$\{\begin{array}{l}{u}_t=D_u{L}^{\delta ,\beta }u+a-(b+1)u+u^{2}v,\\ v_t=D_v{L}^{\delta ,\beta }v+bu-u^{2}v,\end{array}$$

$$D_u=0.0625,D_v=0.12,a=3,\text{and}b=11.$$

$$u(0,x)$$

$$v(0,x)$$

$$L^{\delta ,\beta }u(x)\approx {\mathbb{A}}_{D,0}^{\beta }u(x)=a_{0}u(x)+\sum _{j=1}^r{a}_j[u(x+j\Delta x)+u(x-j\Delta x)].$$

$$\begin{array}{rl}{\mathbb{A}}_{D,0}^{\beta }{\varphi }_k(x)& ={\varphi }_k(x)\left(a_0+\sum _{j=1}^r{a}_j\left[e^{2\pi ikj\Delta x/L}+e^{-2\pi ikj\Delta x/L}\right]\right)\\ & ={\varphi }_k(x)\left(a_0+2\sum _{j=1}^r{a}_j\cos (2\pi kj\Delta x/L)\right).\end{array}$$

$$u_{tt}=L^{\delta ,\beta }u,(x,t)\in (0,20)\times (0,40)$$

$$u(0,x)=\exp (-(x-10{)}^8),u_t(0,x)=0.$$

Full text

Keywords: Nonlocal equations, nonlocal Laplacian, nonlocal operators, peridynamics, nonlocal diffusion, Brusselator, nonlocal wave equation, Fourier multipliers, eigenvalues.

1 Introduction

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

[TABLE]

where $B_{\delta}(x)$ is a ball in ${\mathbb{R}^{n}}$ and the scaling constant $c^{\delta,\beta}$ is defined by

[TABLE]

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 $n$ dimensions. We pursue the Fourier multipliers’ approach introduced in [3]. The Fourier multipliers of $L^{\delta,\beta}$ are given by the integral representation

[TABLE]

for $\beta<n+2$. A key observation in our approach is that this $n$-dimensional integral in (3) is realized through a 1D formula for which the multipliers satisfy

[TABLE]

where $r=\|\nu\|$. Noting that the hypergeometric function ${}_{2}F_{3}$ is smooth in its argument, it follows from (4) that evaluating the multipliers in $\mathbb{R}^{n}$ reduces to computing this 1D smooth function. Numerical methods for computing hypergeometric functions such as ${}_{2}F_{3}$ are available in the literature, see for example [7]. We adopt the implementation of ${}_{2}F_{3}$ 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 $L^{\delta,\beta}$ in Section 3.2.

Utilizing the computations of the eigenvalues of $L^{\delta,\beta}$, 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 $m$ of $L^{\delta,\beta}$ are defined through the Fourier transform by

[TABLE]

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

Theorem 1**.**

Let $n\geq 1$, $\delta>0$ and $\beta<n+2$. The Fourier multipliers can be written in the form

[TABLE]

The smoothness of the ${}_{2}F_{3}$ function in its parameters, allows us to extend formula  (6) for all $\beta\in\mathbb{R}\setminus\{n+4,n+6,n+8,\ldots\}$. Therefore, we utilize  (6) to extend the definition of the multipliers $m(\nu)$ and, consequently, the definition of the operator $L^{\delta,\beta}$ acting on the class of Schwartz functions through the inverse Fourier transform in (5) to the case when $\beta\geq n+2$ (with $\beta\neq n+4,n+6,n+8,\ldots$).

A characterization of the asymptotic behavior of the multipliers, and in particular that of the eigenvalues of $L^{\delta,\beta}$, is given by the following result.

Theorem 2**.**

Let $n\geq 1$, $\delta>0$ and $\beta\in\mathbb{R}\setminus\{n+2,n+4,n+6,\ldots\}$. Then, as $\|\nu\|\to\infty$,

[TABLE]

where $\gamma$ is Euler’s constant and $\psi$ is the digamma function.

Near $\nu=0$, moreover, (6) can be used to approximate the multipliers as $m(\nu)=-\|\nu\|^{2}+O(\|\nu\|^{4})$. Figure 1 provides a comparison, in the case $n=2$, of the multipliers $m(\nu)$ for various choices of $\beta$ and $\delta$. Observe that when $\beta\approx n+2=4$ or $\delta\approx 0$, 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 $L^{\delta,\beta}$ can be used to develop pseudo-spectral solvers for partial integro-differential equations involving $L^{\delta,\beta}$ on an $n$-dimensional (flat) torus. In particular, suppose we consider $L^{\delta,\beta}$ as an operator on the periodic torus

[TABLE]

For any $\alpha\in\mathbb{Z}^{n}$, define

[TABLE]

The functions $\{\phi_{\alpha}\}_{\alpha\in\mathbb{Z}^{n}}$ form a complete set in $L^{2}(\mathbf{T})$. Moreover,

[TABLE]

implying that $\phi_{\alpha}$ is an eigenfunction of $L^{\delta,\beta}$ with eigenvalue $m(\nu_{\alpha})$. Thus, (6) provides an alternative to the methods established in [11] for computing the eigenvalues of $L^{\delta,\beta}$ on $\mathbf{T}$; any method for accurately evaluating the general hypergeometric function ${}_{2}F_{3}$ can be used to evaluate the eigenvalues in any spatial dimension $n$ and for any choices of $\beta\in\mathbb{R}\setminus\{n+4,n+6,n+8,\ldots\}$ and $\delta>0$.

3.1 Efficiency and accuracy using Python’s mpmath library

One implementation of the ${}_{2}F_{3}$ 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 $n$, $\beta$ and $\delta$, the values of $m(\nu)$ with $\|\nu\|$ at $1000$ equispaced values in the interval $[1,318\pi]$. (The interval was chosen in order to avoid relying too heavily on special values of ${}_{2}F_{3}$ 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 $1000$ multipliers using (6). We also report the relative error at $\|\nu\|=m(318\pi)$ 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 $L^{\delta,\beta}$ 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 $m(\nu)$ is radially symmetric, only depending on the magnitude $r=\|\nu\|$ as described in (4). Suppose we wish to evaluate $m$ for many values of $r$ in an interval $[0,K]$. Through the change of variables

[TABLE]

we can transform $m$ into a periodic function of $\theta$. Moreover, since $m(r)$ is an analytic function, so is $m(\theta)$, so $m(\theta)$ has a Fourier series

[TABLE]

with rapidly decaying Fourier coefficients $\hat{m}_{k}$. Figures 2 and 3 show two examples.

The technique for accelerating the evaluation of $m(r)$ proceeds as follows. First, we sample $m(\theta)$ at the points $\theta=2\pi j/N$ for $j=0,1,\ldots,N-1$. The parameter $N$ should be chosen sufficiently large that the coefficients $\hat{m}_{k}$ are very small for $|k|>N/2$. Now let $M>N$ be an even integer. $m(\theta)$ can be resampled with high accuracy on the points $\theta=2\pi j/M$ for $j=0,1,\ldots,M-1$ using zero padding. (One simply uses an FFT on the $N$ sample points, pads the $N$ coefficients with $M-N$ 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 $m(z)$ at the points $r=\frac{K}{2}(1+\cos(2\pi j/M))$ for $j=0,1,\ldots,M/2$, the Chebyshev nodes of the second kind in the interval $[0,K]$. If $M$ is chosen sufficiently large, the cubic spline interpolation of these points will provide an accurate approximation of $m(r)$ 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 $N$-point grid, the FFT-based interpolation from the $N$-point grid to the $M$-point grid, and the construction of the cubic spline interpolant on the $M$-point grid), the “avg interp time” row reports the average time required to evaluate $m(\nu)$ using the cubic interpolation after the interpolant is constructed, and the “avg direct time” row reports the average time required to evaluate $m(\nu)$ using (6) directly. The row titled “max error” reports the maximum value of $|x_{t}-x_{a}|/(1+|x_{t}|)$ where $x_{a}$ is the interpolated value and $x_{t}$ is the value from the formula in (6) at 10,000 uniformly spaced points in $[0,K]$. Again, this study was performed on a 2.7GHz laptop computer using mpmath to compute the values of ${}_{2}F_{3}$. 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 $L^{\delta,\beta}$ 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 $\mathbf{T}$, the peridynamic heat and wave equations take the form,

[TABLE]

and

[TABLE]

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

[TABLE]

where $\lambda_{\alpha}=m(\nu_{\alpha})$ and $\hat{u}_{0,\alpha}$ are the Fourier coefficients for $u_{0}$. Similarly, the solution to (9) can be written as

[TABLE]

where $\lambda_{\alpha}=m(\nu_{\alpha})$ (as before), and $\hat{u}_{0,\alpha}$ and $\hat{v}_{0,\alpha}$ are the Fourier coefficients for $u_{0}$ and $v_{0}$ respectively. Linear nonlocal problems like these can be solved semi-analytically on $\mathbf{T}$ 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 $\mathbf{T}=[-10,10]^{2}$ with initial data $u_{0}(x,y)=\exp(-x^{8}-y^{8})$. 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 $(t,x,y)$-space (with the positive $t$ 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 $800\times 800$ points, which is sufficient to capture the Fourier spectrum of $u_{0}$ to machine precision. The Fourier coefficients for $u_{0}$ were approximated by FFT. The temporal domain $[0,15]$ was discretized using 250 equally spaced points. For each choice of operator and each time level, (10) was used to compute the solution $u(t,x,y)$ in the frequency domain and transformed back to the spatial domain by inverse FFT. To aid in visualization, values of $u<0.1$ are omitted from the 3D plots; the outer “shell” of each plot corresponds to the level surface $u(t,x,y)=0.1$. A wedge running parallel to the $t$ axis has been removed in order to show values of $u$ near $(x,y)=(0,0)$.

Note that, as could be inferred from Figure 1, when $\delta\approx 0$, the close agreement between $m(\nu)$ 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 $\delta$ is relatively large, the dependence on $\beta$ becomes more pronounced. For $\beta<n+2=4$, the solution values spread more slowly than the classical heat equation. In fact, for $\beta<n=2$, the solution does not appear to spread at all. On the other hand, when $\beta>n+2$, 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

[TABLE]

The figure shows how the level lines $u=0.1$ change as a function of time. Again, the results for $\delta\approx 0$ are quite similar to the results from the classical heat equation. When $\delta$ is larger, though, one again sees how the choice of parameter $\beta$ affects the diffusion phenomenon as compared to classical diffusion: values of $\beta<n+2=4$ slow the diffusion while values of $\beta>n+2=4$ 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 $\mathbf{T}=[-48,48]^{2}$ with Gaussian initial conditions $u_{0}(x)=\exp(-\|x\|^{2})$, $v_{0}(x)=0$. Once again, this function is periodic to machine precision. Figure 6 shows the solutions $u(t,x)$ at time $t=10$ 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 $t=10$ was then produced using (11) and the FFT. As in Section 4.1.1, we observe that when $\delta\approx 0$, the behavior is quite similar to the classical wave equation while the cases with larger $\delta$ exhibit more sensitivity to the parameter $\beta$. Again, $\beta=n+2=4$ appears as a threshold with $\beta<4$ displaying slower wave propagation than the classical case and $\beta>4$ 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 $L^{\delta,\beta}$ in the periodic setting allows the development of fast and accurate pseudo-spectral (or Fourier collocation) numerical methods. Essentially, whenever we need to apply $L^{\delta,\beta}$ 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 $L^{\delta,\beta}$ 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

[TABLE]

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

[TABLE]

These equations were solved on the periodic spatial interval $[0,20]$ for times in the interval $[0,40]$. The intial conditions chosen were

[TABLE]

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 $\Delta t=1.9\Delta x^{2}$ was chosen empirically from numerical tests. The periodic interval $[0,20]$ was discretized using $1600$ points, resulting in a temporal discretization using $134,737$ 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 $\beta$ 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 $\Delta x$ to improve accuracy, one finds that, in order for $\delta$ to remain fixed, the number of discretization points (or finite element cells) used in computing $L^{\delta,\beta}$ at a point must grow like $\delta/\Delta x$, leading to increased computation complexity as $\Delta x$ goes to [math]. To avoid this blow-up in computational cost, a common technique is to fix the ratio $\delta/\Delta x$ as $\Delta x\to 0$. The problem with this technique can be seen in Figure 8.

In this figure, we consider the finite difference operator $\mathbb{A}^{\beta}_{D,0}$ found in [22], using the scaling $\delta=r\Delta x$ with $r=3$. This operator approximates $L^{\delta,\beta}u(x)$ on an equispaced grid using a finite difference formula of the form

[TABLE]

When treated as a finite difference operator on a periodic interval, the trigonometric polynomials are eigenfunctions of the discrete operator. To see this, let $\phi_{k}(x)=e^{2\pi ikx/L}$. Then

[TABLE]

The left and center plots in Figure 8 show the eigenvalues $\lambda_{k}$ (upper dashed curve) of $\mathbb{A}^{\beta}_{D,0}$ associated with $1\leq k\leq N$ for the cases $N=100$ and $N=10000$ respectively. In these plots, $L=1$ and $\Delta x=1/(2N+1)$. Shown for comparison are the true eigenvalues of $L^{\delta,\beta}$ (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 $\delta/\Delta x$ ratio is never accurate for high frequency modes, exactly where the eigenvalues of $L^{\delta,\beta}$ 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 $\delta$ and widen the finite difference stencil as $\Delta x\to 0$. 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 $L^{\delta,\beta}$ for $L=1$, $\delta=0.1$ and $\beta=1/3$ along with the corresponding eigenvalues of the finite difference approximations for stencils with radii $r=6$, $r=36$ and $r=216$. Evidently, correctly capturing the spectrum of $L^{\delta,\beta}$ requires finite difference stencils spanning hundreds of points. Such stencils are impractical in numerical solvers.

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

[TABLE]

Note that at $x=0$ and $x=20$, 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 $\mathbb{A}^{\beta}_{D,0}$ operator with $\delta=3\Delta x$, $\beta=1/3$. The bottom row shows the corresponding solutions for a pseudo-spectral $N=2000$ grid points. The dashed curve in all plots shows the true solution evaluated using the FFT and (11) on a grid with $N=8000$ 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 $\Delta x$ and $\delta$ 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 $\delta$, for example, the solutions shown in Figure 10, which demonstrates that, for large $\delta$ very large $r$ 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 $N$ for the spectral method, and the value of $N$ and $r$ for the finite difference method. Also displayed are the associated maximum absolute errors at time $t=40$ (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 $10^{-5}$ with only $N=2000$ grid points in under 4 seconds. On the other hand, the finite differences approach with fixed $r$ reaches an error just under $10^{-1}$ using $N=8000$ points and 22 seconds. Worse, in the case of fixed $\delta=5$, the finite differences method with $N=8000$ points and $r=400$ obtains an error just under $10^{-1}$ 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 $\delta$ with the spatial discretization size (thus making the method accurate only when $L^{\delta,\beta}\approx\Delta$) 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 $n$ dimensions, without the need for integration, as shown in Section 3. In addition, the fact that the nonlocality $\delta$ only appears as a parameter in the argument of the function ${}_{2}F_{3}$ in (6), leads to the uncoupling of $\delta$ and the grid size $\Delta x$. This enables the development of efficient spectral methods in which $\delta$ does not scale with $\Delta x$ as demonstrated in Section 4. Moreover, extending the definition of the operator $L^{\delta,\beta}$, through extending the definition of the multipliers, to all $\beta\in\mathbb{R}\setminus\{n+4,n+6,n+8,\ldots\}$ 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 ($\beta\geq n+2$) or nonlocal ($\beta<n+2$) 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 $\delta$, the parameter $\beta$ is another locality/nonlocality parameter for which smaller values of $\beta$ 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 $\beta$ has a critical value at $\beta=n+2$ corresponding to classical diffusion, with subdiffusion occurring when $\beta<n+2$, and superdiffusion occurring when $\beta>n+2$. This observation combined with the smoothness of the multipliers’ representation (6) in $\beta$ provide justification to extending $L^{\delta,\beta}$ to larger values of $\beta$. Lastly, these observations provide a way to interpret and set apart the roles of these two nonlocality parameters with $\beta$ as the parameter that determines the type of phenomena in a given model while $\delta$ affects the scaling in the given model.