[PPT] - Boundary Approximations for Semi-Lagrangian Schemes Applied to PowerPoint Presentation

SLIDE 1

Overview Boundary treatment Analysis Linear solvers Conclusions

Boundary Approximations for Semi-Lagrangian Schemes Applied to Hamilton-Jacobi-Bellman Equations

Christoph Reisinger⋆

Joint work with Julen Rotaetxe Arto⋆

⋆Mathematical Institute

University of Oxford

Linz, 23 November 2016

SLIDE 2

Overview Boundary treatment Analysis Linear solvers Conclusions

MOTIVATION

Background:

◮ Monotone schemes have “wide stencils” (Motzkin and Wasow, 1952). ◮ They “overstep” domain boundaries. ◮ Standard iterative (multilevel) schemes perform poorly.

We propose a(n)

◮ boundary treatment with same local accuracy as in interior; ◮ unconditionally stable and monotone implicit scheme; ◮ analysis of a convergent aggregation-based multigrid scheme.

SLIDE 3

Overview Boundary treatment Analysis Linear solvers Conclusions

HJB EQUATIONS

Consider ut − inf

α∈A {Lα[u](t, x) + f α(t, x)} = 0,

(t, x) ∈ (0, T] × Ω, u(0, x) = g(x), x ∈ ¯ Ω, u(t, x) = ψ(x), (t, x) ∈ (0, T] × ∂Ω, where Ω is a bounded domain, ¯ Ω := Ω ∪ ∂Ω ⊆ Rd, A is a compact set, and Lα[u](t, x) = tr[aα(t, x)D2u(t, x)] + bα(t, x)Du(t, x), aα = 1

2σασα,T with σα ∈ Rd×P, bα, f α, g and ψ take values in Rd, R, R, and R.

SLIDE 4

Overview Boundary treatment Analysis Linear solvers Conclusions

SEMI-LAGRANGIAN SCHEMES

Define the Linear Interpolation Semi-Lagrangian (LISL) scheme by

Lα

∆x[I∆xφ](t, x) := M p=1 (I∆xφ)(t,x+yα,+

p

(t,x))−2(I∆xφ)(t,x)+(I∆xφ)(t,x+yα,−

p

(t,x)) 2∆x

,

as in Debrabant and Jakobsen (2013), where I is bilinear interpolation and yα,±

p

= ± √ ∆xσα

p for p ≤ P,

yα,±

P+1 = ∆xbα,

and M = P + 1. We consider

◮ an explicit Euler scheme:

Vn+1 − Vn ∆t − inf

α∈A Lα ∆x[I∆xVn] = 0;

◮ an implicit Euler scheme:

Vn+1 − Vn ∆t − inf

α∈A Lα ∆x[I∆xVn+1] = 0.

In practice, we replace A by A with |A| = Nα.

SLIDE 5

Overview Boundary treatment Analysis Linear solvers Conclusions

WIDE STENCILS AND BOUNDARIES

The stencil “oversteps” in a layer of width k = √ ∆x around the boundary.

Figure: Truncation and extrapolation of the stencil for an elliptical domain and a

mesh made of square cells.

Ω ∂Ω

There are two situations where interpolation at x + yα,±

p

(t, x) is not possible:

A. x + yα,±

p

(t, x) / ∈ ¯ Ω (bottom left in Fig.);

B. x + yα,±

p

(t, x) ∈ ¯ Ω, but “its” element has vertices outside ¯ Ω (top right).

SLIDE 6

Overview Boundary treatment Analysis Linear solvers Conclusions

DIFFERENT TREATMENTS

Ω ∂Ω

◮ Constant extrapolation in stencil direction. ◮ Linear extrapolation in stencil direction. ◮ Stencil truncation: Consider instead

ˆ Lα

∆x[I∆xφ](t, x) := M

p=1

Aα

p (I∆xφ)(t, x + ˆ

yα,+

p

(t, x)) − (Aα

p + Bα p )φ(t, x) + Bα p (I∆xφ)(t, x + ˆ

yα,−

p

(t, x)) 2∆x , where ˆ yα,±

p

= µα,±

p

yα,±

p

, where µα,±

p

(t, x) = min

µ ≥ 0 : x + µyα,±

p

(t, x) ∈ ∂Ω

.

