# A fractional model for anomalous diffusion with increased variability. Analysis, algorithms and applications to interface problems

Fractional equations have become the model of choice in several applications where heterogeneities at the microstructure result in anomalous diffusive behavior at the macroscale. In this work we introduce a new fractional operator characterized by a doubly-variable fractional order and possibly truncated interactions. Under certain conditions on the model parameters and on the regularity of the fractional order we show that the corresponding Poisson problem is well-posed. We also introduce a finite element discretization and describe an efficient implementation of the finite-element matrix assembly in the case of piecewise constant fractional order. Through several numerical tests, we illustrate the improved descriptive power of this new operator across media interfaces. Furthermore, we present one-dimensional and two-dimensional h-convergence results that show that the variable-order model has the same convergence behavior as the constant-order model.

## Authors

• 17 publications
• 5 publications
• ### Finite element approximation of fractional Neumann problems

In this paper we consider approximations of Neumann problems for the int...
08/13/2020 ∙ by Francisco M. Bersetche, et al. ∙ 0

• ### Linear and Nonlinear Fractional Diffusion

This paper surveys recent analytical and numerical research on linear pr...
06/10/2019 ∙ by Juan Pablo Borthagaray, et al. ∙ 0

• ### Fractional Operators Applied to Geophysical Electromagnetics

A growing body of applied mathematics literature in recent years has foc...
02/13/2019 ∙ by Chester J. Weiss, et al. ∙ 0

• ### Strong convergence of some Euler-type schemes for the finite element discretization of time-fractional SPDE driven by standard and fractional Brownian motion

In this work, we provide the first strong convergence result of numerica...
11/07/2020 ∙ by Aurelien Junior Noupelah, et al. ∙ 0

• ### Double exponential quadrature for fractional diffusion

We introduce a novel discretization technique for both elliptic and para...
12/10/2020 ∙ by Alexander Rieder, et al. ∙ 0

• ### A stencil scaling approach for accelerating matrix-free finite element implementations

We present a novel approach to fast on-the-fly low order finite element ...
09/20/2017 ∙ by Simon Bauer, et al. ∙ 0

• ### Space-Fractional Diffusion with Variable Order and Diffusivity: Discretization and Direct Solution Strategies

We consider the multidimensional space-fractional diffusion equations wi...
08/29/2021 ∙ by Hasnaa Alzahrani, et al. ∙ 0

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 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, several scientific and engineering applications exhibit hierarchical features that cannot be described by classical models. As an example we mention applications in continuum mechanics

[ha2011characteristics, littlewood2010simulation, ouchi2015fully][Bates1999, Burkovska2020CH, Chen_nonlocalmodels], corrosion [Rokkam2019], turbulence [Bakunin2008, Pang2020, Schekochihin2008], and geoscience [Benson2000, Meerschaert2006, Schumer2001, Schumer2003, WeissBloemenEtAl2020_FractionalOperatorsAppliedToGeophysicalElectromagnetics]. In this work we consider nonlocal operators of fractional type, i.e. integral operators characterized by singular and non-integrable kernels, that correspond to differential operators of fractional orders, as opposed to the integer order in the PDE case.

Note that fractional differential operators are almost as old as their integer counterpart [Gorenflo1997]; however, their usability has increased during the last decades thanks to progress in computational capabilities and to a better understanding of their descriptive power. The simplest form of a fractional operator is the fractional Laplacian; in its integral form, its action on a scalar function is defined as [Gorenflo1997]

 (−Δ)su(x)=Cn,sp.v.∫Rnu(x)−u(y)|x−y|n+2sdy, (1)

where is the fractional order, the spatial dimension, a constant and where p.v. indicates the principal value. It follows from (1) that in fractional modeling the state of a system at a point depends on 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 [DElia2020, DElia2013, Du2012, Pang2020] of the following form [DElia2017]:

 −Lu(x)=∫Bδ(x)(u(x)γ(x,y)−u(y)γ(y,x))dy. (2)

