Use of Artificial Boundary Conditions for Schr odinger equation by - - PowerPoint PPT Presentation

Use of Artificial Boundary Conditions for Schr¨

• dinger equation

by

Christophe Besse Collaborators: X. Antoine, S. Descombes, P. Klein, J. Szeftel,

• F. Xing

Institut de Math´ ematiques de Toulouse, Universit´ e Toulouse 3, CNRS

Colloque Couplages Num´ eriques

INSTITUT de MATHEMATIQUES

de TOULOUSE

Outline of the talk

1 Introduction 2 Derivation of TBC for the linear Wave Eq. 3 Derivation of ABC for the 1D Schr¨

• dinger Eq.

4 Numerical approximation of ABCs 5 Schwarz Waveform Relaxation method

2016/09 – Nice

Introduction

Typical equations

The Schr¨

• dinger Eq. in Rd

(S)          i∂tψ + ∆ψ + V (x, t, ψ) ψ = 0, (x, t) ∈ Rd × [0; T] lim

|x|→+∞ ψ(x, t) = 0,

t ∈ [0; T] ψ(x, 0) = ψ0(x), x ∈ Rd ψ(x, t): wave function, complex V can be → real potential, V = V (x, t) ∈ C∞(Rd × R+, R) → nonlinear functional, V = f(ψ) Ex: V = q|ψ|2 → a combination: V = V (x) + f(ψ) ψ0 compact support in Ω

2016/09 – Nice

Introduction

Typical equations

The Wave Eq. in R

(W)          ∂2

t u − ∂2 xu = 0,

(x, t) ∈ R × [0; T] lim

|x|→+∞ u(x, t) = 0,

t ∈ [0; T] u(x, 0) = u0(x), ∂tu(x, 0) = u1(x) x ∈ R u(x, t): real function, u0, u1 compact support in Ω

2016/09 – Nice

Domain trunction

Problem : Mesh an unbounded domain (here in 1D) Rd × [0; T]

2016/09 – Nice

Domain trunction