and Aα

p ≡ Aα p (t, x) and Bα p ≡ Bα p (t, x), such that

ˆ Lα

∆x a consistent approximation as ∆x → 0.

SLIDE 7

Overview Boundary treatment Analysis Linear solvers Conclusions

THE TEST PROBLEM

Problem A (see Section 9.3 from Debrabant and Jacobsen (2013)). It has exact solution u(t, x1, x2) = 3 2 − t

sin x1 sin x2,

and coefficients and control set are given by f α = 1 2 − t

sin x1 sin x2 +

3 2 − t cos2 x1 sin2 x2 + sin2 x1 cos2 x2+ − 2 sin(x1 + x2) cos(x1 + x2) cos x1 cos x2

,

bα = α, σα = √ 2 sin(x1 + x2) cos(x1 + x2)

,

A = {α ∈ R2 : α2

1 + α2 2 = 1}.

◮ The solution is periodic, and D. & J. (2013) use periodic boundary

conditions in space for (t, x1, x2) ∈ [0, T] × [−π, π]2 with T = 1

2.

◮ Here, we use Dirichlet conditions.

SLIDE 8

Overview Boundary treatment Analysis Linear solvers Conclusions

TRUNCATED STENCILS

Figure: Problem A:

Left the stencil over a 11 × 11 grid and Nα = 10 equally spaced points in A. Right the histograms of the displacement from the central node for 641 × 641 grid.

4
3
2
1

1 2 3 4 x1

4
3
2
1

1 2 3 4 x2 Stencil A = 2.29 B = 1.28 C = 1.00 A = 2.14 B = 1.26 C = 1.00 A = 1.00 B = 1.00 C = 1.00 A = 1.00 B = 1.00 C = 1.00 A = 1.18 B = 1.71 C = 1.00 A = 1.26 B = 2.14 C = 1.00

(a) The finite difference weights are

A ≡ Aα

2,1(x), B ≡ Bα 2,1(x) and

C ≡ (µα

2,2(x))−1, Nα = 10.

(b) The radius of the stencil in σα is 14.27

for this grid, given by σα2

√ ∆x =

640/π.

SLIDE 9

Overview Boundary treatment Analysis Linear solvers Conclusions

CONSTANT AND LINEAR EXTRAPOLATION

Table: Explicit Euler LISL method with Nα = 40 for Problem A. (a) Constant extrapolation: L∞ error over [−π, π)2

Nx ∆t = T ∆t = ∆x

4

∆t = ∆x

3 2

∆t = ∆x2 error rate error rate error rate error rate 41 1.36e+00

3.68e-01
3.72e-01
3.65e-01
81

1.89e+00

0.48

2.61e-01 0.49 2.62e-01 0.51 2.60e-01 0.49 161 2.67e+00

0.49

1.80e-01 0.54 1.80e-01 0.54 1.80e-01 0.53 321 3.77e+00

0.50

1.27e-01 0.51 1.27e-01 0.51 1.27e-01 0.51 641 5.34e+00

0.50

9.18e-02 0.47 9.18e-02 0.47 9.18e-02 0.46

(b) Linear extrapolation: L∞ error over [−π, π)2

Nx ∆t = T ∆t = ∆x

4

∆t = ∆x

3 2

∆t = ∆x2 error rate error rate error rate error rate 41 1.59e-01

1.04e-01
1.05e-01
1.03e-01
81

8.15e-02 0.96 5.25e-02 0.99 5.26e-02 1.00 5.22e-02 0.98 161 4.28e-02 0.93 5.62e-01

3.42

5.63e-01

3.42

5.58e-01

3.42

321 2.75e-02 0.64 4.41e+03

12.94

6.00e+03

13.38

8.00e+03

13.81

641 1.85e-02 0.57 2.77e+20

55.80

2.70e+20

55.32

1.37e+21

57.25

SLIDE 10

Overview Boundary treatment Analysis Linear solvers Conclusions

STENCIL TRUNCATION ON TWO DOMAINS

Table: Explicit Euler LISL method for Problem A. (a) L∞ error over [−π, π)2

Nx ∆t = T ∆t = ∆x

4

∆t = ∆x

3 2

∆t = ∆x2 error rate error rate error rate error rate 41 1.42e-01

4.39e-02
4.39e-02
4.36e-02
81