Here, interactions are limited to a Euclidean ball of radius , often referred to as horizon or interaction radius. The kernel is a modeling choice and determines regularity properties of the solution. Note that for and for the fractional-type kernel the nonlocal operator in (2) is equivalent to the fractional Laplacian in (1). Also, it has been shown in [DElia2013] that for that choice of and solutions corresponding to the nonlocal operator (2) converge to the ones corresponding to the fractional operator (1) as (see [DElia2020] for more convergence results and for a detailed classification of these operators and relationships between them).

In recent years, with the purpose of increasing the descriptive power of fractional operators, new models characterized by a variable fractional order have been introduced for both space- and time-fractional differential operators [Razminia2012, AntilRautenberg2019_SobolevSpacesWithNon, Zheng2020] and several discretization methods have been designed [Chen2015, Zeng2015, Zheng2020optimal-order, Zhuang2009, SchneiderReichmannEtAl2010_WaveletSolutionVariableOrderPseudodifferentialEquations]. However, the analysis of variable-order models is still in its infancy, with [Felsinger2015] and [Schilling2015] being perhaps the only relevant works that address theoretical questions such as well-posedness for space-fractional differential operators with variable order

. The improved descriptive power of variable-order operators has been demonstrated in some works on parameter estimation

[Pang2020, Pang2019fPINNs, Pang2017discovery, WANG]

; here, via machine-learning algorithms or more standard optimization techniques, the authors estimate the variable fractional power for operators of the form

 −Lu(x)=Cn,sp.v.∫Bδ(x)u(x)−u(y)|x−y|n+2s(x)dy, (3)

where is either infinite or finite.

Models such as the one in (3) are convenient to describe anomalous diffusion in case of heterogeneous materials or media, where different regions may be characterized by different diffusion rates, i.e. different values of , see Figure 1 for two-dimensional illustrations where is a piece-wise constant function defined over the domain. In other words, these models can describe physical interfaces by simply tuning the function . Modeling nonlocal interfaces is a nontrivial task due to the unknown nature of the nonlocal interactions across the interfaces. Only a few works in the literature have addressed this problem; among them, we mention [Alali2015] for nonlocal models in mechanics and [Capodaglio2020] for nonlocal diffusion models with integrable kernels.

In this work, we tackle the nonlocal interface problem for fractional-type operators using a generalization of the nonlocal operator in (3). Specifically, we add variability to the fractional order and the kernel function itself. Our main contributions are:

• The introduction of a new variable-order fractional operator characterized by and . This choice enables modeling of a much broader set of interface behaviors and, for symmetric , symmetrizes the kernel, making analysis and implementation a much easier task. In fact, note that in (3) the kernel is no longer symmetric, requiring more sophisticated quadrature rules for the integration (regardless of the discretization technique of choice).

• The analysis of well-posedness of the Poisson problem associated with the new operator in the general nonsymmetric case.

• The design of an efficient finite-element matrix assembly technique for the practical case of piecewise constant 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 the new variable-order fractional operator and the corresponding strong and weak forms of the nonsymmetric diffusion problem. We also describe properties of the associated energy norm and space. In Section 3 we report the main result of the paper that proves via Fredholm alternative the well-posedness of the diffusion problem. In Section 4 we introduce a finite element discretization and propose efficient quadrature rules for the integration of the nonsymmetric kernel. In Section 5 we report several one- and two-dimensional numerical tests that illustrate both the improved variability of our model and the accuracy of the numerical scheme. Finally, in Section 6 we summarize our findings.

## 2 A new variable-order fractional model

In this section we first introduce the notation that will be used throughout the paper and recall state-of-the-art results on nonsymmetric kernels and fractional kernels with variable order. Then, we introduce a new variable-order fractional kernel and discuss the properties of the associated norm and energy space.

### 2.1 Preliminaries

The purpose of this section is two-fold. First, we recall the formulation of a symmetric nonlocal diffusion problem with finite interaction radius following the theory presented in [Du2012, Du2013]; then, we recall the theory presented in [Felsinger2015] on variable-order fractional operators (with infinite interaction radius). Our new model, presented in the next subsection, is a combination of the two.