Problem : Mesh an unbounded domain (here in 1D) Rd × [0; T] Truncation R × [0; T] − → ΩT := ]xℓ, xr[×[0; T] Introduction of a fictitious boundary ΣT := Σ × [0; T] with Σ := ∂Ω = {xl, xr} ⇒ BC on the boundaryΣT

2016/09 – Nice

Domain trunction

Problem : Mesh an unbounded domain (here in 1D) Rd × [0; T] Truncation R × [0; T] − → ΩT := ]xℓ, xr[×[0; T] Introduction of a fictitious boundary ΣT := Σ × [0; T] with Σ := ∂Ω = {xl, xr} ⇒ BC on the boundaryΣT Expression of the boundary condtion with the help of the Dirichlet-to-Neumann map: ∂nψ + iΛ+ψ = 0,

• n ΣT .

2016/09 – Nice

Domain trunction

Problem : Mesh an unbounded domain (here in 1D) Rd × [0; T] Truncation R × [0; T] − → ΩT := ]xℓ, xr[×[0; T] Introduction of a fictitious boundary ΣT := Σ × [0; T] with Σ := ∂Ω = {xl, xr} ⇒ BC on the boundaryΣT Expression of the boundary condtion with the help of the Dirichlet-to-Neumann map: ∂nψ + iΛ+ψ = 0,

• n ΣT .

2016/09 – Nice

What happen if we do not take BC into account

Potential V (x) = x Initial datum: gaussian function ψ0(x) = e−x2+10ix Domain: Ω = [−5; 15]

−5 5 10 15 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 x | u(x,t) |

Evolution

• f

|ψ| w.r.t time ∆t = 0.2

t

Exact Solution

x t −5 5 10 15 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

x

2016/09 – Nice

What happen if we do not take BC into account

Potential V (x) = x Initial datum: gaussian function ψ0(x) = e−x2+10ix Domain: Ω = [−5; 15]

−5 5 10 15 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 x | u(x,t) |

Evolution

• f

|ψ| w.r.t time ∆t = 0.2 Homogeneous Dirichlet Boundary Conditions: ψ|Σ = 0 ⇒ Parasistic reflexion

t

Approximated numerical solution

x t −5 5 10 15 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

x

2016/09 – Nice

Domain decomposition

2016/09 – Nice

Domain decomposition

Global Domain Decomposition in (here 3) subdomains. Iterative procedure: Schwarz waveform relaxation method

2016/09 – Nice

Domain decomposition

Bj = ∂nj + S Informations exchanges from domain 2 → domain 1        L uk+1

1

= 0, (t, x) ∈ (0, T) × (a1, b1), uk+1

1

(0, x) = ψ0(x), x ∈ (a1, b1), ∂n1uk+1

1

= 0, x = a1, B1uk+1

1

= B1uk

2,

x = b1.

2016/09 – Nice

Domain decomposition

Bj = ∂nj + S Informations exchanges from domain 1 → domain 2 and domain 3 → domain 2        L uk+1

2

= 0, (t, x) ∈ (0, T) × (a2, b2), uk+1

2

(0, x) = ψ0(x), x ∈ (a2, b2), B2uk+1

2

= B2uk

1,

x = a2, B2uk+1

2

= B2uk

3,

x = b2.

2016/09 – Nice

Domain decomposition

Bj = ∂nj + S Informations exchanges from domain 2 → domain 3        L uk+1

3

= 0, (t, x) ∈ (0, T) × (a3, b3), uk+1

3

(0, x) = ψ0(x), x ∈ (a3, b3), B3uk+1

3

= B3uk

2,

x = a3, ∂n3uk+1

3

= 0, x = a3.

2016/09 – Nice

Domain decomposition

Examples of operators S for Schr¨

• dinger equation

S = −ip, p real parameter (Halpern, Szeftel ’09-10) S =

• −i∂t − V (a, bj), (Halpern, Szeftel ’09-10)

S: Dirichlet-to-Neumann map iΛ+ (X. Antoine, E. Lorin, A.D. Bandrauk & F. Xing, CB)

2016/09 – Nice

Outline of the talk

1 Introduction 2 Derivation of TBC for the linear Wave Eq. 3 Derivation of ABC for the 1D Schr¨

• dinger Eq.

4 Numerical approximation of ABCs 5 Schwarz Waveform Relaxation method

2016/09 – Nice

TBC for Wave Eq.

Homogeneous Wave Equation in 1D

(P)    ∂2

t ψ − ∂2 xψ = 0 ,

(x, t) ∈ Rx × [0; T] , ψ(x, 0) = ψ0(x) , x ∈ Rx , ∂tψ(x, 0) = ψ1(x) , x ∈ Rx . Hypothesis : supp(ψ0,1) ⊂ B(0, R). Solution ψ(x, t) = 1 2(ψ0(x + t) + ψ0(x − t)) + 1 2 x+t

x−t

ψ1(y)dy. Simply, ψ(x, t) = ϕ1(x + t) + ϕ2(x − t) , with 2ϕ1(x) = ψ0(x) + x

−∞

ψ1(y)dy , 2ϕ2(x) = ψ0(x) − x

−∞

ψ1(y)dy .

2016/09 – Nice

TBC for Wave Eq.

R −R t x L −L supp(ϕ(x + t)) supp(ψ(x − t))

∀t > 0 and ∀|x| > R, the supports are disjoint. x = −L x = L ψ(x, t) = ϕ1(x + t) ψ(x, t) = ϕ2(x − t) ∂xψ = ∂tψ ∂xψ = −∂tψ Then, on x = ±L, the Neumann datum is expressed in function of the Dirichlet one ∂nψ + ∂tψ = 0 where n denotes the outwardly unit normal vector to Ω = (−L, L).

2016/09 – Nice

TBC for Wave Eq.

The problem (P) is thus transformed in (Papp) (Papp)        ∂2

t ψa − ∂2 xψa = 0 ,

(x, t) ∈ Ω × [0; T] , ψa(x, 0) = ψ0(x) , x ∈ Ω , ∂tψa(x, 0) = ψ1(x) , x ∈ Ω , ∂nψa + ∂tψa = 0 , (x, t) ∈ Γ × [0; T] , where Γ = {−L, L}. ψ|Ω(x, t) = ψa(x, t) , (x, t) ∈ Ω × [0; T] . Ths BCs do not perturb the solution: we call them Transparent Boundary Conditions (TBC) Remark: ∂2

t − ∂2 x = (∂t − ∂x)(∂t + ∂x).

2016/09 – Nice

TBC for Wave Eq.

Other way of derivation: Laplace transform L (u)(x, ω) = ˆ u(x, ω) = ∞ u(x, t)e−ωtdt with frequency ω = σ + iτ, σ > 0 L (∂tu)(x, ω) = ωˆ u(x, ω) − u(x, 0) , L (∂2

t u)(x, ω)

= ω2ˆ u(x, ω) − ωu(x, 0) − ∂tu(x, 0) .

1 Transmission problem

Interior Problem      (∂2

t − ∂2 x)v = 0,

x ∈ Ω, t > 0, ∂xv = ∂xw, x ∈ Γ, t > 0, v(x, 0) = ψ0(x), x ∈ Ω, ∂tv(x, 0) = ψ1(x), x ∈ Ω, Exterior problem        (∂t − ∂2

x)w = 0,

x ∈ Ω, t > 0, w(x, t) = v(x, t), x = ±L, t > 0, w(x, 0) = 0, x ∈ Ω, ∂tw(x, 0) = 0, x ∈ Ω.

2016/09 – Nice

TBC for Wave Eq.

2 Laplace tranform w.r.t time for x > L : ω2 ˆ

w − ∂2

x ˆ

w = 0. Solution: ˆ w(x, ω) = A+(ω)eω x + A−(ω)e−ω x

3 How to select the outgoing wave: solution of finite energy −

→ A+ = 0. ˆ w(x, ω) = e−ω (x−L)L (w(L, ·))(ω). Taking derivative ∂x ˆ w(x, ω)|x=L = −ω ˆ w(x, ω)|x=L

4 Inverse Laplace tranform

∂xw(x, t)|x=L = −∂tw(x, t)|x=L.

5 Thus, we obtain the TBC for v

∂nv + ∂tv = 0

2016/09 – Nice

Outline of the talk

1 Introduction 2 Derivation of TBC for the linear Wave Eq. 3 Derivation of ABC for the 1D Schr¨

• dinger Eq.

4 Numerical approximation of ABCs 5 Schwarz Waveform Relaxation method

2016/09 – Nice

The 1D Eq. without potential

1 Splitting between interior and exterior problems

Interior Problem      (i∂t + ∂2

x)v = 0,

x ∈ Ω, t > 0, ∂xv = ∂xw, x ∈ Σ, t > 0, v(x, 0) = ψ0(x), x ∈ Ω. Exterior problem            (i∂t + ∂2

x)w = 0,

x ∈ Ω, t > 0, w(x, t) = v(x, t), x = xl,r, t > 0, lim

|x|→+∞ w(x, t) = 0,

t > 0, w(x, 0) = 0, x ∈ Ω.

• problem

interior problem left exterior problem right exterior

(x,t)

x x

L R

• utput:

input:

L

(x ,t)

L x(x ,t)

Dirichlet data Neumann data

v

v w 2 Laplace transform w.r.t time on Ωr

i∂tw + ∂2

xw = 0

− → iω ˆ w + ∂2

x ˆ

w = 0 Solution: ˆ w(x, ω) = A+(ω)e

+

√−iω x + A−(ω)e− + √−iω x

with the determination of the square root s.t. Re( + √·) > 0.

2016/09 – Nice

The 1D Eq. without potential

3 Argument: w ∈ L2(Ωr) to select the outgoing wave ⇒ A+ = 0

ˆ w(x, ω) = e− +

√−iω (x−xr)L(w(xr, ·))(ω),

then, taking derivative and using continuity ∂x ˆ w(x, ω)|x=xr = −

+

√ −iω ˆ w(x, ω)|x=xr = −e−iπ/4ω ˆ w(x, ω)|x=xr √ω

4 Inverse Laplace transform

∂xw(x, t)|x=xr = −e−iπ/4∂t 1 √π t w(x, s)|x=xr √t − s ds

• = −e−iπ/4∂1/2

t

w(xr, t)

5 Exact boundary condition on xr: ∂xv(x, t)|x=xr = −e−iπ/4∂1/2

t

v(xr, t) where ∂1/2

t

f(x, t) = 1 √π ∂t t f(x, s)|x=xr √t − s ds fractional derivative

• perator of order 1/2

2016/09 – Nice

The 1D Eq. without potential

TBC on Σ : ∂nv + e−iπ/4∂1/2

t

v = 0,

• n Σ × [0; T].

Remark: ∂2

x + i∂t = (∂n + i√i∂t)(∂n − i√i∂t)

The problem (S) is thus transformed in (Sapp)

The Schr¨

• dinger Eq. in Ω

(Sapp)        i∂tψ + ∆ψ = 0, (x, t) ∈ ΩT ∂nψ + e−iπ/4∂1/2

t

ψ = 0,

• n ΣT ,

ψ(x, 0) = ψ0(x), x ∈ ΩT

2016/09 – Nice

Extension to simple cases

V = 0

Expression of the exact Dirichlet-to-Neumann operator ∂nψ + e−iπ/4∂1/2

t

ψ = 0,

• n ΣT .

(ABC0)

V = Vℓ,r constant at the exterior of Ω

∂nψ + e−iπ/4 eitVℓ,r∂1/2

t

• e−itVℓ,rψ
• = 0,
• n ΣT .

V = V (t) : Gauge change

We set: v(x, t) = ψ(x, t)e−iV(t) with V(t) = t V (s) ds. Then i∂tψ = (i∂tv − V (t) v) eiV(t), where v is solution to the equation without potential ⇒ Exact condition ABC0 for v : ∂nv + e−iπ/4∂1/2

t

v = 0,

• n ΣT

⇒ We come back to ψ : ∂nψ + e−iπ/4 eiV(t)∂1/2

t

• e−iV(t)ψ
• = 0,
• n ΣT

2016/09 – Nice

Case V = V (x, t)

If V = V (x, t), then the Laplace transform can not be used anymore. One have to introduce a new tool: pseudodifferential calculus. A pseudodifferential operator P(x, t, ∂t) is described by its total symbol p(x, t, τ) in the Fourier space (τ is the covariable of t) P(x, t, ∂t) u(x, t) = F −1

t

• p(x, t, τ) ˆ

u(x, τ)

• =
• R

p(x, t, τ) ˆ u(x, τ)eitτdτ Notations: P = Op(p) , p(x, t, τ) = σ(P(x, t, ∂t))

2016/09 – Nice

Case V = V (x, t)

Examples

The fractional operators ∂1/2

t

and Iα/2

t

∂1/2

t

f(t) = 1 √π ∂t t f(s) √t − s ds Iα/2

t

f(t) = 1 Γ(α/2) t (t − s)α/2−1f(s) ds Nonlocal w.r.t time convolution operator Operator ∂t ∂1/2

t

I1/2

t

It ↓ ↓ ↓ ↓ Symbol iτ e−iπ/4√−τ eiπ/4 √−τ 1 iτ Class OPS1 OPS1/2 OPS−1/2 OPS−1

2016/09 – Nice

Artificial boundary conditions V = V (x, t)

Factorization of the operator L (Nirenberg)

General Schr¨

• dinger operator

L = i∂t + ∂2

x + A ∂x + B = (∂x + iΛ−)(∂x + iΛ+) + R

Λ± ∈ OPS1/2 is a ΨDO of order 1/2 in time and 0 in space with σ(Λ+) = λ+ ∼

+∞

• j=0

λ+

1/2−j/2 = λ+ 1/2 + λ+ 0 + λ+ −1/2 + λ+ −1 + · · ·

and λ+

1/2−j/2 homogeneous of order 1/2 − j/2.

Artificial condition : ∂nw + iΛ+w = 0

∂nw + i

+∞

• j=0

Op

• λ+

1/2−j/2

• w = 0,
• n ΣT

Approximated condition of order M:

∂nwM + i

M−1

• j=0

Op

• λ+

1/2−j/2

• wM = 0,
• n ΣT

2016/09 – Nice

Artificial boundary conditions V = V (x, t)

Strategy 1: For reason of symmetry and to get adequate estimates, the ABC of 4th order ABC4

1 is

∂nu + e−iπ/4eiV(x,t)∂1/2

t

• e−iV(x,t)u
• (ABC2

1)

+i sg(∂nV )

• |∂nV |

2 eiV(x,t)It

• |∂nV |

2 e−iV(x,t)u

• = 0

Strategy 2: two choices for principal symbol. The ABCs of 4th order ABC4

0,2

are ∂nu + i √ i∂t + V u (ABC2

2)

+ i 4∂nV (i∂t + V )−1u = 0

• r

∂nu + e−iπ/4∂1/2

t

u (ABC2

0)

+ eiπ/4 V 2 I1/2

t

u (ABC3

0) − i∂nV

4 Itu = 0 (ABC0) are actually high frequency approximation of (ABC2).

2016/09 – Nice

Nonlinear case V = V (x, t, u)

“Simply” replace V (x, t) in the previous ABCs by V = V (x, t, u). We have to solve the equation      i∂tψ + ∂2

xψ + V (x, t, ψ) ψ = 0,

in ΩT ψ(·, 0) = ψ0, in Ω ∂nψ + Sψ = 0,

• n ΣT

Potential strategy l = 0

Order 2 : S = S2

0 = e−i π

4 ∂1/2

t

, Order 3 : S = S3

0 = S2 0 − ei π

4 V (t, ·, ψ)

2 I1/2

t

, Order 4 : S = S4

0 = S3 0 − i∂nV (t, ·, ψ)

4 It,

2016/09 – Nice

Nonlinear case V = V (x, t, u)

Gauge change strategy l = 1

Order 2 : S = S2

1 = e−i π

4 eiV(t,x)∂1/2

t

(e−iV(t,x)·), Order 4 : S = S4

1 = S2 1 − isgn(∂nV )

• |∂nV (t, ·)|

2 eiV(t,x)It(

• |∂nV (t, ·)|

2 e−iV(t,x)·), with V(t, x) = t V (s, x, ψ)ds.

Pad´ e approximant strategy l = 2

Order 2 : S = S2

2 = Op(−i

• −τ + V (t, ·, ψ)),

Order 4 : S = S4

2 = S2 2 + 1

4Op( ∂nV −τ + V (t, ·, ψ)), 2D Case: generalization with use of boundary curvature.

2016/09 – Nice

Outline of the talk

1 Introduction 2 Derivation of TBC for the linear Wave Eq. 3 Derivation of ABC for the 1D Schr¨

• dinger Eq.

4 Numerical approximation of ABCs 5 Schwarz Waveform Relaxation method

2016/09 – Nice

Numerical schemes

Let ∆t = T/N be the time step and let us set tn = n∆t and un stands for an approximation of u(tn). Time approximation Semi-discrete Crank-Nicolson symmetrical scheme i un+1 − un ∆t + ∂2

x(un+1 + un

2 ) + V n+1 + V n 2 un+1 + un 2 = 0, for n = 0, . . . , N − 1. Space approximation Finite differences or Finite Element Method ABCM

0,1: discrete convolutions

ABCM

2 : rational approximation of the square root (Pad´

e)

2016/09 – Nice

Numerical schemes for ABCM

0,1

Main tool (from signal processing): the Z-transform

Definition/Properties

1 Definition: Z{xn} = ˆ

x(z) :=

∞

• n=0

xn z−n, |z| > R.

2 InverseZ-transform: xn =

1 2i π

• C

ˆ x(z) zn dz, with C = C(0, S) and S > R.

3 Relation with Laplace transform if un = u(tn) with u(0) = 0 then

δt Z{u(tn)}(esδt) = L(u)(s) + O(δt2), Re(s) ≥ M > 0.

4 Time shift: Z{un−1} = z−1Z{un}(z), Z{un+1} = zZ{un}(z) − u0

Apply Z-transform to CN scheme on exterior domain (u0 = 0) Z

• i un+1 − un

∆t + ∂2

x(un+1 + un

2 )

• = iz − 1

∆t ˆ u + z + 1 2 ∂2

xˆ

u = 0.

2016/09 – Nice

Numerical schemes for ABCM

0,1

Identify the solution of the ODE: ∂xˆ u(z, x) = −e−iπ/4 2/∆t

• 1 − z−1

1 + z−1 ˆ u and we have to compare to the TBC ∂xu + e−iπ/4∂1/2

t

u = 0 ⇐ ⇒ ∂xun + e−iπ/4Z−1

• 2/∆t
• 1 − z−1

1 + z−1

• ∗n un

Power series expansion: (1 − z−1)1/2 =

∞

• p=0

cpz−p, (1 + z−1)k/3

∞

• p=0

dpz−p with cp+1 = −

1 2 − (p − 1)

p cp, p ≥ 0, c0 = 1 and dp+1 = −

1 2 + (p − 1)

p dp, p ≥ 0, d0 = 1

2016/09 – Nice

Numerical schemes

Numerical scheme for ABCM

0,1: discrete convolutions

Approximations of ∂1/2

t

, I1/2

t

and It in agreement with the Crank-Nicolson scheme ⇒ trapezoidal formula [Schmidt-Yevick (97), Antoine-Besse (03)]

∂1/2

t

f(tn) ≈

• 2

∆t

n

• k=0

βn−kfk I1/2

t

f(tn) ≈

• ∆t

2

n

• k=0

αn−kfk It f(tn) ≈ ∆t 2

n

• k=0

γn−kfk

    

(α0, α1, α2, ...) = (1, 1, 1 2 , 1 2 , 3 8 , 3 8 , ...) βk = (−1)kαk, ∀k ≥ 0, (γ0, γ1, γ2, ...) = (1, 2, 2, ...)

2016/09 – Nice

Numerical schemes

Numerical scheme for ABCM

2

The square root is approximated by Pad´ e approximants of order m: √z ≈ Rm(z) =

m

• k=0

am

k − m

• k=1

am

k dm k

z + dm

k

, with am

0 = 0

, am

k =

1 m cos2

• (2k+1)π

4m

• ,

dm

k = tan2

(2k + 1)π 4m

• .

For the conditions ABCM

2

: √ i∂t + V ❀ Rm(i∂t + V ) ⇒ √ i∂t + V ≈ m

• k=0

am

k

• u −

m

• k=1

am

k dm k (i∂t + V + dm k )−1u

• ϕk

Lindmann’s trick (85) : introduction of auxiliary functions

i∂tϕk + (V + dm

k )ϕk = u,

pour 1 ≤ k ≤ m, in x = xl,r, with ϕk(x, 0) = 0.

2016/09 – Nice

Outline of the talk

1 Introduction 2 Derivation of TBC for the linear Wave Eq. 3 Derivation of ABC for the 1D Schr¨

• dinger Eq.

4 Numerical approximation of ABCs 5 Schwarz Waveform Relaxation method

2016/09 – Nice

SWR method

SWR = Schwarz Waveform Relaxation method Aim: to solve the equation (S)        ∂tu + L u = f (x, t) ∈ Ω × (0; T) BCs on ∂Ω (ex: u=g), t ∈ (0; T) u(x, 0) = u0(x), x ∈ Ω L = −i∆ Subdivision Ω =

N

• j=1

Ωj −

N

• j=1

Ωj Ω = [a, b] = [a0, b0]

x t T a b Ω1 b2 a2 Ω2 Ω3 b3 a3 a1 b1

2016/09 – Nice

SWR method

In 1D: Ωj = [aj, bj], Γj,j−1 = {aj} and Γj,j+1 = {bj}.

SWR method, k ≥ 0

(S)        ∂tuk

i + L uk i = f,

(x, t) ∈ Ωi × (0; T) Biuk

i = Biuk−1 j

, (x, t) ∈ Γi,j × (0; T), j = i ± 1 u(x, 0) = u0(x), x ∈ Ωi If i = 1 or i = N, Γ1,0 = {a0} and ΓN,N+1 = {b0}, then, the BC is uk

0 = g, (x, t) ∈ {a0} × (0; T),

uk

N = g, (x, t) ∈ {b0} × (0; T).

If k = 0, u0

i = 0 (or random) on Γi,j.

2016/09 – Nice

Transmission Conditions

Many choices Biuk

i = uk i : Dirichlet

Biuk

i = ∂nijuk i : Neumann

Biuk

i = ∂nijuk i + p uk i , p ∈ C: Robin

Biuk

i = ∂nijuk i + Siuk i

x t T Ωi ai bi Γi,i+1 ni,i+1 ni+1,i Γi+1,i Ωi+1 bi+1 ai

2016/09 – Nice

Principle of the method

a1 b1 Ω1

• a2

b2 Ω2

• One defines for each subdomain

Domain 1 ∂tuk

1 + L uk 1 = f in Ω1 × (0, T)

uk

1(·, 0) = u0 in Ω1

uk

1(b1, ·) = uk−1 2

(b1, ·) on Γ1 × (0, T) Domain 2 ∂tuk

2 + L uk 2 = f in Ω2 × (0, T)

uk

2(·, 0) = u0 in Ω2

uk

2(a2, ·) = uk−1 1

(a2, ·) on Γ2 × (0, T)

Classical lgorithm

u0

1(b1, t) = 0, u0 2(a2, t) = 0, t ∈ (0, T)

k = 1 While

• (uk

1 − uk−1 1

)(b1, T)2

L2 + (uk 2 − uk−1 2

)(a2, T)2

L2

• > ε do

To solve the systems in Ω1 and Ω2, t ∈ (0, T) Exchange informations by transmission conditions. k = k + 1 No convergence without overlap.

2016/09 – Nice

Principle of the method

Optimal Schwarz method Halpern, Szeftel ’09, (V = Cte) a1 b1 Ω1

• a2

b2 Ω2

• Domain 1

(i∂t + ∂2

x + V )uk 1 = 0 in Ω1 × (0, T)

uk

1(·, 0) = u0 in Ω1

(∂x + S1)uk

1(b1, ·)

= (∂x + S1)uk−1

2

(b1, ·) in (0, T) Domain 2 (i∂t + ∂2

x + V )uk 2 = 0 in Ω2 × (0, T)

uk

2(·, 0) = u0 in Ω2

(−∂x + S2)uk

2(a2, ·)

= (−∂x + S2)uk−1

1

(a2, ·) in (0, T)

Convergence in 2 iterations if and only if Sj = √−i∂t − V For non constant V : optimal operators non available

2016/09 – Nice

Principle of the method

Quasi optimal Schwartz method Halpern, Szeftel ’09, (non constant V ) a1 b1 Ω1

• a2

b2 Ω2

• Domain 1

(i∂t + ∂2

x + V )uk 1 = 0 in Ω1 × (0, T)

uk

1(·, 0) = u0 in Ω1

(∂x + S1)uk

1(b1, ·)

= (∂x + S1)uk−1

2

(b1, ·) in (0, T) Domain 2 (i∂t + ∂2

x + V )uk 2 = 0 in Ω2 × (0, T)

uk

2(·, 0) = u0 in Ω2

(−∂x + S2)uk

2(a2, ·)

= (−∂x + S2)uk−1

1

(a2, ·) in (0, T)

S1 =

• −i∂t − V (b1)

S2 =

• −i∂t − V (a2)

The algorithm converges

2016/09 – Nice

Principle of the method

Complex Robin algorithm Halpern, Szeftel ’09; Approximation of order 0 of

• ptimal algo.

a1 b1 Ω1

• a2

b2 Ω2

• Domain 1

(i∂t + ∂2

x + V )uk 1 = 0 in Ω1 × (0, T)

uk

1(·, 0) = u0 in Ω1

(∂x + S1)uk

1(b1, ·)

= (∂x + S1)uk−1

2

(b1, ·) in (0, T) Domain 2 (i∂t + ∂2

x + V )uk 2 = 0 in Ω2 × (0, T)

uk

2(·, 0) = u0 in Ω2

(−∂x + S2)uk

2(a2, ·)

= (−∂x + S2)uk−1

1

(a2, ·) in (0, T)

Sj = −ip Convergence for any p > 0 Optimization w.r.t p > 0 to improve the convergence

2016/09 – Nice

Possible improvements

Use ABC to build transmission conditions Biuk

i = ∂nijuk i + Siuk i , with S

defined in Part 2. No overlap To improve convergence Switch to nonlinear equation Switch to higher dimensions We work on Ω0 = (a0, b0).    L u := (i∂t + ∂xx + V (x))u = 0, (t, x) ∈ (0, T) × Ω0, u(0, x) = u0(x), x ∈ Ω0, ∂nu(t, x) = 0, x = a0, b0,

2016/09 – Nice

Discretization

Time discretization: Crank-Nicolson, ∆t = T/NT 2ivk

j,n

∆t + ∂xxvk

j,n + Wjvk j,n = 2iuk j,n−1

∆t , where vk

j,n = (uk j,n + uk j,n−1)/2 and Wj,n = (Vj,n + Vj,n−1)/2.

Space discretization: Finite elements → weak formulation; φ ∈ H1((aj, bj)) 2i ∆t bj

aj

vk

j,nφdx −

bj

aj

∂xvk

j,n∂xφ +

bj

aj

Wj,nvk

j,nφdx

+ ∂xvk

j,n(bj)φ(bj) − ∂xvk j,n(aj)φ(aj) = 2i

∆t bj

aj

uk

j,n−1φdx,

Discrete transmission conditions for j = 2, 3, ..., N − 1, ∂njvk

j,n + Svk j,n = ∂njvk−1 j−1,n + Svk−1 j−1,n,

x = aj, ∂njvk

j,n + Svk j,n = ∂njvk−1 j+1,n + Svk−1 j+1,n,

x = bj, where S if the discretization of S. Choice for the talk: S = e−iπ/4∂1/2

t

and Svk

j,n = e−iπ/4

• 2

∆t

n

• s=0

˜ βn−svk

j,s,

with (˜ β0, ˜ β1, ˜ β2, ˜ β3, ˜ β4, ˜ β5, . . . ) = (1, 1, 1

2, 1 2, 1·3 2·4, 1·3 2·4, ...).

2016/09 – Nice

Fluxes at interfaces

Definition: Fluxes at left lk

j,n and right interfaces rk j,n

j = 2, 3, ..., N − 1, n = 1, ..., NT lk

j,n = ∂njvk j,n(aj) + Svk j,n(aj),

rk

j,n = ∂njvk j,n(bj) + Svk j,n(bj),

The transmission conditions can be written thanks to fluxes j = 2, 3, ..., N − 1, lk+1

j,n = −rk j−1,n + 2Svk j−1,n(bj−1),

rk+1

j,n = −lk j+1,n + 2Svk j+1,n(aj+1),

2016/09 – Nice

Linear systems

One discretizes Ωj with Nj nodes + P1 FE approx = ⇒ N linear systems (Aj,n − Bj,n)vk

j,n = 2i

∆tMjuk

j,n−1 + bk j,n − QT j

lk

j,n

rk

j,n

• ,,

with Aj,n =

2i ∆tMj − Sj + Mj,Wn, and boundary terms Bj,n =            Bl

j,n

· · · · · · . . . ... . . . . . . · · · Br

j,n

           ∈ CNj ×Nj , bk

j,n =

         bk,l

j,n

. . . bk,r

j,n

         ∈ CNj , Qj =

• 1

· · · · · · 1

• ∈ C2×Nj .

with Qj,l = 1 · · · ∈ CNj , Qj,r = · · · 1 ∈ CNj .

2016/09 – Nice

Interface problem

For S = S2

0, let’s consider two time steps and two subdomains

Domain Ω1 Time step 1 (A1,1 − B1,1)vk

1,1 = 2i

∆tM1uk

1,0 + bk,r 1,1 − QT 1,rrk,T 1,1

Time step 2 (A1,2 − B1,2)vk

1,2 = 2i

∆tM1uk

1,1 + bk,r 1,2 − QT 1,rrk,T 1,2

But, uj,1 = 2vj,1 − uj,0 which leads to (A1 − β0Q1,r)vk

1,1

= 2i ∆tM1uk

1,0 + cβ1Q1,rvk 1,0 − QT 1,rrk,T 1,1

(A1 − β0Q1,r)vk

1,2 − 4i

∆tM1vk

1,1 = − 2i

∆tM1uk

1,0 + cQ1,r(β2vk 1,0 + β1vk 1,1) − QT 1,rrk,T 1,2

• r again
• A1

− 4i

∆tM1

A1 vk

1,1

vk

1,2

• − c

β0Q1,r β1Q1,r β0Q1,r vk

1,1

vk

1,2

• =

2i ∆t M1uk

1,0

−M1uk

1,0

• + c

β1Q1,rvk

1,0

β2Q1,rvk

1,0

• −

QT

1,rrk,T 1,1

QT

1,rrk,T 1,2

• thus

(A1 − B1)vk

1 = F1 + F1 − QT 1,rgk 1

2016/09 – Nice

Interface problem

a1 b1 r1

• a2

b2 l2

• Let us define two interface vectors

gk

1 =

rk

1,1

rk

1,2

• gk

2 =

lk

2,1

lk

2,2

• We have

rk+1

1

= −lk

2 + 2 ¯

Svk

2(b1)

and lk+1

2

= −rk

1 + 2 ¯

Svk

1(a2)

⇓ ⇓ gk+1

1

= −gk

2 + 2Q2,lB2vk 2 + ¯

F2 gk+1

2

= −gk

1 + 2Q1,rB1vk 1 + ¯

F1. But (A1 − B1)vk

1 = F1 + F1 − QT 1,rgk 1 and (A2 − B2)vk 2 = F2 + F2 − QT 2,lgk 2.

Therefore gk+1

1

= −gk

2 − 2Q2,lB2(A2 − B2)−1QT 2,lgk 2 + d1

gk+1

2

= −gk

1 − 2Q1,rB1(A1 − B1)−1QT 1,rgk 1 + d2

gk+1 = Lgk + d.

2016/09 – Nice

Interface problem

gk+1 = Lgk + d. L =

• −I − 2Q2,lB2(A2 − B2)−1QT

2,l

−I − 2Q1,rB1(A1 − B1)−1QT

1,r

• ,

d = 2 Q2,lB2(A2 − B2)−1(F2 + F2) + Q2,lF2 Q1,rB1(A1 − B1)−1(F1 + F1) + Q1,rF1

• .

2016/09 – Nice

Interface problem: general case

Let interface vectors gk

1 = (rk 1,1, rk 1,2, ..., rk 1,NT )T ∈ CNT , gk N = (lk N,1, lk N,2, ..., lk N,NT )T ∈ CNT ,

gk

j = (lk j,1, lk j,2, ..., lk j,NT , rk j,1, rk j,2, ..., rk j,NT )T ∈ C2NT , j = 2, 3, ..., N − 1,

and the global interface vector gk =

• gk

1

gk

2

· · · gk

N

T ∈ C(2N−2)×NT . We still have gk+1 = Lgk + d, where L is a block matrix

L =                 X2,1 X2,2 X1,4 X3,1 X3,2 X2,3 X2,4 · · · X3,3 X3,4 XN−1,1 XN−1,2 · · · XN,1 XN−1,3 XN−1,4                 ,

with Xj,p ∈ CNT ×NT , j = 1, 2, ..., N, p = 1, 2, 3, 4 and d is a vector

d = d1,r d2,l d2,r · · · dN,l T ∈ C(2N−2)×NT , dj,l, dj,r ∈ CNT .

2016/09 – Nice

General interface problems

The interface problem reads gk+1 = Lgk + d = Rgk, with gk = (gk,T

1

, gk,T

2

, ..., gk,T

N ).

= ⇒ fixed point method to solve (I − L)g = d. We can change by other iterative methods, such as Krylov. One defines the application of I − L to vector g by (I − L)g = I − Rg + d. L is easy to compute

2016/09 – Nice

Use of L

Classical algorithm V1

1 Initialize the iterations by g0

j , j = 1, 2, ..., N.

2 Solve Schr¨

• dinger on each subdomain with gk

j . Costly part

3 Exchange values at interfaces and compute gk+1

j

.

4 Do again steps 2 and 3 until error gk+1 − gk is small.

Reinterpretation of algorithm V1

1 Build d = R · 0 explicitely 2 Define the application of I − L to vector g 3 Solve the linear system (I − L)g = d by fixed point method 4 Solve the Schr¨

• dinger equation on each subdomain for each time step using

the BCs obtained at step 3

5 Repeat steps until error gk+1 − gk small.

We can think to change step 3: solve by Krylov method

2016/09 – Nice

Use of L

New algorithm

1 Build L and b explicitely. 2 Initialize the iterations by g0

j .

3 Solve gk+1 = Lgk + b while gk+1 − gk ≥ ε.

⇐ ⇒ solve (I − L)g = b by fixed point or any Krylov method

4 Solve Schr¨

• dinger equation on each subdomain using the boundary conditions

gk

j .

2016/09 – Nice

Results

Tests made on cluster with 92 nodes (8 cores/Cpu, 16 cores/node, Intel Sandy Bridge E5-2670, 32G/node): 1472 cores. (a0, b0) = (−21, 21), T = 0.5, ∆t = 10−3, ∆x = 10−5, V = −x2, u0 = e−(x+10)2+20i(x+10). Transmission condition: S2

N = 2 103 101 10−1 10−3 10−5 10−7 10−9 10−11 Absolute Residual Number of iterations Fixed point GMRES BiCGStab 1 2 3 103 101 10−1 10−3 10−5 10−7 10−9 10−11 Absolute Residual Number of iterations Fixed point GMRES BiCGStab N = 100 20 40 60 80 100

2016/09 – Nice

Results

Table : Computation times (seconds) of the classical and the new algorithms with the fixed point method and the Krylov methods used on the interface problem. The potential is V = −x2 and the mesh is ∆t = 0.001, ∆x = 10−5.

N 2 10 100 500 1000 T ref 403.56 Fixed point T cls 773.07 2937.77 359.30 284.78

• T new

773.72 178.30 18.19 4.76

• GMRES

T cls 771.82 1888.41 332.61 245.51 267.17 T new 777.42 169.50 19.13 5.91 6.55 BiCGStab T cls 774.19 2760.11 679.72 799.09 845.65 T new 774.44 177.02 18.18 6.83 7.12 −: The maximum number of iterations is fixed to 2000. The residual at iteration 2000 is 2.37 × 10−9 and does not reach the convergence tolerance.

2016/09 – Nice

Comparison of transmission conditions, N = 500, V = −x2

Fixed point Gmres Bicgstab Strategy Niter TLd Ttotal Niter TLd Ttotal Niter TLd Ttotal SM S2 357 0.775 4.68 1023 2.883 6.91 368 1.646 5.51 S3 337 0.734 4.62 977 2.620 6.55 345 1.831 5.77 S4 337 0.733 4.65 978 2.681 6.54 350 1.739 5.73 SM

1

S2

1

341 0.745 4.62 1010 2.364 6.20 353 2.102 6.00 S4

1

340 0.743 4.63 1023 3.454 7.19 351 2.225 6.06 SM

2

S2,20

2

• 1240

3.368 7.34 440 2.626 6.64 S2,50

2

• 997

2.320 6.30 352 2.240 6.16 S2,100

2

336 0.735 4.62 998 3.055 7.03 333 1.603 5.62 S4,20

2

• 1216

3.349 7.31 464 2.044 6.05 S4,50

2

• 1043

3.907 7.85 336 1.756 5.63 S4,100

2

336 0.733 4.60 1024 2.424 6.35 334 1.989 5.95 Robin∗ 1690 3.628 7.52 1060 3.000 6.80 318 1.41 5.32

∗: Robin parameters are p = 45 (fixed point), p = 19 (Gmres) and p = 6 (Bicgstab).

• : the algorithm does not converge before 2000 iterations.

TLd : time to solve the interface problem

2016/09 – Nice

Results

Table : Number of iterations and computation times for some different final times T = 0.5, 1.0, 1.5, 2.0 with N = 500. The potential is V = −x2 and the mesh is ∆t = 0.001, ∆x = 10−5.

T 0.5 1.0 1.5 2.0 Niter 313 593 879 1165 T new 5.91 15.57 34.23 71.2

2016/09 – Nice

Results

If we look for V = V (x, t) or nonlinearities, we still have gk+1 = Rgk = gk − (I − R)gk but now R can be nonlinear. One uses the operator Llin as a precondionner and we solve gk+1 = gk + P −1(I − R)gk with P = I − Llin.

2016/09 – Nice

Results for V = 5tx

Table : Number of iterations required and computation times of the classical algorithm and the preconditioned algorithm. The potential is V = 5tx.

N 10 100 500 1000 ∆t = 0.001, ∆x = 10−5 Nnopc 17 71 349 695 Npc 17 32 31 35 T ref 6496.3 Tnopc 10123.1 3217.0 2466.5 2238.0 Tpc 10128.9 1432.7 250.0 170.7

2016/09 – Nice

Results for V = 5tx

0.99998 1.0 1.00002 Real part

• 0.0002

0.0 0.0002 Imaginary part

I − Lh, n P −1(I − Lh, n)

0.99998 1.0 1.00002 Real part

• 0.0002

0.0 0.0002 Imaginary part

I − Lh, n P −1(I − Lh, n)

Figure : Eigenvalues of I − L and P −1(I − L) for N = 10 (left) and N = 100 (right) for S2

0 transmission condition. The potential is V = 5tx, T = 0.01 and the mesh is

∆t = 0.001, ∆x = 10−3.

2016/09 – Nice

Results for V = |u|2

N 10 100 500 1000 Nnopc 12 71 349 694 Npc 11 22 25 26 T ref 3200.8 Tnopc 2582.3 1332.2 1248.0 1129.7 Tpc 2446.7 408.2 117.6 83.8

Table : Number of iterations and computational time in seconds for algorithms without and with preconditionner, S2

0, V = |u|2, ∆t = 0.001, ∆x = 10−5.

2016/09 – Nice

2D case

The same computations can be made in 2D if we consider subdomains like in the following figure Therefore, no cross-points. We present here results for a computation on GP i∂tu + 1

2∆u − V (x, y)u − β|u|2u + ω · Lzu = 0,

u(0, x, y) = u0(x, y), where (t, x, y) ∈ (0, T) × R2 and Lz = −i(x∂y − y∂x). The potential is V (x, y) = 1 2(γ2

xx2 + γ2 yy2), γx, γy ∈ R.

2016/09 – Nice

BEC

Ω = [−16, 16]2, β = 1000, ω = 0.9. The initial datum is built with GPELab with 1024 nodes in each directions (upscaling) and γx = γy = 1. Dynamical test case: γx = 1.05, γy = 0.95, ∆t = 10−4 and N = 32 subdomains and Robin transmission condition

2016/09 – Nice

Conclusion

Derivation of transparent and absorbing boundary conditions for Schr¨

• dinger

equations Use the DtN map as transmission conditions for DDM Extension: application to other dispersive models in the context of water waves modeling

KdV BBM Green Nagdi

Third order derivative: the number of transmission conditions for left and right boundaries is different.

2016/09 – Nice