Finite element approximation of the parabolic fractional obstacle problemEO has been supported in part by NSF grant DMS-1411808 and by CONICYT through project Anillo ACT1106. AJS is partially supported by NSF grant DMS-1418784.

# Finite element approximation of the parabolic fractional obstacle problem1

## Abstract

We study a discretization technique for the parabolic fractional obstacle problem in bounded domains. The fractional Laplacian is realized as the Dirichlet-to-Neumann map for a nonuniformly elliptic equation posed on a semi-infinite cylinder, which recasts our problem as a quasi-stationary elliptic variational inequality with a dynamic boundary condition. The rapid decay of the solution suggests a truncation that is suitable for numerical approximation. We discretize the truncation with a backward Euler scheme in time and, for space, we use first-degree tensor product finite elements. We present an error analysis based on different smoothness assumptions.

o

bstacle problem, thin obstacles, free boundaries, finite elements, fractional diffusion, anisotropic elements.

{AMS}

35J70, 35R11, 35R35, 49M15, 49M20, 49M25, 65N12, 65N30, 65N50

## 1 Introduction

In this work we shall be interested in the design and analysis of a finite element approximation of the so-called parabolic fractional obstacle problem. Let be an open and bounded subset of with . Given , an obstacle , an initial datum and a forcing term , the parabolic fractional obstacle problem asks for a function that satisfies the evolution variational inequality:

 min{\rm dt\textupu+(−Δ)s\textupu−\textupf,\textupu−ψ}=0 (1.1)

and . Here, denotes the fractional powers of the Laplace operator, supplemented with homogeneous Dirichlet boundary conditions, which for convenience we will simply call the fractional Laplacian. We refer to §2.1 for a precise definition. We must immediately remark that although our exposition is for the fractional Laplacian, our techniques and results are equally applicable to fractional powers of a symmetric and uniformly elliptic second order differential operator , supplemented with homogeneous Dirichlet boundary conditions: , with symmetric and positive definite and . The only caveat is that, at the time of this writing, no regularity results are available for (1.1) with the fractional Laplacian replaced by .

The study of numerical techniques for nonlocal problems is a rapidly growing field of research. Fractional diffusion has received a great deal of attention in diverse areas of science and engineering such as mechanics [atanackovic2014fractional], biophysics [bio], turbulence [wow], image processing [GH:14], peridynamics [HB:10] and nonlocal electrostatics [ICH]. In particular, the study of constrained minimization problems such as the parabolic fractional obstacle problem (1.1) has received considerable attention. This type of problems arises, for instance, in financial mathematics as a pricing model for American options. The function u represents the rational price of a non perpetual American option where the assets’ prices are modeled by a Lévy process, and the payoff function is ; see [MR2064019, NTZ:10, MR1424787].

Based on the Caffarelli-Silvestre extension [CS:07] in previous work we provided a comprehensive analysis of the discretization of the linear elliptic case [NOS], evolution equations with fractional diffusion and Caputo fractional time derivative [NOS3] and the elliptic fractional obstacle problem [NOS4]. In this work we proceed in our research program and show the flexibility of the ideas developed in [NOS] by studying the parabolic fractional obstacle problem (1.1). To the best of our knowledge, this is the first work that addresses the numerical approximation of this problem.

Our presentation is organized as follows. The notation and functional setting is described in Section 2, where we also briefly describe, in §2.1, the definition of the fractional Laplacian, its localization via the Caffarelli-Silvestre extension (§2.2) and the well-posedness of the fractional parabolic obstacle problem; see §2.3. The numerical analysis of problem (1.1) begins in Section 3, where we discuss a domain truncation that allows us, in subsequent sections, to consider a space discretization using first-degree tensor product finite elements. The time discretization and its error analysis is described in Section 4. The space discretization and its analysis is the content of Section 5: we provide an error analysis with minimal (§5.2) and maximal (§5.3) regularities. This analysis relies on the construction and approximation properties of a positivity preserving interpolant. For we construct an interpolant with the requisite properties in Section 6.

## 2 Notation and preliminaries

In this work is a convex bounded and open subset of () with polyhedral boundary. Our ideas are equally applicable to domains with curved boundaries, but the exposition becomes rather cumbersome and so we prefer to avoid it. We will follow the notation of [NOS] and define the semi-infinite cylinder and its lateral boundary by , . For we define the truncated cylinder and its lateral boundary . We also define the Dirichlet boundary . Since we will be dealing with objects defined on and , it will be convenient to distinguish the -dimension. For , we denote

 x=(x1,⋯,xd,xd+1)=(x′,xd+1)=(x′,y),x′∈Rd,y∈R.