Let be an open bounded domain. We define its interaction domain as the set of points outside that interact with points inside , i.e.

 ΩI:={y∈Rn∖Ω:|x−y|≤δ,forx∈Ω}, (4)

where , referred to as interaction radius, determines the extent of the nonlocal interactions. Also, let be a nonnegative, symmetric, kernel function, not necessarily integrable. The diffusion operator introduced in [Du2013] is defined as

 −Lu(x)=p.v.∫Ω∪ΩI(u(x)−u(y))γ(x,y)dy. (5)

From now on we remove the explicit reference to the principal value when the operator is used within an equation. The strong form of a truncated, symmetric, diffusion problem is then given by: for , and , find such that

 ∫Ω∪ΩI(u(x)−u(y))γ(x,y)dy=f(x) x∈Ω (6) u(x)=g(x) x∈ΩI.

The well-posedness of problem (6) was studied in [Du2012]. Extensions of this model to the nonsymmetric case have been proposed in [Duadvection], for integrable kernels, and in [DElia2017] for non-integrable kernels by using the operator (2) introduced in the previous section. This operator is often referred to as nonsymmetric nonlocal diffusion operator or, more appropriately, nonlocal convection-diffusion operator as it introduces a nonlocal convection term. Since the non-symmetry of our problem is not related to convection, but only to the presence of an interface that affects the way a quantity diffuses, in this work, we do not employ (2); instead, we use the operator defined in (5) with a nonsymmetric kernel .

Among the state-of-the-art formulations of nonlocal diffusion problems with nonsymmetric kernels, we mention the one introduced in [Felsinger2015] for integrable and non-integrable kernels with infinite support. Here the authors consider the nonlocal operator (5) where the kernel is a nonnegative measurable function, possibly nonsymmetric. Among the kernels studied in [Felsinger2015], we recall the following nonsymmetric variable-order fractional kernel, whose analysis is relevant for the theory presented in this paper:

 γ0(x,y)=ϕ0(x)|x−y|n+2s(x), (7)

where the non-symmetry with respect to is due to the terms and . The latter is a normalization functions with the same role as in the standard fractional Laplacian in (1) [Schilling2015]. Note that the choice of corresponds to the variable-order version of the standard fractional Laplacian operator introduced in (3). Furthermore, when and are constant functions, such operator is a multiple of the standard fractional Laplacian. The well-posedness of problem (6) for nonsymmetric kernels was analyzed in [Felsinger2015] in a very general setting and in [Schilling2015] for kernels of the form (7).

### 2.2 The variable-order model

We the purpose of increasing the descriptive power of the operator , we consider a fractional-type kernel with added variability. We define the new variable-order nonsymmetric fractional-type kernel as follows:

 γ(x,y) =ϕ(x,y)|x−y|n+2s(x,y)X|x−y|≤δ,x,y∈Ω∪ΩI, (8)

with constants , and such that

 0

Thus, the resulting nonlocal operator is given by the following expression

 −Lu(x)=p.v.∫(Ω∪ΩI)∩Bδ(x)(u(x)−u(y))ϕ(x,y)|x−y|n+2s(x,y)dy. (11)
###### Remark 1

As opposed to the kernel model presented in [Schilling2015], the function does not have a normalizing 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, the dependence on and describes changes in material properties across material interfaces.

###### Remark 2

For a specific choice of , the operator defined in (11) corresponds to a variable-order, tempered fractional Laplacian. Specifically, we let

 ϕ(x,y)=exp{−λ|x−y|}, (12)

where is the so-called tempering parameter, whose effect is to fasten the decay of the fractional kernel. Since the overall effect of tempering is similar to the one of truncation[Olson2020], in our numerical tests we do not study the tempered case and we limit ourselves to piecewise constant definitions of for .

###### Remark 3

When the kernel is symmetric. This case is much easier to deal with both in terms of analysis and computation. In fact, in the symmetric case, the theory introduced in [Du2012] can be readily applied and guarantees well-posedness of problem (6). Furthermore, numerical integration of a symmetric function poses minor challenges compared with the nonsymmetric case, which requires the employment of more sophisticated quadrature rule. As we show in our numerical tests, a symmetric still guarantees improved model variability across interfaces, compared to the single-variable model in (7).

