Sparse Monte Carlo method for nonlocal diffusion problems

TL;DR

This paper introduces a sparse Monte Carlo combined with discontinuous Galerkin method for efficiently solving nonlocal diffusion equations, proving convergence and analyzing error sources, with numerical validation.

Contribution

It presents a novel sparse sampling approach for nonlocal diffusion problems, reducing computational cost while maintaining accuracy, with rigorous convergence analysis.

Findings

Sparse Monte Carlo method reduces the number of discretization points needed.

Convergence of the proposed numerical scheme is proven and error estimates are provided.

Numerical experiments confirm the theoretical convergence rates.

Abstract

A class of evolution equations with nonlocal diffusion is considered in this work. These are integro-differential equations arising as models of propagation phenomena in continuum media with nonlocal interactions including neural tissue, porous media flow, peridynamics, models with fractional diffusion, as well as continuum limits of interacting dynamical systems. The principal challenge of numerical integration of nonlocal systems stems from the lack of spatial regularity of the data and solutions intrinsic to nonlocal models. To overcome this problem we propose a semidiscrete numerical scheme based on the combination of sparse Monte Carlo and discontinuous Galerkin methods. An important feature of our method is sparsity. Sparse sampling of points in the Monte Carlo approximation of the nonlocal term allows to use fewer discretization points without compromising the accuracy. We prove convergence of the numerical method and estimate the rate of convergence. There are two principal ingredients in the error of the numerical method related to the use of Monte Calro and Galerkin approximations respectively. We analyze both errors. Two representative examples of discontinuous kernels are presented. The first example features a kernel with a singularity, while the kernel in the second example experiences jump discontinuity. We show how the information about the singularity in the former case and the geometry of the discontinuity set in the latter translate into the rate of convergence of the numerical procedure. In addition, we illustrate the rate of convergence estimate with a numerical example of an initial value problem, for which an explicit analytic solution is available. Numerical results are consistent with analytical estimates.

Equations

$${\partial }_{t}u(t,x)$$

$$u(0,x)$$

$$\left|D(u_1)-D(u_2)\right|\le L_D|u_1-u_2|\text{mbox}and\left|f(u_1,x,t)-f(u_2,x,t)\right|\le L_f|u_1-u_2|$$

$$\underset{u\in \mathbb{R}}{\sup }|D(u)|\le 1.$$

$$\left|D(u_1,v_1)-D(u_2,v_2)\right|\le L_D\left(|u_1-u_2|+|v_1-v_2|\right)\forall u_1,u_2,v_1,v_2\in \mathbb{R}.$$

$$\max \left\{{\mathrm{esssup}}_{x\in Q}\int |W(x,y)|dy,{\mathrm{esssup}}_{y\in Q}\int |W(x,y)|dx\right\}\le W_1,$$

$${\partial }_{t}u(t,x)=f(u,x,t)+\int W^{+}(x,y)D\left(u(t,y)-u(t,x)\right)dy-\int W^{-}(x,y)D\left(u(t,y)-u(t,x)\right)dy,$$

$$W\ge 0.$$

$${\dot{u}}_{n,\bar{i}}$$

$$u_{n,\bar{i}}(0)$$

$$Q_{n,\bar{i}}=\left[\frac{i_1-1}{n},\frac{i_1}{n}\right)\times \left[\frac{i_2-1}{n},\frac{i_2}{n}\right)\times \cdots \times \left[\frac{i_d-1}{n},\frac{i_d}{n}\right),$$

$$g_{n,\bar{i}}=n^d{\int }_{Q_{n,\bar{i}}}g(x)dx,f_{n,\bar{i}}(u,t)=n^d{\int }_{Q_{n,\bar{i}}}f(u,x,t)dx,\bar{i}:=(i_1,i_2,\dots ,i_d)\in [n{]}^d.$$

$${\alpha }_n=n^{-d\gamma },0\le \gamma <1,$$

$$W_{n,\bar{i}\bar{j}}=n^{2d}{\int }_{Q_{n,\bar{i}}\times Q_{n,\bar{j}}}W(x,y)dxdy,$$

$$P(a_{n,\bar{i}\bar{j}}=1)={\alpha }_n{W}_{n,\bar{i}\bar{j}},\bar{i},\bar{j}\in [n{]}^d.$$

$${\tilde{W}}_n(x,y):={\alpha }_n^{-1}\wedge W(x,y)\text{mbox}and{W}_{n,\bar{i}\bar{j}}=n^{2d}{\int }_{Q_{n,\bar{i}}\times Q_{n,\bar{j}}}\tilde{W}(x,y{)}_{n}dxdy,$$

$$P(a_{n,\bar{i}\bar{j}}=1)={\alpha }_n{W}_{n,\bar{i},\bar{j}},\bar{i},\bar{j}\in [n{]}^d.$$

$$\max \left\{{\mathrm{esssup}}_{x\in Q}\int W^k(x,y)dy,{\mathrm{esssup}}_{y\in Q}\int W^k(x,y)dx\right\}\le {\bar{W}}_k,k\in [4].$$

$$\begin{array}{rl}\underset{t\in [0,T]}{\sup }\Vert u(t,\cdot )-u_n(t,\cdot ){\Vert }_{L^2(Q)}& \le C(\Vert g-g_n{\Vert }_{L^2(Q)}+\underset{u\in \mathbb{R},t\in [0,T]}{\sup }\Vert f(u,\cdot ,t)-f_n(u,\cdot ,t){\Vert }_{L^2(Q)}\\ & +\Vert {\tilde{W}}_n-W{\Vert }_{L^2(Q^2)}+\Vert {\tilde{W}}_n-P_n{\tilde{W}}_n{\Vert }_{L^2(Q^2)}+n^{-d(1-\gamma -\delta )/2}),\end{array}$$