Whenever is a normed space we denote by its norm and by its dual. For normed spaces and we write to indicate continuous embedding. We will follow standard notation for function spaces [Adams, Tartar]. In addition, for an open set , , if is a weight and we denote the Lebesgue space of -integrable functions with respect to the measure by ; see [HKM, Kufner80, Turesson]. Similar notation will be used for weighted Sobolev spaces. If and , we consider as a function of with values in a Banach space , . For we will say that if the mapping is in . We also introduce the space of -valued functions of bounded variation [Brezis, Definition A.2]

 VarXg:=supP{J∑j=1∥g(rj)−g(rj−1)∥X}<∞,

where the supremum is taken over all partitions of the time interval . We recall that if , then at every point there exists the right limit [Brezis, Lemma A.1].

Let be the number of time steps. We define the uniform time step as and we set , . Given a function , we denote and . For any sequence , we define the piecewise constant interpolant by

 ¯Wτ(t)=Wk+1t∈(tk,tk+1],k=0,…,K−1.

We also define the piecewise linear interpolant by

 ^Wτ(t)=t−tkτWk+1+tk+1−tτWkt∈[tk,tk+1],k=0,…,K−1.

The first order backward difference operator is defined by . We note that for all and . Finally, we also notice that, for any sequence and we have

and .

The relation indicates that for a constant that does not depend on either or , but it might depend on the problem data. The value of might change at each occurrence.

### 2.1 The fractional Laplacian

For a bounded domain there are several ways, not necessarily equivalent, to define the fractional Laplacian; see [NOS] for a discussion. As in [NOS] we will adopt that based on spectral theory [BS]. Namely, since is an unbounded, positive and closed operator with dense domain and its inverse is compact, there is a countable collection of eigenpairs such that is an orthonormal basis of and an orthogonal basis of . If

 w=∑l∈Nwlφl,wl=∫Ωwφl\rm dx′,

then, for any , we define , As it is well known, the theory of Hilbert scales presented in [Lions, Chapter 1] shows that , i.e., the real interpolation between and . Consequently, the definition of can be extended by density to the space . If, for , we denote by the dual space of , then is an isomorphism.

### 2.2 The Caffarelli-Silvestre extension problem

The Caffarelli-Silvestre result [CS:07, NOS], requires us to deal with a nonuniformly elliptic equation. With this in mind, we define the weighted Sobolev space

 \raisebox{6.9pt}{\tiny∘} \kern-10.3ptH1L(yα,C)={w∈H1(yα,C):w=0 on ∂LC}.

Since , belongs to the Muckenhoupt class ; see [Javier, Turesson]. Consequently, is a Hilbert space, and smooth functions are dense in (cf. [Turesson, Proposition 2.1.2, Corollary 2.1.6]).

As [NOS, (2.21)] shows, the following weighted Poincaré inequality holds:

 ∥w∥L2(yα,C)≲∥∇v∥L2(yα,C),∀w∈ \raisebox{6.9pt}{\tiny∘} \kern-1% 0.3ptH1L(yα,C). (2.1)

Then, the seminorm of is equivalent to the norm in . For we denote by the trace of onto . We recall ([NOS, Prop. 2.5])

 trΩ \raisebox{6.9pt}{\tiny∘} \kern-10.3ptH1L(yα,C)=Hs(Ω),∥trΩw∥Hs(Ω)≤CtrΩ∥w∥\,\raisebox{4.7pt}{\tiny∘} % \kern-9.3ptH1L(yα,C). (2.2)

The seminal work of Caffarelli and Silvestre [CS:07, NOS] showed that the operator can be realized as the Dirichlet-to-Neumann map for a nonuniformly elliptic boundary value problem. Namely, if solves

 −∇⋅(yα∇U)=0 in C,U=0 on ∂LC,∂ανU=\textupdsf on Ω×{0}, (2.3)

where , and is a normalization constant, then solves

 (−Δ)su=f. (2.4)

The reader is referred to [NOS] for a detailed and thorough exposition on how the groundbreaking identity given by (2.3) can be used to design and analyze an efficient finite element approximation of solutions to (2.4).

### 2.3 The parabolic fractional obstacle problem

Given an obstacle that satisfies and on , let

 \mathpzcK(Ω)={w∈Hs(Ω):w(x′)≥ψ(x′) a.e. x′∈Ω}

