Numerical Approximation of the Integral Fractional Laplacian Wenyu - - PowerPoint PPT Presentation

Numerical Approximation of the Integral Fractional Laplacian

Wenyu Lei Department of Mathematics Texas A&M University Joint work with Andrea Bonito and Joseph Pasciak July 26, 2018 • deal.II Workshop • SISSA

Introduction Algorithm Implementation Results Application Conclusion

Outline Introduction Numerical Algorithm Implementation in deal.II Results Application: an obstacle problem for a class of integro-differential operators Conclusion

Approximation of fractional Laplaican Wenyu Lei

Introduction Algorithm Implementation Results Application Conclusion

Outline Introduction Numerical Algorithm Implementation in deal.II Results Application: an obstacle problem for a class of integro-differential operators Conclusion

Approximation of fractional Laplaican Wenyu Lei

Introduction Algorithm Implementation Results Application Conclusion

Integral Fractional Laplacian (−∆)s

• Let s ∈ (0, 1) and η : Rd → R be in the Schwarz class.
• Definition via integral representation:

((−∆)sη)(x) = cd,sP.V.

• Rd

η(x) − η(y) |x − y|d+2s dy, where P.V. stands for the principle value and cd,s is a normalization constant.

• An equivalent definition via Fourier transform is given by

F((−∆)sη)(ζ) = |ζ|2sF(η)(ζ). (∗) Here F is the Fourier transform.

• (∗) defines the unbounded operator (−∆)s on L2(Rd) with domain of definition

D((−∆)s) := {f ∈ L2(Rd) : |ζ|2sF(f) ∈ L2(Rd)} ⊂ Hs(Rd).

• Applications: s-stable L´

evy process, electroconvection and the surface quasigeostrophic models [I. Held et al., 1995].

Approximation of fractional Laplaican Wenyu Lei

Introduction Algorithm Implementation Results Application Conclusion

Boundary Value Problem

• Ω: open, bounded domain in Rd with Lipschitz boundary;
• f ∈ L2(Ω);
• We consider the stationary problem

(−∆)s u|Ω = f, in Ω. Here u means zero extension of u.

• It is a nonlocal problem.
• The test space:
• Hs(Ω) := {f ∈ L2(Ω) :

fHs(Rd) < ∞}.

• For any η, θ ∈

Hs(Ω), define the bilinear form a(η, θ) := ((−∆)s/2 η, (−∆)s/2 θ)L2(Rd) =

• Rd |ζ|2sF(

η)F( θ) dζ.

• Therefore, the weak formulation is: find u ∈

Hs(Ω) satisfying a(u, v) =

• Ω

fv dx, for all v ∈ Hs(Ω).

Approximation of fractional Laplaican Wenyu Lei

Introduction Algorithm Implementation Results Application Conclusion

Finite Element Approximations of the Boundary Value Problem

Finite element approaches:

[Acosta & Borthagaray, 2017]: Directly approximate each entry in the stiffness matrix using the integral form a(u, φ) = cd,s 2