We define the spaces

 H(Ω∪ΩI;γ) :={v:Ω∪ΩI→R:v∈L2(Ω∪ΩI),|||v|||γ<∞}, ˜Hsδ(Ω) :=HΩ(Ω∪ΩI;γ):={v∈H(Ω∪ΩI;γ):v|ΩI=0},

where the semi-norm is defined as

 |||v|||2γ:=∬(Ω∪ΩI)2(v(x)−v(y))2γs(x,y)dxdy.

Here, is the symmetric part of the kernel. We recall that the symmetric and anti-symmetric parts of the kernel are defined as follows:

 γs(x,y) =12(ϕ(x,y)|x−y|n+2s(x,y)+ϕ(y,x)|x−y|n+2s(y,x))X|x−y|≤δ, (13) γa(x,y) =12(ϕ(x,y)|x−y|n+2s(x,y)−ϕ(y,x)|x−y|n+2s(y,x))X|x−y|≤δ. (14)

The following bounds hold trivially:

 1|x−y|d+2s– ≤1|x−y|d+2s(x,y)≤1|x−y|d+2¯s if |x−y|≤1, (15) 1|x−y|d+2¯s ≤1|x−y|d+2s(x,y)≤1|x−y|d+2s– if |x−y|≥1, (16)

In what follows, we will also frequently make use of the following lower bound on the horizon:

 δ–:=min{δ,1}. (17)

The following Lemma shows that the semi-norm above satisfies a Poincaré inequality and, hence, equips with a norm.

###### Lemma 1

For all , the following Poincaré inequality holds:

 Cp|||u|||γ≥||u||L2(Ω∪ΩI).

###### Proof

Assumptions (9) and inequality (15) imply

 |||v|||2γ =∬(Ω∪ΩI)2,|x−y|≤δ–(v(x)−v(y))2γs(x,y)dxdy +∬(Ω∪ΩI)2,|x−y|≥δ–(v(x)−v(y))2γs(x,y)dxdy ≥C∬(Ω∪ΩI)2,|x−y|≤δ–(v(x)−v(y))21|x−y|n+2s–dxdy.

As a consequence, is lower bounded by the semi-norm of a -truncated fractional Sobolev space with fractional order , for which the Poincaré inequality holds, see [Du2012, Lemma 4.3].

The weak form. Let for simplicity; then, we can write the weak form of problem (6) as

 ∫Ω−Luvdx=∬(Ω∪ΩI)2(u(x)−u(y))γ(x,y)v(x)dydx=∫Ωfvdx, (18)