be the convex set of admissible functions, and let be its indicator function, which, since is closed and convex, is a non-smooth lower-semi-continuous convex function. We consider the energy

 J(ϕ)=12∥ϕ∥2Hs(Ω)+\textupInd\mathpzcK(Ω)(ϕ).

With this notation, the parabolic fractional obstacle problem can be understood as the gradient flow for , or an evolution equation for a maximal monotone operator: Given an initial datum and a function , find u such that and it solves the differential inclusion

 \rm dt\textupu(t)+∂J(\textupu(t))∋\textupf(t)a.e. t∈(0,T). (2.5)

Problem (2.5) can be equivalently understood as an evolution variational inequality: Find u such that and for a.e.  and all

 (\rm dt\textupu(t),\textupu(t)−ϕ)L2(Ω)+⟨(−Δ)s\textupu(t),\textupu(t)−ϕ⟩≤(\textupf(t),\textupu(t)−ϕ)L2(Ω), (2.6)

and . Here and in what follows denotes the inner product of and the duality pairing between and .

From these formulations existence and uniqueness of solutions and a priori estimates can be easily obtained with standard techniques on maximal monotone operators [Brezis]. For instance, if u solves (2.6), then it satisfies the energy estimate

 ∥\textupu∥2L∞(0,T;L2(Ω))+∥\textupu∥2L2(0,T;Hs(Ω))≲D2,

where we denoted

 D2=D2(\textupu0,\textupf,ψ)=∥\textupu0∥2L2(Ω)+∥\textupf∥2L2(0,T;L2(Ω))+∥ψ∥2Hs(Ω). (2.7)

If and , then there exists a unique strong solution, that is which is locally absolutely continuous in and satisfies (2.6) at almost every point. In addition, we have that and the mapping is locally absolutely continuous in , which implies that , and that the following estimate holds [Brezis, Theorem 3.6]

 ∥\rm dt\textupu(t)∥2L2(Ω)+\rm dt∥\textupu(t)∥2Hs(Ω)=(\textupf,\rm d% t\textupu(t))L2(Ω)a.e. t∈(0,T).

If, moreover, and , then [Brezis, Proposition 3.3]. Finally, if we have

 12\rm dt∥\rm dt\textupu(t)∥2L2(Ω)≤(\textupft,\rm dt\textupu(t))L2(Ω),

in the distributional sense.

Let us now use the Caffarelli–Silvestre extension detailed in §2.1 to write an obstacle problem that is equivalent to (2.6). To do this, we define the set

 \mathpzcK(C)={w∈ \raisebox{6.9pt}{\tiny∘} % \kern-10.3ptH1L(yα,C):trΩw(x′)≥ψ(x′) a.e. x′∈Ω}.

Problem (2.6) can then be equivalently stated as: Find such that for a.e.  and every

 (trΩ\rm dtU(t),trΩ(U(t)−ϕ))L2(Ω)+a(U(t),U(t)−ϕ)≤(\textupf,trΩ(U(t)−ϕ))L2(Ω), (2.8)

with . Here the bilinear form is defined by

 a(w,ϕ)=1\textupds∫Cyα∇w∇ϕ\rm dx′\rm dy,∀w,ϕ∈ % \raisebox{6.9pt}{\tiny∘} \kern-10.3ptH1L(yα,C).

The description of the functional setting for problem (2.8), together with existence and uniqueness results, follow the analysis developed for problem (2.5); see [Brezis, NSV]. In particular, we have the energy inequality

 ∥trΩU∥2L∞(0,T;L2(Ω))+∥U∥2L2(0,T;\,\raisebox{4.7pt}{\tiny∘}% \kern-9.3ptH1L(yα,C))≲D2, (2.9)

where is defined in (2.7). We shall be more specific on the smoothness of the data and the consequences on the regularity of the solution when we perform the discretization and its analysis. Let us now contempt ourselves with mentioning that, provided is sufficiently smooth, the following complementarity system holds:

 Z:=∂ανU+\textupdstrΩ\rm dtU−\textupds\textupf≥0,trΩU−ψ≥0,Z(trΩU−ψ)=0.

## 3 Truncation

The variational inequality (2.8) is posed on a infinite domain and, consequently, it cannot be directly approximated with finite element-like techniques. A first step towards the discretization is to truncate the domain to a bounded cylinder and study the effect of this truncation. We begin with a result that shows the exponential decay of the solution to (2.8); compare with [NOS, Proposition 3.1], [NOS4, Lemma 4.8] and [NOS3, Proposition 4.1].