1.04e-01 0.45 2.12e-02 1.05 2.11e-02 1.06 2.11e-02 1.05 161 7.36e-02 0.50 1.10e-02 0.94 1.10e-02 0.94 1.10e-02 0.94 321 5.28e-02 0.48 1.34e+23

83.33

5.77e-03 0.93 5.76e-03 0.93 641 3.77e-02 0.48 5.07e+89

221.17

3.10e-03 0.90 3.10e-03 0.89

(b) L∞ error over [ −π

8 , 15π 8 )2

Nx ∆t = T ∆t = ∆x

4

∆t = ∆x

3 2

∆t = ∆x2 error rate error rate error rate error rate 41 1.55e-01

4.71e-02
4.76e-02
4.67e-02
81

1.12e-01 0.47 1.57e+05

21.67

7.90e+05

23.98

2.11e-02 1.15 161 8.04e-02 0.47 1.02e+33

92.39

1.30e+35

97.06

1.10e-02 0.94 321 5.80e-02 0.47 6.73e+103

235.26

5.96e+138

344.35

5.76e-03 0.93 641 4.22e-02 0.46 8.17e+276

574.97

NaN NaN 3.10e-03 0.89

SLIDE 11

Overview Boundary treatment Analysis Linear solvers Conclusions

EMPIRICAL FINDINGS AND CONSISTENCY

◮ Constant extrapolation loses accuracy near boundary. ◮ Linear extrapolation is unstable. ◮ Stencil truncation has stricter CFL condition than in interior. ◮ CFL condition depends on boundary.

The local LISL truncation error is

◮

O(∆x) if neither ˆ yα,+

p

nor ˆ yα,−

p

verstep;

◮

O(∆x1/2) if either ˆ yα,+

p

r ˆ

yα,−

p

verstep;

◮

O(1) if both ˆ yα,+

p

and ˆ yα,−

p

verstep (cf. Max Jensen’s talk).

(NB: The error is actually O(∆x) in the last case if the exact Dirichlet data are used.)

SLIDE 12

Overview Boundary treatment Analysis Linear solvers Conclusions

LOCAL BOUNDARY REFINEMENT

Define (with N(x) the corners of the element x is in) Ω(1)

∆x

:= {x ∈ Ω∆x : N(x ± σ √ ∆x) ⊂ Ω∆x}, Ω(2)

∆x

:=

x∈Ω(1)

∆x

N(x ± σ √ ∆x) ∩

Ω∆x\Ω(1)

∆x

,

Ω(3)

∆x

:= Ω∆x\(Ω(1)

∆x ∪ Ω(2) ∆x).

Now refine Ω(i)

∆x with mesh size ∆x(i) and step k(i),

∆x(1) = O(∆x3/2), k(1) = O(∆x), ∆x(2) = O(∆x3/2), k(2) = √ ∆x, ∆x(3) = ∆x, k(3) = √ ∆x. Then the local error is O(∆x) everywhere, and the total # of points O(|Ω∆x|). Practically not necessary so not done in the following computations.

SLIDE 13

Overview Boundary treatment Analysis Linear solvers Conclusions

STABILITY

Recall ˆ Lα

∆x[I∆xφ](t, x) := M

p=1

Aα

p (I∆xφ)(t, x + ˆ

yα,+

p

(t, x)) − (Aα

p + Bα p )φ(t, x) + Bα p (I∆xφ)(t, x + ˆ

yα,−

p

(t, x)) 2∆x . The θ-scheme (θ = 0 explicit, θ = 1 implicit) is monotone if (1 − θ)∆tn  

M

p=1

Aα

p + Bα p

2∆x   ≤ 1. This implies that for 0 ≤ θ < 1 the scheme is monotone (− → ℓ∞−stable) if

◮

∆t ≤ C∆x and neither ˆ yα,+

p

nor ˆ yα,−

p

verstep;

◮

∆t ≤ C∆x3/2 and either ˆ yα,+

p

r ˆ

yα,−

p

verstep;

◮

∆t ≤ C∆x2 and both ˆ yα,+

p

and ˆ yα,−

p

verstep.

SLIDE 14

Overview Boundary treatment Analysis Linear solvers Conclusions

IMPLICIT SCHEME

