Numerical approximations for a nonlocal evolution equation

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. SIAM J. NUMER. ANAL.c 2011 Society for Industrial and Applied Mathematics Vol. 49, No. 5, pp. 2103–2123 NUMERICAL APPROXIMATIONS FOR A NONLOCAL EVOLUTION EQUATION∗ MAYTE P´ EREZ-LLANOS†AND JULIO D. ROSSI‡ Abstract. In this paper we study numerical approximations of continuous solutions to the nonlocal p-Laplacian type diffusion equation, ut(t, x)=ΩJ(x−y)|u(t, y)−u(t, x)|p−2(u(t, y)− u(t, x)) dy. First, we find that a semidiscretization in space of this problem gives rise to an ODE system whose solutions converge uniformly to the continuous one as the mesh size goes to zero. Moreover, the semidiscrete approximation shares some properties of the continuous problem: it preserves the total mass and the solution converges to the mean value of the initial condition as t goes to infinity. Next, we also discretize the time variable and present a totally discrete method which also enjoys the above mentioned properties. In addition, we investigate the limit as pgoes to infinity in these approximations and obtain a discrete model for the evolution of a sandpile. Finally, we present some numerical experiments that illustrate our results. Key words. numerical approximations, nonlocal diffusion, p-Laplacian, Neumann boundary conditions, sandpiles AMS subject classifications. 65R20, 65N40, 45G10 DOI. 10.1137/110823559 1. Introduction. Our main goal in this paper is to approximate numerically a nonlocal nonlinear diffusion problem, involving the nonlocal p-Laplacian operator (with homogeneous Neumann boundary conditions). More precisely, we deal with the problem (1.1) ⎧ ⎨ ⎩ ut(t, x)=Ω J(x−y)|u(t, y)−u(t, x)|p−2(u(t, y)−u(t, x)) dy, x ∈Ω,t>0, u(x, 0) = u0(x), being J:Rd→Ra nonnegative continuous radial function with compact support satisfying J(0) >0. We also assume that RdJ(x)dx = 1 to simplify our arguments (although this condition is not necessary to prove our results). We take 1 ≤p<+∞ and Ω ⊂Rda bounded domain. Existence and uniqueness of a strong solution to (1.1) (see section 2 for a precise definition) can be found in [3], [4]. Nonlocal evolution equations of the form (1.2) ut(x, t)=J∗u−u(x, t)=Rd J(x−y)(u(y,t)−u(x, t)) dy have many applications in modelling diffusion processes; see the recent book [4] and the references [1], [5], [6], [10], [11], [12], [14], [15], [16], [24], [25], [31], [32], [33]. They are even used in the treatment of images; see [9, 26]. ∗Received by the editors February 7, 2011; accepted for publication (in revised form) August 9, 2011; published electronically October 18, 2011. http://www.siam.org/journals/sinum/49-5/82355.html †Departamento de Matem´atica, Universidad Aut´onoma de Madrid, Campus Cantoblanco, 28049, Madrid, Spain ([email protected]). This author’s research is partly supported by DGICYT grant MTM2010-18128 (Spain). ‡Departamento de An´alisis Matem´atico, Universidad de Alicante, Ap. correos 99, 03080 Alicante, Spain, on leave from Departamento de Matem´atica, FCEyN UBA (1428) Buenos Aires, Argentina ([email protected]). This author’s research is partially supported by UBA X066 and CONICET (Argentina). 2103 Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. 2104 MAYTE P´ EREZ-LLANOS AND JULIO D. ROSSI As stated in [24], if u(x, t) represents the density of a single population at the point xat time t,andJ(x−y) is considered as the probability distribution of jumping from location yto location x, then the convolution (J∗u)(x, t)=RdJ(y−x)u(y,t)dy is the rate at which individuals are arriving to position xfrom any other place, while −u(x, t)=−RdJ(y−x)u(x, t)dy is the rate at which they are leaving location xto travel to any other site. Under these considerations and in the absence of external or internal sources, the density usatisfies (1.2). Equation (1.2) is known as a nonlocal diffusion equation since the diffusion of the density uat a point xand time tdepends on all the values of uin a neighborhood of x, through the convolution term J∗u. In addition, nonlocal equations have been used recently in peridynamics, a model for fractures in elastic materials; see, for example, [17], [28], [29], [30]. For references that include the analysis of numerical approximations for linear diffusion models, we refer to [13] and [18]. Let us mention that in [13] discontinuous solutions are also considered. Here we analyze only continuous solutions but we deal with a nonlinear diffusion problem. The nonlocal diffusion equations share many properties with the corresponding local problems. In [3] the authors found that problem (1.1) is the nonlocal analogous problem to the well-known local p-Laplacian evolution equation ut= div(|∇u|p−2∇u), for p>1 (while the extreme case, p= 1, corresponds to the total variation flow), with homogeneous Neumann boundary conditions. Indeed, if one rescales the kernel Jproperly, solutions to (1.1) converge to solutions of the local problem. Note that, since we are integrating in Ω, we are imposing that diffusion takes place only in Ω. There is no flux of individuals across the boundary. Hence, we are dealing here with the analogue of the Neumann boundary condition. Finally, concerning the large time behavior, solutions to (1.1) converge to the mean value of the initial condition, as it happens for the local problem. Moreover, in [2] the limit as p→∞is considered. It is given by a nonlocal model for the formation and growth of a sandpile, analogous to the local model described in [20, 21, 22]. Now, let us describe our results in the one-dimensional case; that is, Ω = [0,L], since the extension to several space dimensions is straightforward. Assume that the support of Jis the interval [−S, S]. Let us consider a partition x1,...,x N+1 of the interval [0,L] of size h=max(xi−xi−1). This partition is not necessarily uniform but we assume that Nh =O(1). Then, the numerical approximation to problem (1.1) solves the following system of ODEs, for each node xj; (1.3) ⎧ ⎨ ⎩ u j(t)=  i∈Aj J(xi−xj)hi|ui(t)−uj(t)|p−2(ui(t)−uj(t)),t>0, uj(0) = u0(xj) for every j=1,...,N+1, being Aj={isuch that |xi−xj|≤S}and hi=xi−xi−1. In order to obtain a nontrivial scheme we assume that hS. Next, we discretize the time variable using the explicit Euler method and we obtain (1.4) ⎧ ⎪ ⎨ ⎪ ⎩ uk+1 j−uk j τk = i∈Aj J(xi−xj)hi|uk i−uk j|p−2(uk i−uk j),k>0, u0 j=u0(xj). Note that for these nonlocal models the explicit Euler scheme is well suited, since it does not need a restrictive stability constraint for the time step, as it happens for Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. NUMERICAL APPROXIMATIONS FOR A NONLOCAL EQUATION 2105 local problems. This is related to the lack of regularizing effect in nonlocal problems. See Remark 2.10. Our results concerning the semidiscrete and totally discrete approximations read as follows (see sections 2.1 and 2.2 for precise statements). Solutions to the numerical scheme presented here converge uniformly to the continuous solution as the mesh size h(and the time step τ) goes to zero. Moreover, the numerical approximations share some properties of the continuous problem: there is a comparison principle, they preserve the symmetry and the total mass of the initial data and they converge to the mean value of the initial condition as tgoes to infinity. We remark that our results also hold when we deal with approximations in a multidimensional domain. The proofs are similar to the one-dimensional case and are left to the reader. See the short paragraph at the end of section 2. We also study the limit as p→∞and we prove that solutions to the semidiscrete scheme (1.3) converge to a nonlocal evolution problem, that can be regarded as a semidiscrete approximation of a model for the evolution of a sandpile; see [2]. In the local sandpile models it is assumed that the maximum slope of a sandpile is one (otherwise the sand configuration becomes unstable); see [20, 21, 22]. In the nonlocal sandpile model described in [2] the same restriction on the slope is assumed but with some freedom at short distances (distances shorter than the size of the support of J), that can measure irregularities of the sand grains. The rest of the paper is organized as follows: in section 2 we study the semidiscrete and the totally discrete approximations in one space dimension, in section 3 we take the limit as p→∞, and finally in section 4 we show some numerical experiments that illustrate our results. 2. Semidiscrete and totally discrete numerical methods. In this section we develop a fully discrete numerical method to approximate problem (1.1), restricting ourselves to one space dimension; that is, Ω = [0,L], and supp(J)=[−S, S]. We begin by describing our space discretization. 2.1. Semidiscrete scheme. We perform a discretization of the space variable and show the convergence of the semidiscrete scheme and some asymptotic properties of the approximations. To this end we introduce some preliminary notation. Let us consider a partition (not necessarily uniform) 0 = x1<x 2<··· <x N+1 =Lof the interval [0,L]of size h(h=max(xi−xi−1)) and the corresponding standard piecewise linear finite element space Vh. Let us denote by {φj}1≤j≤N+1 the usual Lagrange basis of Vh.We assume that Nh =O(1) (this will be used to show convergence of the scheme) and that hS(this is needed in order to obtain a nontrivial scheme). We define the semidiscrete approximation, (2.1) uh(x, t)= N+1  j=1 uj(t)φj(x), where uj(t) is the solution to the following system of ODEs: (2.2) ⎧ ⎨ ⎩ u j(t)=  i∈Aj J(xi−xj)hi|ui(t)−uj(t)|p−2(ui(t)−uj(t)),t>0, uj(0) = u0(xj) Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. 2106 MAYTE P´ EREZ-LLANOS AND JULIO D. ROSSI for every j=1,...,N+1, being Aj={isuch that |xi−xj|≤S}and hi=xi−xi−1. Note that local existence and uniqueness of solutions to this system is straightforward from the fact that the right-hand side of the equation is locally Lipschitz. Also note that we can consider a discontinuous, piecewise constant approximation in space, but, since solutions to (1.1) with continuous initial datum are continuous, we prefer to use continuous approximations. We first show that the semidiscrete solutions converge to the continuous solution. To this end we prove that our solutions satisfy a comparison principle. Definition 2.1. We say that Uis a supersolution of problem (2.2) if each of its components satisfies ⎧ ⎨ ⎩ u j(t)≥ i∈Aj J(xi−xj)hi|ui(t)−uj(t)|p−2(ui(t)−uj(t)),t>0, uj(0) ≥u0(xj). Analogously, Uis a subsolution of problem (2.2) if it satisfies the previous problem with the reverse inequalities. Lemma 2.2. Let Uand Ube a supersolution and a subsolution of problem (2.2), respectively. Then, calling Uthe solution to (2.2), U(t)≥U(t)≥U(t). Here, by U(t)≥U(t)≥U(t)we mean that uj(t)≥uj(t)≥uj(t)for all j=1,...,N+ 1and for all t≥0. Proof. By an approximation procedure we restrict ourselves to consider strict inequalities for the supersolution. Indeed, we can take ui(t)+δt+δ(δ>0) for every i=1,...,N + 1 as a strict supersolution, and take limit as δ→0 at the end. We show that U(t)>U(t) by contradiction. Let us assume that there exists a first time t0and a node j0such that uj0(t0)=uj0(t0)=a;thenwehave 0≥u j0(t0)−u j0(t0) > i∈Aj hiJ(xi−xj0)|ui(t0)−a|p−2(ui(t0)−a)−|ui(t0)−a|p−2(ui(t0)−a) ≥0, a contradiction. The last inequality follows from ui(t0)≥ui(t0) for every i= 1,...,N + 1 and the fact that |a|p−2a−|b|p−2b≥0 if and only if a≥b. The inequality U(t)≥U(t) can be handled in a similar way. Note that constants, that is, ui≡cfor every iare solutions to the equation. Hence we have that maxiui(0) ≥U(t)≥miniui(0) and then we conclude that solutions are uniformly bounded and globally defined. For completeness, following [4], let us state what we understand by a solution to the continuous problem. Existence, uniqueness, a general contraction principle, and a comparison principle for solutions can be found in Chapter 6 of [4]. Definition 2.3. Asolutionof(1.1) in [0,T]is a function u∈C([0,T]×Ω) ∩W1,1((0,T); C(Ω)), which satisfies u(x, 0) = u0(x)for x∈Ωand ut(x, t)=Ω J(x−y)|u(y,t)−u(x, t)|p−2(u(y,t)−u(x, t)) dy in Ω×(0,T). Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. NUMERICAL APPROXIMATIONS FOR A NONLOCAL EQUATION 2107 Now we are ready to prove the convergence result. Theorem 2.4. Let u∈C([0,T]×Ω) ∩W1,1((0,T); C(Ω)), be a positive solution to (1.1) and let uhbe the numerical approximation defined by (2.1). Then there exists a constant K, such that for every h small enough it holds that max 0≤t≤Tmax x∈[0,L]|u(x, t)−uh(x, t)|≤Kh. Proof. Let us denote by εj(t)=uj(t)−vj(t) the error vector, where vj(t)= u(xj,t). It is easy to check that this vector satisfies ε j= i∈Aj J(xi−xj)hi|ui−uj|p−2(ui−uj)−|vi−vj|p−2(vi−vj) + i∈Aj J(xi−xj)hi|vi−vj|p−2(vi−vj) −Ω J(y−xj)|u(y)−u(xj)|p−2(u(y)−u(xj)) dy ≤ i∈Aj J(xi−xj)hi(p−1)(ηi,j )p−2(εi−εj)+O(h), being ηi,j an intermediate value between ui−ujand vi−vj, resulting from applying the mean value theorem to the first term. The second term is O(h), since each of the approximations to the integrals are O(h2) and we add them taking into account the nodes laying on the supp(J). Therefore, the error vector solves the following problem: (2.3) ⎧ ⎨ ⎩ ε j(t)=  i∈Aj J(xi−xj)hi(p−1)(ηi,j )p−2(εi−εj)+Cjh, t > 0, εj(0) = 0. Thus, if we consider ωj(t)=ω(t)=Cht, with C=max k,x,t Ck, for every j= 1,...,N +1,ωis a supersolution of problem (2.3), then indeed ω j(t)=Ch ≥ i∈Aj J(xi−xj)hi(p−1)(ηi,j )p−2(ωi−ωj)+Cjh. It can be shown in a similar way to Lemma 2.2, that this problem satisfies a comparison principle. Consequently, εj(t)≤ω(t)≤Kh, for every t≤Tand every j=1,...,N+ 1. Here K=CT. Repeating this argument for the vector −εj(t) we conclude that |εj(t)|≤Kh for every j=1,...,N +1, which finishes the proof. Now, let us show that some properties of the continuous problem are preserved by the numerical approximations. First, let us prove that, as it happens for the continuous problem, if the initial datum is symmetric, the numerical scheme preserves this property. Lemma 2.5. Let u0(x)be a symmetric initial datum in [−L, L]; that is, u0(x)= u0(−x), and consider a symmetric partition xi−N−1≤i≤N+1,of this interval, i.e., xi=−x−i. Then, the solution to problem (2.2) is symmetric, that is, it satisfies ui(t)=u−i(t). Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. 2108 MAYTE P´ EREZ-LLANOS AND JULIO D. ROSSI Proof. Let us define ωj(t)=u−j(t). Then, ωjsatisfies ⎧ ⎪ ⎨ ⎪ ⎩ ω j(t)=  i∈Aj J(xi−x−j)hi|ui(t)−ωj(t)|p−2(ui(t)−ωj(t)),t>0, ωj(0) = u0(x−j)=u0(xj) for every j=−1−N,...,N +1,beingAj={isuch that |xi−x−j|≤S}.But note that Aj=−Aj,whereAj={isuch that |xi−xj|≤S}since the partition is symmetric. Taking into account that Jis symmetric, we have that the previous equation can be written as follows: ω j(t)=  i∈Aj J(x−i−x−j)hi|ui(t)−ωj(t)|p−2(ui(t)−ωj(t)) = i∈Aj J(xi−xj)hi|ui(t)−ωj(t)|p−2(ui(t)−ωj(t)). We conclude the proof by uniqueness of the solutions to the ODE system. Now we show that the solutions of the discrete nonlocal problem converge as t→∞to the mean value of the initial condition, as happens for the continuous problem. To carry out this task, we prove the following Poincar´e-type inequality. Lemma 2.6. Let p≥1.Then, (2.4) Ip=I(J, d, p) = inf {v∈RN+1:N+1 i=1 hivi=0} 1 2N i,j=1 hihjJ(xi−xj)|vi−vj|p N+1 j=1 hj|vj|p is strictly positive. Moreover, (2.5) Ip N+1  j=1 hj vj−1 L N+1  i=1 hivi p ≤1 2 N+1  i,j=1 hihjJ(xi−xj)|vi−vj|p for every v∈RN+1. Proof. To show (2.4) it suffices to see that for each v∈RN+1 there exists a constant c>0 such that (2.6) N+1  j=1 hj|vj|p≤c⎛ ⎝ N+1  i,j=1 hihjJ(xi−xj)|vi−vj|p+ N+1  i=1 hivi⎞ ⎠. Arguing by contradiction, if (2.6) does not hold, then for every n∈Nthere exists vn∈RN+1 with N+1 i=1 hi|(vn)i|p= 1, satisfying 1>n⎛ ⎝ N+1  i,j=1 hihjJ(xi−xj)|(vn)i−(vn)j|p+ N+1  i=1 hi(vn)i⎞ ⎠. Hence lim n→∞ N+1  i=1 hi(vn)i=0 and lim n→∞ N+1  i,j=1 hihjJ(xi−xj)|(vn)i−(vn)j|p=0. Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. NUMERICAL APPROXIMATIONS FOR A NONLOCAL EQUATION 2109 This implies that lim n→∞ |(vn)i−(vn)j|p=0. This, together with limn→∞ N+1 i=1 hi(vn)i=0,contradicts the fact that N+1  i=1 hi|(vn)i|p=1. Note that (2.5) follows immediately from (2.4), for vj=ωj−1 LN+1 i=1 hiωi, for any ω∈RN+1. We are now ready to study the asymptotic behavior of the discrete solutions. Theorem 2.7. Let p≥1and uh∈Vhbe the solution to problem (2.2) corresponding to an initial datum, u0. Then, there exists a positive constant C, independent of t, such that (2.7) N+1  j=1 hj|uj(t)−u0|p≤CN+1 j=1 hj|uj(0)|2 t→0,as t→∞, where u0=1 Lihiui(0) <∞. Proof. We define ωj(t)=uj(t)−u0. From (2.2) it is easy to see that ωsatisfies d dt N  j=1 hj|ωj|p =p N+1  j=1 hj|ωj|p−2ωj i∈Aj hiJ(xi−xj)|ωi−ωj|p−2(ωi−ωj) =−p N+1  j=1  i∈Aj hjhiJ(xi−xj)|ωi|p−2ωi−|ωj|p−2ωj|ωi−ωj|p−2(ωi−ωj) +p N+1  j=1  i∈Aj hjhiJ(xi−xj)|ωi|p−2ωi|ωi−ωj|p−2(ωi−ωj) =−p 2 N+1  j=1  i∈Aj hjhiJ(xi−xj)|ωi|p−2ωi−|ωj|p−2ωj|ωi−ωj|p−2(ωi−ωj). Therefore, N+1 j=1 hj|ωj|pis nonincreasing. On the other hand, note that adding the equation in (2.2) top in all of the nodes, we obtain N+1  j=1 hju j(t)= N+1  j=1  i∈Aj hjJ(xi−xj)hi|ui(t)−uj(t)|p−2(ui(t)−uj(t)) =1 2 N+1  j=1  i∈Aj hjJ(xi−xj)hi|ui(t)−uj(t)|p−2(ui(t)−uj(t)) +1 2 N+1  i=1  j∈Aj hjJ(xi−xj)hi|ui(t)−uj(t)|p−2(ui(t)−uj(t)) = 0. Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. 2110 MAYTE P´ EREZ-LLANOS AND JULIO D. ROSSI That is, the following quantity is preserved, 1 LN+1 j=1 hjuj(t)=u0, for every t≥0 and, consequently, 1 LN+1 j=1 hjωj(t)=0,foreveryt≥0. Thus we can apply to this function the Poincare’s inequality (2.5) and we obtain Ip N+1  j=1 hj|ωj|p≤1 2 N+1  i,j=1 hihjJ(xi−xj)|ωi−ωj|p, from which, using the monotonicity of jhj|ωj|p, it easily follows that, (2.8) t N+1  j=1 hj|ωj(t)|p≤t 0 N+1  j=1 hj|ωj(s)|pds ≤Ct 0 N+1  i,j=1 hihjJ(xi−xj)|ωi(s)−ωj(s)|pds =Ct 0 N+1  i,j=1 hihjJ(xi−xj)|ui(s)−uj(s)|pds. Now, we multiply (2.2) top by ujand we add in all the nodes, it is easy to check that 1 2 d dt ⎛ ⎝ N+1  j=1 hj|uj(t)|2⎞ ⎠= N+1  i,j=1 J(xi−xj)hihj|ui(t)−uj(t)|p−2(ui(t)−uj(t))uj(t) =− N+1  i,j=1 J(xi−xj)hihj|ui(t)−uj(t)|p + N+1  i,j=1 J(xi−xj)hihj|ui(t)−uj(t)|p−2(ui(t)−uj(t))ui(t) =−1 2 N+1  i,j=1 J(xi−xj)hihj|ui(t)−uj(t)|p. Integrating in time, we deduce that N+1  j=1 hj|uj(t)|2− N+1  j=1 hj|uj(0)|2=−t 0 N+1  i,j=1 hihjJ(xi−xj)|ui(s)−uj(s)|pds. Thus, t 0 N+1  i,j=1 hihjJ(xi−xj)|ui(s)−uj(s)|pds ≤ N+1  j=1 hj|uj(0)|2, which plugged into (2.8) gives the desired conclusion, (2.7). 2.2. A fully discrete scheme. Now we perform the discretization in time, using the explicit Euler method; that is, (2.9) ⎧ ⎪ ⎨ ⎪ ⎩ uk+1 j−uk j τk = i∈Aj J(xi−xj)hi|uk i−uk j|p−2(uk i−uk j),k>0, u0 j=u0(xj), Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. NUMERICAL APPROXIMATIONS FOR A NONLOCAL EQUATION 2111 being Aj={i:|xi−xj|≤S}, for every j=1,...,N +1. WedenotebyUk= (uk 1,...,u k N+1) the vector whose components solve the previous system (2.9). The symmetry property given in Lemma 2.5 follows similarly from reflection and uniqueness for the totally discrete method. However, the comparison principle for problem (2.9) requires a condition on the time step τk, which is not restrictive. In order to show this comparison principle we give the following definitions. Definition 2.8. We say that Ukis a supersolution of problem (2.9) if each of its components satisfies (2.10) ⎧ ⎪ ⎨ ⎪ ⎩ uk+1 j−uk j τk ≥ i∈Aj J(xi−xj)hi|uk i−uk j|p−2(uk i−uk j),k>0, u0 j≥u0(xj). Analogously, Ukis a subsolution of problem (2.9) if it satisfies (2.10) with the reverse inequalities. Proposition 2.9. Let Ukand Ukbe a super- and a subsolution of problem (2.9), respectively. Note that this implies U0≤U0.Ifhis small and the time step satisfies (2.11) τk<1 2(p−1) maxj|uk j|p−2, then Uk≤Ukfor every k>0. Proof. We define Zk=Uk−Ukand call zk jbe the components of Zk.Notethat, as before, by an approximation argument, we can assume in (2.10) strict inequalities. Therefore, Zksatisfies for each j=1,...,N+1 and every k>0 the following system: ⎧ ⎪ ⎨ ⎪ ⎩ zk+1 j−zk j τk > i∈Aj J(xi−xj)hi|uk i−uk j|p−2(uk i−uk j)−|uk i−uk j|p−2(uk i−uk j), z0 j>0. Let us argue by contradiction and suppose that there exists a first time tk+1 and a node xjsuch that zk+1 j≤0 while zk j>0. From the corresponding equation for that node we obtain zk+1 j>z k j+τk i∈Aj J(xi−xj)hi|uk i−uk j|p−2(uk i−uk j)−|uk i−uk j|p−2(uk i−uk j) =zk j+(p−1)τk i∈Aj J(xi−xj)hi|θk i,j |p−2(zk i−zk j), being θk i,j a value between uk i−uk jand uk i−uk j.Letηk j=max i{θk i,j}. Neglecting the positive terms in the previous inequality and taking (2.11) into account, we get zk+1 j>z k j⎛ ⎝1−(p−1)τk|ηk j|p−2 i∈Aj J(xi−xj)hi⎞ ⎠ =zk j1−(p−1)τk|ηk j|p−2(1 + O(h))≥0, which is a contradiction. Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. 2118 MAYTE P´ EREZ-LLANOS AND JULIO D. ROSSI Since any term in the above sum is nonnegative, we conclude that λ+1≥1 2pk hihjJ(xi−xj)|(upk)i−(upk)j|pk for some i, j such that |xi−xj|≤S. Hence, since (upk)i→uiand (upk)j→ujwe obtain |ui−uj|≤1for|xi−xj|<S. Thus, we conclude that U∈K∞. This ends the proof. Note that our numerical schemes can be used to approximate solutions to the nonlocal sandpile model (3.4). In fact, given ε>0, from (3.3), we have that there exists p0large such that sup t∈[0,T] up0(·,t)−u∞(·,t)L2(Ω) <ε/2. Then, for hsmall enough, sup t∈[0,T] up0,h(·,t)−u∞(·,t)L2(Ω) ≤|Ω|1/2sup t∈[0,T] up0,h(·,t)−up0(·,t)L∞(Ω) +sup t∈[0,T] up0(·,t)−u∞(·,t)L2(Ω) ≤ε. 3.2. An explicit solution. In this subsection we present an explicit example that illustrates the behavior of the solutions when p=+∞. We want to find a solution to (3.8) f(t)−U(t)∈∂G∞(U(t)) a.e. t∈(0,T), U(0) = u0(xj)inΩ, where G∞(v)=0if|vi−vj|≤1,for |xi−xj|<S, +∞in other case. In order to satisfy that a function U:[0,T]→ RN+1 is a solution to (3.8) we need to check that (3.9) G∞(v)≥G∞(U)+f−Ut,v−U∀v∈RN+1. To this end we can assume that v∈K∞(otherwise G∞(v)=+∞and then (3.9) becomes trivial). By (3.1), we need to check that U(t)∈K∞and  j (fj(t)−u j(t))(vj−uj(t)) ≤0 for every v∈K∞. Now, to simplify, assume that the support of Jis (−1,1); that is, S=1. Let us consider a recipient Ω = (0,L) with Lan integer greater than 1, a uniform mesh, xj=(j−1)L/N,j=1,...,N + 1, as initial datum take u0= 0 and a source given by f(x, t)=χ[0,1](x). Then the solution is given by uj(t)=t, xj≤1, 0,x j>1, Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. NUMERICAL APPROXIMATIONS FOR A NONLOCAL EQUATION 2119 for times t∈[0,1]. Let us check that this is indeed the solution for t∈[0,1]. It is clear that U(t)∈K∞since |uj(t)−ui(t)|≤1 for every i, j.Moreover,sincefand Ut coincide in [0,L], we have that (3.9) is immediate. For t∈[1,3] we get uj(t)=⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ 1+t−1 2xj∈[0,1), t−1 2xj∈[1,2), 0xj/∈[0,2). Let us check that this is the solution. We have that U(t)∈K∞. Hence we have to verify (3.9). Since Utand fvanish outside [0,2) we have to see that  xj∈[0,1) (1 −u j(t))(vj−uj(t)) +  xi∈[1,2) (−u i(t))(vi−ui(t)) ≤0. Using the explicit form of ujwe get  xj∈[0,1) 1 2vj−1+t−1 2− xi∈[1,2) 1 2vi−t−1 2≤0. That is equivalent to (note that all the terms that involve tcancel)  xj∈[0,1) vj− xi∈[1,2) vi≤ xj∈[0,1) 1, which holds as a consequence of the fact that |vj−vi|≤1forxj∈[0,1) and xi∈[1,2) (we are using here that v∈K∞). In general we have, until the recipient is full, that is, until uj(t)>0 for every t, for any k=1,...,L and for t∈[tk−1,t k) uj(t)= ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ k−1+t−tk−1 kxj∈[0,1), k−2+t−tk−1 kxj∈[1,2), ··· t−tk−1 kxj∈[k−1,k), 0xj/∈[0,k). Here tk=tk−1+kis the first time at which the solution reaches level k;thatis, uh(tk,0) = k. For times even greater, t≥tL=L(L+1)/2 (the recipient is full of sand from that time, that is, the solution is strictly positive for t>t L=L(L+1)/2), the solution turns out to be uj(t)= ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ L+t−tL Lxj∈[0,1), L−1+t−tL Lxj∈[1,2), ··· 1+t−tL Lxj∈[L−1,L). Hence, when the recipient is full the solution grows with speed 1/L uniformly in (0,L). Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. 2120 MAYTE P´ EREZ-LLANOS AND JULIO D. ROSSI −4 −2 0240 10 20 30 0 5 10 time space p=5 Fig. 1.Evolution in time, symmetric datum. 0 1000 2000 3000 4 4.5 5 5.5 6 N=21 N=41 N=61 mean value Fig. 2.Convergence to the mean value of the initial datum as tincreases. 4. Numerical experiments. In this section we perform, using MATLAB (ode15s subroutine), some numerical experiments just to illustrate our general results. First, we show the evolution in time of a numerical solution, taking p=5,N= 101, and as initial datum u0=6|sin(π 3x)|+0.15(3 −|x|)2in the interval [−3,3]. We can appreciate that, since the initial datum is symmetric, the solution remains symmetric for every positive time; see Figure 1. Next, take u0=9−x2,p=5intheinterval[−3,3]. We show the value of the first node (corresponding to x=−3) of the numerical approximation for different values of Nas tincreases; see Figure 2. We can appreciate the convergence towards the mean value of the numerical initial datum as tincreases for different values of N. Moreover, we can appreciate that the discrete mean value increases with Nto 6, which is the continuous mean value of u0=9−x2. The same behavior can be observed for any other node. Now we include a picture comparing the approximation to the semidiscrete solution using ode15s with the totally discrete one. Here we take p=4,u0(x)=(9−x2), Ω=(−3,3), N= 50, τk=0.1, and k= 8000. See Figure 3. In the next picture (Figure 4) we show the evolution in time of a numerical model of sandpiles with p= 50, corresponding to the initial data u0=0inΩ=(0,4) with asourcef=χ[0,1]. Note that the slope of uj(t) decreases as time evolves, as we Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. NUMERICAL APPROXIMATIONS FOR A NONLOCAL EQUATION 2121 −3 −2 −1 0 1 2 3 5.85 5.9 5.95 6 6.05 6.1 totally discrete solution semidiscrete solution Fig. 3.Semidiscrete (approximated with ode15s) and totally discrete solutions for tk=τkk. 0 2 40 5 10 0 1 2 3 4 5 6 p=50 time space Fig. 4.Evolution in time of the numerical model for sandpiles. 0 2 4 6 0 0.5 1 1.5 2 2.5 3 3.5 4 u8(t) u18(t) u28(t) Fig. 5.Different slopes of the value of the solution at some nodes of the numerical model for sandpiles. illustrate in Figure 5. This is a consequence of the fact that the sand that is added by the source has to be distributed in a larger set. Finally, in Figure 6, we show solutions for different values of pand compare them with the explicit formula of the previous section. In this example we take u0=0inΩ=(0,3) with a source f=χ[0,1]. Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. 2122 MAYTE P´ EREZ-LLANOS AND JULIO D. ROSSI 0 1 2 3 1 2 3 4 5 6 p=8 p=30 p=50 0 1 2 3 1 2 3 4 5 6 explicit solution Fig. 6.Convergence as pincreases to the explicit solution. Acknowledgments. Part of this work was performed during a visit of MPL and JDR to Univ. de Valencia, who the authors want to thank for the warm hospitality and the stimulating working atmosphere found there. REFERENCES [1] F.Andreu,J.M.Maz ´ on,J.D.Rossi,andJ.Toledo,The Neumann problem for nonlocal nonlinear diffusion equations, J. Evol. Equ., 8 (2008), pp. 189–215. [2] F.Andreu,J.M.Mazon,J.D.Rossi,andJ.Toledo,The limit as p→∞in a nonlocal p-Laplacian evolution equation: A nonlocal approximation of a model for sandpiles,Calc. Var. Partial Differential Equations, 35 (2009), pp. 279–316. [3] F.Andreu,J.M.Maz ´ on, J. D. Rossi, and J. Toledo,Anonlocalp-Laplacian evolution equation with Neumann boundary conditions, J. Math. Pures Appl. (9), 90 (2008), pp. 201– 227. [4] F. Andreu, J. M. Maz´ on,J.D.Rossi,andJ.Toledo,Nonlocal Diffusion Problems,AMS. Math. Surveys Monogr. 165, AMS, Providence, RI, 2010. [5] P. Bates and A. Chmaj,An integrodifferential model for phase transitions: Stationary solutions in higher dimensions, J. Statist. Phys., 95 (1999), pp. 1119–1139. [6] P. W. Bates, P. C. Fife, X. Ren, and X. Wang,Travelling waves in a convolution model for phase transitions, Arch. Rational Mech. Anal., 138 (1997), pp. 105–136. [7] H. Brezis,Op´erateur Maximaux Monotones et Semi-groupes de Contractions dans les Espaces de Hilbert, North–Holland, Amsterdam, London, 1973. [8] H. Br´ ezis and A. Pazy,Convergence and approximation of semigroups of nonlinear operators in Banach spaces, J. Functional Analysis, 9 (1972), pp. 63–74. [9] A. Buades, B. Coll, and J. M. Morel,Neighborhood filters and PDE’s,Numer.Math.,105 (2006), pp. 1–34. [10] C. Carrillo and P. Fife,Spatial effects in discrete generation population models,J.Math. Biol., 50 (2005), pp. 161–188. [11] E. Chasseigne, M. Chaves, and J. D. Rossi,Asymptotic behavior for nonlocal diffusion equations, J. Math. Pures Appl. (9), 86 (2006), pp. 271–291. [12] X. Chen,Existence, uniqueness and asymptotic stability of travelling waves in nonlocal evolution equations, Adv. Differential Equations, 2 (1997), pp. 125–160. [13] X. Chen and M. Gunzburguer,Continuous and discontinuous finite elements methods for a peridynamics model of mechanics, Comput. Methods Appl. Mech. Engrg., 200 (2011), pp. 1237–1250. [14] C. Cortazar, M. Elgueta, and J. D. Rossi,A nonlocal diffusion equation whose solutions develop a free boundary, Ann. Henri Poincar´e, 6 (2005), pp. 269–281. [15] C. Cortazar, M. Elgueta, J. D. Rossi, and N. Wolanski,Boundary fluxes for nonlocal diffusion, J. Differential Equations, 234 (2007), pp. 360–390. [16] C. Cortazar, M. Elgueta, J. D. Rossi, and N. Wolanski,How to approximate the heat equation with Neumann boundary conditions by nonlocal diffusion problems, Arch. Rational Mech. Anal., 187 (2008), pp. 137–156. Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. NUMERICAL APPROXIMATIONS FOR A NONLOCAL EQUATION 2123 [17] Q. Du, M. Gunzburger, R. Lehoucq, and K. Zhou,A nonlocal vector calculus, nonlocal volume-constrained problems, and nonlocal balance laws,preprint. [18] Q. Du and K. Zhou,Mathematical and numerical analysis of linear peridynamic models with nonlocal boundary conditions, SIAM J. Numer. Anal., 48 (2010), pp. 1759–1780. [19] I. Ekeland and R. Temam,Convex Analysis and Variational Problems, North–Holland, Amsterdam, Oxford, 1972. [20] L. C. Evans, M. Feldman, and R. F. Gariepy,Fast/slow diffusion and collapsing sandpiles, J. Differential Equations, 137 (1997), pp. 166–209. [21] L. C. Evans and F. Rezakhanlou,A stochastic model for growing sandpiles and its continuum limit, Comm. Math. Phys., 197 (1998), pp. 325–345. [22] M. Feldman,Growth of a sandpile around an obstacle, Monge Amp`ere Equation: Applications to Geometry and Optimization (Deerfield Beach, FL, 1997), Contemp. Math. 226, American Mathematical Society, Providence, RI, 1999, pp. 55–78. [23] R. Ferreira, A. de Pablo, and M. P´ erez-Llanos,Numerical blow-up for the p-Laplacian equation with a source, Comput. Methods Appl. Math., 5 (2005), pp. 137–154. [24] P. Fife,Some nonclassical trends in parabolic and parabolic-like evolutions, Trends in Nonlinear Analysis, Springer, Berlin, 2003, pp. 153–191. [25] P. Fife and X. Wang,A convolution model for interfacial motion: The generation and propagation of internal layers in higher space dimensions, Adv. Differential Equations, 3 (1998), pp. 85–110. [26] S. Kindermann, S. Osher, and P. W. Jones,Deblurring and denoising of images by nonlocal functionals, Multiscale Model. Simul., 4 (2005), pp. 1091–1115. [27] U. Mosco,Convergence of convex sets and solutions of variational inequalities, Advances Math., 3 (1969), pp. 510–585. [28] M. L. Parks, R. B. Lehoucq, S. Plimpton, and S. Silling,Implementing peridynamics within a molecular dynamics code, Computer Physics Comm., 179 (2008), pp. 777–783. [29] S. A. Silling,Reformulation of elasticity theory for discontinuities and long-range forces,J. Mech. Phys. Solids, 48 (2000), pp. 175–209. [30] S. A. Silling and R. B. Lehoucq,Convergence of peridynamics to classical elasticity theory, J. Elasticity, 93 (2008), pp. 13–37. [31] L. Silvestre,H¨older estimates for solutions of integro-differential equations like the fractional Laplace, Indiana Univ. Math. J., 55 (2006), pp. 1155–1174. [32] X. Wang,Metastability and stability of patterns in a convolution model for phase transitions, J. Differential Equations, 183 (2002), pp. 434–461. [33] L. Zhang,Existence, uniqueness and exponential stability of traveling wave solutions of some integral differential equations arising from neuronal networks,J.DifferentialEquations, 197 (2004), pp. 162–196. Reproducedwithpermissionofthecopyrightowner.Furtherreproductionprohibitedwithoutpermission.