{lemma}

[exponential decay] If , and , then, for every , we have

 ∥∇U∥L2(0,T;L2(yα,Ω×(\mathpzcY,∞))≲e−\mathpzcY/2D,

where the hidden constant does not depend on neither nor the problem data. {proof} Consider, for a.e. , the function that solves

 ∇⋅(yα∇w(t))=0 in C,w(t)|∂LC=0,trΩw(t)=trΩU(t) on Ω×{0}. (3.1)

Since solves the fractional parabolic obstacle problem (2.8) and problem (3.1) has a unique solution, we immediately conclude that for a.e. , we have . We now apply the decay estimate of [NOS, Proposition 4.1] to problem (3.1) to obtain

 ∥∇w(t)∥L2(yα,Ω×(\mathpzcY,∞))≲e−\mathpzcY/2∥trΩU(t)∥Hs(Ω).

Finally, integrating over time, and invoking the trace estimate (2.2) and the stability estimate (2.9) for problem (2.8) in terms of , and , we arrive at

 ∥∇w∥L2(0,T;L2(yα,Ω×(\mathpzcY,∞))≲e−\mathpzcY/2∥U∥L2(0,T;\,\raisebox{4.% 7pt}{\tiny∘} \kern-9.3ptH1L(yα,C\mathpzcY))≲e−\mathpzcY/2D,

where is defined in (2.7). This concludes the proof.

The exponential decay of Lemma 3 allows us to consider a truncated version of the variational inequality (2.8). To write this problem we define, for , the Sobolev space

 \raisebox{6.9pt}{\tiny∘} \kern-10.3ptH1L(yα,C\mathpzcY)={w∈H1(yα,C\mathpzcY):w|ΓD=0},

the convex set of admissible functions

 \mathpzcK(C\mathpzcY)={w∈ \raisebox{6.9pt}{% \tiny∘} \kern-10.3ptH1L(yα,C\mathpzcY):trΩw(x′)≥ψ(x′) a.e. x′∈Ω}

and the bilinear form

 a\mathpzcY(w,ϕ)=1\textupds∫C\mathpzcYyα∇w∇ϕ\rm dx′\rm d% y∀w,ϕ∈ \raisebox{6.9pt}{\tiny∘} \kern-10.3ptH1L(yα,C\mathpzcY).

With these definitions we consider the following truncated problem: Find such that and for a.e.  and all

 (trΩ\rm dtv(t),trΩ(v(t)−ϕ))L2(Ω)+a\mathpzcY(v(t),v(t)−ϕ)≤(\textupf,trΩ(v(t)−ϕ))L2(Ω). (3.2)

The analysis of this problem follows that of (2.5), developed in §2.3. For brevity, we only present the energy estimate

 ∥trΩv∥2L∞(0,T;L2(Ω))+∥v∥2L2(0,T;\,\raisebox{4.7pt}{\tiny∘} \kern-9.3ptH1L(yα,C\mathpzcY))≲D2. (3.3)

We define , the truncated -harmonic extension operator as follows: if , then solves

 ∇⋅(yα∇W)=0 in C\mathpzcY,W=0 on ∂LC\mathpzcY∪Ω×{\mathpzcY},W=w on Ω×{0}. (3.4)

We recall [NOS, Theorem 2.7] for problem (3.4). If , then

 ∥∇∇x′W∥L2(yα,C\mathpzcY)+∥∂yyW∥L2(yβ,C\mathpzcY)≲∥w∥H1+s(Ω). (3.5)

The following result shows that by considering (3.2) instead of (2.8) we only incur in an exponentially small error; compare with [NOS, Lemma 3.3], [NOS3, Lemma 4.3] and [NOS4, Proposition 4.20].

{proposition}

[exponential error estimate] Let and be the solutions of (2.8) and (3.2), respectively. Then, for , we have

 ∥trΩ(U−v)∥2L∞(0,T;L2(Ω))+∥∇(U−v)∥2L2(0,T;L2(yα,C\mathpzcY))≲e−\mathpzcY/4D2,

where the hidden constant does not depend on neither , , nor the problem data. {proof} By a trivial zero extension we realize that the solution to problem (3.2) belongs to , then we can set in (2.8). We would like to set in (3.2) but, although it satisfies the constraints, it is not an admissible test function, as it does not have a vanishing trace at . For this reason, instead, we set in (3.2), where is the following smooth cutoff function:

 ρ(y)=10≤y≤\mathpzcY2,ρ(y)=2\mathpzcY(\mathpzcY−y)\mathpzcY2

With these choices of test functions we add the ensuing inequalities to obtain

 12\rm dt∥trΩ(U−v)∥2L2(Ω)+a(U,U−v)≤a\mathpzcY(v,ρU−v).

We now notice that

 a(U,U−v)=a\mathpzcY(U,U−v)+∥∇U∥2L2(yα,Ω×(\mathpzcY,∞))≥a\mathpzcY(U,U−v),

so that we obtain

 12\rm dt∥trΩ(U−v)∥2L2(Ω)+a\mathpzcY(U−v,U−v)≤a\mathpzcY(v,(ρ−1)U). (3.6)

It remains then to bound the right hand side of (3.6). A straightforward computation reveals that if we have that , otherwise

 |∇(ρ−1)U|2≤2(4\mathpzcY2U2+|∇U|2),

and thus

 ∥∇(ρ−1)U∥2L2(yα,C\mathpzcY)≤2(4\mathpzcY2∫\mathpzcY\mathpzcY/2∫Ωyα|U|2\rm dx′% \rm dy+∫\mathpzcY\mathpzcY/2∫Ωyα|∇U|2\rm dx′\rm dy).

Invoking a version of the Poincaré inequality (2.1) based on the interval , we conclude We now use this estimate in (3.6) and integrate in time to obtain

 ∥trΩ(U−v)(t)∥2L2(Ω)+∫t0∥∇(U−v)(s)∥2L2(yα,C\mathpzcY)\rm ds≤∥∇v∥L2(0,T;L2(yα,C\mathpzcY))∥∇U∥L2(0,T;L2(yα,Ω×(\mathpzcY/2,∞)))≲e−\mathpzcY/4D2,

where the last inequality follows from Lemma 3, the fact that and the stability estimate (3.3) of (3.2) in terms of . Since is arbitrary this implies the desired estimate.

## 4 Time discretization

We now proceed with the time discretization of problem (1.1). We could directly apply a suitable time discretization scheme to either (2.5), (2.6) or (2.8) and then argue, for instance, by using the results of [NSV]. However, with the exponential convergence result of Proposition 3 at hand we will consider the time discretization of the truncated problem (3.2) based on the implicit Euler method.

The discrete scheme computes the sequence , an approximation to the solution to problem (3.2) at each time step. We initialize the scheme by setting

 trΩV0=\textupu0, (4.1)

and for , let be such that, for every ,

 (trΩdVk+1τ,trΩ(Vk+1−ϕ))L2(Ω) +a\mathpzcY(Vk+1,Vk+1−ϕ) (4.2) ≤(\textupfk+1,trΩ(Vk+1−ϕ))L2(Ω),

where is defined in Section 2 and .

Existence and uniqueness of a solution to (4.2) follows directly from standard arguments on variational inequalities [Brezis, KS]. The approximate solution to problem (1.1) is then defined by the sequence where

 Uτ=trΩVτ. (4.3)
###### Remark \thetheorem (locality)

The main advantage of scheme (4.1)–(4.2) is its local nature, which mimics that of problem (2.8).

Let us now show the stability of the scheme.

{proposition}

[stability] Assume that and , then and uniformly in . {proof} Set in (4.2), where is the -harmonic extension operator introduced in (3.4). Upon denoting we obtain

 (trΩdWk+1τ,trΩWk+1)L2(Ω)+a\mathpzcY(Wk+1,Wk+1)≤(fk+1,trΩWk+1)L2(Ω)+a\mathpzcY(Hαψ,Wk+1).

The Cauchy Schwartz inequality and summation over yields the result.

The error analysis of (4.1)–(4.2) follows from the general theory presented in [NSV]. To present it we introduce the error

 E(v,Vτ)=∥trΩ(v−^Vτ)∥L∞(0,T;L2(Ω))+∥∇(v−¯Vτ)∥L2(0,T;L2(yα,C\mathpzcY)), (4.4)

where and are defined in Section 2. We also define

 E(\textupu,Uτ)=∥\textupu−^Uτ∥L∞(0,T;L2(Ω))+∥\textupu−¯Uτ∥L2(0,T;Hs(Ω)). (4.5)
{corollary}

[error estimates in time I] If and , then the solutions of (3.2) and of (4.1)–(4.2) satisfy the uniform estimate

 E(v,Vτ)≲τ1/2(∥\textupu0∥Hs(Ω)+∥\textupf∥L2(0,T;L