Table: Stencil truncation for explicit and implicit scheme. (a) L∞ error explicit scheme over [−π, π)2

Nx ∆t = T ∆t = ∆x

4

∆t = ∆x

3 2

∆t = ∆x2 error rate error rate error rate error rate 41 1.42e-01

4.39e-02
4.39e-02
4.36e-02
81

1.04e-01 0.45 2.12e-02 1.05 2.11e-02 1.06 2.11e-02 1.05 161 7.36e-02 0.50 1.10e-02 0.94 1.10e-02 0.94 1.10e-02 0.94 321 5.28e-02 0.48 1.34e+23 -83.33 5.77e-03 0.93 5.76e-03 0.93 641 3.77e-02 0.48 5.07e+89 -221.17 3.10e-03 0.90 3.10e-03 0.89

(b) L∞ error implicit scheme over [−π, π)2

Nx ∆t = T ∆t = ∆x

4

∆t = ∆x

3 2

∆t = ∆x2 error rate error rate error rate error rate 41 3.25e-02

4.21e-02
4.17e-02
4.24e-02
81

1.59e-02 1.03 2.08e-02 1.02 2.08e-02 1.01 2.09e-02 1.02 161 8.39e-03 0.92 1.09e-02 0.93 1.09e-02 0.93 1.10e-02 0.93 321 4.38e-03 0.94 5.75e-03 0.93 5.75e-03 0.93 5.76e-03 0.93 641 2.37e-03 0.89 3.09e-03 0.89 3.10e-03 0.89 3.10e-03 0.89

SLIDE 15

Overview Boundary treatment Analysis Linear solvers Conclusions

SOLUTION OF IMPLICIT SYSTEMS

The scheme is of positive type and can be written as max

α∈A

  Bα,n,n

j,j

Un

j −

i=j

Bα,n,n

j,i

Un

i − N

i=1

Bα,n,n−1

j,i

Un−1

i

− Fα,n

j

   = 0, where Bα,n,n

j,j

= 1 + θ∆tn M 2∆x − lα,n

j,j

, Bα,n,n

j,i

= θ∆tn lα,n

j,i , Bα,n,n−1 j,j

= 1, Bα,n,n−1

j,i

= 0, and lα,n

j,i

=

M

p=1

wi(xj + yα,+

p

(tn, xj)) + wi(xj + yα,−

p

(tn, xj)) 2∆x , in matrix form max

α∈A (Aα i X − Fα i ) = 0,

i = 1, . . . , N, where Aα

i,j = Bα,n,n i,j

, Fα

i = Fα,n i

, and X = (Xi) = (Un

i ).

SLIDE 16

Overview Boundary treatment Analysis Linear solvers Conclusions

NONLINEAR AND LINEAR SOLVERS

We solve max

α∈A (Aαi i X − Fαi i ) = 0,

i = 1, . . . , N, by policy iteration (Howard’s algorithm, semi-smooth Newton method,...). For some x(0), iterate for k = 0, 1, 2, ...:

1. α(k+1)

i

∈ arg supαi

Aαi

i X(k) − Fαi i

;
2. solve Aα(k+1)X(k+1) = Fα(k+1).

Convergence, e.g., in Bokanowski, Maroso and Zidani (2009). See also Iain Smears’ talk. In the following, we focus on second step, solution of linear system for LISL.

SLIDE 17

Overview Boundary treatment Analysis Linear solvers Conclusions

TEST PROBLEM AND METHODS USED

We study the LISL scheme on a linear problem, specifically, for the operator ∇σσT∇ with

◮ σ = 2I2:

interpolation on even levels;

◮ σ =

√ 5I2: interpolation on all levels. We test different iterative solvers, including:

◮ Geometric multigrid (GMG), V and W cycles; ◮ Algebraic multigrid (AMG), Ruge and St¨

uben (1987);

◮ Aggregation-based multigrid (AGMG), Notay (2012), implementation

from http://homepages.ulb.ac.be/ynotay/AGMG/

◮ Bi-CG-stab (van der Vorst, 1992) with and without ILU preconditioning

(as in Ma and Forsyth, 2016).

SLIDE 18

Overview Boundary treatment Analysis Linear solvers Conclusions

RESIDUAL REDUCTION FACTORS