$$W_n(x,y)=\sum _{(\bar{i},\bar{j})\in {[n{]}^d}^2}{W}_{n,\bar{i}\bar{j}}{1}_{Q_{n,\bar{i}}\times Q_{n,\bar{j}}}(x,y).$$

$$f_n(u,x,t)=\sum _{\bar{i}\in [n{]}^d}{f}_{n,\bar{i}}(u,t)1_{Q_{n,\bar{i}}}(x).$$

$$\begin{array}{rl}\underset{t\in [0,T]}{\sup }\Vert u(t,\cdot )-u_n(t,\cdot ){\Vert }_{L^2(Q)}& \le C(\Vert g-g_n{\Vert }_{L^2(Q)}+\underset{u\in \mathbb{R},t\in [0,T]}{\sup }\Vert f(u,\cdot ,t)-f_n(u,\cdot ,t){\Vert }_{L^2(Q)}\\ & +\Vert W_n-P_n{W}_n{\Vert }_{L^2(Q^2)}+n^{-d(1-\gamma -\delta )/2})\text{mbox}a.s..\end{array}$$

$${\dot{v}}_{n,\bar{i}}$$

$$v_{n,\bar{i}}(0)$$

$$v_n(t,x)=\sum _{\bar{i}\in [n{]}^d}{v}_{n,\bar{i}}(t)1_{Q_{n,\bar{i}}}(x).$$

$$u_n(t,x)=\sum _{\bar{i}\in [n{]}^d}{u}_{n,\bar{i}}(t)1_{Q_{n,\bar{i}}}(x).$$

$${\partial }_t{v}_n(t,x)$$

$$v_n(0,x)$$

$$\underset{t\in [0,T]}{\sup }\Vert u_n(t,\cdot )-v_n(t,\cdot ){\Vert }_{L^2(Q)}\le C{n}^{-d(1-\gamma -\delta )/2}\text{mbox}a.s.,$$

$$w_n(t,x)=u(t,x)-v_n(t,x).$$

$$\begin{array}{rl}\int {\partial }_t{w}_n(t,x)w_n(t,x)dx& =\int \left(f(u(t,x),x,t)-f(v_n(t,x),x,t)\right)w_n(t,x)dx\\ & +\int \left(f(v_n(t,x),x,t)-f_n(v_n(t,x),x,t)\right)w_n(t,x)dx\\ & +\int \int W(x,y)\left[D(u(t,y)-u(t,x))-D(v_n(t,y)-v_n(t,x))\right]w_n(t,x)dydx\\ & +\int \int \left(W(x,y)-W_n(x,y)\right)D(v_n(t,y)-v_n(t,x))w_n(t,x)dydx.\end{array}$$

$$\left|\int \left(f(u(t,x),x,t)-f(v_n(t,x),x,t)\right)w_n(t,x)dx\right|\le L_f\int w_n(t,x{)}^{2}dx,$$

Full text

We prove convergence of the numerical method and estimate the rate of convergence. There are two principal ingredients in the error of the numerical method related to the use of Monte Calro and Galerkin approximations respectively. We analyze both errors. Two representative examples of discontinuous kernels are presented. The first example features a kernel with a singularity, while the kernel in the second example experiences jump discontinuity. We show how the information about the singularity in the former case and the geometry of the discontinuity set in the latter translate into the rate of convergence of the numerical procedure. In addition, we illustrate the rate of convergence estimate with a numerical example of an initial value problem, for which an explicit analytic solution is available. Numerical results are consistent with analytical estimates.

1 Introduction

We propose a numerical method for the initial value problem (IVP) for a nonlinear heat equation with nonlocal diffusion

[TABLE]

For analytical convenience, we take $Q=[0,1]^{d}$ as a spatial domain. Throughout this paper, when the domain of integration is not specified, it is assumed to be $Q$. Further, $g\in L^{2}(Q)$, $W\in L^{2}(Q^{2})$, $f$ is a bounded measurable function on ${\mathbb{R}}\times Q\times{\mathbb{R}}^{+}$, which is Lipschitz continuous in $u,$ continuous in $t,$ and integrable in $x,$ and $D$ is a Lipschitz continuous function on ${\mathbb{R}}:$

[TABLE]

for all $(x,t)\in Q\times{\mathbb{R}}.$ Throughout this paper we assume

[TABLE]

This assumption may be dropped if an apriori estimate on $\|u\|_{C(0,T;L^{\infty}(Q))}$ for $T>0$ is available. Furthermore, the analysis below applies to models with the interaction function of a more general form $D(u_{1},u_{2})$ provided

[TABLE]

However, we keep $D(u_{1},u_{2}):=D(u_{2}-u_{1})$ to emphasize the connection to diffusion problems.

