Finite element approximation of the elliptic Isaacs equation Abner - - PowerPoint PPT Presentation

Finite element approximation of the elliptic Isaacs equation

Abner J. Salgado

Department of Mathematics University of Tennessee

RICAM November 24, 2016

Joint work with Wujun Zhang (Rutgers U.) Supported by NSF grant DMS-1418784

Outline

Motivation Building blocks Approximation by an integro-differential operator The scheme Rate of convergence The solution scheme Conclusions

Where do elliptic equations come from?

Diffusion and elliptic equations:

• Particles diffuse in a medium Ω.
• When reaching the boundary ∂Ω, say at point y, the particle

releases a certain amount of energy (payoff) ϕ(y).

• What is the average energy released when starting at point x?
• If the medium is isotropic, then this is given by
• ∆u = 0

in Ω u = ϕ

• n ∂Ω.

Where do elliptic equations come from?

Diffusion and elliptic equations:

• If the medium is not isotropic, the diffusion process can be

described by a matrix A:

• A(x) : D2u(x) = 0

in Ω u = ϕ

• n ∂Ω.
• If

λI A ΛI the equation is uniformly elliptic.

The HJB equation

• Let β ∈ B be a set of design parameters.
• We can choose the material to maximize payoff.

   sup

β∈B

[Aβ(x) : D2u(x)] = 0 in Ω u = ϕ

• n ∂Ω.
• Hamilton Jacobi Bellman equation.

The HJB equation

   F[x, D2u] := sup

β∈B

[Aβ(x) : D2u(x) − fβ(x)] = 0 in Ω u = ϕ

• n ∂Ω.
• The operator F is convex.
• Solutions are u ∈ C2,s(Ω), for some s ∈ (0, 1) (Evans 1982;

Krylov 1982). For this the convexity of F is essential.

The Isaacs equation

• Let α ∈ A and β ∈ B be a set of design parameters.
• Two players compete, the first one wants to maximize payoff,

while the second minimize it.    inf

α∈Asup β∈B

[Aα,β(x) : D2u(x)] = 0 in Ω u = ϕ

• n ∂Ω.
• Isaacs equation.

Differential games

• Two players try to control the process

dXs = ℓ(Xs, s, α, β) ds + σ(Xs, s, α, β) dWs, X0 = x and get a certain payoff J(x, t, α, β) = Ex,t T

t

f(s, Xs, α, β) ds + g(XT )

• .
• The first wants to maximize it, while the second wants to

minimize it inf

α∈Asup β∈B

J(x, t, α, β).

Differential games

• Define the value function

u(x, t) = inf

α∈A sup β∈B

J(x, t, α, β) = J(x, t, α⋆, β⋆).

• Bellman/Isaacs optimality principle (⊕ some stochastic calculus)

yields that inf

α∈A sup β∈B

• ∂tu(x, t) + Aα,β : D2u(x, t) + bα,β · ∇u(x, t) − fα,β

= 0 with Aα,β(x, t) = 1 2σ(x, t, α, β)σ(x, t, α, β)⊺, bα,β(x, t) = ℓ(x, t, α, β), fα,β(x, t) = f(x, t, α, β).

• Parabolic Isaacs equation!

Every equation is Isaacs

Let F : Ω × Sd → R be nondecreasing with respect to its second argument (elliptic). As shown by (Evans 1980), we have the following representation: F[x, M] = inf

α∈Sd sup β∈Sd

1 ∂F ∂M (x, (1 − t)α + tβ) : (M − β) dt + F(x, α)

• Set

Aα,β(x) = 1 ∂F ∂M (x, (1 − t)α + tβ) dt, fα,β(x) = Aα,β : α − F(x, α), A = B = Sd. Notice that Aα,β(x) 0 and we have F[x, M] = inf

α∈A sup β∈B

• Aα,β(x) : M − fα,β(x)
• .

The Isaacs equation

   F[x, D2u] := inf

α∈A sup β∈B

• Aα,β(x) : D2u(x) − fα,β(x)
• = 0

in Ω u = ϕ

• n ∂Ω.
• The operator F is not convex.
• What can we say about the regularity of the solution?
• (Nirenberg 1953) showed that, if d = 2, the solution to every

elliptic equation is locally C2,s for some s ∈ (0, 1).

The counterexample of Nadirashvili and friends

• Let d = 5.
• Define

P5(x) = x3

1 + 3

2x1

• x2

3 + x2 4 − 2x2 5 − 2x2 2