Table: Refinement level l for the two-dimensional Laplacian; the length of the stencil m. (a) σ = 2I2

ℓ m GMG V(1,1) GMG W(1,1) AMG AGMG BICGSTAB BICGSTAB with ILU(0) 6 16 0.4193 0.4204 0.0415 0.1015 0.2502 0.0151 7 22 0.2612 0.2666 0.0633 0.1502 0.5158 0.0767 8 32 0.7561 0.7564 0.0981 0.1551 0.5763 0.1234 9 45 0.5076 0.4905 0.1216 0.1858 0.7109 0.2392 10 64 0.8823 0.8841 0.1219 0.2001 0.7621 0.3382

(b) σ =

√ 5I2

ℓ m GMG V(1,1) GMG W(1,1) AMG AGMG BICGSTAB BICGSTAB with ILU(0) 6 17 0.2999 0.2857 0.0403 0.1162 0.3718 0.0262 7 25 0.2565 0.2568 0.0532 0.1367 0.4408 0.0302 8 35 0.4348 0.4300 0.0740 0.1656 0.5938 0.1302 9 50 0.4480 0.4314 0.1124 0.2030 0.6751 0.1746 10 71 0.5547 0.4799 0.1272 0.1992 0.7478 0.3374

SLIDE 19

Overview Boundary treatment Analysis Linear solvers Conclusions

SMOOTHING FACTORS

−4 −2 2 4 −4 −2 2 4 0.1 0.2 0.3 0.4 0.5 θ1 Smoothing factor for Gauss−Seidel high frequencies θ2 f(θ1, θ2)

(c) Smoothing factor standard FD scheme.

−4 −2 2 4 −4 −2 2 4 0.2 0.4 0.6 0.8 1 θ1 Smoothing factor for Gauss−Seidel high frequencies θ2 f(θ1, θ2)

(d) LISL scheme with parameters

((m1 = 9, m2 = 3), (γ1 = 0.5, γ2 = 1)).

Figure: Smoothing factor for high frequencies, i.e. θ ∈ [−π, π]2 \ [−π/2, π/2]2, for the

Gauss-Seidel iteration. The maxima are 0.49 for FD and 0.95 for LISL.

SLIDE 20

Overview Boundary treatment Analysis Linear solvers Conclusions

ONE DIMENSION

◮ Let

m :=

σ

√ ∆x

,

and γ := (m + 1) − σ √ ∆x ,

◮ where m ∈ N denotes the stencil length and ◮ γ ∈ [0, 1] is the interpolation weight, such that ◮ I∆x(φ)(xi +

√ ∆xσ) = γφ(xi + m∆x) + (1 − γ)φ(xi + (m + 1)∆x). Assume that ∆x = 1. Denote the N × N Laplacian matrix

LN :=        2 −1 · · · −1 2 −1 · · · −1 2 · · · . . . . . . . . . ... . . . · · · 2        .

Let now m = 2 and ∆x = 1, then the LISL matrix is

LN,m,γ

SL

:=        2 −γ −1 + γ · · · 2 −γ −1 + γ · · · −γ 2 −γ · · · . . . . . . . . . . . . . . . ... . . . · · · 2        = γLm

N + (1 − γ)Lm+1 N

.

SLIDE 21

Overview Boundary treatment Analysis Linear solvers Conclusions

EIGENVALUES AND EIGENVECTORS IN 1D

5 10 15 20 25 30 35 0.5 1 1.5 2 2.5 3 3.5 4 Eigenvalues λ(L31,5, 0.346

SL

) λ(L6

31)

λ(L5

31)

λ(L31)

(a) Eigenvalues of L31,5,0.346

SL

, L5

31, L6 31, L31.

5 10 15 20 25 30

0.1

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Eigenvectors for eigenvalue 1 V(L31,5, 0.346

SL

), λ1=0.23 V(L6

31), λ1=0.20

V(L5

31), λ1=0.15

V(L31), λ1=0.01

(b) Eigenvectors with the smallest

eigenvalue for L31,5,0.346

SL

, L5

31, L6 31 and L31.

5 10 15 20 25 30

0.6
0.4
0.2

0.2 0.4 0.6 0.8 Eigenvectors for eigenvalue 15 V(L31,5, 0.346

SL

), λ15=1.97 V(L6

31), λ15=2.00