Equation (1.1) is a nonlocal diffusion problem. It arises as a continuum limit of interacting particle systems [24, 17]. Equations of this form are used for modeling population dynamics [31, 7, 30, 2, 6], neural tissue [8], porous media flows [10, 11], and various other biological and physicochemical processes involving anomalous diffusion [32, 4]. The key distinction of the evolution equations with nonlocal diffusion from their classical counterparts is the lack of smoothening property. A priori the solution of (1.1), (1.2) is a square integrable function in $x$ for all $t>0$ [24] and it may not possess much more regularity beyond that, unless the initial data and kernel $W$ are smooth [21, Theorem 3.3]. The lack of smoothness is a serious challenge for constructing numerical schemes for (1.1), (1.2) and for analyzing their convergence. All deterministic quadrature formulas require at least piecewise differentiability for a guaranteed convergence rate. The problem is even more challenging in high dimensional spatial domains. The main idea underlying our approach is to use the Monte Carlo approximation of the nonlocal term in (1.1). We take advantage of the essential feature of the Monte Carlo method: the independence of the convergence rate on the regularity of the integrand. The second key idea is sparsity, whose use is twofold: First, sparse sampling of points in the Monte Carlo method is used to minimize computation without compromising the accuracy. For $W$ with jump discontinuity across Lipschitz hypersurfaces, the use of sparsity is computationally beneficial starting from $d=2$. If the discontinuity set has nontrivial fractal dimension, the sparse Monte Carlo method performs better than its dense counterpart already for $d=1$ (cf. Lemma 4.4). Furthermore, sparsity is the key for extending the Monte Carlo method for models with singular kernels (see § 4.1). Not only does it allow to apply the Monte Carlo method for unbounded functions, it also makes it automatically adaptive: more sample points are selected near the singularities. The combination of these ideas together with the discontinuous Galerkin method yields a numerical scheme for the IVP (1.1), (1.2) that performs well under minimal assumptions on the regularity of $W$ and initial data.

This paper is based on our previous work on convergence of interacting particle systems on convergent graph sequences [21, 22, 18, 24]. Continuum limit is a powerful tool for studying various aspects of network dynamics including existence, stability, and bifurcations of spatiotemporal patterns [33, 23, 25, 26]. Very often the derivation of the continuum limit is based on heuristic considerations and its rigorous mathematical justification is a nontrivial problem. Recently, motivated by the theory of graph limits [20, 19, 5], we proved convergence to the continuum limit for a broad class of dynamical systems on graphs [21]. Importantly, our proof applies to models on random graphs [22] including sparse random graphs [18, 24]. These results prepared the ground for the numerical method proposed in this paper. There is an intimate relation between the problem of the continuum limit for interacting particle systems and numerical integration of nonlocal diffusion models. Given a continuum model (1.1), one can construct the corresponding particle system, approximating (1.1). This idea had been already mentioned in [21], but has never been detailed. Further, recent results for the continuum limit of coupled systems on sparse graphs indicate a strong potential of sparse discretization for numerical integration of nonlocal problems. It is the goal of this paper to present these ideas in detail.

In the next section, we present a discretization of (1.1), which can be viewed as an interacting dynamical system on a sparse random graph. The structure of the graph is determined by the kernel $W$, which defines the asymptotic connectivity of the graph sequence parametrized by the size of the graph. In the theory of graph limits, such functions are called graphons [19], the term we adopt for the reminder of this paper. In Section 3, we prove convergence of the semidiscrete (discrete in space and continuous in time) approximation of (1.1) and turn to estimating the rate of convergence in Sections 4 and 5. There are two main factors contributing to the error of approximation. The first is due to approximating the nonlocal term in (1.1) by a random sum (Monte Carlo method), while the second is due to approximating the kernel and the initial data by piecewise constant functions (discontinuous Galerkin method). The rate of convergence of the sparse Monte Carlo approximation follows from our previous results [24, Theorem 4.1]. Convergence of piecewise constant approximation in the $L^{2}$–norm follows from classical theorems of analysis (cf. the Lebesgue-Besicovitch Theorem [15] or $L^{2}$–Martingale Convergence Theorem [34]). However, neither of these theorems elucidates the rate of convergence. In fact, the example in § 4.2 shows that without additional hypotheses the algebraic convergence may be arbitrarily slow. To this end, we study what determines the rate of convergence of piecewise constant approximations for a square integrable function. For Hölder continuous functions the answer is simple (cf. Lemma 4.1). For discontinuous functions, on the other hand, the answer naturally depends on the type of discontinuity. In Section 4, we consider two examples elucidating this issue. The first example is based on a singular (unbounded) graphon. It shows how the information about the singularity translates into the rate of convergence estimate. Here, we also see how to use sparsity to optimize computation. The second example adapted from [21], on the other hand, presents a bounded graphon with jump discontinuity (§ 4.2). In this case, the convergence rate depends on the geometry of the set of the discontinuity (cf. Lemma 4.4). In the light of these examples, in Section 5, we perform convergence analysis under general assumptions on $W$. In Section 6, we illustrate rate of convergence estimates with a numerical example. Here, we choose a nonlinear problem, which has an explicit solution. This allows us to verify the rate of convergence of the $L^{2}$-error as the discretization step tends to zero. Special attention is paid to the dependence of the convergence rate on sparsity. Finally, in Section 7 we present a proof of a technical Lemma 3.5, which extends the corresponding result in [24] to models in multidimensional domains and affords a wider range of sparsification.

Numerical methods for nonlocal diffusion problems have been subject of intense research recently due to their increased use in modeling [28, 14, 13, 4, 3, 29]. Compared to the existing literature, the contribution of the present work is that our method applies to problems with nonlinear diffusivity as well as to problems with more general form of the interaction function (cf. (1.5)). The main focus of this paper is how to deal with models with low regularity of the data. We believe that the combination of the Monte Carlo and discontinuous Galerkin methods provides an effective tool for numerical integration of nonlocal problems under minimal regularity assumptions.

2 The model and its discretization

