A PDE Approach to Numerical Fractional Diffusion Enrique Ot arola - PowerPoint PPT Presentation

Muckenhoupt weights There is a constant C such that for every a , b ∈ R , with a > b , � b � b 1 1 | y | α d y · | y | − α d y ≤ C b − a b − a a a which means y α belongs to the Muckenhoupt class A 2 . Then ◮ The Hardy-Littlewood maximal operator is continuous on L 2 ( y α , C ). ◮ Singular integral operators are continuous on L 2 ( y α , C ). ◮ L 2 ( y α , C ) ֒ → L 1 loc ( C ). ◮ H 1 ( y α , C ) is Hilbert and C ∞ b ( C ) is dense. ◮ Traces on ∂ L C are well defined.

Weighted Sobolev spaces ◮ Weighted Poincar´ e inequailty: There is a constant C , s.t. � � y α | w | 2 ≤ C y α |∇ w | 2 H 1 ◦ L ( y α , C ) . ∀ w ∈ C C ◦ ◮ Surjective trace operator tr Ω : H 1 L ( y α , C ) → H s (Ω). ◮ Lax-Milgram = ⇒ existence and uniqueness for every f ∈ H − s (Ω). Also �U� ◦ L ( y α , C ) = � u � H s (Ω) = d s � f � H − s (Ω) . H 1 We will discretize the α -harmonic extension! · ( y α ∇U ) = 0  ∇ in C   ◦ H 1 L ( y α , C ) : U ∈ U = 0 on ∂ L C  ∂ ν α U = d s f on Ω × { 0 } 

Weighted Sobolev spaces ◮ Weighted Poincar´ e inequailty: There is a constant C , s.t. � � y α | w | 2 ≤ C y α |∇ w | 2 H 1 ◦ L ( y α , C ) . ∀ w ∈ C C ◦ ◮ Surjective trace operator tr Ω : H 1 L ( y α , C ) → H s (Ω). ◮ Lax-Milgram = ⇒ existence and uniqueness for every f ∈ H − s (Ω). Also �U� ◦ L ( y α , C ) = � u � H s (Ω) = d s � f � H − s (Ω) . H 1 We will discretize the α -harmonic extension! · ( y α ∇U ) = 0  ∇ in C   ◦ H 1 L ( y α , C ) : U ∈ U = 0 on ∂ L C  ∂ ν α U = d s f on Ω × { 0 } 

Separation of Variables • Apply separation of variables (Capella et al. 2011) ∞ ∞ u ( x ′ ) = � u k ϕ k ( x ′ ) = ⇒ U ( x ′ , y ) = � u k ϕ k ( x ′ ) ψ k ( y ) , k =1 k =1 where the functions ψ k solve k + α  ψ ′′ y ψ ′ k − λ k ψ k = 0 , in (0 , ∞ ) ,   ψ k (0) = 1 , y →∞ ψ k ( y ) = 0 . lim   • It turns out that � s �� � ψ k ( y ) = c s λ k y K s ( λ k y ) , K s – modified Bessel function of the second kind. • Since ψ ′ k ( y ) ≈ y − α , U �∈ H 1 ( C ) . y ↓ 0 = ⇒

Outline Motivation The elliptic linear problem case Formulation The Caffarelli-Silvestre extension Discretization Interpolation estimates in weighted spaces Regularity and a priori error estimates Numerical Experiments A posteriori error analysis and adaptivity Space-time fractional parabolic problem The fractional obstacle problem An optimal control problem

Domain truncation The domain C is infinite. We need to consider a truncated problem. Theorem (exponential decay) For every Y > 0 L ( y α , Ω × ( Y , ∞ )) � e −√ λ 1 Y / 2 � f � H − s (Ω) . �U� ◦ H 1 Let v solve · ( y α ∇ v ) = 0  ∇ in C Y = Ω × (0 , Y ) ,   v = 0 on ∂ L C Y ∪ Ω × { Y } ,  ∂ ν α v = d s f on Ω × { 0 } .  Theorem (exponential convergence) For all Y > 0 , L ( y α , C Y ) � e −√ λ 1 Y / 4 � f � H − s (Ω) . �U − v � ◦ H 1

Finite element method I ◮ Continuous solution. V -Hilbert space. Find u s.t. B [ u , v ] = F [ v ] , ∀ v ∈ V . B -continuous and coercive bilinear form, and F -continuous linear functional. ◮ Approximate solution. Let V N be a finite dimensional space. Find U N s.t. B [ U N , V N ] = F [ V N ] , ∀ v ∈ V N .

Finite element method I ◮ Continuous solution. V -Hilbert space. Find u s.t. B [ u , v ] = F [ v ] , ∀ v ∈ V . B -continuous and coercive bilinear form, and F -continuous linear functional. ◮ Approximate solution. Let V N be a finite dimensional space. Find U N s.t. B [ U N , V N ] = F [ V N ] , ∀ v ∈ V N .

Finite element method I ◮ Continuous solution. V -Hilbert space. Find u s.t. B [ u , v ] = F [ v ] , ∀ v ∈ V . B -continuous and coercive bilinear form, and F -continuous linear functional. ◮ Approximate solution. Let V N be a finite dimensional space. Find U N s.t. B [ U N , V N ] = F [ V N ] , ∀ v ∈ V N .