• Rd
• Rd(

u(x) − φ(y))( u(x) − ( φ(y)) 1 |x − y|d+2s dy dx together with a certain mesh setting and special quadrature formulas (boundary element approach). [D’ Elia & Gunzburger, 2013]: Approximate the above integral in a bounded truncated domain. Approximate the bilinear form based on its Dunford-Taylor integral representation.

Approximation of fractional Laplaican Wenyu Lei

Introduction Algorithm Implementation Results Application Conclusion

A Dunford-Taylor Integral Representation of a(·, ·)

An equivalent representation

a(η, θ) = 2 sin(πs) π ∞ t1−2s

• Rd
• (−∆)(I − t2∆)−1

η θ dx dt =: I

Proof

• Parseval’s Theorem:
• Rd
• (−∆)(I − t2∆)−1

η θ dx =

• Rd

|ζ|2 1 + t2|ζ|2 F( η)(ζ)F( θ)(ζ) dζ.

• Apply the Fubini’s Theorem and use the change of variable y = t|ζ|:

I = 2 sin(πs) π

• Rd |ζ|2sF(η)(ζ)F(θ)(ζ)

∞ t1−2s |ζ|2−2s 1 + t2|ζ|2 dt dζ =

• Rd |ζ|2sF(η)(ζ)F(θ)(ζ) dζ.

Note that cs := ∞ y1−2s 1 + y2 dy −1 = 2 sin(πs) π .

Approximation of fractional Laplaican Wenyu Lei

Introduction Algorithm Implementation Results Application Conclusion

Game Plan a(η, θ) = cs ∞ t−1−2s

• Ω
• (−t2∆)(I − t2∆)−1

η

• w(t)

θ dx dt for η, θ ∈ Hs(Ω). Here w(t) = η−(I − t2∆)−1 η := η + v(t) with v(t) satisfying

• Rd v(t)φ dx + t2
• Rd ∇v(t) · ∇φ dx = −
• Ω

ηφ dx, for all φ ∈ H1(Rd),

• Discretize the outer integral with a quadrature spacing k:

ak(η, θ) = Ck

• j

t−1−2s

j

• Ω

w(tj)θ dx.

• Approximate w(t) or (v(t)) in a bounded domain BM(t):

ak,M(η, θ) = Ck

• j

t−1−2s

j

• Ω

wM(tj)θ dx.

• Approximate wM(t) (or vM(t)) using the finite element method:

ak,M

h

(ηh, θh) = Ck

• j

t−1−2s

j

• Ω

wM

h (tj)θh dx.

Approximation of fractional Laplaican Wenyu Lei

Introduction Algorithm Implementation Results Application Conclusion

Outline Introduction Numerical Algorithm Implementation in deal.II Results Application: an obstacle problem for a class of integro-differential operators Conclusion

Approximation of fractional Laplaican Wenyu Lei

Introduction Algorithm Implementation Results Application Conclusion

Sinc Quadrature

• Use the change of variable t = e−y/2;

a(η, θ) = cs 2 ∞

−∞

esy

• Ω

w(t(y))θ dx dy.

• Use a truncated equally spaced quadrature: let k > 0 and N+ and N− are

positive integers chosen to be on the order of 1/k2. Define ak(η, θ) := csk 2

N+

• j=−N−

esyj

• Ω

w(yj)θ dx. with yj = jk.

• Consistency Error:

|a(η, θ) − ak(η, θ)| ≤ Ce−c/kη

Hδ(Ω)θ Hs(Ω).

Here δ ∈ (s, 2 − s].

Approximation of fractional Laplaican Wenyu Lei

Introduction Algorithm Implementation Results Application Conclusion

Domain Truncation

• w(t) =

η + v(t) where v(t) = −(I − t2∆)−1 η.

• Let B be a ball enclosing Ω and assume that diamB = 1, then for a truncation

parameter M, we define the dilated domains BM(t) := {y = (1 + t(1 + M))x : x ∈ B} , t ≥ 1 {y = (2 + M)x : x ∈ B} , t < 1.

• Truncated solution: wM(t) =

η + vM(t) and

• BM (t)

vM(t)φ dx+t2

• BM (t)

∇vM(t)·∇φ dx = −

• Ω

ηφ dx, for all φ ∈ H1

0(BM(t)),

• Truncated bilinear form:

ak,M(η, θ) := csk 2

N+

• j=−N−

esyj

• Ω

wM(yj)θ dx.

• Consistency Error:

|ak(η, θ) − ak,M(η, θ)| ≤ Ce−cMηL2(Ω)θL2(Ω), for all η, θ ∈ Hs(Ω).

Approximation of fractional Laplaican Wenyu Lei

Introduction Algorithm Implementation Results Application Conclusion

Two Finite Element Spaces

• The red part is the grid of the original domain and it is quasi-uniform.
• The red+black part is the grid of the truncated domain. The black part may not

be quasi-uniform for implementation. But we assume the quasi-uniformity for numerical analysis.

• Vh(Ω) ⊂ Vh(BM(t)) := VM

h (t).

Approximation of fractional Laplaican Wenyu Lei

Introduction Algorithm Implementation Results Application Conclusion

Finite Element Approximation

• The discrete solution: wM

h (t) =

η + vM

h (t) and for all φh ∈ VM(t),

• BM (t)

vM

h (t)φh dx + t2

• BM (t)

∇vM

h (t) · ∇φh dx = −

• Ω

ηhφh dx.

• The discrete bilinear form: for ηh, θh ∈ Vh(Ω)

ak,M

h

(ηh, θh) := csk 2

N+

• j=−N−

esyj

• Ω

wM

h (yj)θh dx.

• Consistency Error: let β ∈ (s, 3/2),

|ak,M(ηh, θh) − ak,M

h

(ηh, θh)| ≤ C(1 + ln(h−1))hβ−sηh

Hβ(Ω)θh Hs(Ω)

Approximation of fractional Laplaican Wenyu Lei

Introduction Algorithm Implementation Results Application Conclusion

Discrete Problem

• The discrete problem: find uh ∈ Vh(Ω) satisfying

ak,M

h

(uh, vh) =

• Ω

fvh dx for all vh ∈ Vh(Ω).

• Vh(Ω)-ellipticity: For a fixed h, choose k small enough so that

Ce−c/khs−δ < 1/2. Then, ak,M

h

(ηh, ηh) > 1/2ηh2

• Hs(Ω).

Assume u ∈ Hβ(Ω) with β ∈ (s, 3/2) and α = min(s, 1/2).

Error estimates

Strang’s Lemma: u − uh

Hs(Ω) ≤ C(e−c/k + e−cM + (1 + ln (h−1))hβ−s)u Hβ(Ω).

Duality argument: if the domain is smooth, u − uhL2(Ω) ≤ C ln(h−1)(e−c/k + e−cM + (1 + ln (h−1))hβ−s+α)u

Hβ(Ω).

Approximation of fractional Laplaican Wenyu Lei

Introduction Algorithm Implementation Results Application Conclusion

Outline Introduction Numerical Algorithm Implementation in deal.II Results Application: an obstacle problem for a class of integro-differential operators Conclusion

Approximation of fractional Laplaican Wenyu Lei

Introduction Algorithm Implementation Results Application Conclusion

Non-uniform Mesh Setting: 1d case

• Let Ω be a unit ball.
• Th(Ω) is uniform with the mesh size h.
• We consider an exponential grading in BM(t)\Ω.
• Recall that the radius of BM(t) is 1 + (M + 1)t for t > 1.

1 eh0 e2h0 · · · 1 + (M + 1)t = eN h0 R Here h0 = h log(1 + (M + 1)t)/M and N = ⌈M/h⌉.

Approximation of fractional Laplaican Wenyu Lei

Introduction Algorithm Implementation Results Application Conclusion

Non-uniform Mesh Setting: 2d case coarse grid refine once refine twice three times refinement

• The red part is the grid of Ω.
• In the coarse grid: all vertices of a square have the same radial coordinates.
• The position of the vertices along the radial direction and outside of Ω follows

the same exponential distribution as in the one dimensional case.

• We refine each cell in Ω by connecting the midpoints between opposite edges.
• For the cells outside of Ω, we consider the same refinement in the polar

coordinate system (ln r, θ) with r > 1 and θ ∈ [0, 2π].

• Therefore, we maintain the same number of mesh points for all BM(t).

Approximation of fractional Laplaican Wenyu Lei

Introduction Algorithm Implementation Results Application Conclusion

System Matrix sin (πβ)k π MΩ,h

N+

• i=−N−

esyiR(eyiMh(ti) + Ah(ti))−1Ah(ti)EU = F

• ti = e−yi/2.
• F ∈ Rdim(Vh(Ω)): vector of inner product between f and shape functions.
• U ∈ Rdim(Vh(Ω)): coefficient vector of uh.
• E: zero extension from Rdim(Vh(Ω)) to Rdim(VhM (t)).
• R: restriction to Rdim(Vh(Ω)).
• Mh(ti) and Ah(ti): mass and stiffness matrix in VhM(t).
• MΩ,h: mass matrix in Vh(Ω).

About the evaluation:

• We parallelized the summation using mpi and sum up the results by

Utilities::MPI::sum().

• We apply preconditioned conjugate gradient method with multigrid

preconditioner to compute the evaluation of (eyiMh(ti) + Ah(ti))−1 (step-16).

Approximation of fractional Laplaican Wenyu Lei

Introduction Algorithm Implementation Results Application Conclusion

Preconditioned CG Iteration

• The system matrix is symmetric.
• We use the conjugated gradient method to solve the linear system.
• Condition Number: Ch−2s (because

Hs(Ω) coincides with [L2(Ω), H1

0(Ω)]s).

Reconditioner

• Negative power of the Dirichlet Laplaican: (−∆D)−s.
• In 1d uniform mesh, this can be done by discrete sine transform.
• For a general domain: solve a Caffarelli-Silvestre extension problem [Nochetto,

Otarola, and Salgado, 2015] or use Dunford-Taylor integral approach [Bonito and Pasciak, 2015].

• For nested triangular meshes: use a multilevel preconditioner [Bramble, Pasciak,

and Vassilevski, 2000] BJ =

J

• j=1

h2s

j (

Qj − Qj−1)2, where J is the number of levels and

• Qjv =

dim(Vh(Ω))

• i=1

(v, φi)Ω (1, φi)Ω φi.

Approximation of fractional Laplaican Wenyu Lei

Introduction Algorithm Implementation Results Application Conclusion

Outline Introduction Numerical Algorithm Implementation in deal.II Results Application: an obstacle problem for a class of integro-differential operators Conclusion

Approximation of fractional Laplaican Wenyu Lei

Introduction Algorithm Implementation Results Application Conclusion

One Dimensional Numerical Test

• f = 1, u = c(s)(1 − x2)s.
• Compute u − uhL2(D).

k M

• h = 1/8192
• h = 1/8192
• M = 20
• k = 0.2

Approximation of fractional Laplaican Wenyu Lei

Introduction Algorithm Implementation Results Application Conclusion

One Dimensional Numerical Test

• u = 1 − x2, f(x) = c(s)2F1
• 1/2 + s, s − 1, 1/2, x2

.

• M = 5, k = 0.2.
• Compute u − uhL2(D).

h s = 0.3(1.5) s = 0.4(1.5) s = 0.7(1.3) 1/16 4.51 × 10−4 3.47 × 10−4 9.27 × 10−4 1/32 1.42 × 10−4 1.58 1.02 × 10−4 1.77 4.16 × 10−4 1.16 1/64 4.25 × 10−5 1.63 3.31 × 10−5 1.62 1.80 × 10−4 1.21 1/128 1.34 × 10−5 1.66 1.14 × 10−5 1.54 7.66 × 10−5 1.23 1/256 4.43 × 10−6 1.59 4.06 × 10−6 1.49 3.21 × 10−5 1.25 1/512 1.50 × 10−6 1.56 1.46 × 10−6 1.48 1.33 × 10−5 1.27

Approximation of fractional Laplaican Wenyu Lei

Introduction Algorithm Implementation Results Application Conclusion

Two Dimensional Numerical Test

• f = 1, u = c(s)(1 − |x|2)s, M = 4, k = 0.25.
• Compute u − uhL2(D).

#DOFS s = 0.3 (0.8) s = 0.5 (1.0) s = 0.7 (1.0) 345 2.69 × 10−1

• 1.63 × 10−1
• 1.03 × 10−1
• 1361

1.59 × 10−1 0.7575 9.07 × 10−2 0.8426 5.55 × 10−2 0.8918 5409 9.56 × 10−2 0.7323 5.05 × 10−2 0.8438 2.95 × 10−2 0.9091 21569 5.71 × 10−2 0.7447 2.78 × 10−2 0.8633 1.54 × 10−2 0.9366 86145 3.38 × 10−2 0.7547 1.51 × 10−2 0.8832 7.91 × 10−3 0.9641 344321 1.99 × 10−2 0.7644 8.07 × 10−3 0.9004 3.97 × 10−3 0.9936 s = 0.3 s = 0.7

Approximation of fractional Laplaican Wenyu Lei

Introduction Algorithm Implementation Results Application Conclusion

Outline Introduction Numerical Algorithm Implementation in deal.II Results Application: an obstacle problem for a class of integro-differential operators Conclusion

Approximation of fractional Laplaican Wenyu Lei

Introduction Algorithm Implementation Results Application Conclusion

Obstacle Problem (joint work with A. Bonito and A. Salgado) Motivation: perpetual American options under s-stable L´ evy precesses. min {−σ∆u + β · ∇w + (−∆)su − f, u − χ} = 0, in Ω, u = 0, in Ωc. where

• σ ≥ 0 and β ∈ L∞(Ω)d, div(β) = 0.
• χ is a smooth obstacle such that χ|∂Ω < 0.
• f ∈ L2(Ω).

The solution u represents the rational price of a perpetual American option against the log-prices of the stocks.

Approximation of fractional Laplaican Wenyu Lei

Introduction Algorithm Implementation Results Application Conclusion

Results in 1d

• min(−σu′′ + βu′ + (−∆)s − 1, u − χ) = 0.
• Ω = (−1, 1), and χ(x) = 3 − 6x2.
• Use the primal-dual active set method to solve the discrete obstacle problem.

• 1
• 0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1

• bstcale

• 1
• 0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1

• bstcale

• 1
• 0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1

• bstcale

σ = β = 0. σ = 0, β = 0.5. σ = β = 0.5.

Approximation of fractional Laplaican Wenyu Lei

Introduction Algorithm Implementation Results Application Conclusion

Results in 2d

• min((−∆)s − 1, u − χ) = 0.
• Ω = [−0.5, 0.5]2/[0, 0.52], and χ(x) = 256x1(x1 + 0.5)x2(x2 − 0.5).

s = 0.3 s = 0.5

• 0.5 -0.45 -0.4 -0.35 -0.3 -0.25 -0.2 -0.15 -0.1 -0.05

• bstacle

s = 0.3 x = −0.25

Approximation of fractional Laplaican Wenyu Lei

Introduction Algorithm Implementation Results Application Conclusion

Results in 2d

• min((−∆)s + (−0.5, 0)t · ∇u − f, u − χ) = 0.
• f = 1 if |(x − (0.5, 0)| < 0.5 and f = 0 otherwise
• Ω =unit ball, and χ(x) = 3 − 6|x|2.

• 1
• 0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1

• bstacle (1)

s = 0.5 s = 0.7 x2 = 0

Approximation of fractional Laplaican Wenyu Lei

Introduction Algorithm Implementation Results Application Conclusion

Outline Introduction Numerical Algorithm Implementation in deal.II Results Application: an obstacle problem for a class of integro-differential operators Conclusion

Approximation of fractional Laplaican Wenyu Lei

Introduction Algorithm Implementation Results Application Conclusion

Conclusion

• Approximation to the bilinear form

a(·, ·) Sinc Quadrature − − − − − − − − − → ak(·, ·) Domain Truancation − − − − − − − − − − − − → ak,M(·, ·) Finite Element − − − − − − − − → ak,M

h

(·, ·)

• Strang’s Lemma:

u − uh

Hs(Ω) ≤ C(e−c/k + e−cM + (1 + ln (h−1))hβ−s)u Hβ(Ω).

• Non-uniform mesh setting (1d & 2d):
• This method is based on usual finite element subroutines.
• The method is applicable in high-dimension.
• But for each iteration, we have to solve a series of sub-problems (about

C log(1/h)2).

• Application to obstacle problems.

Figure: The approximated singular solutions for s = 0.3 on the unit disk (left) and on the unit ball

Approximation of fractional Laplaican Wenyu Lei