In this section, we formulate the technical assumptions on the kernel $W$ and describe the numerical scheme for solving the IVP (1.1)-(1.2).

We assume that $W\in L^{2}(Q^{2})$ is subject to the following assumptions:

[TABLE]

Theorem 2.1**.**

Let $W\in L^{2}(Q^{2})$ satisfy (W-1). Then for any $g\in L^{2}(Q)$ and $T>0$ there is a unique solution of the IVP (1.1), (1.2) $u\in C^{1}(0,T;L^{2}(Q))$.

Proof.

The proof is as in [18, Theorem 3.1] with minor adjustments. ∎

Next, we note that the kernel in the nonlocal term may be assumed nonnegative. Indeed, by writing $W=W^{+}-W^{-}$ as the difference of its positive and negative parts, one can rewrite (1.1) as

[TABLE]

where the nonlocal terms splits into the difference of two terms with nonnegative kernels. Thus, without loss of generality, in the remainder of this paper we will assume

[TABLE]

We approximate the IVP (1.1), (1.2) by the following system of ordinary differential equations

[TABLE]

where

[TABLE]

[TABLE]

The semidiscrete system (2.3) can be viewed as a system of interacting particles on a random graph $\Gamma_{n}$ with the node set $[n]^{d}$ and adjacency matrix $(a_{n,\bar{i}\bar{j}})$. The positive sequence

[TABLE]

is used to control the sparsity of $\Gamma_{n}$. The adjacency matrix $(a_{n,\bar{i}\bar{j}})$ is defined as follows. The case $W\in L^{\infty}(Q^{2})$ is slightly different and so we treat it separately. Thus, there are two cases to consider.

(I)

Suppose $W\in L^{\infty}(Q^{2})$. Without loss of generality, we further assume that $0\leq W\leq 1$. Then let

[TABLE]

and

[TABLE]

If $\alpha_{n}\equiv 1(\gamma=0)$ (cf. (2.7)), $\Gamma_{n}$ is a dense W-random graph [20], otherwise $\Gamma_{n}$ is sparse with the mean degree $O(n^{d(1-\gamma)}),\;0\leq\gamma<1.$

(II)

Alternatively, if $W$ is in $L^{2}(Q^{2})$ but not in $L^{\infty}(Q^{2})$ then let

[TABLE]

where $\alpha_{n}$ defined in (2.7) with $\gamma\in(0,1)$. Then

[TABLE]

3 Convergence of the numerical method

In this section, we study convergence of the discrete scheme (2.3), (2.4). We first deal with the more general case of unbounded graphon $W$ and then specialize the result for $W\in L^{\infty}(Q^{2})$.

The following additional mild assumption on $W$ is used to get a wider range of sparsity. Let nonnegative $W\in L^{4}(Q^{2})$ satisfy

[TABLE]

Theorem 3.1**.**

Suppose nonnegative $W\in L^{4}(Q^{2})$ is subject to (W-1s), $D,$ $f,$ and $g$ are as in (1.1), (1.2). Further, $\alpha_{n}=n^{-d\gamma}$ for some $\gamma\in(0,1)$. Then for arbitrary $0<\delta<1-\gamma,$ we have

[TABLE]

where $C$ is a positive constant independent of $n$, and $P_{n}\tilde{W}_{n}=:W_{n}$ stands for the $L^{2}$-projection of $\tilde{W}_{n}$ onto the finite–dimensional subspace $X_{n}=\operatorname{span}\{{\mathbf{1}}_{Q_{n,\bar{i}}\times Q_{n,\bar{j}}},\,(\bar{i},\bar{j})\in[n]^{2d}\}$:

[TABLE]

and

[TABLE]

Estimate (3.1) holds almost surely (a.s.) with respect to the random graph model.

Remark 3.2*.*

The theorem still holds without (W-1s), i.e., for square integrable $W$ subject to (W-1). In this case, the last term on the right-hand side of (3.1) is replaced by $n^{-d(1/2-\gamma-\delta)},$ $\gamma\in(0,1/2),$ and $\delta<1/2-\gamma.$

The first two terms on the right–hand side of (3.1) correspond to the error of approximation of the initial data $g\in L^{2}(Q)$ and $f(u,x,t)$ by the step functions in $x$. Further, $\|\tilde{W}_{n}-W\|^{2}_{L^{2}(Q^{2})}$ and $\|\tilde{W}_{n}-P_{n}\tilde{W}_{n}\|^{2}_{L^{2}(Q^{2})}$ bound the error of approximation of $W$ by a bounded step function $W_{n}$. Here, the first term $\|\tilde{W}_{n}-W\|^{2}_{L^{2}(Q^{2})}$ is the error of truncating $W$ and the second term $\|\tilde{W}_{n}-P_{n}\tilde{W}_{n}\|^{2}_{L^{2}(Q^{2})}$ is the error of approximation of the truncated function $\tilde{W}_{n}$ by projecting it onto a finite–dimensional subspace. Finally, the last term on the right–hand side of (3.1) is the error of the approximation of the nonlocal term by the random sum in (2.3).

For bounded graphons $W,$ Theorem 3.1 implies the following result.

Corollary 3.3**.**

Let $W\in L^{\infty}(Q^{2})$. Then under the assumptions of Theorem 3.1 we have

[TABLE]

where $C$ is a positive constant independent of $n.$

Remark 3.4*.*