• ,

w(x) = 1 |x|1+δ P5(x), δ ∈ [0, 1)

• We have that w ∈ C1,1−δ( ¯

B1) \ C2(B1).

• There exists a uniformly elliptic Isaacs operator that depends
• nly on M, it is Lipschitz and such that

F[D2w] = 0, in B1 w = P5 on ∂B1.

• Solutions of the Isaacs equation (for d ≥ 5) are not C2.
• Open question: What can we say for d = 3, 4?

Outline

Motivation Building blocks Approximation by an integro-differential operator The scheme Rate of convergence The solution scheme Conclusions

Construction of numerical schemes I

Let Fh[uh] = 0 be an approximation scheme for the elliptic problem F[u] = 0. As shown in (Barles, Souganidis 1991), if

• The scheme is monotone, that is: If uh − vh has a nonnegative

maximum at z, then Fh[uh](z) ≤ Fh[vh](z).

• The scheme is consistent, i.e., if Ih is the interpolation onto the

mesh Fh[Ihw] ≈ F[w] whenever w is smooth. Then the scheme is convergent: uh → u as h → 0.

Construction of numerical schemes II

Question: How do we construct schemes with these properties? Finite Differences:

• Kuo and Trudinger in a series of papers (1990-1992) devised a

program to construct consistent and monotone finite difference schemes. Question: How do we obtain rates of convergence? From (Kuo, Trudinger 1992): So far, such estimates have eluded us and to the best of our knowledge they remain a challenging and important open problem.

Analysis of numerical schemes I

Let F be elliptic and uε be a smooth subsolution such that u − uεL∞(Ω) ≈ εκ1.

• By consistency

Fh[uε] ≈ hκ2,

• By monotonicity

uε − uhL∞(Ω) ≈ hκ2,

• Triangle inequality

u − uhL∞(Ω) ≈ εκ1 + hκ2, and choosing ε we get a rate of convergence.

Analysis of numerical schemes II

Question: How do we construct uε? If the operator F is convex

• (Krylov 1997; 2000) “shake” the coefficients.
• (Barles, Jakobsen 2002) “shake” the coefficients of the scheme.
• (Barles, Jakobsen 2005) relate to a QVI.
• ...

What if F is not convex?

• (Caffarelli, Souganidis 2008) Use the “sup” convolution

u+

ε (x) = sup y∈¯ Ω

• u(y) − 1

2ε|x − y|2

• and very subtle properties of viscosity solutions to get a rate

u − uhL∞(Ω) hσ. The rate σ cannot be determined.

Outline

Motivation Building blocks Approximation by an integro-differential operator The scheme Rate of convergence The solution scheme Conclusions

Approximation: construction

• Introduced in (Caffarelli, Silvestre 2010).

For Ω ⊂ Rd consider    F[u](x) := inf

α∈A sup β∈B

• Aα,β(x) : D2u(x)
• = f(x)

in Ω u = 0

• n ∂Ω,

with λI ≤ Aα,β ≤ ΛI ∀α ∈ A, ∀β ∈ B. Then λ

2I ≤ Aα,β − λ 2I and we can write

F[u](x) = λ 2 ∆u + inf

α∈A sup β∈B

• Aα,β(x) − λ

2 I

• : D2u(x)
• .

Approximation: construction

Let ϕ be a radially symmetric, smooth function, supported in the unit ball, such that

• Rd |y|2ϕ(y) dy = d and

dw(x, y) = w(x + y) − 2w(x) + w(x − y). For ε > 0 define Iα,β

ε

[w] (x) = 1 εd+2 det Mα,β(x)

• Rd dw(x, y)ϕ

1 εMα,β(x)−1y

• dy.

with Mα,β(x) =

• Aα,β(x) − λ

2I

1/2. We approximate F[w](x) by Fε[w](x) := λ 2 ∆w + inf

α∈A sup β∈B

• Iα,β

ε

[w] (x)

• .

Approximation: construction

And consider, instead,

• Fε[uε](x) = f(x)

in Ω uε = 0 in Ωε \ Ω where Ωε =

• y ∈ Rd : dist(y, Ω) < Qε
• with Q =
• 2/λ.

This is a good idea because

Lemma

If w ∈ C2(¯ Ωε) then, as ε → 0, we have Iα,β

ε

[w] (x) →

• Aα,β(x) − λ

2 I

• : D2w(x)

Approximation: properties

Proposition (CS ‘10)