V(L5

31), λ15=1.55

V(L31), λ15=1.80

(c) Eigenvectors for the 15-th eigenvalue

for L31,5,0.346

SL

, L5

31, L6 31 and L31.

5 10 15 20 25 30

0.6
0.4
0.2

0.2 0.4 0.6 0.8 Eigenvectors for eigenvalue 31 V(L31,5, 0.346

SL

), λ31=3.73 V(L6

31), λ31=3.80

V(L5

31), λ31=3.85

V(L31), λ31=3.99

(d) Eigenvectors for the largest

eigenvalue for L31,5,0.346

SL

, L5

31, L6 31 and L31.

SLIDE 22

Overview Boundary treatment Analysis Linear solvers Conclusions

FURTHER TESTS

The grid complexity cG for the total number of variables N., is cG = 1 N1

nlevels

ℓ=1

Nℓ. The algebraic complexity cA is cA = 1 nnz(A1)

nlevels

ℓ=1

nnz(Aℓ).

Table: Grid and algebraic complexities, for the two-dimensional case. (a) Complexities for σ = 2I2

ℓ GMG AMG AGMG cG cA cG cA cG cA 6 1.31 3.66 1.75 1.74 1.24 1.18 7 1.32 2.69 1.76 6.92 1.26 1.26 8 1.33 3.90 1.72 2.05 1.33 1.28 9 1.33 2.75 1.71 11.79 1.25 1.31 10 1.33 3.97 1.70 2.16 1.25 1.23

(b) Complexities for σ =

√ 5I2

GMG AMG AGMG cG cA cG cA cG cA 1.31 2.59 1.67 3.78 1.18 1.10 1.32 2.69 1.74 6.48 1.24 1.22 1.33 2.65 1.71 7.47 1.20 1.22 1.33 2.76 1.69 9.34 1.22 1.30 1.33 2.67 1.64 7.39 1.32 1.45

SLIDE 23

Overview Boundary treatment Analysis Linear solvers Conclusions

PROPERTIES OF LISL MATRIX

Notay (2012) proves convergence of AGMG for matrices with non-negative row and column sum.

◮ Let A be the LISL discretization matrix for a given time step, ◮ on an equispaced Cartesian grid Ω∆x of Ω ⊂ Rd with ∆x > 0, and ◮ a control vector (αi), αi ∈ A, associated with mesh points (xi).

The row sum is non-negative by construction. Assume further, for β ∈ (0, 1], η ∈ 1

2, 1

,

σαi

p (s, xi) − σαl p (s, xl)∞ ≤ Lσxi − xlη ∞,

bαi(s, xi) − bαl(s, xl)∞ ≤ Lbxi − xlβ

∞.

Then the column sum of the matrix is non-negative provided ∆t ≤ ∆x (M − 1)(P + 1), where

◮ P is the number of SL directions, usually P = d; ◮ M = 3d for sufficiently small ∆x.

SLIDE 24

Overview Boundary treatment Analysis Linear solvers Conclusions

SOLVER CPU TIMES

103 104 105 106 System size 10-1 100 101 102 103 104 Elapsed time (sec) Solve time vs system size for Problem A

Direct: O(Na), a ≈2.80 RSAMG: O(Na), a ≈1.46 AGMG: O(Na), a ≈1.27

(e) Time on solver for single step.

103 104 105 106 System size 10-2 10-1 100 101 102 103 104 Elapsed time (sec) Solve time vs system size for Problem A

Direct: O(Na), a ≈2.33 RSAMG: O(Na), a ≈1.07 AGMG: O(Na), a ≈1.08

(f) Average time for ∆t = ∆x. Figure: Average number of seconds per time step for solving the linear systems

versus the size of the systems. We use equispaced Cartesian grids in space with 81, 161, 321 and 641 nodes per dimension.

SLIDE 25

Overview Boundary treatment Analysis Linear solvers Conclusions

CONCLUSION

Analyzed LISL schemes on bounded domains from an applications perspective. The necessary non-locality leads to complications.

◮ A modification of the scheme is needed close to the boundary. ◮ It is possible to maintain the same order of accuracy. ◮ More restrictive CFL conditions. ◮ Standard multigrid schemes struggle. AGMG performs best. ◮ Open question:

◮ provable convergence orders.