Error analysis of some nonlocal diffusion discretization schemes
EError analysis of some nonlocal diffusion discretization schemes ∗ Gonzalo Galiano † Abstract
We study two numerical approximations of solutions of nonlocal diffusion evolution prob-lems which are inspired in algorithms for computing the bilateral denoising filtering of animage, and which are based on functional rearrangements and on the Fourier transform.Apart from the usual time-space discretization, these algorithms also use the discretizationof the range of the solution (quantization). We show that the discrete approximations con-verge to the continuous solution in suitable functional spaces, and provide error estimates.Finally, we present some numerical experiments illustrating the performance of the algo-rithms, specially focusing in the execution time.
Keywords:
Numerical approximation, nonlocal diffusion, functional rearrangement, Fouriertransform, error analysis.
Let
T > ⊂ R d ( d ≥
1) be an open and bounded set with Lipschitz continuous boundary,and consider the following problem: Find u : [0 , T ] × Ω → R + such that ∂ t u ( t, x ) = (cid:90) Ω J ( x , y ) A ( u ( t, y ) − u ( t, x )) d y + f ( t, x , u ( t, x )) , (1) u (0 , x ) = u ( x ) , (2)for ( t, x ) ∈ Q T = (0 , T ] × Ω, and for some u : Ω → R . The functions J and A are called,respectively, the spatial and the range kernels.Nonlocal diffusion problems of the type (1)-(2) have many applications in modelling diffusionprocesses in which the interactions among particles occur at a distance, in contrast to localdiffusion modelled by partial differential equations (PDE).They appear in a number of fields, among which image processing [6, 19] and populationdynamics [22, 21] are specially relevant. For the latter, the usual choice of the monotone rangekernel A ( s ) = | s | p − s , with p ≥
1, gives rise to the so called nonlocal p –Laplacian due to somesimilarities to the usual evolution p –Laplacian PDE. For the former, non-monotone but Lipschitzcontinuous range kernels like the Gaussian function are more frequently used. In both of them,the spatial kernel is assumed to be, at least, an integrable non-negative symmetric function.In addition to the well-posedness theory, see the monograph by Andreu et al. [2] for themonotone case and [17] for the Lipschitz continuous case, the numerical discretization of thistype of problems has also been investigated, being most of the analyzed methods related to the ∗ Supported by Spanish MCI Project MTM2017-87162-P. † Department of Mathematics, University of Oviedo, Spain ( [email protected] ). a r X i v : . [ m a t h . NA ] F e b . Galianousual techniques for PDE. For instance, P´erez-Llano and Rossi [24] analyze a discrete spatialdecomposition in terms of Lagrange basis while Tiang and Du [28] study, in addition, thestandard finite element approach. In [9], Du et al. consider a discontinuous Galerkin method.See the recent monographs by D’Elia et al. for a review on finite element based methods [8] andother methods [7] that also apply to non-integrable spatial kernels.In this article we study some discretization techniques which were introduced in the imageprocessing literature without performing a rigorous mathematical analysis. The purpose ofthese techniques is the fast computation of the approximations, being its accuracy relegated asof secondary importance. As stated by Weiss [30], the fundamental property that concerns usis the runtime per pixel, as a function of the spatial kernel support.The first method we analyze is based on the work of Weiss [30], later improved by Porikli [25]and generalized in [14, 15]. The main idea is changing the space of integration in (1) by usingthe coarea formula. In the simplest discrete situation, the problem is to compute efficiently thefilter F u [ i ] = (cid:88) j ∈ J J [ i − j ] A ([ u [ j ] − u [ i ]) , (3)where J ⊂ Z is the set of pixels and u : J → Q = { , , . . . , } is an intensity image. When J is a box kernel (taking a constant positive value inside a small square and zero outside) thenthe filter may be computed using the local histogram , h i , of u in the box centered at each pixel i ∈ J , F u [ i ] = (cid:88) k ∈ Q h i ( k ) A ( k − u [ i ]) , (4)that is, changing the space of integration from the pixel space to the image value space. In thecontinuous case, and for general integrable spatial kernels, this change involves in a natural waythe coarea formula, leading to the consideration of some functional rearrangements involvingthe image (decreasing rearrangement) and the spatial kernel (relative rearrangement). We givethe details in the Appendix.The second method was introduced by Durand and Dorsey [10] and later computationallyimproved by Yang et al. [31, 32]. Considering the collection of convolution filters F k u [ i ] = (cid:88) j ∈ J J [ i − j ] A ([ u [ j ] − k ) , for k ∈ Q, (5)one observes that F u [ i ] = F k u [ i ] if u [ i ] = k . The convolution filters may be efficiently computedby means of the fast Fourier transform. Moreover, a natural and fast approximation to the exactfilter may be defined by calculating the convolution filters only for a small subset of Q and thendefining the approximation through linear interpolation.For image processing tasks, specially for image denoising, filters involving expressions like(3) are known as bilateral filters [33, 27, 29, 3, 11, 5, 19] and are mostly applied in a singlediscrete time step, although other tasks like image segmentation and saliency detection mayrequire more steps [14, 16, 18]. For other applications, these algorithms may result interestingwhen accuracy is secondary to execution time, e.g. for real time computations.The error analysis we perform in this article fills a gap, already pointed out by Durand andDorsey [10], with the aim of establishing a rigorous mathematical setting for these algorithms,2rror analysis of some nonlocal diffusion discretization schemesusually formulated in the discrete case and for concrete spatial and range kernels, mostly thebox and the Gaussian kernels. With a well-posedness theory at hand [2, 17], we extend the ideasleading to the filters (4) and (5) to the continuous evolution case of the form (1)-(2), which apartof adding a reaction term, it also allows the consideration of general spatial and range kernels.More concretely, we assume the following hypothesis on the data (H) :1. The spatial kernel J ∈ BV ( R d × R d ) is symmetric and non-negative. We also assume (cid:90) Ω × Ω J ( x , y ) d x d y = 1 . (6)2. The range kernel A ∈ C ,α ( R ) with α ∈ (0 ,
1) satisfies A ( − s ) = − A ( s ) , A ( s ) s ≥ s ∈ R , and either A ∈ W , ∞ loc ( R ) or A is non-decreasing in R .3. The reaction function f ∈ (cid:0) L ∞ ( Q T ) ∩ BV ( Q T ) (cid:1) × W , ∞ loc ( R ) satisfies | f ( · , · , s ) | ≤ C f (1 + | s | ) in Q T , for s ∈ R , and for some constant C f > .
4. The initial datum u ∈ L ∞ (Ω) ∩ BV (Ω).Under the set of assumptions (H), the existence and uniqueness of a solution, u , of (1)-(2) suchthat u ∈ W , ∞ (0 , T ; L ∞ (Ω)) ∩ C ([0 , T ]; L ∞ (Ω) ∩ BV (Ω)), is proven in [17, Theorem 1]. Inparticular, it is shown that (cid:107) u (cid:107) L ∞ ( Q T ) ≤ C M , (7)where C M only depends on T , (cid:107) u (cid:107) L ∞ (Ω) and C f . In addition, if u ≥ f ( · , · , ≥ Q T (8)then the solution is non-negative in Q T . In the rest of the paper, we assume (8), for simplicity.Let us make some observations. On one hand, the well-posedness of problem (1)-(2) is provenin [17] under more general conditions than those stated in (H). For instance, condition (6) maybe removed. We may also include a time-space dependence on the range kernel, or a nonlocalreaction term in place or in addition to the local reaction term, f . With the aim of keepingsome notational simplicity we have opted to deal with a simpler model, being the techniquesdescribed in this article easy to extend to the general case.On the other hand, if the spatial regularity of J , f and u is improved, for instance to a W ,p Sobolev space, with p ∈ (1 , ∞ ], then the corresponding spatial regularity of the solutionis improved to the same space. Notice that one of the most remarkable properties of nonlocaldiffusion evolution problems, when compared to their local diffusion versions, is the lack of aregularizing effect on the solutions. This is, the solutions are not more regular than the data.Finally, observe that the important limit case of the nonlocal 1–Laplacian equation is notcovered by the assumptions (H). However, a proof of its well-posedness may be found in [2].In contrast, (H) allows to handle the case A ( s ) = | s | p ( s ) − s , this is, the nonlocal version of the3. Galianoproblem ∂ t u − div (cid:0) |∇ u | p ( ∇ u ) − ∇ u (cid:1) = f , already employed in the literature [4], but for whichthe well-posedness is an open problem.The contents of the rest of the article is the following. In Section 2 we introduce the gen-eral discretization scheme, based on an explicit time Euler method and a piecewise constantapproximation in space. Since the only difference among the methods is in the quadrature ofthe integral term in (1), we can give in this section the error estimates for the other terms. InSection 3 we provide the analysis of the three quadrature methods considered in the article: thepointwise method ( Ptw ), the method based in functional rearrangements ( RR ) and the methodbased on the Fourier transform ( FFT ). Finally, in Section 4 we perform some numerical exper-iments illustrating the performance of each method for different choices of the space and rangekernels and for a variety of discretization parameters.
The proof of existence of solutions in [17] suggests the use of an explicit scheme to discretize(1)-(2) due to the absence of a stability constraint for the time step in terms of the spatial meshsize. This is also pointed out in [24].Thus, all the fully discrete schemes considered in this article are based on the followingsemi-discrete scheme for problem (1)-(2): Set τ = T /N , for some N ∈ N and t n = τ n for n = 0 , , . . . , N . Then, for x ∈ Ω, and n = 0 , , . . . , N −
1, we set u ( x ) = u ( x ) and u n +1 ( x ) = u n ( x ) + τ (cid:0) A ( u n )( x ) + f n ( x , u n ( x )) (cid:1) , where f n is some approximation of f ( t, · , · ) for t ∈ [ t n , t n +1 ) and where, for v ∈ L ∞ (Ω) and x ∈ Ω, we define A ( v )( x ) = (cid:90) Ω J ( x , y ) A ( v ( y ) − v ( x )) d y . (9)We assume that Ω ⊂ R d is the hyper-cube Ω = (0 , L ) d , with L ∈ R and consider the uniformmesh of Ω of size h , M h (Ω) = { x j } j ∈ J , given by x j = (cid:0) j − (cid:1) h, for j ∈ J = { ( j , . . . , j d ) , j i = 1 , . . . , J Ω } , (10)where ∈ R d is given as = (1 , . . . ,
1) and J Ω = L/h is assumed to be integer. Thus, thecollection of hypercubes of volume | Ω j | = h d centered at x j Ω j = { x ∈ Ω : max k =1 ,...,d | x k − ( x j ) k | < h/ } induces a partition of Ω.Let v h be a generic (G) piecewise constant approximation of v ∈ L (Ω) defined by v G h ( x ) = v G h [ j ] if x ∈ Ω j , (11)where, abusing on notation, we use the same symbol, v G h , for the function v G h : Ω → R , whoseevaluation at x ∈ Ω is denoted by v G h ( x ), and for the constant value this function takes in4rror analysis of some nonlocal diffusion discretization schemesΩ j , which is denoted by v G h [ j ]. The value v G h [ j ] will be specified for each discretization method.Similarly, assume the following generic forms for the approximations of A and f : A G h ( v G h )( x ) = A G h ( v G h )[ j ] if x ∈ Ω j , and f G τh ( t, x , s ) = f G n j ( s ) (12)if ( t, x ) ∈ [ t n , t n +1 ) × Ω j . Then, we consider the discrete scheme given by, for j ∈ J and n = 0 , . . . , N − u [ j ] = u G h [ j ] ,u n +1 [ j ] = u n [ j ] + τ (cid:0) A G h ( u τh ( t n , · ))[ j ] + f G n j ( u n [ j ]) (cid:1) , (13)where u τh is the piecewise constant interpolator of the collection { u n [ j ] } j ∈ J , defined as u τh ( t, x ) = u n [ j ] if ( t, x ) ∈ [ t n , t n +1 ) × Ω j . We also introduce here, for future reference, the corresponding time piecewise linear interpolator u τh ( t, x ) = u n +1 [ j ] + t j +1 − tτ ( u n [ j ]) − u n +1 [ j ]) if ( t, x ) ∈ [ t n , t n +1 ) × Ω j , and fix the initial datum approximations as u τh (0 , · ) = u τh (0 , · ) = u G h in Ω.Observe that, for notational consistency, we should write u G τh and u n G instead of u τh and u n . For the sake of clearness, we drop this index from the approximating functions in the hopethat the context will avoid any confusion. We shall deduce, for each discretization method, error bounds in the usual norm of L ∞ (0 , T ; L (Ω)).Subtracting the equations satisfied by u and u τh we get, for ( t, x ) ∈ [ t n , t n +1 ) × Ω j , ∂ t (cid:0) u ( t, x ) − u τh ( t, x ) (cid:1) = A ( u ( t, · ))( x ) − A G h ( u τh ( t, · ))( x )+ f ( t, x , u ( t, x )) − f G τh ( t, x , u τh ( t, x )) . Integrating in Q t and using the triangle inequality, we obtain the estimate (cid:90) Ω | u ( t, x ) − u τh ( t, x ) | d x ≤ (cid:90) Ω | u ( x ) − u G h ( x ) | d x + (cid:90) Ω | u τh ( t, x ) − u τh ( t, x ) | d x + (cid:90) Q t |A ( u ( s, · ))( x ) − A ( u τh ( s, · ))( x ) | d x ds + (cid:90) Q t |A ( u τh ( s, · ))( x ) − A G h ( u τh ( s, · ))( x ) | d x ds + (cid:90) Q t | f ( s, x , u ( s, x )) − f ( s, x , u τh ( s, x )) | d x ds + (cid:90) Q t | f ( s, x , u τh ( s, x )) − f G τh ( s, x , u τh ( s, x )) | d x ds =: I + · · · + I . (14)From the assumptions (H) we easily deducemax { I , I } ≤ C (cid:90) Q t | u ( s, x ) − u τh ( s, x ) | d x ds. (15)The remaining terms of (14) will be estimated in the following section, according to the type ofapproximation considered. 5. Galiano Definition of the approximation.
We start fixing the values of the approximations definedin (11) and (12). We set G = P ( Pointwise ), and v P h [ j ] = 1 h d (cid:90) Ω j v ( x ) d x , for j ∈ J , (16) f P n j ( s ) = 1 τ h d (cid:90) t n +1 t n (cid:90) Ω j f ( t, x , s ) d x dt, for n = 0 , . . . , N, j ∈ J , s ∈ R , (17) A P h ( v P h )[ j ] = h d (cid:88) k ∈ J J P h [ j , k ] A ( v P h [ k ] − v P h [ j ]) , for j ∈ J , with J P h [ j , k ] = 1 h d (cid:90) Ω j (cid:90) Ω k J ( x , y ) d y d x , for j , k ∈ J . (18)Observe that the definitions of A P h and J P h are such that A P h ( v P h )[ j ] = 1 h d (cid:90) Ω j A ( v P h )( x ) d x , if x ∈ Ω j . The following lemma allows to estimate the approximation error between a BV function and itsconstantwise approximation defined in terms of its mean values in Ω j , for j ∈ J . Lemma 1.
Let v ∈ BV (Ω) and consider the function v P h defined by (11) and (16). Then thereexists a constant C > , independent of v and h , such that (cid:107) v P h − v (cid:107) L (Ω) ≤ Ch | Dv | (Ω) , (19) where | Dv | (Ω) is the variation of v in Ω . Similarly, we have (cid:107) J P h − J (cid:107) L (Ω × Ω) ≤ Ch | DJ | (Ω × Ω) , (20) (cid:107) f P τh ( · , · , s ) − f ( · , · , s ) (cid:107) L ( Q T ) ≤ C ( τ + h ) | Df ( · , · , s ) | ( Q T ) , (21) where J P h ( x , y ) = J P h [ j , k ] if x ∈ Ω j and y ∈ Ω k . Proof.
According to Theorem 3.9 of [1], if v ∈ BV (Ω) then there exists a sequence { v ε } ⊂ C ∞ (Ω) such that v ε → v in L (Ω) and (cid:107)∇ v ε (cid:107) L (Ω) → | Dv | (Ω) as ε →
0. Let v εh : Ω → R bedefined by v εh ( x ) = 1 h d (cid:90) Ω j v ε ( x ) d x if x ∈ Ω j . On one hand, we have (cid:90) Ω | v εh ( x ) − v P h ( x ) | d x = (cid:88) j ∈ J (cid:90) Ω j h d (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:90) Ω j ( v ε ( x ) − v ( x )) d x (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ≤ (cid:90) Ω | v ε ( x ) − v ( x ) | d x → ε →
0. On the other hand, since v εh | Ω j is the mean value of v ε in Ω j , Poincar´e’s inequality inthe square Ω j of diameter h gives (cid:107) v εh − v ε (cid:107) L (Ω j ) ≤ Ch (cid:107)∇ v ε (cid:107) L (Ω j ) , with C independent of h and v ε . Summing for j ∈ J , we deduce (cid:107) v εh − v ε (cid:107) L (Ω) ≤ Ch (cid:107)∇ v ε (cid:107) L (Ω) . Thus, we obtain (cid:107) v P h − v (cid:107) L (Ω) ≤ (cid:107) v P h − v εh (cid:107) L (Ω) + (cid:107) v ε − v (cid:107) L (Ω) + Ch (cid:107)∇ v ε (cid:107) L (Ω) , and we deduce the result by taking the limit ε →
0. The estimates (20) and (21) are proven ina similar manner, replacing Ω by the corresponding space product. (cid:50)
Here, and in what follows, C denotes a positive constant, independent of τ , h and otherapproximation parameters, which may change of value among different expressions.The following lemma states that, like the solution u of problem (1)-(2), the sequence ofapproximating functions u τh remains uniformly bounded in L ∞ ( Q T ). The result is a consequenceof the Lipschitz continuity of A and the boundedness of u , which imply the sublinearity of A in[0 , (cid:107) u (cid:107) L ∞ ( Q T ) ]. Lemma 2.
There exists a constant
C > , independent of τ and h , such that (cid:107) u τh (cid:107) L ∞ ( Q T ) ≤ C. (22) Proof.
Let u be the unique solution of problem (1)-(2). Since it satisfies the bound (7), wededuce from (H) the existence of a positive constant C A such that | A ( s ) | ≤ C A s for s ∈ [0 , C M ].We may assume that this property is satisfied for all s ∈ R since, on the contrary, redefining A in s ∈ R \ [0 , C M ] as A ( s ) = C A s , the solution to problem (1)-(2) is still the same function thanfor the original range kernel A . Hence, from now on we assume, without loss of generality, | A ( s ) | ≤ C A | s | for all s ∈ R . Consider the change of unknown w n = e − µt n u n , for some µ > w n satisfies w [ j ] = u [ j ] ,w n +1 [ j ] = e − µτ w n [ j ] + τ e − µt n +1 (cid:0) A P h ( e µt n w τh ( t n , · ))[ j ] + f P n j ( e µt n w n [ j ]) (cid:1) , for j ∈ J . We have |A P h ( e µt n w τh ( t n , · ))[ j ] | ≤ e µt n C A max k ∈ J | w n [ k ] | h d (cid:88) k ∈ J J P h [ j , k ] , | f P n j ( e µt n w n [ j ]) | ≤ e µt n C f (1 + max k ∈ J | w n [ k ] | ) , implying | w n +1 [ j ] | ≤ e − µτ (cid:16) τ C f + (cid:0) τ (2 C A + C f ) (cid:1) max k ∈ J w n [ k ] (cid:17) . Therefore, taking µ > C A + C f and returning to the original unknown this differences inequalityyields the uniform estimate max k ∈ J | u n [ k ] | ≤ C , with C depending only on max k ∈ J | u [ k ] | , T and the constants C A and C f . (cid:50)
7. Galiano
Remark 1.
Since A ( s ) and f ( · , · , s ) are Lipschitz continuous in the interval I = (0 , (cid:107) u (cid:107) L ∞ ( Q T ) ) a sufficient condition for stability of the numerical discretization (13) is given by τ < (cid:107) A (cid:107) W , ∞ ( I ) + (cid:107) f (cid:107) L ∞ ( Q T ) × W , ∞ ( I ) , (23) which is similar to that found by P´erez-Llano and Rossi [24, Proposition 2.9] when A ( s ) = | s | p − s and f = 0 . In fact, (23) ensures that the comparison principle is satisfied in the discrete settingof that problem. Notice that, in contrast to the corresponding local diffusion problem, (23) doesnot depend on the spatial mesh size.The proof of Lemma 2 may be easily adapted to the discretizations introduced in Sections 3.2and 3.3. We, therefore, omit the corresponding proofs. Remaining estimates.
Estimate (19) applied to the initial datum yields I ≤ Ch. (24)By direct computation we get, for ( t, x ) ∈ [ t n , t n +1 ) × Ω j , u τh ( t, x ) − u τh ( t, x ) = ( t − t n ) u n +1 [ j ] − u n [ j ] τ = ( t − t n ) (cid:0) A P h ( u τh ( t, · ))( x ) + f P τh ( t, x , u τh ( t, x )) (cid:1) . Using (22) and that | t − t n | < τ for t ∈ [ t n , t n +1 ), we deduce I ≤ Cτ. (25)For ( s, x ) ∈ [ t n , t n +1 ) × Ω j , we have A ( u τh ( s, · ))( x ) − A P h ( u τh ( s, · ))( x )= (cid:88) k ∈ J A ( u n [ k ] − u n [ j ]) (cid:16) (cid:90) Ω k J ( x , y ) d y − h d J P h [ j , k ] (cid:17) = (cid:88) k ∈ J A ( u n [ k ] − u n [ j ]) (cid:16) (cid:90) Ω k (cid:0) J ( x , y ) − J P h ( x , y ) (cid:1) d y (cid:17) . From (20) and (22), we obtain I ≤ C (cid:90) Ω (cid:90) Ω (cid:12)(cid:12) J ( x , y ) d y − J P h ( x , y ) (cid:12)(cid:12) d y d x ≤ Ch , (26)and, in view of (21), we deduce I ≤ Cτ h. (27)Using estimates (15), (24), (25), (26) and (27) in (14), we finally obtain (cid:90) Ω | u ( t, x ) − u τh ( t, x ) | d x ≤ C ( τ + h ( τ + h ))+ C (cid:90) Q t | u ( s, x ) − u τh ( s, x ) | d x ds, and, for τ, h <
1, Gronwall’s lemma implies (cid:107) u − u τh (cid:107) L ∞ (0 ,T ; L (Ω)) ≤ C ( τ + h ) . In this section, we assume J ∈ L ∞ (Ω × Ω).
Definition of the approximation.
We set G = R ( Rearrangement ). The definitions of u R h and f R n j are the same as (16) and (17). For those of A R h and J R h we need to introduce someadditional tools.Let L Q ∈ N and consider the uniform mesh, { q i } L Q i =0 , of an interval [ C m , C M ] ⊃ [inf( u ) , sup( u )],where u is the solution of problem (1)-(2). We choose q i decreasingly ordered, that is, given by q = C M , q i = q − iq, for i = 1 , . . . , L Q , with q = C M − C m L Q . (28)Similarly, let L R ∈ N and consider the uniform mesh, { r i } L R i =0 , of the interval [inf( J ) , sup( J )].We also take { r i } Ri =0 decreasingly ordered, r = sup( J ) , r j = r − jr, for j = 1 , . . . , L R , with r = sup( J ) − inf( J ) L R . We introduce the decomposition of a real number, ρ , in its quantized and remainder partsrelative to a mesh { s i } L S i =0 , for S = Q, R , given by ρ = Q S ( ρ ) + R S ( ρ ), where Q S ( ρ ) = s m , with m = argmin ≤ i ≤ S | ρ − s i | , and |R S ( ρ ) | ≤ s, (29)being s the size of the mesh.Consider the functional A R ( v )( x ) = (cid:90) Ω Q R ( J ( x , y )) A ( Q Q ( v ( y )) − Q Q ( v ( x ))) d y . We have A ( v )( x ) = A R ( v )( x ) + ε ( x ) , with (30) ε ( x ) = (cid:90) Ω Q R ( J ( x , y )) (cid:16) A ( v ( y ) − v ( x )) − A ( Q Q ( v ( y )) − Q Q ( v ( x ))) (cid:17) d y + (cid:90) Ω (cid:0) J ( x , y ) − Q R ( J ( x , y )) (cid:1) A ( v ( y ) − v ( x )) d y . Since J ∈ L (Ω × Ω) and A is Lipschitz continuous, if v ∈ L ∞ (Ω) we deduce | ε ( x ) | ≤ C ( q + r ) . (31)The quantized functions are simple finite functions, which may be expressed as Q Q ( v ( x )) = L Q (cid:88) i =0 q i X U i ( x ) , Q R ( J ( x , y )) = L R (cid:88) m =0 r m X V m ( x , y ) , (32)where U i and V m are the level sets of Q Q ◦ v and Q R ◦ J . Then, as shown in [15, Theorem 1],we have that if x ∈ U k , for some k = 0 , . . . , Q , then A R ( v )( x ) = L Q (cid:88) i =0 A ( q i − q k ) L R (cid:88) m =0 r m µ im ( x ) , (33)9. Galianowhere µ im ( x ) = meas { y ∈ U i : Q R ( J ( x , y )) = r m } .If v ≡ v h is a piecewise constant function defined as in (16), then the level sets U i introducedin (32) are unions of hypercubes Ω j . Inspired in the expression (33), we define the discretenonlocal diffusion operator as, for j ∈ J such that Ω j ⊂ U k , A R h ( v h )[ j ] = L Q (cid:88) i =0 A ( q i − q k ) L R (cid:88) m =0 r m µ im [ j ] , with µ im [ j ] = 1 h d (cid:90) Ω j µ im ( x ) d x . (34)Finally, J R h is defined like (18), with J replaced by Q R ◦ J . Remaining estimates.
The estimates for I , I and I are the same as in (24), (25) and (27).Thus, it only rests to estimate I . Using (30) and (31), we get |A ( u τh ( s, · ))( x ) − A R h ( u τh ( s, · ))( x ) | ≤ |A R ( u τh ( s, · ))( x ) − A R h ( u τh ( s, · ))( x ) | + C ( q + r ) . (35)Using (33) we get, for ( s, x ) ∈ [ t n , t n +1 ) × Ω j such that Ω j ⊂ U nk = { y ∈ Ω : u n ( y ) = q k } , A R ( u τh ( s, · ))( x ) − A R h ( u τh ( s, · ))( x ) = L Q (cid:88) i =0 A ( q i − q k ) L R (cid:88) m =0 r m ( µ i,nm ( x ) − µ i,nm [ j ]) , (36)where µ i,nm is defined like µ im with U i replaced by U ni , see (34). We have L R (cid:88) m =0 r m µ i,nm [ j ] = 1 h d L R (cid:88) m =0 r m (cid:90) Ω j µ i,nm ( x ) d x = 1 h d (cid:90) Ω j M ( x ) d x , with M ( x ) = (cid:80) L R m =0 r m µ i,nm ( x ). Therefore, L R (cid:88) m =0 r m ( µ i,nm ( x ) − µ i,nm [ j ]) = M ( x ) − h d (cid:90) Ω j M ( x ) d x . Let J ni = { k ∈ J : Ω k ⊂ U ni } so that U ni = ∪ k ∈ J ni Ω k . Then M ( x ) = L R (cid:88) m =0 r m (cid:88) k ∈ J ni meas { y ∈ Ω k : Q R ( J ( x , y )) = r m } = (cid:88) k ∈ J ni (cid:90) Ω k L R (cid:88) m =0 r m { y ∈ Ω k : Q R ( J ( x , y ))= r m } ( y ) d y = (cid:88) k ∈ J ni (cid:90) Ω k Q R ( J ( x , y )) d y = (cid:90) U ni Q R ( J ( x , y )) d y . Since J ∈ BV (Ω × Ω), we deduce that M ∈ BV (Ω) and, therefore, from (36), the boundednessof A and the estimate (19) applied to M we deduce |A R ( u τh ( s, · ))( x ) − A R h ( u τh ( s, · ))( x ) | ≤ Ch, I ≤ C ( h + q + r ) . (37)We conclude like in the previous section: using estimates (15), (24), (25), (27) and (37) in (14),we obtain, for τ, h < (cid:90) Ω | u ( t, x ) − u τh ( t, x ) | d x ≤ C ( τ + h + q + r )+ C (cid:90) Q t | u ( s, x ) − u τh ( s, x ) | d x ds. Gronwall’s lemma implies (cid:107) u − u τh (cid:107) L ∞ (0 ,T ; L (Ω)) ≤ C ( τ + h + q + r ) . For notational simplicity, we relocate the set Ω = (0 , L ) d to Ω F = ( − L/ , L/ d , where weintroduce the mesh M h (Ω F ) = { x j } j ∈ J F similar to (10), but with J F = (cid:26) j = ( j , . . . , j d ) , j i = − J Ω , . . . , J Ω , i = 1 , . . . , d (cid:27) , (38)with J Ω even, for simplicity. In this section, we assume that the spatial kernel is of the convo-lution form J ( x , y ) = w ( x − y ), and replace Assumption ( H ) by:( H ) F w ∈ BV (Ω F ) is even and non-negative. In addition, either (i) w is periodic in Ω F and w ∈ L ∞ (Ω F ), or (ii) w is compactly supported in Ω F .Observe that to have w ( x − y ) well defined for x , y ∈ Ω F , we need to be able to evaluatethis function in B = ( − L, L ) d . However, in ( H ) F , w is only defined in Ω F . We remove thisobstacle by considering, according to the assumptions in ( H ) F , either (i) its periodic extensionto B , or (ii) its extension by zero to B . Correspondingly, we also extend the mesh and set M h ( B ) = { x j } j ∈ J B similar to (10), but with J B = { j = ( j , . . . , j d ) , j i = − J Ω , . . . , J Ω , i = 1 , . . . , d } . Thus, if x ∈ Ω j and y ∈ Ω k , with j , k ∈ J F , we have x − y ∈ Ω j − k , with j − k ∈ J B . Definition of the approximation.
We set G = F ( Fourier ). Consider the functional A F ( v )( x ) = (cid:90) Ω F w ( x − y ) A ( v ( y ) − Q Q ( v ( x ))) d y , with Q Q introduced in (29) for the uniform mesh of the range of u given in (28). We have A v ( x ) = A F ( v )( x ) + ε ( x ) , with (39) ε ( x ) = (cid:90) Ω F w ( x − y ) (cid:16) A ( v ( y ) − v ( x )) − A ( v ( y ) − Q Q ( v ( x ))) (cid:17) d y .
11. GalianoSince w ∈ L (Ω F ) and A is Lipschitz continuous, if v ∈ L ∞ (Ω F ) we deduce | ε ( x ) | ≤ Cq. (40)We also introduce the functional A F ,i ( v )( x ) = (cid:90) Ω F w ( x − y ) A ( v ( y ) − q i ) d y if Q Q ( v ( x )) = q i , for some i = 0 , . . . , Q , which may be rewritten as the convolution A F ,i ( v )( x ) = w ∗ H i ( v )( x ) , with H i ( v )( y ) = A ( v ( y ) − q i ) . (41)If v and w are periodic functions in Ω F then A F ,i may be computed using the Fourier transform,this is, A F ,i ( v ) = F − ( F ( w ) F ( H i ( v ))) . (42)The natural discretization of (42) is given in terms of the discrete Fourier transform, F d . Toemploy this tool, the discrete sequence generated by the scheme (13) must be periodic, propertythat we obtain if the piecewise constant approximations of u , f and w are periodic. Since thespatial domain is an hyper-cube, we may do this by averaging the function values correspondingto opposed faces of the cube. For instance, if Ω ⊂ R , we split Ω F in three disjoint components:the interior Ω h FI = [ − L/ h, L/ − h ] × [ − L/ h, L/ − h ], the corners Ω h FC = (cid:26) L z , z ) : z j ∈ {− , } for j = 1 , (cid:27) , and the sides without the corners , Ω h FS = Ω F \ (Ω h FI ∪ Ω h FC ). Then, we define the periodic constant-wise approximation of u as u F h ( x ) = h (cid:90) Ω j u ( x ) d x if x ∈ Ω j ∩ Ω h FI for some j ∈ J F , h (cid:90) Ω j ∪ Ω j s u ( x ) d x if x ∈ Ω j ∩ Ω h FS for some j ∈ J F , h (cid:90) Ω h FC u ( x ) d x if x ∈ Ω h FC , where j s is the node opposed to j . This is, if j = ( ± , j ) then j s = ( ∓ , j ), and if j = ( j , ± j s = ( j , ∓ j , j (cid:54) = − ,
1. We define in a similar way the periodic approximations of f F τh and w F h . The extension to higher dimensional domains is straightforward.For the nonlocal diffusion operator, we define, for j ∈ J F , A F h ( v h )[ j ] = A F ,ih ( v h )[ j ] if Q Q ( v h [ j ]) = q i , (43)with A F ,ih ( v h )[ j ] = F − d ( F d ( w F h ) F d ( H i ( v h )))[ j ] . Due to the periodicity imposed on the piecewise approximations of the data, we may use thecircular convolution theorem to obtain F − d ( F d ( w F h ) F d ( H i ( v h )))[ j ] = w F h (cid:126) H i ( v h )[ j ] , u , w and f are periodic functions in Ω F , see (42). However, notice that suchperiodicity has not been assumed and that the identity (42) is, in general, not valid. Remaining estimates.
The estimate for I is the same than (27). Thus, it only rests toestimate I and I . We have I = (cid:90) Ω h FI | u − u F h | + (cid:90) Ω \ Ω h FI | u − u F h | . Since the restriction of a BV function to a smooth set is still a BV function, we may use (19)to estimate the first term of the right hand side by Ch . For the second (the border ), we have (cid:90) Ω F \ Ω h FI | u − u F h | ≤ C (cid:107) u (cid:107) L ∞ (Ω F ) Per(Ω F ) h, (44)implying that I ≤ Ch. (45)Using (39) and (40), we get |A ( u τh ( s, · ))( x ) − A F h ( u τh ( s, · ))( x ) | (46) ≤ |A F ( u τh ( s, · ))( x ) − A F h ( u τh ( s, · ))( x ) | + Cq.
Since u τh is periodic, we may use the circular convolution theorem to express A F ,ih as a convo-lution. The continuous counterpart, A F ,i , may be also expressed in convolution form, see (41).Therefore, the estimation of (46) is reduced to estimate the difference w ∗ H i ( u n )( x ) − w F h (cid:126) H i ( u n )( x ) . (47)In the interior , i.e., if x ∈ Ω j ∩ Ω h FI for some j ∈ J F , we have, if Q Q ( u n [ j ]) = q i , w ∗ H i ( u n )( x ) − w F h (cid:126) H i ( u n )( x )= (cid:88) k ∈ J F A ( u n [ k ] − q i ) (cid:16) (cid:90) Ω k w ( x − y ) d y − (cid:90) Ω j − k w ( z ) d z (cid:17) . Using the change of variable z = ˆ x j − y , with ˆ x j = ( j − / h , we get (cid:12)(cid:12)(cid:12) (cid:90) Ω k w ( x − y ) d y − (cid:90) Ω j − k w ( z ) d z (cid:12)(cid:12)(cid:12) ≤ (cid:90) Ω k (cid:12)(cid:12) w ( x − y ) − w (ˆ x j − y ) (cid:12)(cid:12) d y . If w ∈ C ( B ), using that | x − ˆ x j | ≤ Ch , we deduce (cid:88) k ∈ J F (cid:12)(cid:12)(cid:12) (cid:90) Ω k w ( x − y ) d y − (cid:90) Ω j − k w ( z ) d z (cid:12)(cid:12)(cid:12) ≤ Ch (cid:107)∇ w (cid:107) L ( B ) . (48)By density, we obtain the same result for w ∈ BV ( B ), but with (cid:107)∇ w (cid:107) L ( B ) replaced by (cid:107) w (cid:107) BV ( B ) . 13. Galiano Initial data x u Initial data profile x u Initial data profile
Figure 1:
Initial datum and its profiles at x = 0 .
45 and x = 0 . To get a bound of (47) in the border , we use the assumption ( H ) F . If w ∈ L ∞ (Ω) then wereason like in the deduction of (44). Otherwise, if w is compactly supported in Ω then, for h small enough, we have w F h ( x ) = w ( x ) = 0 for x ∈ Ω h FS ∪ Ω h FC . From these observations andrecalling that u τh is bounded in L ∞ ( Q T ), we deduce from (46) and (48) I ≤ C ( h + q ) . (49)We conclude like in the previous sections: using estimates (15), (27), (25), (45) and (49) in(14), we obtain, for τ, h < (cid:90) Ω | u ( t, x ) − u τh ( t, x ) | d x ≤ C ( τ + h + q + r )+ C (cid:90) t (cid:90) Ω | u ( s, x ) − u τh ( s, x ) | d x ds, and Gronwall’s lemma implies, (cid:107) u − u τh (cid:107) L ∞ (0 ,T ; L (Ω)) ≤ C ( τ + h + q ) . Remark 2.
For periodic data, i.e. when the identity (42) is satisfied, we may obtain strongerestimates of (47) by using the results of Epstein [12], which establish error bounds in the ap-proximation of F by F d . In particular, his results are valid for piecewise continuous periodicfunctions. For the next experiments, we always fix the spatial domain as Ω = (0 , × (0 ,
1) and the finaltime as T = 1. The initial datum is given by u ( x , x ) = (cid:88) j =1 2 (cid:89) i =1 J ( x i ; µ ji , σ ji ) , (50)where J ( x ; µ, σ ) is a Gaussian function of mean µ and variance σ , with unitary norm in L ( R ),and µ = 0 . , σ = 0 . , µ = µ = 0 . , σ = 0 . , µ = 0 . , σ = σ = 0 . J Ω ∈ { , , } . The time discretization is taken uniform too, with timestep τ ∈ { . , . , . } . Although in some cases these parameters do not satisfy the sufficientstability condition (23) we always deal with cases in which the schemes are actually convergent.For the spatial kernel we use either a L unitary box spatial kernel of radius r , S ( x ) = 14 r B r ( x ) , J ( x , y ) = S ( x − y ) , (51)where B r ( x ) = ( x − r, x + r ) × ( x − r, x + r ), with r ∈ { . , . , . , . , . } , or a centeredand L unitary truncated Gaussian kernel, this is, G ( x ) = c exp (cid:16) − | x | σ (cid:17) B r ( x ) , J ( x , y ) = G ( x − y ) (52)where c is a normalizing constant, and the radius of B r ( x ) is fixed as r = 6 σ , for the set of data σ ∈ { . , . , . , . , . } .For the range kernel we take either the identity A ( s ) = s or A ( s ) = | s | p − s , for p = 2 (whichis actually the identity) or p = 3. The latter is the range kernel associated to the nonlocal p –Laplacian.The quantization of the range used in the functional rearrangements and the Fourier basedmethods is fixed, respectively, as follows. For the former, we take L Q = 500 levels of quantization,from min j u n [ j ] to max j u n [ j ], where the n –dependence indicates that we requantize in each timestep. For the latter, we take a smaller number of levels: 10 in Experiments 1 and 2, and 150 inExperiment 3, without requantization, that is, as a uniform quantization between min j u [ j ] andmax j u [ j ]. In the Fourier transform based algorithm, after each time step the approximation isdefined as a linear interpolation in terms of the quantized levels: u n +1 [ j ] = 1 q (cid:16) ( q i +1 − u n [ j ]) A F ,ih ( u n )[ j ] + ( u n [ j ] − q i ) A F ,i +1 h ( u n )[ j ] (cid:17) , see (28) and (43).We perform the following experiments, where the notation Ptw , RR , and FFT is used torefer to the pointwise, the functional rearrangements, and the Fourier transform based methods,respectively:1. Box spatial kernel and linear range kernel. This is, computationally, the simplest situation.Moreover, an explicit solution may be calculated. Thus, in addition to a time executioncomparison of the algorithms we may compute their accuracy and rate of convergence.2. Box spatial kernel and power-like range kernel. We implement the p –Laplacian rangekernel A ( s ) = | s | p − s , but actually taking p = 2 so that the problem we solve is the sameas in Experiment 1. However, it is computationally different since the evaluation of therange kernel A ( u ( y ) − u ( x )) is performed via an user defined function, instead of directlyimplemented as u ( y ) − u ( x ). This is clearly more time consuming.3. Gaussian spatial kernel and power-like range kernel. We take the Gaussian space kernel(51) and the range kernel A ( s ) = | s | s . An explicit solution is not at hand. In additionto the comparison among execution times, we compare the approximations to the finest15. Galiano( J Ω = 300, τ = 0 . Ptw method. In this experimentwe introduced a modification of the
FFT algorithm consisting on a zero padding of thesamples to the next power of two. As it is well known [13], the computation of the fastFourier transform is specially efficient in this case.The experiments have been performed in a standard laptop with an i7 processor. The codes arewritten in C ++ and run on a single core. The execution time is measured in seconds, using the clock() library. Although this gives a rough measure, its accuracy is enough for our purposes,given the large differences among the measured execution times of the different methods anddata sets. The codes are not specially optimized from the programming point of view, butwritten literally as described in the previous sections. However, notice that the structure of the FFT method is specially suitable for parallelization since the slides, A F ,ih ( u n ), that define theapproximation may be computed independently, in separated threads. In this case, we may calculate an explicit solution. Let r ∈ (0 , /
4) and consider the spatial kernel(52). Let u ( x ) be compactly supported in ( r, − r ) × ( r, − r ). Then u ( x ) (cid:82) Ω J ( x − y ) d y = u ( x ),and, for i, j = 1 , i (cid:54) = j , the convolution term reduces to I ( x ) = r (cid:82) B r ( x ) u ( y ) d y , whichmay be calculated exactly for our choice of the initial datum, see (50). Fixing the reaction termas f ( t, x , s ) = e − λt (cid:0) (1 − λ ) u ( x ) − I ( x ) (cid:1) , for some λ ∈ R ( λ = 0 . u ( t, x ) = e − λt u ( x )is the solution of problem (1)-(2). Observe that although the initial datum defined in (50) isnot compactly supported in Ω, its values close to the boundary are so small that its discretizedversion is actually compactly supported in Ω.In Figures 2 and 3 we show the execution times and the relative errors, respectively, asfunctions of the spatial kernel radius, corresponding to each method and the correspondingdiscretization parameters. As we shall check in Experiment 2, the linearity of the range kernelis crucial for the Ptw method to provide lower execution times than the other methods whenthe radius of the spatial kernel and the number of nodes are small. However, for large radius the
FFT method has always the best figures, and the RR method outperforms the Ptw method. Itis also interesting to notice the similarity between the slopes of the execution time versus theradius corresponding to the
FFT and RR methods, in contrast to that of the Ptw method.With respect to the relative error, computed as RE ( u τh , u ) = (cid:88) j ∈ J Ω | u τh ( T, x j ) − u ( T, x j ) | (cid:46) (cid:88) j ∈ J Ω | u ( T, x j ) | , the Ptw discretization is, in general, the most accurate. In any case, the other methods alsogive good approximations, having the RR method a threshold accuracy which depends on thenumber of quantization levels fixed for the range kernel. This is specially noticeable for smalltime steps.In this example, and whatever the function u we choose, we observe that, for a given timestep, the error is not specially affected by the reduction of the spatial mesh size, suggesting that16rror analysis of some nonlocal diffusion discretization schemes J = 100, = 0.1 PtwRRFFT J = 100, = 0.01 J = 100, = 0.001 J = 200, = 0.1 J = 200, = 0.01 J = 200, = 0.001 J = 300, = 0.1 J = 300, = 0.01 J = 300, = 0.001 Figure 2:
Experiment 1: Execution times (seconds, log scale) versus radius. Legends are the same forall the plots, see the top-left plot. the nonlocal equation behaves like an ordinary differential equation. Indeed, writing the firstiterations of the time semi-discrete Ptw method we get u = u + τ (cid:0) I − u + e − λτ ((1 − λ ) u − I ) (cid:1) = u + τ (cid:0)(cid:0) (1 − λ )(1 − λτ ) − (cid:1) u + (cid:0) − (1 − λτ ) (cid:1) I (cid:1) + o ( τ )= (1 − λτ ) u + o ( τ ) , where I j ( x ) = r (cid:82) B r ( x ) u j ( y ) d y . Similarly, u = u + τ (cid:0) I − u + e − λτ ((1 − λ ) u − I ) (cid:1) = (1 − λτ ) u + o ( τ ) , so that the resulting sequence is an explicit Euler approximation of the differential equation ∂ t u ( t, x ) = λu ( t, x ), for each x ∈ Ω. Thus, the only source of error due to the spatial dis-cretization (excluding the initial datum discretization), the integrals I j , only contribute to theapproximation error as o ( τ ) and pass, in practice, unnoticed in the numerical experiments.Notice that this property is particular to the choice of the linear range kernel and the sourceterm of this example. It is not a general behavior of the problem. In this experiment we study the increase in execution time due to function evaluation in therange kernel. We modify Experiment 1 by replacing the identity range kernel by A ( s ) = | s | p − s ,for p ≥
1, which defines the p –Laplacian nonlocal evolution equation. We, actually, take p = 2,17. Galiano J = 100, = 0.1 PtwRRFFT J = 100, = 0.01 J = 100, = 0.001 Figure 3:
Experiment 1: Relative errors (log scale) versus radius. We show the case J Ω = 100 nodes.The cases J Ω = 200 ,
300 practically reproduce the same figures. so that the problem is analytically the same as in Experiment 1. However, the codes are designedto compute the general nonlinear case, which includes the evaluation of the function A ( s ), andare therefore slower than the codes of Example 1 that implement this term directly as theidentity, without any call to a user defined function.In Figure 4 we show the execution times corresponding to the different methods and dis-cretization parameters. The corresponding relative errors are the same as in Experiment 1.The behaviors of the RR and the FFT methods is very similar to those of Experiment 1.However, the
Ptw method suffers from the intensive function evaluation and performs up to threeorders of magnitude slower than its counterparts. The reason is that for the RR method thefunction evaluations are allocated in a fixed matrix given by A ( q i − q j ) for i, j = 1 , . . . ,
500 in ourexample, see (34), while for the
FFT method only few slices of the range kernel A F ,ih ( v h )[ j ] = F − d ( F d ( w F h ) F d ( H i ( v h )))[ j ], for i = 1 , . . . ,
10, are computed, see (43). On the contrary, thecomputation of
Ptw in each time step involves two nested loops, one covering all the nodes ofthe domain, j , and another covering the support of the spatial kernel (a box of radius r ), i ,where the evaluation of A ( u [ i ] − u [ j ]) must be accomplished. Therefore, for increasing radiussize the computation becomes much more expensive. In this experiment we take both the spatial and the range kernels as nonlinear functions definedby the Gaussian (52) and the p –Laplacian kernel A ( s ) = | s | p − s , with p = 3. In addition, wedeal with the purely diffusive case, i.e. we set f ≡
0. The stability of the numerical schemes,ensured when condition (23) is satisfied, fails in this case for the time step τ = 0 .
1. Thus, welimit this experiment to the time steps τ = 0 .
01 and τ = 0 . FFT method exhibits a somehowerratic behavior in which the execution time decreases or stays nearly flat when the radius ofthe spatial kernel is increased. See, e.g., the Experiments 1 and 2 with J Ω = 200. As it is wellknown, the efficiency of the fast Fourier transform is specially high when the size of the sampleto be transformed is a power of two or, more in general, of the form [13] 2 a a a a a a ,for integer numbers a , . . . , a . Thus, we incorporated to the set of methods a variant of the FFT method (
FFTO method) consisting on a slight optimization of the fft by oversampling byzero to the next power of two. This optimization is effective when the resulting padded samplesare not much bigger than the original ones.As we may check in Figure 5, only this optimization is capable of outperforming the RR method,specially for large spatial kernel radius and for the finer spatial mesh. Let us mention here that18rror analysis of some nonlocal diffusion discretization schemes J = 100, = 0.1 PtwRRFFT J = 100, = 0.01 J = 100, = 0.001 J = 200, = 0.1 J = 200, = 0.01 J = 200, = 0.001 J = 300, = 0.1 J = 300, = 0.01 J = 300, = 0.001 Figure 4:
Experiment 2: Execution times (seconds, log scale) versus radius. Legends are the same forall the plots, see the top-left plot. the small number (10) of slides (quantized levels) used in the Experiments 1 and 2 for the Fouriertransform based method must be replaced by a larger number (150) in order to keep the errorin similar bounds than the other methods.With respect to the RR method, the main change in this experiment is the computation ofthe discretization of the relative rearrangement J ∗ u given by, see (34), L R (cid:88) m =0 r m µ im [ j ] , (53)where r m are the discrete level lines of the quantized version of J , Q R ( J ), and µ im [ j ] is themeasure of the set of nodes { k : Q Q ( u n [ k ]) = q i , and Q R ( J [ j , k ]) = r m } .In the previous experiments, J has only the two level lines { , / (4 r ) } , making the compu-tation of (53) straightforward. In fact, fast algorithms for this task are easy to implement, see[25]. However, if J is a continuous non-constant function, a choice of the number of discretelevel lines must be done. Just like for the unknown, u , itself.In this experiment we fixed, like in the previous, L Q = 500 levels of discretization for u n (i.e., with requantization en each time step). The number L R of discrete levels of J is fixed asthe maximum number of levels given by the precision (double) used in the codes. This is, thenumber of levels of G ( x k ) for x k ∈ B r ( ). This number is determined by the mesh size andby the radius r = 6 σ , where σ ∈ { . , . , . , . , . } . For instance, for J Ω = 100 and σ = 0 .
01 we get L R = 518, while for J Ω = 300 and σ = 0 .
05 we have L R = 4183. Of course,both L Q and L R are arbitrarily chosen and may be used to find a compromise between accuracyand execution time. 19. Galiano J = 100, = 0.01 PtwRRFFTFFTO J = 200, = 0.01 J = 300, = 0.01 J = 100, = 0.001 J = 200, = 0.001 J = 300, = 0.001 Figure 5:
Experiment 3: Execution times (seconds, log scale) versus radius. Legends are the same forall the plots, see the top-left plot. The results of this experiment are certainly clear, see Figures 5 and 6 . With respect tothe execution time, the
Ptw method is much slower than the other methods, specially for largeradius sizes, and the RR method outperforms the FFT method up to an order of magnitudefor small radius. The only hope for the Fourier transform based method seems to be a furtheroptimization of the sample sizes introduced in the
FFTO implementation.With respect to the error, and taking into account that a exact solution is not known,we have computed the relative differences between the
Ptw approximation obtained for τ =0 . RR and the FFT -approximations (or the
FFTO , since they are very close to each other). In contrast with theprevious experiments, we see a slight improvement of the RR and FFT approximations fordecreasing spatial mesh. Observe that the only reason for improvement of the approximationsbecause the mesh refinement is due to a higher accuracy in the quadrature of the nonlocal term.However, since this quadrature is reduced to a square of size 4 r (cid:28) | Ω | , the improvement is hardto notice. Appendix: Functional rearrangements
Under suitable regularity assumptions on the functions v and g defined in Ω, the coarea formulastates (cid:90) Ω g ( y ) |∇ v ( y ) | d y = (cid:90) ∞−∞ (cid:90) v = s g ( y ) d Γ( y ) ds. (54)Taking g ( y ) = J ( x , y ) A ( v ( y ) − v ( x )) / |∇ v ( y ) | , and assuming v ( x ) ∈ [0 , Q ] for all x ∈ Ω we getfrom (9) and (54), A ( v )( x ) = (cid:90) Q A ( s − v ( x )) (cid:90) v = s J ( x , y ) |∇ v ( y ) | d Γ( y ) ds. (55)In this expression, the inner integral is computed relative to the level sets of v , but independentlyof the ordering or such sets. This is the main motivation to introduce a reformulation of (9) interms of the decreasing and the relative rearrangements.20rror analysis of some nonlocal diffusion discretization schemes J = 100, = 0.01 PtRRFFT J = 200, = 0.01 J = 300, = 0.01 J = 100, = 0.001 J = 200, = 0.001 J = 300, = 0.001 Figure 6:
Experiment 3: Relative differences with respect to the
Ptw approximation with τ = 0 . scale) versus radius. Legends are the same for all the plots, see the top-left plot. The decreasing and the relative rearrangements