Finite element method II ◮ Continuous solution. V -Hilbert space. Find u s.t. B [ u , v ] = F [ v ] , ∀ v ∈ V . B -continuous and coercive bilinear form, and F -continuous linear functional. ◮ Approximate solution. Let V N be a finite dimensional space. Find U N s.t. B [ U N , V N ] = F [ V N ] , ∀ v ∈ V N . ◮ Error estimates. ◮ A priori. Convergence, a rate of convergence, and know the depende of the error on different factors. Typical estimate: � u − U N � V � N − a � u � ≈ h b � u � ◮ A posteriori. Information beyond asymptotics; computable in terms of F and U N . Quality assessment; adaptivity.

Galerkin method: mesh Let T Ω = { K } be triangulation of Ω (simplices or cubes) ◮ T Ω is conforming and shape regular. Let T Y = { T } be a triangulation of C Y into cells of the form T = K × I , K ∈ T Ω , I = ( a , b ) . Why? Natural on the cylinder C Y , deal.ii, and U yy ≈ y − α − 1 as y ≈ 0+ Approximation of singular functions = ⇒ anisotropic elements Shape regularity condition does NOT hold!

Galerkin method: discrete spaces We only require that if T = K × I and T ′ = K ′ × I ′ are neighbors | I | | I ′ | ≃ 1 , so the lengths of I and I ′ are comparable. This weak condition allows us to consider anisotropic meshes Define: W ∈ C 0 ( C Y ) : W | T ∈ P 1 ( T ) , � � V ( T Y ) = W | Γ D = 0 with Γ D = ∂ L C ∪ Ω × { Y } , and W ∈ C 0 (¯ � � U ( T Ω ) = tr Ω V ( T Y ) = Ω) : W | K ∈ P 1 ( K ) , W ∂ Ω = 0