We denote the bilinear form in (18) by

 A(u,v)=∬(Ω∪ΩI)2γ(x,y)v(x)(u(x)−u(y)dydx. (19)

While this form is formally equivalent to the one introduced in [Felsinger2015], the kernel (8) does not belong to class of kernels studied in that paper. In the next section we show how to extend the well-posedness theory developed in [Felsinger2015] and [Schilling2015] to our new operator.

For discretization purposes, it is convenient to further rewrite as

 A(u,v)=12∬(Ω∪ΩI)2(u(x)−u(y))(γ(x,y)v(x)−γ(y,x)v(y))dydx. (20)

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 (18) by adapting the theory developed in the works of Felsinger, Kassmann and Voigt [Felsinger2015] and of Schilling and Wang [Schilling2015]. 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 the kernel and is a generalization of Proposition 3.1 in the paper by Schilling and Wang [Schilling2015].

###### Proposition 1

Assume that

 β(r):=sup|x−y|≤r|s(x,y)−s(y,x)| (21)

satisfies

 ∫10dr (β(r)|logr|)2r1+2¯s<∞. (22)

Moreover, assume that

 α(r):=sup|x−y|≤r|ϕ(x,y)−ϕ(y,x)|

is -Hölder continuous for arbitrarily small. Then there exists a constant such that

 supx∈Ω∪ΩI∫{γs(x,y)≠0}dy γa(x,y)2γs(x,y)≤Aγ<∞. (23)

###### Proof

First, by definition of symmetric and anti-symmetric parts of the kernel in (13), we have that . This and assumptions (9), (10) and inequality (16) imply that there exists a positive constant such that

 ∫{γs(x,y)≠0},|x−y|≥1dy γa(x,y)2γs(x,y)≤∫|x−y|≥1dy γs(x,y)≤C∫∞1r−1−2s–=:C1.

Moreover, due to the definition of ,

 ∫{γs(x,y)≠0},|x−y|∈[δ–,1]dy γa(x,y)2γs(x,y)=0

so that we obtain

 ∫{γs(x,y)≠0},|x−y|≥δ–dy γa(x,y)2γs(x,y)≤C1 (24)

Next, for , we write:

 γa(x,y) +12ϕ(y,x)|x−y|−n(|x−y|−2s(x,y)−|x−y|−2s(y,x)). =:γa,1(x,y)+γa,2(x,y)

Then, by (10) and the Hölder continuity of we have that there exists a positive constant such that

 ∫{γs(x,y)≠0},|x−y|≤δ–dy γa,1(x,y)2γs(x,y) =12∫{γs(x,y)≠0},|x−y|≤δ–dy (ϕ(x,y)−ϕ(y,x))2|x−y|−2n−4s(x,y)ϕ(x,y)|x−y|−n−2s(x,y)+ϕ(y,x)|x−y|−n−2s(y,x) ≤C∫{γs(x,y)≠0},|x−y|≤δ–dy (ϕ(x,y)−ϕ(y,x))2|x−y|n+2s(x,y) ≤C∫{γs(x,y)≠0},|x−y|≤δ–dy α(|x−y|)2|x−y|n+2¯s ≤C∫δ–0dr r2¯s+2εr1+2¯s=:C2. (25)

Here, we have used (15). On the other hand, since for

then, for , we have that

 (|x−y|−2s(x,y)−|x−y|−2s(y,x))2 ≤(log|x−y|)2(s(x,y)−s(y,x))2|x−y|−4min{s(x,y),s(y,x)}.

Hence, by (10) and assumption (22), there exists a positive constant such that

 ∫{γs(x,y)≠0},|x−y|≤δ–dy γa,2(x,y)2γs(x,y) = 12∫{γs(x,y)≠0},|x−y|≤δ–dy ϕ(y,x)2|x−y|−2n(|x−y|−2s(x,y)−|x−y|−2s(y,x))2ϕ(x,y)|x−y|−n−2s(x,y)+ϕ(y,x)|x−y|−n−2s(y,x) ≤ ≤ ≤ C∫δ–0dr (β(r)|logr|)2r1+2¯s=:C3. (26)

Here, we have again used (15). Hence, from (24), (25) and (26), we obtain:

 supx∈Ω∪ΩI∫Ω∪ΩIdy γa(x,y)2γs(x,y) ≤C1+2(C2+C3)=:Aγ.

and the result follows.

The next Gårding inequality follows directly from [Felsinger2015, Lemma 3.1]; we report its proof for completeness.

###### Lemma 2 (Gårding inequality)

There exist constants such that

 A(u,u)≥C|||u|||2γ−c||u||2L2(Ω∪ΩI)∀u∈˜Hsδ(Ω). (27)

###### Proof

By using Young’s inequality for , we obtain

 A(u,u) =∬(Ω∪ΩI)2γ(x,y)u(x)(u(x)−u(y))dydx ≥12∬(Ω∪ΩI)2γs(x,y)|u(x)−u(y)|2dydx−∬(Ω∪ΩI)2|γa(x,y)u(x)(u(x)−u(y))|dydx =12|||u|||2γ−∬(Ω∪ΩI)2∣∣γa(x,y)γs(x,y)1/2γs(x,y)−1/2u(x)(u(x)−u(y))∣∣dydx ≥12|||u|||γ−ε∬(Ω∪ΩI)2γs(x,y)|(u(x)−u(y))|2dydx−14ε∬{γs(x,y)≠0}2γa(x,y)2γs(x,y)u(x)2dydx ≥(12−ε)|||u|||2γ−Aγ4ε||u||2L2(Ω∪ΩI).

Hence, for small enough, we obtain

 A(u,u)≥C|||u|||2γ−c||u||2L2(Ω∪ΩI).

The Poincaré inequality and Proposition 1 yield the continuity of the bilinear form , as illustrated in the following lemma.

###### Lemma 3

For all

 |A(u,v)|≤C(CP,Aγ)|||u|||γ|||v|||γ,

i.e., the bilinear form is continuous on .

###### Proof

By combining the Poincaré inequality and Proposition 1, we have

 |A(u,v)| ≤|||u|||γ|||v|||γ+⎛⎜ ⎜⎝∬(Ω∪ΩI)2γs(x,y)(u(x)−u(y))2dydx⎞⎟ ⎟⎠1/2⎛⎜ ⎜⎝∬{γs(x,y)≠0}2γa(x,y)γs(x,y)v(x)2dydx⎞⎟ ⎟⎠1/2 ≤(1+CPA1/2γ)|||u|||γ|||v|||γ.

Before proceeding with the well-posedness analysis, we prove a result for a perturbation of the weak problem (18). First, let

 γ––(x,y):=1|x−y|n+2s–X|x−y|≤δ–,and˜Hs–δ–(Ω):=HΩ(Ω∪ΩI;γ––).

Then, there exists some such that

 γs(x,y)≥λγ––(x,y),

and hence for all ,

 |||u|||2γ=∬(Ω∪ΩI)2(u(x)−u(y))2γs(x,y)dydx ≥λ∬(Ω∪ΩI)2(u(x)−u(y))2γ––(x,y)dydx=|||u|||2γ––. (28)

Hence, is continuously embedded in . Since is compactly embedded in , and is compactly embedded in , forms a Gelfand triple, and the embedding of into is compact. This fact, and the Gårding inequality (27), directly imply the well-posedness of a perturbation of problem (18) and a bound on its solution, as shown in the following lemma. In what follows, we denote the duality pairing between and its dual in the usual manner, i.e. .

###### Lemma 4

Given a positive constant and a function , the problem

 Find u∈˜Hsδ(Ω) such thatA(u,v)+c(u,v)L2(Ω) =⟨f,v⟩,∀v∈˜Hsδ(Ω),and u=0 in ΩI

is well-posed and its solution is such that , for a positive constant .

The proof is simply based on application of (27) to the perturbed bilinear form, and is hence omitted. Note that (28) also allows us to apply [Felsinger2015, Theorem 4.1], which states that the bilinear form satisfies a weak maximum principle, which we report for completeness.

###### Lemma 5

Let and be defined as in (7) and (19) and let satisfy

 A(u,v)≤0∀v∈˜Hsδ(Ω),

then, . In particular,

 A(u,v)=0∀v∈˜Hsδ(Ω)⟺u=0. (29)

An immediate consequence of Lemma 4 is the fact that the operator mapping , is a compact operator. This fact allows us to use the Fredholm alternative theorem to prove the well-posedness of problem (18).

###### Theorem 1

Given a function , the problem

 Find u∈˜Hsδ(Ω) such thatA(u,v)=⟨f,v⟩,∀v∈˜Hsδ(Ω),and u=0 in ΩI

has a unique solution .

###### Proof

Consider the operator . We have the following identity:

 K−ηI=−ηK(K−1−η−1I)=−ηK(L+(c−η−1)I).

By the Fredholm alternative, either

is an eigenvalue of

, or is a bounded linear operator. We show that the first option leads to a contradiction.

Let and assume that is an eigenvalue of . This implies that for all and, by (29), that . Hence, is not an eigenvalue of . This ends the proof since we can conclude that , and subsequently , are bounded linear operators. In other words, the variational problem associated with the operator is well-posed.

## 4 Finite element discretization and efficient matrix assembly

In this section we briefly introduce the finite element (FE) discretization of problem (18) and provide details regarding numerical integration.

We choose a family of finite-dimensional subspaces , parameterized by the mesh size . We define the Galerkin approximation as the solution of the following problem: find such that

 A(uh,vh)=∫Ωfvhdx,∀vh∈Vh. (30)

As in [Du2012], we consider finite element approximations for the case that both and are polyhedral domains. We partition into finite elements and denote by the diameter of the largest element in the partition. We assume that the partition is shape-regular and quasi-uniform [Brenner2008] as the grid size and we choose the subspace to consist of piecewise polynomials of degree no more than defined with respect to the partition.

We perform numerical integration on pairs of elements, by rewriting the bilinear form (20) as the sum of integrals over the partition, i.e.

 A(uh,vh)=12∑K∑˜K∬K×˜K(uh(x)−uh(y))(γ(x,y)vh(x)−γ(y,x)vh(y))dydx.

The design of quadrature rules for the integrals above depends on the regularity properties of the functions and . If and are constant on a pair of elements , we can use the state-of-the-art quadrature rules described in [AinsworthGlusa2018, AinsworthGlusa2017] designed for the constant-order fractional Laplacian, i.e. for the case . However, when this is not the case, the use of more sophisticated and often more expensive quadrature rules, such as adaptive quadrature rules [PiessensDoncker-KapengaEtAl2012_Quadpack], becomes mandatory to prevent the integration error to be dominant. Conditions for well-posedness derived in Proposition 1 are such that any non-symmetric kernel falls into the latter category, hence requiring higher computational effort. In addition, note that, as for the constant-order fractional Laplacian, when , the integrand function features a singularity at . In this case, the domain of integration must be partitioned into subdomains in order to avoid the singularity.

In the specific case of piecewise constant fractional order, we speed up the assembly of the discrete operator and subsequent matrix-vector multiplications by generalizing the panel-clustering approach of

[AinsworthGlusa2018, AinsworthGlusa2017]

developed for the constant-order fractional Laplacian. Informally speaking, in the standard version of this approach, the unknowns are recursively grouped into nodes of a tree structure, with the root node containing all unknowns. Pairs of clusters correspond to sub-blocks of the discrete operator. If certain conditions on their location in physical space are satisfied, sub-blocks are approximated using Chebyshev interpolation, otherwise they are assembled using quadrature

aaaSee the works by Ainsworth and Glusa[AinsworthGlusa2018, AinsworthGlusa2017] for a detailed description of this approach for constant parameters and infinite horizon.. We generalize this algorithm to finite horizon and piecewise constant coefficients by restricting the approximation of matrix sub-blocks to cluster pairs that are strictly contained 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 of available fractional kernels, we restrict ourselves to a one-dimensional problem where model parameters are constant over two subdomains. This configuration can be considered a one-dimensional counterpart of the two-dimensional configurations displayed in Figure 1. Specifically, we let , , and consider kernels given by (8) that are piecewise constant with respect to the partition . In order for the assumptions of Proposition 1 to be satisfied, such kernels need to be symmetric, i.e. and need to take the form

 ψ(x,y;ψ+,ψ−,ψ±) =⎧⎨⎩ψ+if x,y>0,ψ−if x,y<0,ψ±else.

Note that this excludes the interesting cases, that can be found in the literature, of piecewise constant or that only depend on . Hence, we also consider functions of the following form and compare them with the proposed model.

 η(x;η+,η−,r,κ) =⎧⎪⎨⎪⎩η+if x≥r12(η++η−)+12(η+−η−)arctanκxarctanκrif −r≤x≤rη−if x≤−r.

Here, determines the size of the transition region between the two states and , and the steepness of the transition. Since is Lipschitz continuous, we can use for both and and satisfy the conditions of Proposition 1.

We numerically solve the diffusion problem (6) for homogeneous Dirichlet data and forcing term . In all tests conducted in this section we use continuous piecewise linear FE on a uniform mesh of size .

With the purpose of highlighting how a variable fractional order affects the qualitative behavior of the solution across the interface, we first compare the family of symmetric kernels reported in Table 1 and display the corresponding solutions in Figure 2. We observe that the solutions for constant seem to be differentiable, whereas the cases with piecewise defined are only continuous at .