Let Aα,β(x) = Aα,β, Ω be bounded and convex, f ∈ C0,1(¯ Ω). Then there exists a unique function uε ∈ C2,s(Ω) that solves the approximate equation. Moreover, there is s ∈ (0, 1], such that uεC1,s(Ω) + uεC0,1(¯

Ω) uεL∞(Ω) + fL∞(Ω).

In addition, there is a σ > 0 such that u − uεL∞(Ω) εσfC0,1(¯

Ω).

Corollary

There is a s > 0 such that u ∈ C1,s(Ω) ∩ C0,1(¯ Ω). We will discretize the integro-differential equation!

Outline

Motivation Building blocks Approximation by an integro-differential operator The scheme Rate of convergence The solution scheme Conclusions

Notation

• We follow (Nochetto, Zhang 2014).
• Th is a quasiuniform triangulation of Ω of size h.
• Nh are the interior and N ∂

h boundary nodes, respectively.

• Vh is a P1 continuous finite element space, V0

h = Vh ∩ H1 0(Ω).

• {φz}z∈Nh is the Lagrange nodal basis of V0

h.

• For wh ∈ Vh and z ∈ Nh the discrete Laplacian ∆h is

∆hwh(z) = −

• Ω

φz(x) dx −1

Ω

∇wh(x)∇φz(x) dx.

Description of the scheme

Define, for z ∈ Nh, Fε

h[wh](z) := λ

2 ∆hwh(z) + inf

α∈A sup β∈B

• Iα,β

ε

[wh] (z)

• .

We seek for uε

h ∈ V0 h such that

Fε

h[uε h](z) = fz

∀z ∈ Nh, where fz =

• Ω

φz(x) dx −1

Ω

f(x)φz(x) dx.

Properties of the scheme: monotonicity, existence and uniqueness

Lemma (monotonicity)

Assume the mesh Th is weakly acute, then if vh, wh ∈ Vh with equality at some node z ∈ Nh implies Fε

h[vh](z) ≤ Fε h[wh](z).

Monotonicity implies:

Theorem (existence and uniqueness)

For every h > 0 and ε > 0 the scheme has a unique solution.

• The proof is based on a discrete Perron method.

Properties of the scheme: consistency

• We want to use (a variant of) the framework of (Barles,

Souganidis 1991), but as we know operator consistency is too restrictive for finite elements.

• We will instead use the notion of weak consistency originally

proposed in (Jensen, Smears 2013).

• Introduce the Galerkin projection Gh:
• Ω

∇Ghw(x)∇vh(x) dx =

• Ω

∇w(x)∇vh(x) dx ∀vh ∈ Vh, and recall that it satisfies w − GhwL∞(Ωε) h| log h|wC0,1(¯

Ωε),

from which we obtain |Iα,β

ε

[Ghw] (z) − Iα,β

ε

[w] (z)| h ε2 | log h|wC0,1(¯

Ωε)

∀z ∈ Nh.

Properties of the scheme: consistency

If ωz =

• Ω φz(x) dx, define

Rh,ε[w](z) = 1 ωz

• Ω
• inf

α∈A sup β∈B

• Iα,β

ε

[Ghw] (z)

• − inf

α∈A sup β∈B

• Iα,β

ε

[w] (x)

• φz(x) dx

Lemma (integral operator vs. Galerkin projection)

|Rh,ε[w](z)| h ε2 | log h|wC0,1(¯

Ωε)

∀w ∈ C0,1(¯ Ωε).

• In a sense this provides a weak consistency result.

Outline

Motivation Building blocks Approximation by an integro-differential operator The scheme Rate of convergence The solution scheme Conclusions

Rate of convergence: limitations

• Convergence requires ε h1/2| log h|1/2+t so we operate in this

setting.

• The results of (CS ‘10) require constant “coefficients”, which

from now on we assume. (Turanova 2015)??

• In this setting uε ∈ C1,s(Ω) and there is σ > 0 such that

u − uεL∞(Ω) εσ.

The discrete ABP estimate

• consistency ⊕ stability =

⇒ . . .

• We have the right notion of consistency.
• We need a notion of stability. We use the discrete ABP estimate
• f (NZ 2014).

Theorem (discrete ABP)

If vh ∈ Vh, vh(z) ≥ 0 for all z ∈ N ∂

h and

∆hvh(z) ≤ gz ∀z ∈ Nh, then sup

Ω

v−

h

 

• z∈C−(vh)

|gz|dωz  

1/d

