A fractional model for anomalous diffusion with increased variability. Analysis, algorithms and applications to interface problems
AA FRACTIONAL MODEL FOR ANOMALOUS DIFFUSION WITH INCREASEDVARIABILITY. ANALYSIS, ALGORITHMS AND APPLICATIONS TO INTERFACEPROBLEMS
MARTA D’ELIA ∗ AND
CHRISTIAN GLUSA † Abstract.
Fractional equations have become the model of choice in several applications where heterogeneities at themicrostructure result in anomalous diffusive behavior at the macroscale. In this work we introduce a new fractional operatorcharacterized by a doubly-variable fractional order and possibly truncated interactions. Under certain conditions on the modelparameters and on the regularity of the fractional order we show that the corresponding Poisson problem is well-posed. Wealso introduce a finite element discretization and describe an efficient implementation of the finite-element matrix assembly inthe case of piecewise constant fractional order. Through several numerical tests, we illustrate the improved descriptive power ofthis new operator across media interfaces. Furthermore, we present one-dimensional and two-dimensional h -convergence resultsthat show that the variable-order model has the same convergence behavior as the constant-order model. Key words.
Variable-order fractional operators, anomalous diffusion, subsurface diffusion, interface problems
1. Introduction.
Nonlocal models are becoming a popular alternative to partial differential equations(PDEs) when the latter fail to capture effects such as multiscale and anomalous behavior. In fact, severalscientific and engineering applications exhibit hierarchical features that cannot be described by classicalmodels. As an example we mention applications in continuum mechanics [37, 39, 43], phase transitions[8, 13, 15], corrosion [50], turbulence [7, 44, 51], and geoscience [9, 40, 55, 54, 58]. In this work we considernonlocal operators of fractional type, i.e. integral operators characterized by singular and non-integrablekernels, that correspond to differential operators of fractional orders, as opposed to the integer order in thePDE case.Note that fractional differential operators are almost as old as their integer counterpart [34]; however,their usability has increased during the last decades thanks to progress in computational capabilities and toa better understanding of their descriptive power. The simplest form of a fractional operator is the fractionalLaplacian; in its integral form, its action on a scalar function u is defined as [34](1.1) ( − ∆) s u ( x ) = C n,s p . v . (cid:90) R n u ( x ) − u ( y ) | x − y | n +2 s d y , where s ∈ (0 ,
1) is the fractional order , n the spatial dimension, C n,s a constant and where p.v. indicates theprincipal value. It follows from (1.1) that in fractional modeling the state of a system at a point dependson the value of the state at any other point in the space; in other words, fractional models are nonlocal.Specifically, fractional operators are special instances of more general nonlocal operators [21, 22, 28, 44] ofthe following form [19]:(1.2) − L u ( x ) = (cid:90) B δ ( x ) ( u ( x ) γ ( x , y ) − u ( y ) γ ( y , x )) d y . Here, interactions are limited to a Euclidean ball B δ ( x ) of radius δ , often referred to as horizon or interactionradius. The kernel γ ( x , y ) is a modeling choice and determines regularity properties of the solution. Notethat for δ = ∞ and for the fractional-type kernel γ ( x , y ) = | x − y | − n − s the nonlocal operator in (1.2) isequivalent to the fractional Laplacian in (1.1). Also, it has been shown in [22] that for that choice of δ and γ solutions corresponding to the nonlocal operator (1.2) converge to the ones corresponding to the fractionaloperator (1.1) as δ → ∞ (see [21] for more convergence results and for a detailed classification of theseoperators and relationships between them).In recent years, with the purpose of increasing the descriptive power of fractional operators, new mod-els characterized by a variable fractional order have been introduced for both space- and time-fractionaldifferential operators [49, 6, 62] and several discretization methods have been designed [16, 60, 61, 64, 53].However, the analysis of variable-order models is still in its infancy, with [33] and [52] being perhaps the ∗ Computational Science and Analysis, Sandia National Laboratories, CA, USA, ([email protected]) † Center for Computing Research, Sandia National Laboratories, NM, USA1 a r X i v : . [ m a t h . NA ] J a n Fig. 1.1 . Interface-problem configuration. only relevant works that address theoretical questions such as well-posedness for space-fractional differentialoperators with variable order s = s ( x ). The improved descriptive power of variable-order operators has beendemonstrated in some works on parameter estimation [44, 45, 46, 63]; here, via machine-learning algorithmsor more standard optimization techniques, the authors estimate the variable fractional power for operatorsof the form(1.3) − L u ( x ) = C n,s p . v . (cid:90) B δ ( x ) u ( x ) − u ( y ) | x − y | n +2 s ( x ) d y , where δ is either infinite or finite.Models such as the one in (1.3) are convenient to describe anomalous diffusion in case of heterogeneousmaterials or media, where different regions may be characterized by different diffusion rates, i.e. differentvalues of s , see Figure 1.1 for two-dimensional illustrations where s ( x ) is a piece-wise constant functiondefined over the domain. In other words, these models can describe physical interfaces by simply tuningthe function s ( x ). Modeling nonlocal interfaces is a nontrivial task due to the unknown nature of thenonlocal interactions across the interfaces. Only a few works in the literature have addressed this problem;among them, we mention [3] for nonlocal models in mechanics and [14] for nonlocal diffusion models withintegrable kernels. In this work, we tackle the nonlocal interface problem for fractional-type operators usinga generalization of the nonlocal operator in (1.3). Specifically, we add variability to the fractional order andthe kernel function itself. Our main contributions are: • The introduction of a new variable-order fractional operator characterized by s = s ( x , y ) and δ ∈ (0 , ∞ ]. This choice enables modeling of a much broader set of interface behaviors and, for symmetric s , symmetrizes the kernel, making analysis and implementation a much easier task. In fact, notethat in (1.3) the kernel is no longer symmetric, requiring more sophisticated quadrature rules forthe integration (regardless of the discretization technique of choice). • The analysis of well-posedness of the Poisson problem associated with the new operator in the generalnonsymmetric case. • The design of an efficient finite-element matrix assembly technique for the practical case of piecewiseconstant definition of the fractional order that often appears in, e.g., subsurface applications.
Outline of the paper.
In Section 2 we first recall some state-of-the-art results; then, we introduce thenew variable-order fractional operator and the corresponding strong and weak forms of the nonsymmetricdiffusion problem. We also describe properties of the associated energy norm and space. In Section 3 wereport the main result of the paper that proves via Fredholm alternative the well-posedness of the diffusionproblem. In Section 4 we introduce a finite element discretization and propose efficient quadrature rulesfor the integration of the nonsymmetric kernel. In Section 5 we report several one- and two-dimensionalnumerical tests that illustrate both the improved variability of our model and the accuracy of the numericalscheme. Finally, in Section 6 we summarize our findings. . A new variable-order fractional model. In this section we first introduce the notation thatwill be used throughout the paper and recall state-of-the-art results on nonsymmetric kernels and fractionalkernels with variable order. Then, we introduce a new variable-order fractional kernel and discuss theproperties of the associated norm and energy space.
The purpose of this section is two-fold. First, we recall the formulation of asymmetric nonlocal diffusion problem with finite interaction radius following the theory presented in [28, 31];then, we recall the theory presented in [33] on variable-order fractional operators (with infinite interactionradius). Our new model, presented in the next subsection, is a combination of the two.Let Ω ⊂ R n be an open bounded domain. We define its interaction domain as the set of points outsideΩ that interact with points inside Ω, i.e.(2.1) Ω I := { y ∈ R n \ Ω : | x − y | ≤ δ, for x ∈ Ω } , where δ ∈ (0 , ∞ ], referred to as interaction radius, determines the extent of the nonlocal interactions. Also,let γ ( x , y ) : R n × R n → R + be a nonnegative, symmetric, kernel function, not necessarily integrable. Thediffusion operator introduced in [31] is defined as(2.2) − L u ( x ) = p . v . (cid:90) Ω ∪ Ω I ( u ( x ) − u ( y )) γ ( x , y ) d y . From now on we remove the explicit reference to the principal value when the operator is used within anequation. The strong form of a truncated, symmetric, diffusion problem is then given by: for f : Ω → R ,and g : Ω I → R , find u : Ω ∪ Ω I → R such that(2.3) (cid:90) Ω ∪ Ω I ( u ( x ) − u ( y )) γ ( x , y ) d y = f ( x ) x ∈ Ω u ( x ) = g ( x ) x ∈ Ω I . The well-posedness of problem (2.3) was studied in [28]. Extensions of this model to the nonsymmetric casehave been proposed in [32], for integrable kernels, and in [19] for non-integrable kernels by using the operator(1.2) introduced in the previous section. This operator is often referred to as nonsymmetric nonlocal diffusionoperator or, more appropriately, nonlocal convection-diffusion operator as it introduces a nonlocal convectionterm. Since the non-symmetry of our problem is not related to convection, but only to the presence of aninterface that affects the way a quantity diffuses, in this work, we do not employ (1.2); instead, we use theoperator defined in (2.2) with a nonsymmetric kernel γ .Among the state-of-the-art formulations of nonlocal diffusion problems with nonsymmetric kernels, wemention the one introduced in [33] for integrable and non-integrable kernels with infinite support. Herethe authors consider the nonlocal operator (2.2) where the kernel γ is a nonnegative measurable function,possibly nonsymmetric. Among the kernels studied in [33], we recall the following nonsymmetric variable-order fractional kernel, whose analysis is relevant for the theory presented in this paper:(2.4) γ ( x , y ) = φ ( x ) | x − y | n +2 s ( x ) , where the non-symmetry with respect to ( x , y ) is due to the terms s ( x ) and φ ( x ). The latter is a nor-malization functions with the same role as C n,s in the standard fractional Laplacian in (1.1) [52]. Notethat the choice of γ corresponds to the variable-order version of the standard fractional Laplacian operatorintroduced in (1.3). Furthermore, when s and φ are constant functions, such operator is a multiple of thestandard fractional Laplacian. The well-posedness of problem (2.3) for nonsymmetric kernels was analyzedin [33] in a very general setting and in [52] for kernels of the form (2.4). We the purpose of increasing the descriptive power of the operator L ,we consider a fractional-type kernel with added variability. We define the new variable-order nonsymmetricfractional-type kernel as follows: γ ( x , y ) = φ ( x , y ) | x − y | n +2 s ( x , y ) X | x − y |≤ δ , x , y ∈ Ω ∪ Ω I , (2.5) ith constants s, s , and φ, φ such that0 < s ≤ s ( x , y ) ≤ s < , ∀ x , y ∈ Ω ∪ Ω I , (2.6) 0 < φ ≤ φ ( x , y ) ≤ φ < ∞ ∀ x , y ∈ Ω ∪ Ω I . (2.7)Thus, the resulting nonlocal operator is given by the following expression(2.8) −L u ( x ) = p . v . (cid:90) (Ω ∪ Ω I ) ∩ B δ ( x ) ( u ( x ) − u ( y )) φ ( x , y ) | x − y | n +2 s ( x , y ) d y . Remark φ ( x , y ) does not have anormalizing purpose but depends on the type of phenomenon that we are targeting. As an example, φ could describe a material property, such as diffusivity or permeability in the subsurface; in this context, thedependence on x and y describes changes in material properties across material interfaces. Remark φ , the operator defined in (2.8) corresponds to a variable-order,tempered fractional Laplacian. Specifically, we let(2.9) φ ( x , y ) = exp {− λ | x − y |} , where λ > φ for δ ∈ (0 , ∞ ]. Remark s ( x , y ) = s ( y , x ) the kernel γ is symmetric. This case is much easier to deal withboth in terms of analysis and computation. In fact, in the symmetric case, the theory introduced in [28]can be readily applied and guarantees well-posedness of problem (2.3). Furthermore, numerical integrationof a symmetric function poses minor challenges compared with the nonsymmetric case, which requires theemployment of more sophisticated quadrature rule. As we show in our numerical tests, a symmetric s stillguarantees improved model variability across interfaces, compared to the single-variable model in (2.4).We define the spaces H (Ω ∪ Ω I ; γ ) := (cid:110) v : Ω ∪ Ω I → R : v ∈ L (Ω ∪ Ω I ) , ||| v ||| γ < ∞ (cid:111) , (cid:101) H sδ (Ω) := H Ω (Ω ∪ Ω I ; γ ) := { v ∈ H (Ω ∪ Ω I ; γ ) : v | Ω I = 0 } , where the semi-norm ||| v ||| γ is defined as ||| v ||| γ := (cid:90) (cid:90) (Ω ∪ Ω I ) ( v ( x ) − v ( y )) γ s ( x , y ) d x d y . Here, γ s is the symmetric part of the kernel. We recall that the symmetric and anti-symmetric parts of thekernel γ are defined as follows: γ s ( x , y ) = 12 (cid:18) φ ( x , y ) | x − y | n +2 s ( x , y ) + φ ( y , x ) | x − y | n +2 s ( y , x ) (cid:19) X | x − y |≤ δ , (2.10) γ a ( x , y ) = 12 (cid:18) φ ( x , y ) | x − y | n +2 s ( x , y ) − φ ( y , x ) | x − y | n +2 s ( y , x ) (cid:19) X | x − y |≤ δ . (2.11)The following bounds hold trivially:1 | x − y | d +2 s ≤ | x − y | d +2 s ( x , y ) ≤ | x − y | d +2 s if | x − y | ≤ , (2.12) | x − y | d +2 s ≤ | x − y | d +2 s ( x , y ) ≤ | x − y | d +2 s if | x − y | ≥ , (2.13)In what follows, we will also frequently make use of the following lower bound on the horizon: δ := min { δ, } . (2.14)The following Lemma shows that the semi-norm above satisfies a Poincar´e inequality and, hence, equips (cid:101) H sδ (Ω) with a norm. Lemma
For all u ∈ (cid:101) H sδ (Ω) , the following Poincar´e inequality holds: C p ||| u ||| γ ≥ || u || L (Ω ∪ Ω I ) . Proof.
Assumptions (2.6) and inequality (2.12) imply ||| v ||| γ = (cid:90) (cid:90) (Ω ∪ Ω I ) , | x − y |≤ δ ( v ( x ) − v ( y )) γ s ( x , y ) d x d y + (cid:90) (cid:90) (Ω ∪ Ω I ) , | x − y |≥ δ ( v ( x ) − v ( y )) γ s ( x , y ) d x d y ≥ C (cid:90) (cid:90) (Ω ∪ Ω I ) , | x − y |≤ δ ( v ( x ) − v ( y )) | x − y | n +2 s d x d y . As a consequence, |||·||| γ is lower bounded by the semi-norm of a δ -truncated fractional Sobolev space withfractional order s , for which the Poincar´e inequality holds, see [28, Lemma 4.3]. The weak form.
Let g = 0 for simplicity; then, we can write the weak form of problem (2.3) as(2.15) (cid:90) Ω −L u v d x = (cid:90) (cid:90) (Ω ∪ Ω I ) ( u ( x ) − u ( y )) γ ( x , y ) v ( x ) d y d x = (cid:90) Ω f v d x , We denote the bilinear form in (2.15) by(2.16) A ( u, v ) = (cid:90) (cid:90) (Ω ∪ Ω I ) γ ( x , y ) v ( x )( u ( x ) − u ( y ) d y d x . While this form is formally equivalent to the one introduced in [33], the kernel (2.5) does not belong toclass of kernels studied in that paper. In the next section we show how to extend the well-posedness theorydeveloped in [33] and [52] to our new operator.For discretization purposes, it is convenient to further rewrite A as(2.17) A ( u, v ) = 12 (cid:90) (cid:90) (Ω ∪ Ω I ) ( u ( x ) − u ( y )) ( γ ( x , y ) v ( x ) − γ ( y , x ) v ( y )) d y d x . This form allows for simpler numerical integration and is used in our numerical implementation.
3. Well-posedness analysis.
In this section we prove the well-posedness of problem (2.15) by adaptingthe theory developed in the works of Felsinger, Kassmann and Voigt [33] and of Schilling and Wang [52].We first report the building blocks of our main well-posedness result.The following proposition establishes a relation between the symmetric and anti-symmetric part of thekernel and is a generalization of Proposition 3.1 in the paper by Schilling and Wang [52].
Proposition
Assume that β ( r ) := sup | x − y |≤ r | s ( x , y ) − s ( y , x ) | (3.1) atisfies (cid:90) dr ( β ( r ) | log r | ) r s < ∞ . (3.2) Moreover, assume that α ( r ) := sup | x − y |≤ r | φ ( x , y ) − φ ( y , x ) | is ( s + ε ) -H¨older continuous for ε > arbitrarily small. Then there exists a constant A γ such that sup x ∈ Ω ∪ Ω I (cid:90) { γ s ( x , y ) (cid:54) =0 } d y γ a ( x , y ) γ s ( x , y ) ≤ A γ < ∞ . (3.3) Proof.
First, by definition of symmetric and anti-symmetric parts of the kernel in (2.10), we have that | γ a ( x , y ) | ≤ γ s ( x , y ). This and assumptions (2.6), (2.7) and inequality (2.13) imply that there exists apositive constant C such that (cid:90) { γ s ( x , y ) (cid:54) =0 } , | x − y |≥ d y γ a ( x , y ) γ s ( x , y ) ≤ (cid:90) | x − y |≥ d y γ s ( x , y ) ≤ C (cid:90) ∞ r − − s =: C . Moreover, due to the definition of δ , (cid:90) { γ s ( x , y ) (cid:54) =0 } , | x − y |∈ [ δ, d y γ a ( x , y ) γ s ( x , y ) = 0so that we obtain (cid:90) { γ s ( x , y ) (cid:54) =0 } , | x − y |≥ δ d y γ a ( x , y ) γ s ( x , y ) ≤ C (3.4)Next, for | x − y | ≤ δ , we write: γ a ( x , y ) = 12 ( φ ( x , y ) − φ ( y , x )) | x − y | − n − s ( x , y ) + 12 φ ( y , x ) | x − y | − n (cid:16) | x − y | − s ( x , y ) − | x − y | − s ( y , x ) (cid:17) . =: γ a, ( x , y ) + γ a, ( x , y )Then, by (2.7) and the H¨older continuity of α we have that there exists a positive constant C such that (cid:90) { γ s ( x , y ) (cid:54) =0 } , | x − y |≤ δ d y γ a, ( x , y ) γ s ( x , y )= 12 (cid:90) { γ s ( x , y ) (cid:54) =0 } , | x − y |≤ δ d y ( φ ( x , y ) − φ ( y , x )) | x − y | − n − s ( x , y ) φ ( x , y ) | x − y | − n − s ( x , y ) + φ ( y , x ) | x − y | − n − s ( y , x ) ≤ C (cid:90) { γ s ( x , y ) (cid:54) =0 } , | x − y |≤ δ d y ( φ ( x , y ) − φ ( y , x )) | x − y | n +2 s ( x , y ) ≤ C (cid:90) { γ s ( x , y ) (cid:54) =0 } , | x − y |≤ δ d y α ( | x − y | ) | x − y | n +2 s ≤ C (cid:90) δ dr r s +2 ε r s =: C . (3.5) ere, we have used (2.12). On the other hand, since for a, b > | x − y | − a − | x − y | − b = − (cid:90) ba du | x − y | − u log | x − y | , then, for | x − y | ≤ δ ≤
1, we have that (cid:16) | x − y | − s ( x , y ) − | x − y | − s ( y , x ) (cid:17) ≤ (log | x − y | ) ( s ( x , y ) − s ( y , x )) | x − y | − { s ( x , y ) ,s ( y , x ) } . Hence, by (2.7) and assumption (3.2), there exists a positive constant C such that (cid:90) { γ s ( x , y ) (cid:54) =0 } , | x − y |≤ δ d y γ a, ( x , y ) γ s ( x , y )= 12 (cid:90) { γ s ( x , y ) (cid:54) =0 } , | x − y |≤ δ d y φ ( y , x ) | x − y | − n (cid:16) | x − y | − s ( x , y ) − | x − y | − s ( y , x ) (cid:17) φ ( x , y ) | x − y | − n − s ( x , y ) + φ ( y , x ) | x − y | − n − s ( y , x ) ≤ C (cid:90) { γ s ( x , y ) (cid:54) =0 } , | x − y |≤ δ d y (log | x − y | ) ( s ( x , y ) − s ( y , x )) | x − y | n +4 min { s ( x , y ) ,s ( y , x ) } φ ( x , y ) | x − y | − n − s ( x , y ) + φ ( y , x ) | x − y | − n − s ( y , x ) ≤ C (cid:90) { γ s ( x , y ) (cid:54) =0 } , | x − y |≤ δ d y (log | x − y | ) ( s ( x , y ) − s ( y , x )) | x − y | n +2 s ≤ C (cid:90) δ dr ( β ( r ) | log r | ) r s =: C . (3.6)Here, we have again used (2.12). Hence, from (3.4), (3.5) and (3.6), we obtain:sup x ∈ Ω ∪ Ω I (cid:90) Ω ∪ Ω I d y γ a ( x , y ) γ s ( x , y ) ≤ C + 2( C + C ) =: A γ . and the result follows.The next G˚arding inequality follows directly from [33, Lemma 3.1]; we report its proof for completeness. Lemma
There exist constants
C, c > such that A ( u, u ) ≥ C ||| u ||| γ − c || u || L (Ω ∪ Ω I ) ∀ u ∈ (cid:101) H sδ (Ω) . (3.7) Proof.
By using Young’s inequality for ε >
0, we obtain A ( u, u ) = (cid:90) (cid:90) (Ω ∪ Ω I ) γ ( x , y ) u ( x )( u ( x ) − u ( y )) d y d x ≥ (cid:90) (cid:90) (Ω ∪ Ω I ) γ s ( x , y ) | u ( x ) − u ( y ) | d y d x − (cid:90) (cid:90) (Ω ∪ Ω I ) | γ a ( x , y ) u ( x )( u ( x ) − u ( y )) | d y d x = 12 ||| u ||| γ − (cid:90) (cid:90) (Ω ∪ Ω I ) (cid:12)(cid:12)(cid:12) γ a ( x , y ) γ s ( x , y ) / γ s ( x , y ) − / u ( x )( u ( x ) − u ( y )) (cid:12)(cid:12)(cid:12) d y d x ≥ ||| u ||| γ − ε (cid:90) (cid:90) (Ω ∪ Ω I ) γ s ( x , y ) | ( u ( x ) − u ( y )) | d y d x − ε (cid:90) (cid:90) { γ s ( x , y ) (cid:54) =0 } γ a ( x , y ) γ s ( x , y ) u ( x ) d y d x ≥ (cid:18) − ε (cid:19) ||| u ||| γ − A γ ε || u || L (Ω ∪ Ω I ) . Hence, for ε small enough, we obtain A ( u, u ) ≥ C ||| u ||| γ − c || u || L (Ω ∪ Ω I ) . he Poincar´e inequality and Proposition 3.1 yield the continuity of the bilinear form A , as illustratedin the following lemma. Lemma
For all u, v ∈ (cid:101) H sδ (Ω) |A ( u, v ) | ≤ C ( C P , A γ ) ||| u ||| γ ||| v ||| γ , i.e., the bilinear form A ( · , · ) is continuous on (cid:101) H sδ (Ω) × (cid:101) H sδ (Ω) .Proof. By combining the Poincar´e inequality and Proposition 3.1, we have |A ( u, v ) | = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:90) (cid:90) (Ω ∪ Ω I ) ( γ s ( x , y ) + γ a ( x , y )) ( u ( x ) − u ( y )) v ( x ) d y d x (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ≤ ||| u ||| γ ||| v ||| γ + (cid:90) (cid:90) (Ω ∪ Ω I ) γ s ( x , y )( u ( x ) − u ( y )) d y d x / (cid:90) (cid:90) { γ s ( x , y ) (cid:54) =0 } γ a ( x , y ) γ s ( x , y ) v ( x ) d y d x / ≤ ||| u ||| γ ||| v ||| γ + A / γ ||| u ||| γ || v || L (Ω ∪ Ω I ) ≤ (1 + C P A / γ ) ||| u ||| γ ||| v ||| γ . Before proceeding with the well-posedness analysis, we prove a result for a perturbation of the weakproblem (2.15). First, let γ ( x , y ) := 1 | x − y | n +2 s X | x − y |≤ δ , and (cid:101) H sδ (Ω) := H Ω (Ω ∪ Ω I ; γ ) . Then, there exists some λ > γ s ( x , y ) ≥ λγ ( x , y ) , and hence for all u ∈ (cid:101) H sδ (Ω), ||| u ||| γ = (cid:90) (cid:90) (Ω ∪ Ω I ) ( u ( x ) − u ( y )) γ s ( x , y ) d y d x ≥ λ (cid:90) (cid:90) (Ω ∪ Ω I ) ( u ( x ) − u ( y )) γ ( x , y ) d y d x = ||| u ||| γ . (3.8)Hence, (cid:101) H sδ (Ω) is continuously embedded in (cid:101) H sδ (Ω). Since (cid:101) H sδ (Ω) is compactly embedded in L (Ω), and L (Ω)is compactly embedded in (cid:101) H sδ (Ω) ∗ , (cid:16) (cid:101) H sδ (Ω) , L (Ω) , (cid:101) H sδ (Ω) ∗ (cid:17) forms a Gelfand triple, and the embedding of (cid:101) H sδ (Ω) into (cid:101) H sδ (Ω) ∗ is compact. This fact, and the G˚arding inequality (3.7), directly imply the well-posednessof a perturbation of problem (2.15) and a bound on its solution, as shown in the following lemma. In whatfollows, we denote the duality pairing between (cid:101) H sδ (Ω) and its dual in the usual manner, i.e. (cid:104)· , ·(cid:105) . Lemma
Given a positive constant c and a function f ∈ (cid:101) H sδ (Ω) ∗ , the problemFind u ∈ (cid:101) H sδ (Ω) such that A ( u, v ) + c ( u, v ) L (Ω) = (cid:104) f, v (cid:105) , ∀ v ∈ (cid:101) H sδ (Ω) , and u = 0 in Ω I is well-posed and its solution is such that ||| u ||| γ ≤ C || f || (cid:101) H sδ (Ω) ∗ , for a positive constant C . The proof is simply based on application of (3.7) to the perturbed bilinear form, and is hence omitted.Note that (3.8) also allows us to apply [33, Theorem 4.1], which states that the bilinear form A satisfies aweak maximum principle, which we report for completeness. emma Let γ and A be defined as in (2.4) and (2.16) and let u ∈ (cid:101) H sδ (Ω) satisfy A ( u, v ) ≤ ∀ v ∈ (cid:101) H sδ (Ω) , then, sup Ω u ≤ . In particular, (3.9) A ( u, v ) = 0 ∀ v ∈ (cid:101) H sδ (Ω) ⇐⇒ u = 0 . An immediate consequence of Lemma 3.4 is the fact that the operator K := ( L + cI ) − mapping (cid:101) H sδ (Ω) ∗ → (cid:101) H sδ (Ω), is a compact operator. This fact allows us to use the Fredholm alternative theorem toprove the well-posedness of problem (2.15). Theorem
Given a function f ∈ (cid:101) H sδ (Ω) ∗ , the problemFind u ∈ (cid:101) H sδ (Ω) such that A ( u, v ) = (cid:104) f, v (cid:105) , ∀ v ∈ (cid:101) H sδ (Ω) , and u = 0 in Ω I has a unique solution u ∈ (cid:101) H sδ (Ω) .Proof. Consider the operator
K − ηI . We have the following identity: K − ηI = − η K ( K − − η − I ) = − η K ( L + ( c − η − ) I ) . By the Fredholm alternative, either η is an eigenvalue of K , or ( K − ηI ) − is a bounded linear operator. Weshow that the first option leads to a contradiction.Let η = 1 /c and assume that η is an eigenvalue of K . This implies that A ( u, v ) = 0 for all v ∈ (cid:101) H sδ (Ω)and, by (3.9), that u = 0. Hence, η is not an eigenvalue of K . This ends the proof since we can concludethat ( K − ηI ) − , and subsequently ( − η K − ηI ) − K = L − , are bounded linear operators. In other words,the variational problem associated with the operator L is well-posed.
4. Finite element discretization and efficient matrix assembly.
In this section we briefly in-troduce the finite element (FE) discretization of problem (2.15) and provide details regarding numericalintegration.We choose a family of finite-dimensional subspaces V h ⊂ (cid:101) H sδ (Ω), parameterized by the mesh size h . Wedefine the Galerkin approximation u h ∈ V h as the solution of the following problem: find u h ∈ V h such that(4.1) A ( u h , v h ) = (cid:90) Ω f v h d x , ∀ v h ∈ V h . As in [28], we consider finite element approximations for the case that both Ω and Ω I are polyhedral domains.We partition Ω into finite elements and denote by h the diameter of the largest element in the partition.We assume that the partition is shape-regular and quasi-uniform [10] as the grid size h → V h to consist of piecewise polynomials of degree no more than m defined with respect to thepartition.We perform numerical integration on pairs of elements, by rewriting the bilinear form (2.17) as the sumof integrals over the partition, i.e. A ( u h , v h ) = 12 (cid:88) K (cid:88) (cid:101) K (cid:90) (cid:90) K × (cid:101) K ( u h ( x ) − u h ( y )) ( γ ( x , y ) v h ( x ) − γ ( y , x ) v h ( y )) d y d x . The design of quadrature rules for the integrals above depends on the regularity properties of thefunctions s and φ . If s and φ are constant on a pair of elements K × (cid:101) K , we can use the state-of-the-art quadrature rules described in [2, 1] designed for the constant-order fractional Laplacian, i.e. for thecase δ = ∞ . However, when this is not the case, the use of more sophisticated and often more expensivequadrature rules, such as adaptive quadrature rules [48], becomes mandatory to prevent the integration errorto be dominant. Conditions for well-posedness derived in Proposition 3.1 are such that any non-symmetrickernel falls into the latter category, hence requiring higher computational effort. In addition, note that, as or the constant-order fractional Laplacian, when K ∩ (cid:101) K (cid:54) = ∅ , the integrand function features a singularityat x = y . In this case, the domain of integration must be partitioned into subdomains in order to avoid thesingularity.In the specific case of piecewise constant fractional order, we speed up the assembly of the discreteoperator and subsequent matrix-vector multiplications by generalizing the panel-clustering approach of [2, 1]developed for the constant-order fractional Laplacian. Informally speaking, in the standard version of thisapproach, the unknowns are recursively grouped into nodes of a tree structure, with the root node containingall unknowns. Pairs of clusters correspond to sub-blocks of the discrete operator. If certain conditions ontheir location in physical space are satisfied, sub-blocks are approximated using Chebyshev interpolation,otherwise they are assembled using quadrature a . We generalize this algorithm to finite horizon and piecewiseconstant coefficients by restricting the approximation of matrix sub-blocks to cluster pairs that are strictlycontained within the horizon δ of each other and whose interaction reduces to a constant fractional kernel.
5. Numerical tests.5.1. Comparison of different types of kernels.
In order to highlight the different behaviors ofavailable fractional kernels, we restrict ourselves to a one-dimensional problem where model parametersare constant over two subdomains. This configuration can be considered a one-dimensional counterpartof the two-dimensional configurations displayed in Figure 1.1. Specifically, we let Ω = ( − , δ = 1,and consider kernels given by (2.5) that are piecewise constant with respect to the partition Ω ∪ Ω I =((Ω ∪ Ω I ) ∩ R − ) ∪ ((Ω ∪ Ω I ) ∩ R + ). In order for the assumptions of Proposition 3.1 to be satisfied, suchkernels need to be symmetric, i.e. s and φ need to take the form ψ ( x, y ; ψ + , ψ − , ψ ± ) = ψ + if x, y > ,ψ − if x, y < ,ψ ± else . Note that this excludes the interesting cases, that can be found in the literature, of piecewise constant φ or s that only depend on x . Hence, we also consider functions of the following form and compare them withthe proposed model. η ( x ; η + , η − , r, κ ) = η + if x ≥ r ( η + + η − ) + ( η + − η − ) arctan κ x arctan κr if − r ≤ x ≤ rη − if x ≤ − r. Here, r determines the size of the transition region between the two states η + and η − , and κ the steepnessof the transition. Since η is Lipschitz continuous, we can use η for both s and φ and satisfy the conditionsof Proposition 3.1.We numerically solve the diffusion problem (2.3) for homogeneous Dirichlet data g ( x ) = 0 and forcingterm f ( x ) = X | x − . | < . . In all tests conducted in this section we use continuous piecewise linear FE on auniform mesh of size h = 2 − .With the purpose of highlighting how a variable fractional order affects the qualitative behavior of thesolution across the interface, we first compare the family of symmetric kernels reported in Table 5.1 anddisplay the corresponding solutions in Figure 5.1. We observe that the solutions for constant s seem to bedifferentiable, whereas the cases with piecewise defined s are only continuous at x = 0.In a second experiment, we analyze how a variation of the fractional order of interactions betweenthe two sub-domains affects the regularity of the solution. We choose s = ψ ( · ; 0 . , . , α ) for α ∈{ . , . , . , . , . } and φ = 1 and report the corresponding solutions in Figure 5.2; we observe that thechoice of the fractional order for the interactions impacts the entity of the jump of the derivative at x = 0.In a third experiment, we analyze the impact of the parameter φ . We keep s = 0 .
75 fixed and vary the in-teraction between the two sub-domains by setting φ = ψ ( · ; 0 . , . , β ) for β ∈ { . , . , . , . , , . } .Results are reported in Figure 5.3; we observe that while β affects the steepness of the transition from x < x >
0, the derivative of the solution does not display any jumps. dentifier s φ const, const 0 .
75 1const, sym 0 . ψ ( · ; 1 , , . ψ ( · ; 0 . , . , .
5) 1sym, sym ψ ( · ; 0 . , . , . ψ ( · ; 1 , , . Table 5.1
Symmetric kernels − . − . − . − .
25 0 .
00 0 .
25 0 .
50 0 .
75 1 . . . . . . . . .
07 const, constconst, symsym, constsym, sym
Fig. 5.1 . Solutions for symmetric kernels with s and φ given in Table 5.1. − . − . − . − .
25 0 .
00 0 .
25 0 .
50 0 .
75 1 . . . . . . . . α = 0 . α = 0 . α = 0 . α = 0 . α = 0 . Fig. 5.2 . Solutions for symmetric kernels with s = ψ ( · ; 0 . , . , α ) for α ∈ { . , . , . , . , . } and φ = 1 . . − . − . − .
25 0 .
00 0 .
25 0 .
50 0 .
75 1 . . . . . . . β = 1 . β = 1 β = 0 . β = 0 . β = 0 . β = 0 . Fig. 5.3 . Solutions for symmetric kernels with s = 0 . and φ = ψ ( · ; 0 . , . , β ) for β ∈ { . , . , . , . , , . } . identifier s φ const, const 0 .
75 1const, nonsym 0 . η ( · ; 1 , . , . , η ( · ; 0 . , . , . , η ( · ; 0 . , . , . , η ( · ; 1 , . , . , Table 5.2
Non-symmetric kernels − . − . − . − .
25 0 .
00 0 .
25 0 .
50 0 .
75 1 . . . . . . . . .
07 const, constconst, nonsymnonsym, constnonsym, nonsym
Fig. 5.4 . Solutions for non-symmetric kernels with s and φ given in Table 5.2. − − − h − − − problem: constant – s=leftRightFractionalOrder(ll=0.25,rr=0.75,lr=0.5,rl=0.5,sym=1) L numerical h . ˜ H s numerical h . Fig. 5.5 . Convergence in L and energy norm for a one-dimensional example. With the purpose of determining the impact of a nonsymmetric kernel on the qualitative behavior ofthe solutions, in a forth set of experiments, we consider non-symmetric kernels, as given in Table 5.2, andreport the corresponding solutions in Figure 5.4. These results do not highlight a significant change in thequalitative behavior of u across the interface.We conclude that the restriction to constant fractional orders s limits the types of solution behaviorthat can be recovered, even when allowing for a variable coefficient function φ . On the other hand, usingnonsymmetric kernels does not appear to lead to significantly different solution behavior. Given the signifi-cant drawbacks in terms of implementation difficulties posed by nonsymmetric kernels, we therefore restrictourselves to symmetric kernels with constant φ in what follows. We consider the convergence of solutions withrespect to the mesh size h . In light of the results presented in the previous section, we only focus onsymmetric kernel functions γ . Since, in general, no closed-form analytic solutions of (2.3) are known, wecompute errors with respect to a reference solution u h obtained on a fine mesh. We compute the energy and L norms as follows: (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) u h − u h (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) γ = A ( u h − u h , u h − u h ) = (cid:90) Ω f ( u h − u h ) d x , (cid:12)(cid:12)(cid:12)(cid:12) u h − u h (cid:12)(cid:12)(cid:12)(cid:12) L (Ω) = (cid:90) Ω ( u h − u h ) d x . Firstly, we consider the one-dimensional case where Ω = ( − ,
1) and δ = 1, with parameters s ( x , y ) = ψ ( x , y ; 0 . , . , . φ ≡ f ≡
1. We display the obtained errors in Figure 5.5. We observe that (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) u h − u h (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) γ ∼ h / , whereas (cid:12)(cid:12)(cid:12)(cid:12) u h − u h (cid:12)(cid:12)(cid:12)(cid:12) L (Ω) ∼ h / . The apparent speedup for smaller values of the meshsize h are due to the fact that we are comparing against u h , and not against the unknown exact solution u . It is known[11] that, for constant fractional kernels, the energy norm error converges as O ( h / − ε ) forany ε >
0, whereas the L -error converges as O ( h min { , / s }− ε ). Our numerical results suggest that thisresult can be generalized to O ( h / − ε ) convergence in energy norm and O ( h min { , / s }− ε ) convergence in L -norm for the case of variable fractional order. Here, s is defined as in (2.6).In a second convergence experiment, we consider the two-dimensional case of Ω = ( − , , and δ = 1 / a See the works by Ainsworth and Glusa[2, 1] for a detailed description of this approach for constant parameters and infinitehorizon. 13 . − . − . . . . . − . − . . . . . . . . . Fig. 5.6 . Solution for constant forcing, finite horizon and fractional order (5.1) with four layers. with φ ≡ f ≡ s with four layers: s ( x , y ) = 12 ( σ ( x ) + σ ( y )) with σ ( z ) = / z < − / , / − / ≤ z < , / ≤ z < / , / / ≤ z. (5.1)This definition is such that the resulting configuration mimics the realistic setting displayed in Figure 1.1,right. The reference solution for h = 0 . n ≈ ,
000 is displayed in Figure 5.6. As expected, it can beobserved that the solution behavior is more diffusive for layers with larger fractional order. Convergenceresults are reported in Figure 5.7, where we display the energy and L -norm of the discretization error withrespect to a reference solution. Similar to the one-dimensional case above, the observed rates suggest thatthe convergence results for constant kernels can be generalized to the variable coefficient case. We consider a configuration that mimics therealistic setting of Figure 1.1, left. Specifically, we consider the two-dimensional case of Ω = B ( ) andinfinite horizon δ = ∞ . We set φ ≡ s by inclusions, i.e. s ( x , y ) = .
25 if | x i | , | y i | ∈ [0 . , . , i = 1 , , .
75 if | x i | , | y i | (cid:54)∈ [0 . , . , i = 1 , ,α else , (5.2)for α ∈ { . , . , . } . We enforce the homogeneous Dirichlet condition on Ω I = R \ Ω using Gauss’stheorem, see [2]. Note that this shortcut can be applied because, for a given x ∈ Ω, s ( x , y ) and s ( y , x ) areconstant for all y ∈ Ω I . Results are reported in Figure 5.8; there, we observe that the solution behavioris less diffusive in the inclusion regions than in the surrounding domain. Also, as in the one-dimensionaltest cases, we observe that the choice of the fractional order for the interactions impacts the entity of thejump of the derivative at the interface between the inclusions and the rest of the domain. This fact clearlyemphasizes the need of identification methods for variable-order fractional models. In particular, in this case,it would be sufficient to simply learn the interaction order α to better predict the behavior of the solutionacross the interface. − h − − problem: constant – s=layersFractionalOrder(numLayers=4) L numerical h . ˜ H s numerical h . Fig. 5.7 . Convergence in L and energy norm for a two-dimensional example with four material layers.
6. Concluding remarks.
This work introduces a new variable-order fractional model that, comparedwith state-of-the-art alternatives, features improved variability. The latter is obtained by allowing thefractional order s to be a function of both x and y , hence allowing an agile treatment of changes in theproperties of the underlying physical system. This new model finds practical use in presence of materialheterogeneities and, in particular, in presence of abrupt changes in material properties, i.e. of physicalinterfaces.Our theoretical results guarantee the feasibility of the proposed operator by showing that, under certainconditions on the model parameters, the associated diffusion problem is well-posed. In fact, we point outthat by proving coercivity of the time-independent problem, we automatically guarantee well-posedness ofthe associated parabolic problem. On the other hand, our computational tests not only do they illustrate theimproved descriptive power of the proposed model, but they also show that solutions to doubly-variable orderequation converge with the same rates as the ones corresponding to the minimum value, s , of the fractionalorder. Furthermore, our implementation of the finite-element matrix assembly proves to be efficient inpresence of piecewise constant definitions of s ( x , y ).Clearly, the problem of finding the fractional order profile that best fits a physical system remains openand raises many challenges. However, our proposed model and the corresponding numerical tests show thatthe identification problem can be reduced to determining a handful of parameters, namely, the fractionalorders in each subregion and the interaction orders. Furthermore, the dependence on both x and y allowsus to consider symmetric kernels, simplifying significantly the finite element implementation and allowingfor more efficient assembly algorithms. Acknowledgments.Financial disclosure.
MD and CG are supported by Sandia National Laboratories (SNL), SNL is amultimission laboratory managed and operated by National Technology and Engineering Solutions of Sandia,LLC., a wholly owned subsidiary of Honeywell International, Inc., for the U.S. Department of Energys Na-tional Nuclear Security Administration contract number DE-NA0003525. This work was supported throughthe Sandia National Laboratories Laboratory-directed Research and Development (LDRD) program, project218318 and by the U.S. Department of Energy, Office of Advanced Scientific Computing Research under theCollaboratory on Mathematics and Physics-Informed Learning Machines for Multiscale and MultiphysicsProblems (PhILMs) project. This paper describes objective technical results and analysis. Any subjectiveviews or opinions that might be expressed in the paper do not necessarily represent the views of the U.S.Department of Energy or the United States Government. SAND Number: SAND2021-0842 O. . − . − . . . . . − . − . . . . . . . . . . . . − . − . − . . . . . − . − . . . . . . . . . − . − . − . . . . . − . − . . . . . . . . . Fig. 5.8 . Solution for constant forcing, infinite horizon and fractional order (5.2) with inclusions. From top to bottom: α ∈ { . , . , . } . onflict of interest. The authors declare no potential conflict of interests.
Supporting information.
None reported.
REFERENCES[1] M. Ainsworth and C. Glusa. Aspects of an adaptive finite element method for the fractional Laplacian: A priori and aposteriori error estimates, efficient implementation and multigrid solver.
Computer Methods in Applied Mechanicsand Engineering , 327:4–35, 2017.[2] M. Ainsworth and C. Glusa. Towards an efficient finite element method for the integral fractional Laplacian on polygonaldomains. In
Contemporary Computational Mathematics-A Celebration of the 80th Birthday of Ian Sloan , pages17–57. Springer, 2018.[3] B. Alali and M. Gunzburger. Peridynamics and material interfaces.
Journal of Elasticity , 120(2):225–248, 2015.[4] H. Antil, E. Otarola, and A.J. Salgado. Optimization with respect to order in a fractional diffusion model: Analysis,approximation and algorithmic aspects.
Journal of Scientific Computing , 77:204 – 224, 2018.[5] H. Antil and M. Warma. Optimal control of fractional semilinear PDEs.
ESAIM Control Optimisation and Calculus ofVariations , 2019. To appear.[6] Harbir Antil and Carlos N Rautenberg. Sobolev spaces with non-muckenhoupt weights, fractional elliptic operators, andapplications.
SIAM Journal on Mathematical Analysis , 51(3):2479–2503, 2019.[7] O. Bakunin.
Turbulence and Diffusion: Scaling versus Equations . Springer-Verlag, New York, 2008.[8] P. W. Bates and A. Chmaj. An integrodifferential model for phase transitions: stationary solutions in higher spacedimensions.
J. Statist. Phys. , 95:1119–1139, 1999.[9] D. A. Benson, S. W. Wheatcraft, and M. M. Meerschaert. Application of a fractional advection-dispersion equation.
Water Resources Research , 36(6):1403–1412, 2000.[10] S. Brenner and R. Scott.
The Mathematical Theory of Finite Element Methods . Springer-Verlag, New York, 2008.[11] O. Burkovska and M. Gunzburger. Regularity analyses and approximation of nonlocal variational equality and inequalityproblems.
Journal of Mathematical Analysis and Applications , 478(2):1027 – 1048, 2019.[12] Olena Burkovska and Max Gunzburger. Affine Approximation of Parametrized Kernels and Model Order Reduction forNonlocal and Fractional Laplace Models.
SIAM Journal on Numerical Analysis , 58(3):1469–1494, 2020.[13] Olena Burkovska and Max Gunzburger. On a nonlocal Cahn–Hilliard model permitting sharp interfaces. arXiv:2004.14379,2020.[14] G. Capodaglio, M. D’Elia, P. Bochev, and M. Gunzburger. An energy-based coupling approach to nonlocal interfaceproblems.
Computers and Fluids , 2019. To appear.[15] C. K. Chen and P. C. Fife. Nonlocal models of phase transitions in solids.
Advances in Mathematical Sciences andApplications , 10(2):821–849, 2000.[16] Y.-M. Chen, Y.-Q. Wei, D.-Y. Liu, and H. Yu. Numerical solution for a class of nonlinear variable order fractionaldifferential equations with legendre wavelets.
Applied Mathematics Letters , 46:83–88, 2015.[17] C. Cortazar, M. Elgueta, J. Rossi, and N. Wolanski. How to approximate the heat equation with Neumann boundaryconditions by nonlocal diffusion problems.
Archive for Rational Mechanics and Analysis , 187:137–156, 2008.[18] M. D’Elia, Q. Du, C. Glusa, M. Gunzburger, X. Tian, and Z. Zhou. Numerical methods for nonlocal and fractionalmodels.
Acta Numerica , 2020. To appear.[19] M. D’Elia, Q. Du, M. Gunzburger, and R. Lehoucq. Nonlocal convection-diffusion problems on bounded domains andfinite-range jump processes.
Computational Methods in Applied Mathematics , 29:71–103, 2017.[20] M. D’Elia, C. Glusa, and E. Ot´arola. A priori error estimates for the optimal control of the integral fractional Laplacian.
SIAM Journal on Control and Optimization , 57:2775–2798, 2019.[21] M. D’Elia, M. Gulian, H. Olson, and G. E. Karniadakis. A unified theory of fractional, nonlocal, and weighted nonlocalvector calculus. arXiv:2005.07686, 2020.[22] M. D’Elia and M. Gunzburger. The fractional Laplacian operator on bounded domains as a special case of the nonlocaldiffusion operator.
Computers and Mathematics with applications , 66:1245–1260, 2013.[23] M. D’Elia and M. Gunzburger. Optimal distributed control of nonlocal steady diffusion problems.
SIAM Journal onControl and Optimization , 55:667–696, 2014.[24] M. D’Elia and M. Gunzburger. Identification of the diffusion parameter in nonlocal steady diffusion problems.
AppliedMathematics and Optimization , 73:227–249, 2016.[25] M. D’Elia, M. Gunzburger, and C. Vollman. A cookbook for finite element methods for nonlocal problems, includingquadrature rule choices and the use of approximate neighborhoods. arXiv:2005.10775, 2020.[26] M. D’Elia, J.-C. De los Reyes, and A. Miniguano Trujillo. Bilevel parameter optimization for nonlocal image denoisingmodels. arXiv:1912.02347, 2019.[27] M. D’Elia, X. Tian, and Y. Yu. A physically-consistent, flexible and efficient strategy to convert local boundary conditionsinto nonlocal volume constraints.
Accepted for publication in SIAM Journal of Scientific Computing , 2020.[28] Q. Du, M. Gunzburger, R. Lehoucq, and K. Zhou. Analysis and approximation of nonlocal diffusion problems with volumeconstraints.
SIAM Review , 54(4):667–696, 2012.[29] Q. Du, L. Ju, Xiao Li, and Z. Qiao. Maximum principle preserving exponential time differencing schemes for the nonlocalallen–cahn equation.
SIAM Journal on Numerical Analysis , 57(2):875–898, 2019.[30] Q. Du and J. Yang. Asymptotically Compatible Fourier Spectral Approximations of Nonlocal Allen–Cahn Equations.
SIAM Journal on Numerical Analysis , 54(3):1899–1919, 2016.[31] Qiang Du, Max Gunzburger, R. B. Lehoucq, and Kun Zhou. A nonlocal vector calculus, nonlocal volume-constrained17roblems, and nonlocal balance laws.
Mathematical Models and Methods in Applied Sciences , 23:493–540, 2013.[32] Qiang Du, Zhan Huang, and Richard B. Lehoucq. Nonlocal convection-diffusion volume-constrained problems and jumpprocesses.
Discrete & Continuous Dynamical Systems - B , 19:373, 2014.[33] M. Felsinger, M. Kassmann, and P. Voigt. The Dirichlet problem for nonlocal operators.
Mathematische Zeitschrift ,279:779–809, 2015.[34] R. Gorenflo and F. Mainardi. Fractional calculus: integral and differential equations of fractional order.
Fractals andFractional Calculus in Continuum Mechanics , pages 223–276, 1997.[35] G. Grubb. Fractional Laplacians on domains, a development of H¨ormander’s theory of µ -transmission pseudodifferentialoperators. Adv. Math. , 268:478–528, 2015.[36] M. Gulian, M. Raissi, P. Perdikaris, and G. E. Karniadakis. Machine learning of space-fractional differential equations.
SIAM Journal on Scientific Computing , 41(4):A2485–A2509, 2019.[37] Y. D. Ha and F. Bobaru. Characteristics of dynamic brittle fracture captured with peridynamics.
Engineering FractureMechanics , 78(6):1156–1168, 2011.[38] A. Lischke, G. Pang, M. Gulian, F. Song, C. Glusa, X. Zheng, Z. Mao, W. Cai, M. M. Meerschaert, M. Ainsworth,and G. E. Karniadakis. What is the fractional Laplacian? A comparative review with new results.
Journal ofComputational Physics , 404(109009), 2020.[39] D. J. Littlewood. Simulation of dynamic fracture using peridynamics, finite element modeling, and contact. In
ASME2010 International Mechanical Engineering Congress and Exposition , pages 209–217. American Society of MechanicalEngineers Digital Collection, 2010.[40] M. M. Meerschaert, J. Mortensen, and S. W. Wheatcraft. Fractional vector calculus for fractional advection–dispersion.
Physica A: Statistical Mechanics and its Applications , 367:181–190, 2006.[41] T. Mengesha and Q. Du. Analysis of a scalar nonlocal peridynamic model with a sign changing kernel.
Discrete &Continuous Dynamical Systems-B , 18(5):1415–1437, 2013.[42] H. A. Olson, M. Gulian, and M. D’Elia. The tempered fractional laplacian as a special case of the nonlocal laplaceoperator. pages 111–126. Computer Science Research Institute Summer Proceedings 2020, A.A. Rushdi and M.L.Parks, eds., Sandia National Laboratories, 2020. Technical Report SAND2020-12580R.[43] H. Ouchi, A. Katiyar, J. York, J. T. Foster, and M. M. Sharma. A fully coupled porous flow and geomechanics model forfluid driven cracks: a peridynamics approach.
Computational Mechanics , 55(3):561–576, 2015.[44] G. Pang, M. D’Elia, M. Parks, and G. E. Karniadakis. nPINNs: nonlocal Physics-Informed Neural Networks for aparametrized nonlocal universal Laplacian operator. Algorithms and Applications. arXiv:2004.04276, 2020.[45] G. Pang, L. Lu, and G. E. Karniadakis. fPINNs: Fractional physics-informed neural networks.
SIAM Journal on ScientificComputing , 41:A2603–A2626, 2019.[46] G. Pang, P. Perdikaris, W. Cai, and G. E. Karniadakis. Discovering variable fractional orders of advection–dispersionequations from field data using multi-fidelity Bayesian optimization.
Journal of Computational Physics , 348:694 –714, 2017.[47] M. Pasetto.
Enhanced Meshfree Methods for Numerical Solution of Local and Nonlocal Theories of Solid Mechanics . PhDthesis, UC San Diego, 2019.[48] Robert Piessens, Elise de Doncker-Kapenga, Christoph W ¨Uberhuber, and David K Kahaner.
Quadpack: a subroutinepackage for automatic integration , volume 1. Springer Science & Business Media, 2012.[49] A. Razminia, A. F. Dizaji, and V. J. Majd. Solution existence for non-autonomous variable-order fractional differentialequations.
Mathematical and Computer Modelling , 55(3-4):1106–1117, 2012.[50] S. Rokkam, M. Gunzburger, M. Brothers, N. Phan, and K. Goel. A nonlocal peridynamics modeling approach for corrosiondamage and crack propagation.
Theoretical and Applied Fracture Mechanics , 101:373–387, 2019.[51] A. A. Schekochihin, S. C. Cowley, and T. A. Yousef. MHD turbulence: Nonlocal, anisotropic, nonuniversal? In
IUTAMSymposium on computational physics and new perspectives in turbulence , pages 347–354. Springer, 2008.[52] R. L. Schilling and J. Wang. Lower bounded semi-dirichlet forms associated with l´evy type operators. In
FestschriftMasatoshi Fukushima: In Honor of Masatoshi Fukushima’s Sanju , pages 507–526. World Scientific, 2015.[53] R Schneider, O Reichmann, and Christoph Schwab. Wavelet solution of variable order pseudodifferential equations.
Calcolo , 47(2):65–101, 2010.[54] R. Schumer, D. A. Benson, M. M. Meerschaert, and B. Baeumer. Multiscaling fractional advection-dispersion equationsand their solutions.
Water Resources Research , 39(1):1022–1032, 2003.[55] R. Schumer, D.A. Benson, M.M. Meerschaert, and S.W. Wheatcraft. Eulerian derivation of the fractional advection-dispersion equation.
Journal of Contaminant Hydrology , 48:69–88, 2001.[56] S. A Silling and E. Askari. A meshfree method based on the peridynamic model of solid mechanics.
Computers &structures , 83(17-18):1526–1535, 2005.[57] H. Wang, K. Wang, and T. Sircar. A direct O ( N log N ) finite difference method for fractional diffusion equations. Journalof Computational Physics , 229(21):8095–8104, 2010.[58] Chester J Weiss, Bart G van Bloemen Waanders, and Harbir Antil. Fractional operators applied to geophysical electro-magnetics.
Geophysical Journal International , 220(2):1242–1259, 2020.[59] H. You, Y. Yu, N. Trask, M. Gulian, and M. D’Elia. Data-driven learning of robust nonlocal physics from high-fidelitysynthetic data. arXiv:2005.10076, 2020.[60] F. Zeng, Z. Zhang, and G. E. Karniadakis. A generalized spectral collocation method with tunable accuracy for variable-order fractional differential equations.
SIAM Journal on Scientific Computing , 37(6):A2710–A2732, 2015.[61] X. Zheng and H. Wang. An optimal-order numerical approximation to variable-order space-fractional diffusion equationson uniform or graded meshes.
SIAM Journal on Numerical Analysis , 58(1):330–352, 2020.[62] X. Zheng and H. Wang. Wellposedness and regularity of a variable-order space-time fractional diffusion equation.
Analysis nd Applications , 18(04):615–638, 2020.[63] Xiangcheng Zheng, Yiqun Li, Jin Cheng, and Hong Wang. Inverting the variable fractional order in a variable-orderspace-fractional diffusion equation with variable diffusivity: analysis and simulation. Journal of Inverse and Ill-posedProblems , page 000010151520190040, 10 Nov. 2020.[64] P. H. Zhuang, F. W. Liu, V. Anh, and I. Turner. Numerical methods for the variable-order fractional advection- diffusionequation with a nonlinear source term.