From (3.2) one can see how to use sparsity to optimize computation. Already for $d=1,$ if the largest of the two errors of approximation of $g$ and $W$ by step functions is $O(n^{-\kappa})$ with $1/2<\kappa<1$ (cf. §4.2) and the nonlinearity $f(\cdot)$ does not depend on $x,$ then taking $\gamma=1-2\kappa$ one can use sparse discretization without compromising the accuracy. This has obvious computational advantages over dense random and, moreover, deterministic spatial discretization schemes, e.g., Galerkin method. Sparse random discretization is even more efficient when $d>1$.

The proof of Theorem 3.1 modulo a few minor details proceeds as the proof of convergence to the continuum limit in [22, 24]. First, the solution of the IVP (2.3), (2.4) is compared to that of the IVP for the averaged equation:

[TABLE]

Then the solution of the averaged problem is compared to the solution of the IVP (1.1), (1.2). It is convenient to view the solution of the averaged problem as a function on ${\mathbb{R}}^{+}\times Q$:

[TABLE]

Likewise, we interpret the solution of the discrete problem (2.3), (2.4) as a function on ${\mathbb{R}}^{+}\times Q:$

[TABLE]

We recast the IVP (3.3), (3.4) as follows

[TABLE]

The first step of the proof of convergence of the numerical scheme (2.3), (2.4) is accomplished in the following lemma.

Lemma 3.5**.**

Let nonnegative $W\in L^{4}(Q^{2})$ subject to (W-1s), and $\alpha_{n}=n^{-d\gamma},\;\gamma\in(0,1)$ (cf. (4.24)). Then for any $T>0$ for solutions of (2.3) and (3.3) subject to the same initial conditions, we have

[TABLE]

where $0<\delta<1-\gamma,$ and positive constant $C$ independent of $n$.

The proof of the lemma is technical and is relegated to Section 7. The result still holds without for square integrable functions the additional assumption (W-1s) albeit for a narrower range of $\gamma\in(0,0.5)$ (cf. [24, Theorem 4.1]).

Proof of Theorem 3.1.

Denote the difference between the solutions of the original IVP (1.1), (1.2) and the averaged IVP (3.7), (3.8)

[TABLE]

By subtracting (3.3) from (1.1), multiplying the resultant equation by $w_{n}$, and integrating over $Q$, we obtain

[TABLE]

Using Lipschitz continuity of $f(u,x,t)$ in $u$ and an elementary case of the Young’s inequality, we obtain

[TABLE]

[TABLE]

where $\|\cdot\|$ stands for the $L^{2}(Q)$-norm. Recall that $D$ is bounded by $1$ (cf. (1.4)). Using this bound and the Young’s inequality, we obtain

[TABLE]

Finally, using Lipschitz continuity of $D$ and Young’s inequality, we estimate

[TABLE]

where we used Fubini theorem and (W-1) in the last line.

By combining (3.11)-(3.15), we arrive at

[TABLE]

where $L=1+2L_{f}+L_{D}(3W_{1}+W_{2}).$

By Gronwall’s inequality, we have

[TABLE]

∎

4 Two examples

The error of approximation of the nonlocal term by a random sum, the last term on the right–hand side of (3.1), is known explicitly. Next in importance is the error of approximation of the square integrable graphon $W$ by the step function $W_{n}.$ This error depends on the regularity of the graphon. In this section, we consider two representative examples of $W$: a singular graphon (§ 4.1) and a bounded graphon with jump discontinuities (§ 4.2). Motivated by these examples in the next section, we will analyze the rate of convergence estimates under general assumptions on graphon $W$.

We will begin with the following estimate for Hölder continuous functions. To this end, $\phi\in L^{p}(Q),\;p\geq 1,$ and

[TABLE]

Lemma 4.1**.**

Suppose $\phi\in L^{p}(Q),\;p\geq 1,$ is a Hölder continuous function

[TABLE]

Then

[TABLE]

Here and below, $h:=n^{-1}$.

Proof.

Using Jensen’s inequality and (4.18), we have

[TABLE]

∎

Remark 4.2*.*

Below, we will freely apply Lemma 4.1 to functions on $Q$ and on $Q^{2}$. The latter are clearly covered by the lemma by taking $Q:=Q^{2}$.

4.1 A singular graphon

Consider the problem of approximation by step functions of the singular kernel graphon

[TABLE]

where $0<\lambda<d/2$.

Lemma 4.3**.**

For $\gamma\in(0,1/2)$ and $0<\lambda<d/2$ we have

[TABLE]

Proof.

Below, we will use the following change of variables $(x,y)=T(u,v)$ for $(x,y)$ and $(u,v)$ from ${\mathbb{R}}^{2d},$ defined by

[TABLE] 2. 2.

Let $\alpha_{n}=n^{-d\gamma}$ and recall that $\tilde{W}_{n}=h^{-d\gamma}\wedge W$. Denote $\tilde{Q}=\left\{(x,y)\in Q^{2}:\;|x-y|^{-\lambda}\geq h^{-d\gamma}\right\}$. Further,

[TABLE]

where we used (4.22) followed by the change to polar coordinates.

Thus,

[TABLE] 3. 3.

Next we turn to estimating $\left\|\tilde{W}_{n}-W_{n}\right\|_{L^{2}(Q^{2})}$. Since the truncated function $\tilde{W}_{n}$ is Lipschitz continuous on $Q^{2}$, by Lemma 4.1,

[TABLE]

It remains to estimate the Lipschitz constant $L(\tilde{W}_{n})\leq\operatorname{ess\,sup}_{Q^{2}}|\nabla\tilde{W}_{n}|.$ On $Q^{2}-\tilde{Q}$,

[TABLE]