, where C−(vh) is the nodal contact set.

The discrete ABP estimate

The error equation

• We know that

u − uεL∞(Ω) εσ.

• By approximation

uε − GhuεL∞(Ω) h| log h|.

• It is “only” left to estimate

Ghuε − uε

hL∞(Ω) ❆!

We need to find an equation for Ghuε − uε

h!

The error equation

• Set eh = Ghuε − uε

h.

• The error satisfies

λ 2 ∆eh(z) + inf

α∈A sup β∈B

• Iα,β

ε

[Ghuε] (z)

• −

inf

α∈A sup β∈B

• Iα,β

ε

[uε

h] (z)

• = Rh,ε[uε](z).
• It turns out that at the contact set C−(eh):

inf

α∈A sup β∈B

• Iα,β

ε

[Ghuε] (z)

• − inf

α∈A sup β∈B

• Iα,β

ε

[uε

h] (z)

• ≤ 0

so that using ABP sup

Ω

e−

h

 

• z∈C−(eh)

|Rh,ε[uε](z)|dωz  

1/d

.

Error estimates

Gluing things together we obtain:

Theorem (rate of convergence)

If u ∈ C1,s(Ω) ∩ C0,1(¯ Ω) solves the Isaacs equation and uε

h ∈ V0 h is

the solution to our scheme. Choose ε h1/2| log h|1/2+t, then there is σ > 0 such that u − uε

hL∞(Ω) εσ + h

ε2 | log h|.

Corollary

If we know σ, set ε ≈ (h| log h|)

1 σ+2 and the error estimate

becomes u − uε

hL∞(Ω) (h| log h|)

σ σ+2 .

• The rate is O(h1/3| log h|1/3) at best!

Outline

Motivation Building blocks Approximation by an integro-differential operator The scheme Rate of convergence The solution scheme Conclusions

The solution scheme: difficulties

• The operator Fε

h is not convex/concave.

• As shown in (Bokanowski, Maroso, Zidani 2009) semi-smooth

Newton methods may not converge!

• We adapt the algorithm of (BMZ ‘09) to our setting, which is

an extension/generalization of Howard’s algorithm (policy iterations), but is not a Newton-like method.

The solution scheme: description

• For vh ∈ V0

h and z ∈ Nh set α(vh, z) ∈ A to be

inf

α∈A sup β∈B

• Iα,β

ε

[vh] (z)

• = sup

β∈B

• Iα(vh,z),β

ε

[vh](z)

• .
• Set α(vh) = {α(vh, z) : z ∈ Nh}.
• For vh ∈ V0

h define the operator

Fα(vh)

h,ε

[wh](z) := λ 2 ∆hwh(z) + sup

β∈B

• Iα(vh,z),β

ε

[wh](z)

• .

The solution scheme: description

The solution scheme: description

• Initialization: Choose w−1

h

∈ V0

h and set α0 = α(w−1 h ).

• Iteration: For k ≥ 0 find wk

h ∈ V0 h that solves

Fαk

h,ε[wk h](z) = fz

∀z ∈ Nh. (IT)

• Convergence test: If

Fε

h[wk h](z) = fz

∀z ∈ Nh stop, else set αk+1 = α(wk

h).

Fact

The iterative step (IT) can be rewritten as: find wk

h ∈ V0 h such that

sup

β∈B

• I

α(wk−1

h

,z),β ε

[wk

h](z)

• = fz.

Discrete HJB equation! It can be solved by semi-smooth Newton methods.

The solution scheme: properties

The functions Fα

h,ε satisfy:

Lemma (monotonicity and comparison)

For every α ∈ A#Nh the operator Fα

h,ε is monotone and satisfies:

if vh, wh ∈ V0

h and

Fα

h,ε[vh] ≤ Fα h,ε[wh]

then vh ≥ wh. With this we obtain:

Theorem (convergence)

If #A < ∞, then the sequence {wk

h}k≥0 ⊂ V0 h converges, in a

finite number of steps, to uε

h.

Outline

Motivation Building blocks Approximation by an integro-differential operator The scheme Rate of convergence The solution scheme Conclusions

Conclusions

• Two scale approximation of the Isaacs equation.
• Existence of discrete solutions via a discrete Perron method.
• Convergence to the solution provided ε h1/2| log h|1/2.
• For constant coefficients, we provide an algebraic rate of

convergence: u − uε

hL∞(Ω) (h| log h|)

σ σ+2 .

• Solution scheme with exact and inexact solves.