Galerkin method: discrete problem Galerkin method for the extension: Find V T Y ∈ V ( T Y ) such that � y α ∇ V T Y ∇ W = d s � f , tr Ω W � H − s (Ω) , H s (Ω) , ∀ W ∈ V ( T Y ) C Y Define: U T Ω = tr Ω V T Y ∈ U ( T Ω ) A trace estimate and C` ea’s Lemma imply quasi-best approximation: � u − U T Ω � H s (Ω) � � v − V T Y � ◦ H 1 L ( y α , C Y ) = W ∈ V ( T Y ) � v − W � ◦ inf H 1 L ( y α , C Y ) Setting W = Π v ∈ V ( T Y ) = ⇒ interpolation analysis!

Outline Motivation The elliptic linear problem case Formulation The Caffarelli-Silvestre extension Discretization Interpolation estimates in weighted spaces Regularity and a priori error estimates Numerical Experiments A posteriori error analysis and adaptivity Space-time fractional parabolic problem The fractional obstacle problem An optimal control problem

The averaged Taylor polynomial Consider ω ∈ A 2 ( R N ) and φ ∈ L 2 ( ω, D ), with D ⊂ R N . Given a node z of the mesh, we define Given m ∈ N , we define 1 � � α ! D α φ ( x )( y − x ) α ψ z ( x ) d x . Q m z φ ( y ) = | α |≤ m A weighted Poincar´ e inequailty yields � φ − Q 0 z φ � L 2 ( ω, S z ) � diam ( S z ) �∇ φ � L 2 ( ω, S z ) , which, via an induction argument, allows us to derive | φ − Q m z φ | H k ( ω, S z ) � diam ( S z ) m +1 − k | φ | H m +1 ( ω, S z ) , k = 0 , 1 , ..., m +1 Extension to the weighted case! Simple argument!

The averaged Taylor polynomial Consider ω ∈ A 2 ( R N ) and φ ∈ L 2 ( ω, D ), with D ⊂ R N . Given a node z of the mesh, we define Given m ∈ N , we define 1 � � α ! D α φ ( x )( y − x ) α ψ z ( x ) d x . Q m z φ ( y ) = | α |≤ m A weighted Poincar´ e inequailty yields � φ − Q 0 z φ � L 2 ( ω, S z ) � diam ( S z ) �∇ φ � L 2 ( ω, S z ) , which, via an induction argument, allows us to derive | φ − Q m z φ | H k ( ω, S z ) � diam ( S z ) m +1 − k | φ | H m +1 ( ω, S z ) , k = 0 , 1 , ..., m +1 Extension to the weighted case! Simple argument!

The quasi-interpolant operator We introduce an averaged interpolation operator Π ´ a la Duran Lombardi, 2005 (Sobolev 1950, Dupont Scott 1980) Π φ ( z ) = Q m z φ ( z ) . Notice that: ◮ This is defined for all polynomial degree m and any element shape (simplices or rectangles). ◮ We do not go back to the reference element — This is important for anisotropic estimates. The mesh is rectangular and Cartesian. If R and S are neighbors h i R / h i S � 1 , i = 1 , N .

The quasi-interpolant operator We introduce an averaged interpolation operator Π ´ a la Duran Lombardi, 2005 (Sobolev 1950, Dupont Scott 1980) Π φ ( z ) = Q m z φ ( z ) . Notice that: ◮ This is defined for all polynomial degree m and any element shape (simplices or rectangles). ◮ We do not go back to the reference element — This is important for anisotropic estimates. The mesh is rectangular and Cartesian. If R and S are neighbors h i R / h i S � 1 , i = 1 , N .

The quasi-interpolant operator We introduce an averaged interpolation operator Π ´ a la Duran Lombardi, 2005 (Sobolev 1950, Dupont Scott 1980) Π φ ( z ) = Q m z φ ( z ) . Notice that: ◮ This is defined for all polynomial degree m and any element shape (simplices or rectangles). ◮ We do not go back to the reference element — This is important for anisotropic estimates. The mesh is rectangular and Cartesian. If R and S are neighbors h i R / h i S � 1 , i = 1 , N .

Error estimates on rectangles Theorem If φ ∈ H 1 ( ω, S R ) N � h i � φ − Π φ � L 2 ( ω, R ) � R � ∂ i φ � L 2 ( ω, S R ) . i =1 If φ ∈ H 2 ( ω, S R ) N � h i � ∂ j ( φ − Π φ ) � L 2 ( ω, R ) � R � ∂ i ∂ j φ � L 2 ( ω, S R ) , i =1 N � R h j h i � φ − Π φ � L 2 ( ω, R ) � R � ∂ i ∂ j φ � L 2 ( ω, S R ) . i , j =1

Error estimates on rectangles Theorem If ω ∈ A p ( R N ) , and φ ∈ ✘✘✘✘✘ H 1 ( ω, S R ) W 1 p ( ω, S R ) N � h i � φ − Π φ � L p ( ω, R ) � R � ∂ i φ � L p ( ω, S R ) . i =1 If φ ∈ ✘✘✘✘✘ H 2 ( ω, S R ) W 2 p ( ω, S R ) N � h i � ∂ j ( φ − Π φ ) � L p ( ω, R ) � R � ∂ i ∂ j φ � L p ( ω, S R ) , i =1 N � R h j h i � φ − Π φ � L p ( ω, R ) � R � ∂ i ∂ j φ � L p ( ω, S R ) . i , j =1 Estimates on simplicial elements, different metrics and applications in RHN, EO, AJS. Piecewise polynomial interpolation in Muckenhoupt weighted Sobolev spaces and applications, Numer. Math. (2015)

Error estimates on rectangles Theorem If ω ∈ A p ( R N ) , and φ ∈ ✘✘✘✘✘ H 1 ( ω, S R ) W 1 p ( ω, S R ) N � h i � φ − Π φ � L p ( ω, R ) � R � ∂ i φ � L p ( ω, S R ) . i =1 If φ ∈ ✘✘✘✘✘ H 2 ( ω, S R ) W 2 p ( ω, S R ) N � h i � ∂ j ( φ − Π φ ) � L p ( ω, R ) � R � ∂ i ∂ j φ � L p ( ω, S R ) , i =1 N � R h j h i � φ − Π φ � L p ( ω, R ) � R � ∂ i ∂ j φ � L p ( ω, S R ) . i , j =1 Estimates on simplicial elements, different metrics and applications in RHN, EO, AJS. Piecewise polynomial interpolation in Muckenhoupt weighted Sobolev spaces and applications, Numer. Math. (2015)

Outline Motivation The elliptic linear problem case Formulation The Caffarelli-Silvestre extension Discretization Interpolation estimates in weighted spaces Regularity and a priori error estimates Numerical Experiments A posteriori error analysis and adaptivity Space-time fractional parabolic problem The fractional obstacle problem An optimal control problem

Regularity of Extension U Using properties of Bessel functions we obtain ψ ′′ k ( y ) ≈ y − α − 1 , ∈ H 2 ( C , y α ) . y ↓ 0 = ⇒ U / But Theorem (regularity of the extension) If f ∈ H 1 − s (Ω) and Ω is C 1 , 1 or a convex polygon � ∆ x ′ U� 2 L 2 ( C , y α ) + � ∂ y ∇ x ′ U� 2 L 2 ( C , y α ) = d s � f � 2 H 1 − s (Ω) If β > 1 + 2 α � ∂ yy U� L 2 ( C , y β ) � � f � L 2 (Ω) Anisotropic estimates compensate singular behavior!

Error Estimates: Quasi-uniform Meshes On uniform meshes h T ≈ h K ≈ h I for all T ∈ T Y , then Theorem (error estimates) The following estimate holds for all ǫ > 0 �∇ ( v − V T Y ) � L 2 ( C Y , y α ) � h K � ∂ y ∇ x ′ v � L 2 ( C , y α ) + h s − ǫ � ∂ yy v � L 2 ( C , y β ) I � h s − ǫ � f � H 1 − s (Ω) Consequently, � u − U T Ω � H s (Ω) � h s − ǫ � f � H 1 − s (Ω) . • This is suboptimal in terms of order (only order s − ǫ ) • It cannot be improved as numerical experimentation reveals!

Numerical Experiment: Quasi-uniform Mesh Let Ω = (0 , 1) and f = π 2 s sin( π x ), then U = 2 1 − s π s sin( π x ′ ) y s K s ( π y ) Γ( s ) If s = 0 . 2, then −0.3 10 H 1 (y α ) error C (DOF’s) −0.1 −0.4 10 −0.5 10 Error −0.6 10 −0.7 10 −0.8 10 0 1 2 3 4 5 10 10 10 10 10 10 Degrees of Freedom (DOF’s) The energy error behaves like DOFS − 0 . 1 ≈ h 0 . 2 , as predicted!

Error Estimates: Graded Meshes We use the principle of error equilibration. Mesh on (0 , Y ) � j � γ y j = Y , j = 0 , M , γ > 1 M k ( y ) ≈ y − α − 1 = ψ ′′ ⇒ energy equidistribution for γ > 3 / (1 − α ). Theorem (error estimates) If f ∈ H 1 − s (Ω) and Y ≈ | log N | , � u − U T Ω � H s (Ω) = �∇ ( U − V T Y ) � L 2 ( C , y α ) 1 � | log N | s N − n +1 � f � H 1 − s (Ω) , or equivalently � u − U T Ω � H s (Ω) � | log N Ω | s N − 1 / n � u � H 1+ s (Ω) . Ω

Outline Motivation The elliptic linear problem case Formulation The Caffarelli-Silvestre extension Discretization Interpolation estimates in weighted spaces Regularity and a priori error estimates Numerical Experiments A posteriori error analysis and adaptivity Space-time fractional parabolic problem The fractional obstacle problem An optimal control problem

Numerical Experiments: Meshes for Circle and s = 0 . 3 Set Ω = D (0 , 1) ⊂ R 2 , Figure : Uniform mesh in x ′ and anisotropic mesh in y

Experimental Rates for Circle and s = 0 . 3 and s = 0 . 7 Set Ω = D (0 , 1) ⊂ R 2 , f = j 2 s 1 , 1 J 1 ( j 1 , 1 r )( A 1 , 1 cos( θ ) + B 1 , 1 sin( θ )). With graded meshes: DOFs − 1 / 3 s = 0 . 3 s = 0 . 7 0 10 Error −1 10 2 4 6 10 10 10 Degrees of Freedom (DOFs) The experimental convergence rate − 1 / 3 is optimal! RHN, EO, AJS. A PDE approach to fractional diffusion in general domains: a priori error analysis , Found. Comput. Math. (2014).

Outline Motivation The elliptic linear problem case Formulation The Caffarelli-Silvestre extension Discretization Interpolation estimates in weighted spaces Regularity and a priori error estimates Numerical Experiments A posteriori error analysis and adaptivity Space-time fractional parabolic problem The fractional obstacle problem An optimal control problem

Adaptivity Adaptivity is motivated by • Computational efficiency: extra n + 1-dimension. • The a priori theory requires ◮ regularity of the datum: f ∈ H 1 − s (Ω). ◮ regularity of the domain: Ω is C 1 , 1 or a convex polygon. • If one of these conditions is violated, the solution U may have singularities in Ω and exhibit fractional regularity. • Quasi-uniform refinement of Ω would not result in an efficient solution technique. • We need an adaptive loop.

An adaptive loop Our adaptive loop is almost standard: SOLVE → ESTIMATE → MARK → REFINE with ◮ SOLVE : Finds V T Y , the Galerkin solution. ◮ ESTIMATE : Compute E z ′ for every node z ′ ∈ Ω. ◮ MARK : For θ ∈ (0 , 1) choose a minimal subset of nodes M : � E 2 z ′ ≥ θ 2 E 2 T . z ′ ∈M ◮ REFINE : Given M : 1. ∀ z ′ ∈ M refine the cells K ∋ z ′ to get ˜ T Ω . 2. Create an anisotropic mesh { ˜ I } with M so that grading holds. T Y = ˜ ˜ T Ω × { ˜ 3. The refined mesh is I } .

Adaptivity ◮ One of the main ingredients of our adaptive loop is an a posteriori error estimator. Despite of what might be claimed, the theory of a posteriori error estimation on anisotropic discretizations is still in its infancy. We propose an error estimator based on solving local problem on stars:

An ideal error estimator Define w ∈ H 1 ( y α , C z ′ ) : w = 0 on ∂ C z ′ \ Ω × { 0 } � � W ( C z ′ ) = . For z ′ ∈ Ω a node, we define the ideal estimator ζ z ′ ∈ W ( C z ′ ): � � y α ∇ ζ z ′ ∇ ψ = d s � f , tr Ω ψ � H − s (Ω) × H s (Ω) − y α ∇ V ∇ ψ C z ′ C z ′ for all ψ ∈ W ( C z ′ ), and � 1 / 2 �� ˜ E 2 ˜ E T Y = , E z ′ = �∇ ζ z ′ � L 2 ( y α , C z ′ ) . z ′ z ′ Theorem (ideal estimator) We have �∇ e � L 2 ( y α , C Y ) � ˜ E T Y and, for every node z ′ ∈ Ω ˜ E z ′ ≤ �∇ e � L 2 ( y α , C z ′ ) .

An ideal error estimator Define w ∈ H 1 ( y α , C z ′ ) : w = 0 on ∂ C z ′ \ Ω × { 0 } � � W ( C z ′ ) = . For z ′ ∈ Ω a node, we define the ideal estimator ζ z ′ ∈ W ( C z ′ ): � � y α ∇ ζ z ′ ∇ ψ = d s � f , tr Ω ψ � H − s (Ω) × H s (Ω) − y α ∇ V ∇ ψ C z ′ C z ′ for all ψ ∈ W ( C z ′ ), and � 1 / 2 �� ˜ E 2 ˜ E T Y = , E z ′ = �∇ ζ z ′ � L 2 ( y α , C z ′ ) . z ′ z ′ Theorem (ideal estimator) We have �∇ e � L 2 ( y α , C Y ) � ˜ E T Y and, for every node z ′ ∈ Ω ˜ E z ′ ≤ �∇ e � L 2 ( y α , C z ′ ) .

Local Problems on Stars • Discretization based on P 2 : Discrete space W ( C z ′ ). Then, we define � E 2 y α |∇ W z ′ | 2 , E 2 � E 2 z ′ := T Ω := z ′ . C z ′ z ′ 1 � • Define the data oscillation. If f z ′ | K = K f then | K | osc T Ω ( f ) 2 = osc z ′ ( f ) 2 = d s h 2 s � osc z ′ ( f ) 2 , z ′ � f − f z ′ � 2 L 2 ( S z ′ ) z ′ Theorem (computable estimator) E 2 � �∇ ( v − V T Y ) � 2 L 2 ( y α , C Y ) � E 2 + osc ( y α , V T Y , f , C z ′ ) 2 • Is this enough for convergence and optimality?

L-shaped domain with incompatible data ◮ Ω is the standard L-shaped domain. f = 1. For s < 1 2 the data is not compatible with the problem. ◮ The nature of the singularity of the solution is not known for this problem. DOFs − 1 DOFs − 1 0 3 0 3 10 10 s = 0 . 2 s = 0 . 2 s = 0 . 4 s = 0 . 4 s = 0 . 6 s = 0 . 6 s = 0 . 8 s = 0 . 8 Estimator Error −1 10 −1 10 2 3 4 2 3 4 10 10 10 10 10 10 Degrees of Freedom (DOFs) Degrees of Freedom (DOFs)

L-shaped domain with incompatible data Some meshes: s = 0 . 2 s = 0 . 8

A posteriori error analysis and adaptivity LC, RHN, EO, AJS: A PDE approach to fractional diffusion in general domains: a posteriori error analysis . J. Comput. Phys. (2015). ◮ Question: Is there any theory on anisotropic error estimators? (Cohen Mirebeau 2010-2012) (Petrushev 2007-2009)? ◮ A posteriori error estimators, convergence of AFEM, convergence rates for AFEM are still open questions.

Outline Motivation The elliptic linear problem case Space-time fractional parabolic problem Formulation Localization Discretization The fractional obstacle problem An optimal control problem Conclusions and future work

Outline Motivation The elliptic linear problem case Space-time fractional parabolic problem Formulation Localization Discretization The fractional obstacle problem An optimal control problem Conclusions and future work

Space-time fractional parabolic problem Let T > 0 be some positive time. Given f : Ω → R and u 0 : Ω → R the problem reads: Find u such that ∂ γ t u + ( − ∆) s u = f in Ω × (0 , T ] u | t =0 = u 0 in Ω . Here γ ∈ (0 , 1]. For γ = 1 this is the usual time derivative, if γ < 1 we consider the Caputo derivative � t 1 ∂ r u ( x , r ) ∂ γ t u ( x , t ) = ( t − r ) γ d r . Γ(1 − γ ) 0 Nonlocality in space and time! We will overcome the nonlocality in space using the Caffarelli-Silvestre extension.

Outline Motivation The elliptic linear problem case Space-time fractional parabolic problem Formulation Localization Discretization The fractional obstacle problem An optimal control problem Conclusions and future work

Extended evolution problem The Caffarelli-Silvestre extension turns our problem into a quasi stationary elliptic problem with dynamic boundary condition · ( y α ∇U ) = 0 ,  −∇ in C , t ∈ (0 , T ) ,    U = 0 , on ∂ L C , t ∈ (0 , T ) ,    t U + ∂ U d s ∂ γ ∂ν α = d s f , on Ω × { 0 } , t ∈ (0 , T ) ,      U = u 0 , on Ω × { 0 } , t = 0 .  Connection: u = tr Ω U , α = 1 − 2 s . Nonlocality just in time! Weak formulation: seek U ∈ V such that for a.e. t ∈ (0 , T ), � tr Ω ∂ γ � t U , tr Ω φ � H − s (Ω) × H s (Ω) + a ( w , φ ) = � f , tr Ω φ � H − s (Ω) × H s (Ω) , tr Ω U (0) = u 0 ◦ H 1 L ( C , y α ), where for all φ ∈ a ( w , φ ) = 1 � y α ∇ w · ∇ φ. d s C

Extended evolution problem The Caffarelli-Silvestre extension turns our problem into a quasi stationary elliptic problem with dynamic boundary condition · ( y α ∇U ) = 0 ,  −∇ in C , t ∈ (0 , T ) ,    U = 0 , on ∂ L C , t ∈ (0 , T ) ,    t U + ∂ U d s ∂ γ ∂ν α = d s f , on Ω × { 0 } , t ∈ (0 , T ) ,      U = u 0 , on Ω × { 0 } , t = 0 .  Connection: u = tr Ω U , α = 1 − 2 s . Nonlocality just in time! Weak formulation: seek U ∈ V such that for a.e. t ∈ (0 , T ), � tr Ω ∂ γ � t U , tr Ω φ � H − s (Ω) × H s (Ω) + a ( w , φ ) = � f , tr Ω φ � H − s (Ω) × H s (Ω) , tr Ω U (0) = u 0 ◦ H 1 L ( C , y α ), where for all φ ∈ a ( w , φ ) = 1 � y α ∇ w · ∇ φ. d s C

Truncation C is infinite, but we have exponential decay. Theorem Let γ ∈ (0 , 1] and s ∈ (0 , 1) . If Y > 1 then �∇U� L 2 (0 , T ; L 2 (Ω × ( Y , ∞ ) , y α )) � e −√ λ 1 / 2 ◮ This allows us to consider a truncated problem. ◮ In doing so we commit only an exponentially small error L 2 (Ω) + �∇ ( U − v ) � L 2 (0 , T ; L 2 ( C Y , y α )) � e −√ λ 1 Y I 1 − γ � tr Ω ( U − v ) � 2 where I σ is the Riemann Liouville fractional integral of order σ .

Outline Motivation The elliptic linear problem case Space-time fractional parabolic problem Formulation Localization Discretization The fractional obstacle problem An optimal control problem Conclusions and future work

Time discretization for γ = 1 Time step τ = T / K . Compute V τ = { V k } K ◦ H 1 L ( y α , C ), k =0 ⊂ where V k denotes an approximation at each time step. For γ = 1, we consider backward Euler ◮ We initialize by setting tr Ω V 0 = u 0 . ◮ For k = 0 , . . . , K − 1, we find V k +1 ∈ ◦ H 1 L ( y α , C ) solution of ( tr Ω ∂ V k +1 , tr Ω W ) L 2 (Ω) + a ( V k +1 , W ) = � f k +1 , tr Ω W � H − s (Ω) × H s (Ω) , L ( C , y α ), where f k +1 = f ( t k +1 ). ◦ H 1 for all W ∈ ◮ Unconditional stability: � tr Ω V τ � 2 ℓ ∞ ( L 2 (Ω)) + � V τ � 2 L ( y α , C )) � � u 0 � 2 L 2 (Ω) + � f τ � 2 ℓ 2 ( H − s (Ω)) . ◦ H 1 ℓ 2 (

Time discretization for γ ∈ (0 , 1) For γ ∈ (0 , 1), we consider the so-called L1 scheme � t k +1 1 ∂ r u ( x , r ) ∂ γ t u ( x , t k +1 ) = ( t k +1 − r ) γ d r Γ(1 − γ ) 0 k 1 u ( x , t k +1 − j ) − u ( x , t k − j ) � ≈ a j τ γ Γ(2 − γ ) j =0 = D γ u ( x ) k +1 where a j = ( j + 1) 1 − γ − j 1 − γ . For γ ∈ (0 , 1), the scheme reads ◮ We initialize by setting tr Ω V 0 = u 0 . ◮ For k = 0 , . . . , K − 1, we find V k +1 ∈ ◦ H 1 L ( C , y α ) solution of ( tr Ω D γ V k +1 , tr Ω W ) L 2 (Ω) + a ( V k +1 , W ) = � f k +1 , tr Ω W � H − s (Ω) × H s (Ω) .

Time discretization for γ ∈ (0 , 1). Stability The lack of fractional integration by parts makes it difficult to obtain energy estimates. We obtain new semidiscrete energy estimates for the L1 scheme Theorem (stability) I 1 − γ � tr Ω V τ � 2 L 2 (Ω) + � V τ � 2 L ( C , y α )) ≤ I 1 − γ � u 0 � 2 L 2 (Ω) + � f τ � 2 ℓ 2 ( H − s (Ω)) , ◦ H 1 ℓ 2 ( Since these are uniform in τ and the scheme is consistent 1 we derive a novel continuous energy estimate Theorem (energy estimates) I 1 − γ � tr Ω U� 2 L 2 (Ω) + �U� 2 L ( C , y α )) ≤ I 1 − γ � u 0 � 2 L 2 (Ω) + � f τ � 2 ℓ 2 ( H − s (Ω)) . ◦ H 1 ℓ 2 ( 1 see next slide

Time discretization for γ ∈ (0 , 1). Consistency ◮ The literature analyzes the L1 scheme assuming smoothness of the solution u ∈ C 2 ([0 , T ] , H − s (Ω)). ◮ However, in general, this assumption is not valid! ◮ We showed that ∂ t u ∈ L log L (0 , T , H − s (Ω)) and ∂ tt u ∈ L 2 ( t σ , (0 , T )) , for σ > 3 − 2 γ . These are valid under realistic assumptions on f and u 0 .

Time discretization for γ ∈ (0 , 1). Consistency ◮ Using these new regularity estimates we can provide an analysis of the L1 scheme. ◮ Since t u ( x , t k +1 ) = D γ u ( x ) k +1 + r τ ∂ γ γ and the remainder satisfies � r τ γ � H − s (Ω) � τ θ � � � u 0 � H 2 s (Ω) + � f � H 2 (0 , T ; H − s (Ω)) , where θ < 1 2 . ◮ Key result: u t ∈ L log L (0 , T ; H − s (Ω)). Hardy and Littlewood yields I 1 − γ : L log L (0 , T ) → L 1 γ (0 , T ) boundedly.

Error estimates for fully discrete schemes Discretization in time and space: stability + consistency yield ◮ Error estimates for U : s ∈ (0 , 1) and γ ∈ (0 , 1) 1 [ I 1 − γ � tr Ω ( v τ − V τ 2 � τ θ + | log N | 2 s N − (1+ s ) T Y ) � L 2 (Ω) ( T )] n +1 − 1 � v τ − V τ L ( C Y , y α )) � τ θ + | log N | s N n +1 . T Y � ℓ 2 ( ◦ H 1 ◮ Error estimates for u : s ∈ (0 , 1) and γ ∈ (0 , 1) 1 − (1+ s ) [ I 1 − γ � u τ − U τ � L 2 (Ω) ( T )] 2 � τ θ + | log N | 2 s N n +1 � u τ − U τ � ℓ 2 ( H s (Ω)) � τ θ + | log N | s N − 1 n +1 , where θ < 1 2 . RHN, EO, AJS. A PDE approach to space-time fractional parabolic problems . SIAM J. Numer. Analysis. 2014 (submitted).

Outline Motivation The elliptic linear problem case Space-time fractional parabolic problem The fractional obstacle problem Motivation Formulation Truncation Error estimates An optimal control problem Conclusions and future work

Outline Motivation The elliptic linear problem case Space-time fractional parabolic problem The fractional obstacle problem Motivation Formulation Truncation Error estimates An optimal control problem Conclusions and future work

Classical obstacle problem ◮ Consider a surface given by the graph of a function u . ◮ u solves ∆ u = 0 for fixed boundary data (elastic membrane). ◮ Let us now slide an obstacle from below. The surface must stay above it.

Classical obstacle problem ◮ Consider a surface given by the graph of a function u . ◮ u solves ∆ u = 0 for fixed boundary data (elastic membrane). ◮ Let us now slide an obstacle from below. The surface must stay above it.

Classical obstacle problem ◮ Consider a surface given by the graph of a function u . ◮ u solves ∆ u = 0 for fixed boundary data (elastic membrane). ◮ Let us now slide an obstacle from below. The surface must stay above it.

Classical obstacle problem ◮ Consider a surface given by the graph of a function u . ◮ u solves ∆ u = 0 for fixed boundary data (elastic membrane). ◮ Let us now slide an obstacle from below. The surface must stay above it. ◮ For a given obstacle ψ , we obtain a function u ≥ ψ , that will try to be as harmonic as possible.

Classical obstacle problem ◮ ∆ u = 0 when u > ψ , since there u is free to move. ◮ ∆ u ≤ 0 everwhere, since the surface pushes down. ◮ u ≥ ψ . ◮ Complementarity system: λ = − ∆ u ≥ 0 , u − ψ ≥ 0 , ∆ u ( u − ψ ) = 0 a.e. in Ω .

Motivation for the fractional obstacle problem ◮ Consider τ E ( ψ ( X x u = sup τ )) , where X x τ is a purely jump process starting at x and τ denotes any stopping time. ◮ Then u ≥ ψ, Lu ≥ 0 , Lu = 0 if u > ψ, where the operator L is � Lu ( x ) = P.V. ( u ( x ) − u ( x + u )) K ( y ) . ◮ Natural example: K ( y ) = | y | − ( n +2 s ) with s ∈ (0 , 1) gives ( − ∆) s u = 0 where u > ψ, ( − ∆) s u ≥ 0 everywhere , u ≥ ψ.

Outline Motivation The elliptic linear problem case Space-time fractional parabolic problem The fractional obstacle problem Motivation Formulation Truncation Error estimates An optimal control problem Conclusions and future work

The fractional obstacle problem ◮ Given f ∈ H − s (Ω) and an obstacle ψ ∈ H s (Ω) ∩ C (¯ Ω) satisfying ψ ≤ 0 on ∂ Ω: � ( − ∆) s u , u − w � ≤ � f , u − w � u ∈ K : ∀ w ∈ K . ◮ K := { w ∈ H s (Ω) : w ≥ ψ a.e. in Ω } . ◮ Nonlinear and nonlocal problem since ( − ∆) s ! ◮ We use Caffarelli-Silvestre extension! In fact, the study of the regularites properties of the fractional obstacle problem motivated the Caffarelli-Silvestre extension.

Thin obstacle problem ◮ We convert the fractional obstacle problem in a thin obstacle problem. ◮ The restriction U > ψ only applies when y = 0 (thin obstacle).

Outline Motivation The elliptic linear problem case Space-time fractional parabolic problem The fractional obstacle problem Motivation Formulation Truncation Error estimates An optimal control problem Conclusions and future work

Thin obstacle problem ◮ Truncation of the cylinder: �∇ ( U − V ) � L 2 ( y α , C Y ) � e −√ λ 1 Y / 8 � � � ψ � H s (Ω) + � f � H − s (Ω) . ◮ To derive an error estimate the following regularity results are fundamental: ◮ u ∈ C 1 ,α for α < s by Silvestre (2007). ◮ Optimal regularity: u ∈ C 1 , s by Cafarelli, Salsa and Silvestre (2008). ◮ ∂ α ν U ( · , 0) ∈ C 0 , 1 − s (Ω). ◮ Optimal regularity by Allen, Lindgren, and Petrosyan (2014) s ≤ 1 2 ⇒ V ∈ C 0 , 2 s ( C Y ) and s > 1 2 ⇒ V ∈ C 1 , 2 s − 1 ( C Y ).

Outline Motivation The elliptic linear problem case Space-time fractional parabolic problem The fractional obstacle problem Motivation Formulation Truncation Error estimates An optimal control problem Conclusions and future work

Thin obstacle problem ◮ Nearly optimal error estimate: L ( y α , C ) ≤ C | log N | s N − 1 / ( n +1) , �U − V T Y � ◦ H 1 where C depends on the H¨ older moduli of smoothness of U and V , � f � H − s (Ω) and � ψ � H s (Ω) . ◮ Same techniques can be applied for the Signori or thin obstacle problem. RHN, EO, AJS. Convergence rates for the obstacle problem: classical, thin and fractional , Phil. Trans. R. Soc. A (2015).

Outline Motivation The elliptic linear problem case Space-time fractional parabolic problem The fractional obstacle problem An optimal control problem Formulation Localization Discretization Conclusions and future work

Motivation: Cardiac Microstructure ◮ The heart has its own internal electrical system that controls the rate and rhythm of heartbeat. ◮ Heartbeat produces an electrical signal that spreads from the top to the bottom: it causes the heart to contract and pump blood. ◮ Problems with this electrical system cause arrhythmia! ◮ Implantable cardioverter defibrillator (ICD): monitors the heart rhythm. ◮ If an irregular rhythm is detected, it will use low-energy electrical pulses to restore a normal rhythm. ◮ Fundamental modeling to understand the propagation of electrical excitation is: ∂ t u − ∆ u = f

Motivation: Cardiac Microstructure ◮ This conventional model neglects the highly complex, heterogeneous nature of the underlying tissues. ◮ Bueno-Orovio, Kay, Grau, Rodriguez, and Burrage (2014): ∂ t u + ( − ∆) s u = f .

Outline Motivation The elliptic linear problem case Space-time fractional parabolic problem The fractional obstacle problem An optimal control problem Formulation Localization Discretization Conclusions and future work

Problem Formulation Define J (u , z) = 1 L 2 (Ω) + λ 2 � u − u d � 2 2 � z � 2 L 2 (Ω) . We are interested in the optimal control problem: min J (u , z) subject to the non-local state equation L s u = z in Ω , u = 0 on ∂ Ω , and the control constraints a.e. x ′ ∈ Ω } . z ∈ Z ad := { w ∈ L 2 (Ω) : a( x ′ ) ≤ w( x ′ ) ≤ b( x ′ ) Here, L w = −∇ · x ′ ( A ∇ x ′ w ) + cw .

Outline Motivation The elliptic linear problem case Space-time fractional parabolic problem The fractional obstacle problem An optimal control problem Formulation Localization Discretization Conclusions and future work

An equivalent control problem The Caffarelli-Silvestre result allows us to rewrite our control problem as follows: min J ( tr Ω U , z) = 1 L 2 (Ω) + λ 2 � tr Ω U − u d � 2 2 � z � 2 L 2 (Ω) subject to the linear and local state equation 1 � ◦ y α ∇U · ∇ φ = � z , tr Ω φ � H − s (Ω) , H s (Ω) , H 1 L ( y α , C ) , ∀ φ ∈ d s C and the control constraints a.e. x ′ ∈ Ω } . z ∈ Z ad := { w ∈ L 2 (Ω) : a( x ′ ) ≤ w( x ′ ) ≤ b( x ′ ) z , ¯ Existence and uniquess of an optimal pair (¯ U ) follows standard arguments.

Outline Motivation The elliptic linear problem case Space-time fractional parabolic problem The fractional obstacle problem An optimal control problem Formulation Localization Discretization Conclusions and future work

A truncated control problem min J ( tr Ω v , r) = 1 L 2 (Ω) + λ 2 � tr Ω v − u d � 2 2 � r � 2 L 2 (Ω) , subject to the truncated state equation 1 � y α ∇ v · ∇ φ = � r , tr Ω φ � H − s (Ω) × H s (Ω) , H 1 ◦ L ( y α , C Y ) , ∀ φ ∈ d s C Y and the control constraints r ∈ Z ad . First order necessary and sufficient optimality conditions: ◦ H 1 L ( y α , C Y ) solution of state equation ,  ¯ v = ¯ v (¯ r) ∈   H 1 ◦ L ( y α , C Y ) solution of adjoint equation , ¯ p = ¯ p (¯ r) ∈   ¯ r ∈ Z ad , ( tr Ω ¯ p + λ ¯ r , r − ¯ r) L 2 (Ω) ≥ 0 ∀ r ∈ Z ad . Exponential convergence: For every Y ≥ 1, we have z � L 2 (Ω) � e −√ λ 1 Y / 4 � � � ¯ r − ¯ � ¯ r � L 2 (Ω) + � u d � L 2 (Ω) , � ¯ � L 2 ( y α , C ) � e −√ λ 1 Y / 4 ( � ¯ � �∇ U (¯ z) − ¯ v (¯ r) r � L 2 (Ω) + � u d � L 2 (Ω) ) .

A truncated control problem min J ( tr Ω v , r) = 1 L 2 (Ω) + λ 2 � tr Ω v − u d � 2 2 � r � 2 L 2 (Ω) , subject to the truncated state equation 1 � y α ∇ v · ∇ φ = � r , tr Ω φ � H − s (Ω) × H s (Ω) , H 1 ◦ L ( y α , C Y ) , ∀ φ ∈ d s C Y and the control constraints r ∈ Z ad . First order necessary and sufficient optimality conditions: ◦ H 1 L ( y α , C Y ) solution of state equation ,  ¯ v = ¯ v (¯ r) ∈   H 1 ◦ L ( y α , C Y ) solution of adjoint equation , ¯ p = ¯ p (¯ r) ∈   ¯ r ∈ Z ad , ( tr Ω ¯ p + λ ¯ r , r − ¯ r) L 2 (Ω) ≥ 0 ∀ r ∈ Z ad . Exponential convergence: For every Y ≥ 1, we have z � L 2 (Ω) � e −√ λ 1 Y / 4 � � � ¯ r − ¯ � ¯ r � L 2 (Ω) + � u d � L 2 (Ω) , � ¯ � L 2 ( y α , C ) � e −√ λ 1 Y / 4 ( � ¯ � �∇ U (¯ z) − ¯ v (¯ r) r � L 2 (Ω) + � u d � L 2 (Ω) ) .