The gradient approaches its greatest value as $|x-y|\searrow h^{d\gamma\over\lambda}$. Thus,

[TABLE]

and

[TABLE] 4. 4.

The statement of the lemma follows (4.24) and (4.26) and the triangle inequality.

∎

Next we choose $\gamma$ to optimize the rate of convergence in (4.21). By setting the two exponents of $h$ on the right–hand side of (4.21) equal, we see that the rate is optimal for

[TABLE]

With this choice of $\gamma,$

[TABLE]

To optimize the rate of convergence of the numerical scheme (2.3), (2.4), one has to choose $\gamma\in(0,1)$ to maximize the smallest of the following three exponents

[TABLE]

where the last exponent comes from the error of the Monte Carlo approximation (cf. (3.1)).

4.2 $\{0,1\}$-valued functions

The following example is adapted from [21]. It shows how jump discontinuities affect the rate of convergence of approximation by piecewise constant functions. The accuracy of approximation depends on the geometry of the hypersurface of discontinuity, more precisely, on its fractal dimension.

Let $Q^{+}$ be a closed subset of $Q$ and consider

[TABLE]

Denote by $\partial Q^{+}$ the boundary of $Q$ and recall the upper box-counting dimension of $\partial Q^{+}$

[TABLE]

where $N_{h}(\partial Q^{+})$ stands for the number of $Q_{n,\bar{i}},\;\bar{i},\bar{j}\in[n]^{d},$ having nonempty intersection with $\partial Q^{+}$ (cf. [16]).

Lemma 4.4**.**

[TABLE]

for some positive $C$ independent on $n$.

Proof.

As in (4.19), we have

[TABLE]

Note that the only nonzero terms in the sum on the right–hand side of (4.30) are the integrals over $Q_{n,\bar{i}}\times Q_{n,\bar{j}}$’s having nonempty intersection with $\partial Q^{+}.$ Thus,

[TABLE]

where we used (4.28). ∎

Remark 4.5*.*

Note that as $\beta\to d-0$ the rate of convergence in (4.29) can be made arbitrarily low.

5 The rate of convergence of the numerical method

5.1 Approximation by step functions

In this section, we address the rate of convergence of the Galerkin component of the numerical scheme (2.3), (2.4). Specifically, we study the error of the approximation of the graphon $W\in L^{2}(Q^{2})$ and the initial data $g\in L^{2}(Q)$ by step functions.

We will need an $L^{p}$-modulus of continuity of function on a unit $d$-cube $Q=[0,1]^{d}$. In fact, we only need the $L^{2}$-modulus of continuity, but present the analysis in the more general $L^{p}$-setting, since this does not require any extra effort. For functions on the real line, the definition of the $L^{p}$-modulus of continuity can be found in [1, 12]. Here, we present a suitable adaptation of this definition for the problem at hand.

Definition 5.1**.**

For $\phi\in L^{p}(Q)$ we define the $L^{p}$-modulus of continuity

[TABLE]

where $|\xi|_{\infty}:=\max_{i\in[d]}|\xi_{i}|$ and $Q_{\xi}=\{x\in Q:\;x+\xi\in Q\}$.

For $\alpha\in(0,1],$ we define a generalized Lipschitz space111 Below, we will freely apply the definitions and various estimates established for functions on $Q$ to functions on $Q^{2},$ for which they are trivially valid by setting $d:=2d$. In particular, the definitions of the modulus of continuity and the corresponding Lipschitz spaces obviously translate to functions on $Q^{2}$.

[TABLE]

Clearly $\operatorname{Lip}\left(\alpha,L^{p}(Q)\right)$ contains $\alpha$-Hölder continuous functions. However, Lipschitz spaces are much larger than Hölder spaces. For instance, $\operatorname{Lip}\left(1/p,L^{p}(Q)\right)$ contains discontinuous functions.

Below, we express the error of approximation of $\phi\in L^{p}(Q)$ by a step function through $\omega_{p}(\phi,h)$. The analysis works out a little cleaner for dyadic discretization of $Q$, which will be assumed for the remainder of this section. Thus, we approximate $\phi\in L^{p}(Q)$ by a piecewise constant function

[TABLE]

where $\phi_{Q_{2^{m},\bar{i}}}$ stands for the mean value of $\phi$ on $Q_{2^{m},\bar{i}}$

[TABLE]

Lemma 5.2**.**

For $\phi\in\operatorname{Lip}\left(\alpha,L^{p}(Q)\right)$, we have

[TABLE]

where $C$ is independent of $m$.

Proof.

Fix $m\in{\mathbb{N}}$ and denote $h:=2^{-m}$. To simplify notation, throughout the proof we drop $2^{m}$ in the subscript of $x_{2^{m},i}$ and $Q_{2^{m},\bar{i}}$.

We write

[TABLE]

where

[TABLE]

[TABLE]

and

[TABLE]

Rewrite (5.3)

[TABLE]

By subtracting (5.5) from (5.6) we have

[TABLE]

Further,

[TABLE]

Integrating both sides of (5.7) over $Q$ and using Jensen’s inequality, we continue

[TABLE]

Thus,

[TABLE]

Since $\phi\in\operatorname{Lip}\left(\alpha,L^{p}(Q)\right)$, we have

[TABLE]

where $C$ depends on $\phi,d,$ and $p$ but not $m$.

Let $m\in{\mathbb{N}}$ be arbitrary but fixed. For any integer $M>m$ we have

[TABLE]

By passing $M$ to infinity in (5.10), we get (5.4). ∎

5.2 The rate of convergence

We now can combine Theorem 3.1 and Lemma 5.4 to estimate the convergence rate for (2.3), (2.4). For the model with a bounded graphon $W$ (cf. (I), Section 2) we have the following theorem.

Theorem 5.3**.**

Suppose that in addition to the assumptions of Theorem 3.1, for some $\alpha_{i}\in(0,1],\;i\in[3],$ $g\in\operatorname{Lip}\left(\alpha_{1},L^{2}(Q)\right),$ $W\in\operatorname{Lip}\left(\alpha_{2},L^{2}(Q^{2})\right)\cap L^{\infty}(Q^{2}),$ and $f(u,\cdot,t)\in\operatorname{Lip}\left(\alpha_{3},L^{2}(Q)\right)$ uniformly for $(u,t)\in{\mathbb{R}}\times[0,T],$ i.e.,

[TABLE]

where $C>0$ is independent of $(u,t).$

Then

[TABLE]

where $C$ is independent of $n$.

If $W\in L^{2}(Q^{2})$ has singularities then the convergence rate may also depend on the accuracy of approximation of $W$ by the truncated function $\tilde{W}_{n}$. We do not estimate the truncation error for a general $W\in L^{2}(Q^{2})$. For an example of how this error can be estimated for a given graphon in practice, we refer to the example in § 4.1.

6 Numerical example

In this section, we illustrate convergence analysis in the previous sections with a numerical example. To this end, we consider an IVP for the continuum Kuramoto model with nonlocal nearest–neighbor coupling [27]:

[TABLE]

where $u(t,x)\in{\mathbb{T}},\;{\mathbb{T}}={\mathbb{R}}/2\pi{\mathbb{Z}},$ stands for the phase of the oscillator at $x\in[0,1]$, $\omega\in{\mathbb{R}}$ is its intrinsic frequency. Function $K,$ describing the connectivity of the network, is first defined on $[0,1/2)$ by

[TABLE]

and then extended as a $1$–periodic function on ${\mathbb{R}}$. The initial condition

[TABLE]

is called a $q$–twisted state (Figure 1a). For $\omega=0$, $u^{(q)}$ is a stationary solution of (6.12). Thus,

[TABLE]

solves the IVP (6.12), (6.13). We use the explicit solution (6.16) to compute the error of the numerical integration of (6.12), (6.13).

To estimate the rate of convergence of the numerical scheme (2.3), (2.4) we use the following values of parameters: $r=0.2,$ $\omega=0.5,$ and $q=3.$ For these parameter values, the travelling wave solution (6.16) is unstable. We integrated (6.12) numerically for $t\in[0,1],$ using the fourth order Runge–Kutta method with the time step $10^{-2}.$ Note that the error of the Runge–Kutta method, i.e., of the discretization in time is significantly smaller than of that of the discretizing in space (c.f. (2.3), (2.4)). We integrated (6.12) numerically for different values of $\gamma$ and for $n\in\{128,256\}$. For each pair $(\gamma,n)$ we repeated the numerical experiment $200$ times and computed the mean value of the error of numerical integration (compared to the exact solution (6.16)). The mean errors $\bar{e}_{\gamma,128}$ and $\bar{e}_{\gamma,256}$ computed for $n=128$ and $n=256$ respectively are used to determine the convergence rate:

[TABLE]

The results of this numerical experiment are shown in Figure 1b. Our main goal was to verity the dependence of the convergence rate on sparsity controlled by $\gamma$. The pixel pictures for the adjacency matrices of random graphs corresponding to the nonlocal nearest–neighbor coupling for $n=512$ and different values of $\gamma$ are shown in Figure 2. The plot in Figure 1b shows a clear linear relation between the exponent $\alpha_{\gamma}$ and $\gamma$. The numerical rates plotted by blue stars are slightly lower the theoretical rates $(1-\gamma)/2$ plotted in red. Overall numerical rates show a good fit with the analytical estimate.

7 Proof of Lemma 3.5

In this section, we prove Lemma 3.5. The proof follows the lines of the proof of Theorem 4.1 in [24], which covers $\gamma\in(0,1/2)$ for $d=1$. Extension to the multidimensional case $d>1$ is straightforward. Lemmas 7.3 and 7.4 adapted from [9] allow to extend the range of $\gamma$ to $(0,1)$. The reader not interested in the extended range of $\gamma$ may find a simpler proof in [24] easier to follow. For those interested in the full range of $\gamma,$ below we present the following proof of Lemma 3.5.

Theorem 7.1**.**

Let nonnegative $W\in L^{4}(Q^{2})$ satisfy

[TABLE]

and

[TABLE]

Then for solutions of (2.3) and (3.3) subject to the same initial conditions and arbitrary $0<\epsilon<1/2$, we have

[TABLE]

for arbitrary $T>0$ and positive constant $C$ independent of $n$. In particular, for $\alpha_{n}=n^{-d\gamma},\;\gamma\in(0,1)$, we have

[TABLE]

where $0<\delta<1-\gamma$ can be taken arbitrarily small.

We precede the proof of Theorem 7.1 with several auxiliary estimates.

Lemma 7.2**.**

From (W-1s) it follows

[TABLE]

Proof.

We prove (7.4) assuming that nonnegative $W$ is in $L^{2}(Q^{2}),$ but not in $L^{\infty}(Q^{2})$. In this case, $W_{n,\bar{i}\bar{j}}$ are defined by (4.24). For arbitrary $\bar{k}\in[n]^{d}$ and $n\in{\mathbb{N}},$ we have

[TABLE]

where we used Jensen’s inequality in the second line and (W-1s) in the last line. Thus,

[TABLE]

The bound for $\sup_{n\in{\mathbb{N}}}\max_{\bar{j}\in[n]^{d}}n^{-d}\sum_{\bar{i}\in[n]^{d}}W^{k}_{n,\bar{i}\bar{j}}$ is proved similarly. ∎

Lemma 7.3**.**

For $K\geq 2\bar{W}_{1}$, we have

[TABLE]

In particular, with probability $1$ there exists $n_{0}\in{\mathbb{N}}$ such that

[TABLE]

for all $n\geq n_{0}$.

For the next lemma, we will need the following notation

[TABLE]

and $Z_{n}=(Z_{n,\bar{i}},\;\bar{i}\in[n]^{d})$.

Lemma 7.4**.**

For arbitrary $\epsilon>0$, we have

[TABLE]

where $C$ is a positive constant independent of $n$ and

[TABLE]

Proof of Theorem 7.1.

Recall that $f(u,x,t)$ and $D$ are Lipschitz continuous function in $u$ with Lipschitz constants $L_{f}$ and $L_{D}$ respectively.

Further, $a_{n,ij},$ are Bernoulli random variables

[TABLE]

Denote $\psi_{n,\bar{i}}:=v_{n,\bar{i}}-u_{n,\bar{i}}.$ By subtracting (2.3) from (3.3), multiplying the result by $n^{-d}\psi_{n,\bar{i}},$ and summing over $\bar{i}\in[n]^{d}$, we obtain

[TABLE]

where $\|\cdot\|^{2}_{2,n^{d}}$ is the discrete $L^{2}$-norm (cf. (7.12)).

Using Lipschitz continuity of $f$ in $u$, we have

[TABLE]

Using Lipschitz continuity of $D$ and the triangle inequality, we have

[TABLE]

Using Lemma 7.3 and (7.4), we obtain

[TABLE]

Similarly,

[TABLE]

By plugging (7.17) and (7.18) into (7.16), we have

[TABLE]

It remains to bound $I_{2}$:

[TABLE]

The combination of (7.14), (7.15), (7.19) and (7.20) yields

[TABLE]

where $L=L_{f}+L_{D}\left(2K+{3\over 2}\bar{W}_{1}+{1\over 2}\bar{W}_{2}\right)+{1\over 2}.$

Using the Gronwall’s inequality and Lemma 7.4, we have

[TABLE]

∎

Proof of Lemma 7.3.

Let

[TABLE]

Note that for fixed $\bar{i}\in[n]^{d},$ $\{\xi_{n,\bar{i}\bar{j}},\,\bar{j}\in[n]^{d}\}$ are mean zero independent RVs. Further, using the definition of $\xi_{n,\bar{i}\bar{j}},$ it is straightforward to bound

[TABLE]

From (7.25), we have

[TABLE]

Using Bernstein’s inequality and the union bound, we have

[TABLE]

Finally, the combination of (7.23) and (7.27) yields

[TABLE]

This proves (7.13). By Borel-Cantelli Lemma, (7.7) follows. ∎

Proof of Lemma 7.4.

Recall (7.8)-(7.10) and rewrite

[TABLE]

where

[TABLE]

By (7.1), one can choose a sequence $\delta_{n}\searrow 0$ such that

[TABLE]

Specifically, let

[TABLE]

and define events

[TABLE]

[TABLE]

Clearly,

[TABLE]

We want to show that $\operatorname{\mathbb{P}}\left(\Omega_{n}\;\mbox{infinitely often}\right)=0$. By Borel-Cantelli Lemma, it is sufficient to show that

[TABLE]

From Lemma 7.3, we know that $\sum_{n\geq 1}\operatorname{\mathbb{P}}(A_{n})<\infty$ for $K>2\bar{W}_{1}.$ In the remainder of the proof, we show that $\sum_{n\geq 1}\operatorname{\mathbb{P}}(\Omega_{n}\cap A_{n}^{c})$ is convergent.

Applying the exponential Markov inequality to $\operatorname{\mathbb{P}}(\Omega_{n}\left|A_{n}^{c}\right.)$, from $\operatorname{\mathbb{P}}(\Omega_{n}\cap A_{n}^{c})\leq\operatorname{\mathbb{P}}(\Omega_{n}\left|A_{n}^{c}\right.)$ and (7.32), we have

[TABLE]

Using the independence of $\eta_{n,\bar{i}\bar{k}\bar{l}}$ in $\bar{i}\in[n]^{d}$, we have

[TABLE]

Using

[TABLE]

and the Cauchy-Schwartz inequality, we bound the right–hand side of (7.36) as follows

[TABLE]

From (7.10), (7.29), and under $A_{n}^{c}$ (cf. (7.33)), we have

[TABLE]

Further,

[TABLE]

Using (7.10), we estimate sum of the fourth moments of $\eta_{n,\bar{i}\bar{j}}$

[TABLE]

where we also use (7.4). Similarly,

[TABLE]

By combining (7.39)-(7.41), we obtain

[TABLE]

By plugging (7.38) and (7.42) into (7.37), we obtain

[TABLE]

Using this bound on the right–hand side of (7.36), we further obtain

[TABLE]

Using (7.44), from (7.35) we obtain

[TABLE]

Furthermore, using (7.31) it is straightforward to check that

[TABLE]

The statement of the lemma then follows from (7.32)-(7.34) via Borel-Cantelli Lemma. ∎

Acknowledgements. This work was supported in part by the NSF grant DMS 1715161.