Stable Outgoing Wave Filters for Anisotropic Waves Chris Stucchio - - PowerPoint PPT Presentation

▶
stable outgoing wave filters for anisotropic waves
SMART_READER_LITE
LIVE PREVIEW

Stable Outgoing Wave Filters for Anisotropic Waves Chris Stucchio - - PowerPoint PPT Presentation

Mar 25, 2023 399 likes •1.26k views

Stable Outgoing Wave Filters for Anisotropic Waves Chris Stucchio Courant Institute of Mathematical Sciences Joint work with Avy Soffer. Supported by NSF and Bevier Fellowship. Monday, July 7, 2008 1 Linear Waves Linear wave equation: u

slide-1
SLIDE 1

Chris Stucchio Courant Institute of Mathematical Sciences Joint work with Avy Soffer. Supported by NSF and Bevier Fellowship.

Stable Outgoing Wave Filters for Anisotropic Waves

1 Monday, July 7, 2008

slide-2
SLIDE 2

Linear Waves

  • Linear wave equation:
  • ut(x, t)

= H u(, t) H(i∇) = −H†(i∇)

2 Monday, July 7, 2008

slide-3
SLIDE 3

Linear Waves

  • Linear wave equation:
  • Schrodinger equation:
  • ut(x, t)

= H u(, t) H(i∇) = −H†(i∇) H = i∆ u(x, t) = ψ(x, t)

3 Monday, July 7, 2008

slide-4
SLIDE 4

Linear Waves

  • Linear wave equation:
  • Maxwell’s equation
  • ut(x, t)

= H u(, t) H(i∇) = −H†(i∇) H = −µ−1/2∇ × ǫ−1/2 ǫ−1/2∇ × µ−1/2

  • u(x, t)

= (√µ H, √ǫ E)

4 Monday, July 7, 2008

slide-5
SLIDE 5

Linear Waves

  • Linear wave equation:
  • Linearized Euler equation:
  • ut(x, t)

= H u(, t) H(i∇) = −H†(i∇) H =   M∂x1 −∂x1 −∂x2 −∂x1 M∂x1 −∂x2 M∂x1   (x, y) = (p(x, t), vx(x, t), vy(x, t))

5 Monday, July 7, 2008

slide-6
SLIDE 6

Linear Waves

  • Linear wave equation:
  • Relativistic Schrodinger Equation
  • ut(x, t)

= H u(, t) H(i∇) = −H†(i∇) H =

  • −∆ + m2 − m

u(x, t) = ψ(x, t)

6 Monday, July 7, 2008

slide-7
SLIDE 7

Linear Waves

  • Linear wave equation:
  • Linear part of Benjamin-Ono equation:
  • ut(x, t)

= H u(, t) H(i∇) = −H†(i∇) H = |∂x|∂x u(x, t) = h(x, t)

7 Monday, July 7, 2008

slide-8
SLIDE 8

Numerical Solution

  • Finite Differences
  • Finite Elements
  • Spectral methods

I’ll stay agnostic FFT spectral methods rock.

8 Monday, July 7, 2008

Mention Fundamental Complexity

slide-9
SLIDE 9

Numerical Solution

  • Sample spacing:
  • Fundamental complexity of timestepping on
  • Solution on requires careful choice of boundary conditions.

Memory = O((Lkmax)N) Complexity = O((Tmax/δt)(Lkmax)N) [−L, L]N

RN

δx ≤ O(2π/kmax)

9 Monday, July 7, 2008

slide-10
SLIDE 10

Outgoing Waves are a Problem

10 Monday, July 7, 2008

slide-11
SLIDE 11

The Problem

1D Schrodinger Equation

11 Monday, July 7, 2008

slide-12
SLIDE 12

The Problem

Anisotropic Maxwell Incorrect Boundaries

12 Monday, July 7, 2008

Quality is a MAC issue, not an accuracy one. Looks better on Linux.

slide-13
SLIDE 13

Possible Solution: Exact NRBC

  • Dirichlet-to-Neumann boundaries: impose exact non-reflecting boundary

conditions, constructed from Green’s function to free wave.

  • Nonlocal in time, nonlocal on boundary
  • Internal solver restricted (no Fourier spectral methods)
  • Geometry restricted
  • Majda-Engquist, Bayliss-Turkell, Hagstrom, Greengard, Grote, ...

13 Monday, July 7, 2008

slide-14
SLIDE 14

Possible Solution: Perfectly Matched Layers

  • Extend with absorbing layer
  • Dissipation inside layer
  • Must be Perfectly Matched to

avoid reflection at the interface.

  • Equivalent to complex scaling

Region

  • f

Interest PML

14 Monday, July 7, 2008

slide-15
SLIDE 15

Possible Solution: Perfectly Matched Layers

  • Complex scaling for Wave equation:
  • PML (Conjugate Operator) for general linear waves:

H → ezAHe−zA A = x · i∇ + i∇ · x H → ezAHe−zA A = x · vg(i∇) + vg(i∇) · x

15 Monday, July 7, 2008

slide-16
SLIDE 16

Complex scaling is easy

  • Change coordinates
  • Make layer perfectly matched
  • Stable if

A = x · i∇ + i∇ · x ezA = Dilation(z) k1vg,1(k) ≥ 0

Picture from Becache, Fauqueux, Joly, JCP 188 (2003) 399–433.

16 Monday, July 7, 2008

slide-17
SLIDE 17

PML Instability

  • PML unstable for some

anisotropic waves (Becache, Fauqueux, Joly, 2003).

Pictures from Becache, Fauqueux, Joly, JCP 188 (2003) 399–433.

17 Monday, July 7, 2008

Picture and graph are for difgerent equations

slide-18
SLIDE 18

Conjugate operators are hard

A = x · vg(i∇) + vg(i∇) · x ezA = ?

18 Monday, July 7, 2008

If stability condition not satisfied, complex scaling shifts spectrum the wrong way.

slide-19
SLIDE 19
  • Identify outgoing waves
  • Filter them off
  • Nothing hits the boundary

Phase Space Filters

19 Monday, July 7, 2008

slide-20
SLIDE 20

Outgoing waves

20 Monday, July 7, 2008

slide-21
SLIDE 21

Outgoing Waves, Schrodinger Equation

  • 1D Schrodinger Equation
  • Center of mass at ,

x = vt width = σ + t/σ ψ0(x) = eivx √σ exp −x2 2σ2

  • ψ(x, t)

= eivx

  • σ + it/σ

exp −(x − vt)2 2σ2(1 + it/σ)

  • 21

Monday, July 7, 2008

slide-22
SLIDE 22

Outgoing Waves, Schrodinger Equation

  • Outgoing wave
  • Incoming wave

ψ0(x) = e+ivxe−(x−L)2/σ2 ψ0(x) = e−ivxe−(x−L)2/σ2 Trajectory = L + vt Trajectory = L − vt

22 Monday, July 7, 2008

slide-23
SLIDE 23

Outgoing Waves, Schrodinger Equation

  • Outgoing wave
  • Incoming wave

ψ0(x) = e+ivxe−(x−L)2/σ2 ψ0(x) = e−ivxe−(x−L)2/σ2 Trajectory = L + vt Trajectory = L − vt

23 Monday, July 7, 2008

slide-24
SLIDE 24

Outgoing Waves, Schrodinger Equation

  • Outgoing wave
  • Incoming wave

ψ0(x) = e+ivxe−(x−L)2/σ2 ψ0(x) = e−ivxe−(x−L)2/σ2 Trajectory = L + vt Trajectory = L − vt

24 Monday, July 7, 2008

slide-25
SLIDE 25

Outgoing Waves, Schrodinger Equation

  • Mixed wave:

ψ0(x) = e−ivxe−(x−L)2/σ2+e+ivxe−(x−L)2/σ2

25 Monday, July 7, 2008

slide-26
SLIDE 26
  • Mixed wave:

Outgoing Waves, Schrodinger Equation

ψ0(x) = e−ivxe−(x−L)2/σ2

Incoming wave Outgoing wave

26 Monday, July 7, 2008

slide-27
SLIDE 27
  • Mixed wave:

Outgoing Waves, Schrodinger Equation

ψ0(x) = e−ivxe−(x−L)2/σ2

Incoming wave Outgoing wave

26 Monday, July 7, 2008

slide-28
SLIDE 28
  • Mixed wave:

Outgoing Waves, Schrodinger Equation

ψ0(x) = e−ivxe−(x−L)2/σ2

Incoming wave Problem solved!

27 Monday, July 7, 2008

slide-29
SLIDE 29

It really is that easy

  • Windowed Fourier Transform:
  • Outgoing waves:

The Wavelet Transform, Time Frequency Localization and Signal Analysis, Ingrid Daubechies, IEEE Trans. Info. Theory, Vol 36 5 1990

ψ(x) =

  • a∈Z
  • b∈Z

ψa,beibk0xg(x − ax0) g(x) = e−x2/σ2 ax0 > L bk0 > σ−1

28 Monday, July 7, 2008

slide-30
SLIDE 30

Quantum Phase Space

  • Quantum phase space is set of

points , x a position and k a frequency.

  • Heisenberg Uncertainty

principle: localizing on region of volume causes error .

  • A function is localized near a

point if it is localized in position near and it’s Fourier transform is localized near .

O(2π ln(ǫ))

ǫ

(x0, k0) x0 k0 (x, k) ∈ RN × RN

29 Monday, July 7, 2008

slide-31
SLIDE 31
  • Ambiguous waves
  • Spreads in both directions

Outgoing Waves, Schrodinger Equation

ψ0(x) = ei0xe−(x−L)2/σ2

Issue can be resolved.

30 Monday, July 7, 2008

slide-32
SLIDE 32

Phase space filters

31 Monday, July 7, 2008

slide-33
SLIDE 33

Phase Space Filters

Outgoing Waves:

ax0 > L bk0 > σ−1

32 Monday, July 7, 2008

slide-34
SLIDE 34

Phase Space Filters

Outgoing Waves:

ax0 > L bk0 > σ−1

32 Monday, July 7, 2008

slide-35
SLIDE 35

Phase Space Filters

Outgoing Waves:

ax0 > L bk0 > σ−1

32 Monday, July 7, 2008

slide-36
SLIDE 36

A simpler version

  • We want rightward moving

waves near the boundary.

  • Extend computational domain
  • Localize in boundary layer

[−L − w, L + w]N 1(x > L)1(k > kmin)1(x > L) = O+

33 Monday, July 7, 2008

Step function built out of erfc(x) functions, as smooth as possible without overflowing.

slide-37
SLIDE 37

How does it work?

  • Take wave comprised of incoming and outgoing waves, plus interior waves.

ψ0(x) ≈ eivxg(x − L − 1) + e−ivxg(x − L − 1) + interior waves

34 Monday, July 7, 2008

slide-38
SLIDE 38

How does it work?

  • Take wave comprised of incoming and outgoing waves, plus interior waves.

ψ0(x) ≈ eivxg(x − L − 1) + e−ivxg(x − L − 1) + interior waves

Our target

35 Monday, July 7, 2008

slide-39
SLIDE 39

How does it work?

  • Take wave comprised of incoming and outgoing waves, plus interior waves.
  • We don’t care about interior waves

1(x > L)ψ0(x) ≈ eivxg(x − L − 1) + e−ivxg(x − L − 1) + ψ0(x) ≈ eivxg(x − L − 1) + e−ivxg(x − L − 1) + interior waves

36 Monday, July 7, 2008

slide-40
SLIDE 40

How does it work?

  • Take wave comprised of incoming and outgoing waves, plus interior waves.
  • We don’t care about interior waves
  • Or incoming waves

1(x > L)ψ0(x) ≈ eivxg(x − L − 1) + e−ivxg(x − L − 1) + ψ0(x) ≈ eivxg(x − L − 1) + e−ivxg(x − L − 1) + interior waves 1(k > kmin)1(x > L)ψ0(x) ≈ eivxg(x − L − 1) + 0

37 Monday, July 7, 2008

slide-41
SLIDE 41

How does it work?

  • Or incoming waves
  • Symmetry is always good (for stability, etc):

1(k > kmin)1(x > L)ψ0(x) ≈ eivxg(x − L − 1) + 0 1(x > L)1(k > kmin)1(x > L)ψ0(x) ≈ eivxg(x − L − 1)

38 Monday, July 7, 2008

slide-42
SLIDE 42

How does it work?

  • Or incoming waves
  • Symmetry is always good (for stability, etc):
  • Operator localizes outgoing waves, and lets us remove them:

1(k > kmin)1(x > L)ψ0(x) ≈ eivxg(x − L − 1) + 0 1(x > L)1(k > kmin)1(x > L)ψ0(x) ≈ eivxg(x − L − 1) ψ0(x) − O+ψ0(x) = 0 + e−ivxg(x − L) + Interior Waves O+

39 Monday, July 7, 2008

slide-43
SLIDE 43

let let on domain for n = 1 to :

  • utput

Propagation Algorithm

[−L − w, L + w]N

Tmax/Ts u(x) := u0(x) Ts := O(w/3vmax ln(ǫ)) u(x) ←

all sides

(1 − O+)

  • u(x)

u(x) = u(x, nTs) u(x) ← ei∆Tsu(x)

40 Monday, July 7, 2008

slide-44
SLIDE 44

let let on domain for n = 1 to :

  • utput

Propagation Algorithm

[−L − w, L + w]N

Tmax/Ts u(x) := u0(x) Ts := O(w/3vmax ln(ǫ)) u(x) ←

all sides

(1 − O+)

  • u(x)

u(x) = u(x, nTs) u(x) ← ei∆Tsu(x) Not enough time to travel distance w

41 Monday, July 7, 2008

slide-45
SLIDE 45

let let on domain for n = 1 to :

  • utput

Propagation Algorithm

[−L − w, L + w]N

Tmax/Ts u(x) := u0(x) Ts := O(w/3vmax ln(ǫ)) u(x) ←

all sides

(1 − O+)

  • u(x)

u(x) = u(x, nTs) u(x) ← ei∆Tsu(x) Not enough time to travel distance w Propagate any way you like.

42 Monday, July 7, 2008

slide-46
SLIDE 46

let let on domain for n = 1 to :

  • utput

Propagation Algorithm

[−L − w, L + w]N

Tmax/Ts u(x) := u0(x) Ts := O(w/3vmax ln(ǫ)) u(x) ←

all sides

(1 − O+)

  • u(x)

u(x) = u(x, nTs) u(x) ← ei∆Tsu(x) Not enough time to travel distance w Propagate any way you like. (Nothing reached the boundary yet.)

43 Monday, July 7, 2008

slide-47
SLIDE 47

let let on domain for n = 1 to :

  • utput

Propagation Algorithm

[−L − w, L + w]N

Tmax/Ts u(x) := u0(x) Ts := O(w/3vmax ln(ǫ)) u(x) ←

all sides

(1 − O+)

  • u(x)

u(x) = u(x, nTs) u(x) ← ei∆Tsu(x) Not enough time to travel distance w Propagate any way you like. (Nothing reached the boundary yet.) Filter outgoing waves about to reach the boundary.

44 Monday, July 7, 2008

slide-48
SLIDE 48

let let on domain for n = 1 to :

  • utput

Propagation Algorithm

[−L − w, L + w]N

Tmax/Ts u(x) := u0(x) Ts := O(w/3vmax ln(ǫ)) u(x) ←

all sides

(1 − O+)

  • u(x)

u(x) = u(x, nTs) u(x) ← ei∆Tsu(x) Not enough time to travel distance w Propagate any way you like. (Nothing reached the boundary yet.) Filter outgoing waves about to reach the boundary. Next propagation step is accurate: waves which would have reached boundary were filtered.

45 Monday, July 7, 2008

slide-49
SLIDE 49

Phase Space Filtering, Schrodinger Equation

  • is “blurring” operator in frequency domain
  • Characteristic distance of “blurring” (in k domain)

1(x1 > L)

[

  • 1(x1 > L)f](k) ≈ (...)e−k2/w2 ⋆ ˆ

f(k) kmin = O(ln(ǫ−1)/w) O+ = 1(x1 > L)1(k > kmin)1(x1 > L)

46 Monday, July 7, 2008

slide-50
SLIDE 50

Schrodinger Equation

Results

47 Monday, July 7, 2008

slide-51
SLIDE 51
  • Measured error as

function of frequency of initial data.

  • Errors are large for low

frequencies, small for high.

  • By increasing width of

buffer, one reduce errors for low frequencies.

Schrodinger equation: Error vs Frequency

48 Monday, July 7, 2008

slide-52
SLIDE 52

Phase Space Filters for Vector Systems

49 Monday, July 7, 2008

slide-53
SLIDE 53

Vector Systems

  • Linear wave equation:
  • ut(x, t)

= H u(, t) H(i∇) = −H†(i∇) H =   H11(k) . . . H1N(k) . . . . . . . . . −H1N(k) . . . HNN(k)  

50 Monday, July 7, 2008

slide-54
SLIDE 54

Wavepackets

  • Not a 1-way wavepacket:
  • Will split into N different wavepackets.

u0(x) =   eikxg(x) . . .  

51 Monday, July 7, 2008

slide-55
SLIDE 55

Wavepackets

  • Diagonalize hamiltonian to find dispersion relation
  • For each frequency, is skew adjoint matrix. Can always do this.
  • Plane Waves:

H = D†   iω1(k) . . . . . . iωj(k) . . . . . . iωM(k)   D H h(x, t) =   d1,1(k) . . . d1,N(k)   ei(kx−ω1(k)t)

52 Monday, July 7, 2008

slide-56
SLIDE 56

Wavepackets

  • Localize a plane wave:

u0(x) =   d11(k0) . . . d1N(k0)   eik0xg(x)

53 Monday, July 7, 2008

slide-57
SLIDE 57

Wavepackets

  • Localize a plane wave:
  • Wavepacket propagation:
  • (Fourier transform and do stationary phase)

u0(x) =   d11(k0) . . . d1N(k0)   eik0xg(x) u(x, t) =   d11(k0) . . . d1N(k0)   ei(k0x−ω1(k0)t)[eDtg](x − ∇kω1(k0)t)

54 Monday, July 7, 2008

slide-58
SLIDE 58

Wavepackets

  • Localize a plane wave:
  • Wavepacket propagation:
  • (Fourier transform and do stationary phase)

u0(x) =   d11(k0) . . . d1N(k0)   eik0xg(x) u(x, t) =   d11(k0) . . . d1N(k0)   ei(k0x−ω1(k0)t)[eDtg](x − ∇kω1(k0)t)

Translation

55 Monday, July 7, 2008

slide-59
SLIDE 59

Wavepackets

  • Localize a plane wave:
  • Wavepacket propagation:
  • (Fourier transform and do stationary phase)

u0(x) =   d11(k0) . . . d1N(k0)   eik0xg(x) u(x, t) =   d11(k0) . . . d1N(k0)   ei(k0x−ω1(k0)t)[eDtg](x − ∇kω1(k0)t)

Dispersion Translation

56 Monday, July 7, 2008

slide-60
SLIDE 60

Wavepackets

  • Envelope obeys Schrodinger like equation:
  • is the Hessian of the dispersion relation.

Hω1(k0)

  • [eDtg](k)

= exp((ωq(k) − ω1(k0) − [∇k(ωq)](k0)(k − k0))t)ˆ g(k) ≈ e(k−k0)[Hω1(k0)](k−k0)tˆ g(k)

Hessian is Quadratic Differential

  • perator, like Laplacian.

57 Monday, July 7, 2008

slide-61
SLIDE 61

Wavepackets

  • Schrodinger:
  • Vector system:

u0(x) =   d11(k0) . . . d1N(k0)   eik0xg(x) ψ0(x) = eikxe−x2/σ2 position = kt position = [∇kω1(k0)]t = vg(k0)t

58 Monday, July 7, 2008

slide-62
SLIDE 62

Phase space filtering for vector systems

  • We want rightward moving

waves near the boundary.

  • Extend computational domain
  • Localize in boundary layer

[−L − w, L + w]N

59 Monday, July 7, 2008

Step function built out of erfc(x) functions, as smooth as possible without overflowing.

slide-63
SLIDE 63

Phase space filtering for vector systems

  • Project onto rightward-moving group velocities
  • Un-diagonalize to project onto rightward moving waves
  • Localize:

P(k) =   1(∇kω1(k) · e1 > 0) . . . ... . . . . . . 1(∇kωM(k) · e1 > 0   D†P(k)D O+ = 1(x1 > L)D†P(k)D1(x1 > L)

60 Monday, July 7, 2008

slide-64
SLIDE 64

let let on domain for n = 1 to :

  • utput

Propagation Algorithm

[−L − w, L + w]N

Tmax/Ts u(x) := u0(x) u(x) ←

all sides

(1 − O+)

  • u(x)

u(x) = u(x, nTs) Ts = O(w/3vmax ln(ǫ)), u(x) ← eiHTsu(x)

61 Monday, July 7, 2008

slide-65
SLIDE 65

Numerical Results, Anisotropic Waves

62 Monday, July 7, 2008

slide-66
SLIDE 66

Maxwell’s Equations in Birefringent Medium

  • In a birefringent medium, Maxwell’s equations take the form
  • The wavefield is defined as
  • Assume is a scalar, and assume
  • Then with

H = −µ−1/2∇ × ǫ−1/2 ǫ−1/2∇ × µ−1/2

  • u(x, t) = (√µ

H, √ǫ E)T µ ǫ =   1 b b 1 c   f = (1/2)( √ 1 + b + √ 1 − b), g = (1/2)(− √ 1 + b + √ 1 − b), ωj=1,2(k) = (−1)1+jic−1|k| ωj=3,4(k) = (−1)1+ji

  • (f 2 + g2)(k2

1 + k2 2) − 4fgk1k2

ωj=5,6(k) = 0.

63 Monday, July 7, 2008

slide-67
SLIDE 67

Maxwell’s Equations, mode

Birefringent medium, b=0.25

TMz

64 Monday, July 7, 2008

slide-68
SLIDE 68

5 10 15 20

Frequency

10

  • 5

10

  • 4

10

  • 3

10

  • 2

10

  • 1

10

Relative and Absolute errors

Errors for Maxwell’s Equations

t d

1

L

1

L

| u | / | u ! u | p u s

t d

"

L

"

L

| u | / | u ! u | p u s

t d

2

L

2

L

| u | / | u ! u | p u s

Maxwell’s Equations: Error vs Frequency

65 Monday, July 7, 2008

slide-69
SLIDE 69

Linearized Euler Equations

  • Euler equations, linearized about jet flow:
  • The dispersion relations are:

H =   M∂x1 −∂x1 −∂x2 −∂x1 M∂x1 −∂x2 M∂x1   ω1(k) = Mk1 + |k|, ω2(k) = Mk1 − |k|, ω3(k) = Mk1

66 Monday, July 7, 2008

slide-70
SLIDE 70

Euler Equations

Results, M=0.5

67 Monday, July 7, 2008

slide-71
SLIDE 71

Linearized Quasi-Geostrophic Equations

  • Quasi-geostrophic equations, midlatitude:
  • is a streamfunction:
  • Geostrophic balance: Coriolis force = horizontal pressure gradient
  • Anisotropic and non-local

V = Mean wind, F ∼ (earth’s rotation)2 g ˜ β = FV + β, β = R cos(φ) ψ

  • v = ∇⊥ψ

H = V ∂x − ˜ β(−∆ + F)−1∂x

68 Monday, July 7, 2008

\phi = latitude. Mention that Tom Hagstrom prompted this example.

slide-72
SLIDE 72

Dispersion Relations

  • Complicated dispersion

relation, not quite hyperbolic

  • PML unstable in y direction for

k0 < 0

69 Monday, July 7, 2008

slide-73
SLIDE 73

Dispersion Relations

  • Complicated dispersion

relation, not quite hyperbolic

  • PML unstable in y direction for
  • PML unstable in x direction on

irregular region

k0 < 0

Unstable Stable

70 Monday, July 7, 2008

Tom Hagstrom can’t fix this by completing the square.

slide-74
SLIDE 74

Quasi-Geostrophic Equations

71 Monday, July 7, 2008

slide-75
SLIDE 75

5 10 15 20

Frequency

10

  • 5

10

  • 4

10

  • 3

10

  • 2

10

  • 1

10

Error

Relative Errors

Maxwell Euler Quasi-geostrophic

  • Measured error as

function of frequency of initial data.

  • Errors are large for low

frequencies, small for high.

  • By increasing width of

buffer, one reduce errors for low frequencies.

Errors

72 Monday, July 7, 2008

Higher error for quasi-geostrophic due to nonlocality. Can be fixed.

slide-76
SLIDE 76

Stability

73 Monday, July 7, 2008

slide-77
SLIDE 77

Stability of Phase Space Filtering

  • Operator is self-adjoint, and .
  • Implies filtering is dissipative:
  • Propagation operator has norm 1:
  • Numerical solution is strongly stable:

O+ σ(O+) ⊆ [0, 1] σ(1 − O+) ⊆ 1 − [0, 1] = [0, 1] ||eHTs

all sides

(1 − O+)

  • || ≤ ||eHTs||
  • all sides

||(1 − O+)|| ≤ 1

  • all sides

1 = 1

||u(x, t)|| ≤ ||u0(x)||

74 Monday, July 7, 2008

slide-78
SLIDE 78

PML

Phase Space Filter

75 Monday, July 7, 2008

slide-79
SLIDE 79

Low Frequencies

76 Monday, July 7, 2008

slide-80
SLIDE 80

The Low Frequency Problem

  • Heisenberg Uncertainty principle limits phase space filters for low

frequencies.

  • Filter width :
  • PML has similar issues: low frequencies dissipate over long distances.
  • Dirichlet-to-Neumann immune to this problem in homogeneous case. In

inhomogeneous case, Dirichlet-to-Neumann built using approximations valid

  • nly for high frequencies (Pseudo/Paradifferential calculus, see Szeftel).

w = O(ln(ǫ)/kmin)

Memory = O((kmax/kmin)N)

77 Monday, July 7, 2008

slide-81
SLIDE 81

50 100 150 200 Position 5 10 15 20 25 30 35 Frequency Filter 0 Filter 1 Filter 2 ) 2 ( 8 9 . = x ! )

1

2 ( 8 9 . = x ! )

2

2 ( 8 9 . = x ! Interaction Region

Filter Regions in Phase Space

Multiscale Solution

  • Narrow filter for high frequency.
  • Use filter with double the width

to filter low frequencies; cut sampling rate in half.

  • Filter width w = O(ln(ǫ)/kmin)

Memory = O(ln(ǫ) log2(kmax/kmin))

78 Monday, July 7, 2008

slide-82
SLIDE 82

The Low Frequency Problem: Resolution

  • Implemented for 1-dimensional

Schrodinger equation

  • Cost:
  • If is unknown, cost is:
  • Works for long range potential/

inhomogeneity.

kmin

O(log2(kmax/kmin)) O

  • Tmax log2
  • Tmax

vk≈0 L

  • 79

Monday, July 7, 2008

slide-83
SLIDE 83

The Low Frequency Problem: Resolution

  • Implemented for 1-dimensional

Schrodinger equation

  • Cost:
  • If is unknown, cost is:
  • Works for long range potential/

inhomogeneity.

kmin

O(log2(kmax/kmin)) O(log2(Tmax))

80 Monday, July 7, 2008

slide-84
SLIDE 84

The Low Frequency Problem: Resolution

  • Implemented for 1-dimensional

Schrodinger equation

  • Cost:
  • If is unknown, cost is:
  • Works for long range potential/

inhomogeneity.

kmin

O(log2(kmax/kmin)) O(log2(Tmax))

80 Monday, July 7, 2008

slide-85
SLIDE 85

Conclusion

  • Phase space filtering a new method of filtering outgoing waves.
  • Works for anisotropic, inhomogeneous and even non-local waves.
  • Stable and accurate: confirmed by rigorous theorem and numerical tests.
  • [1] Open Boundaries for the Nonlinear Schrodinger Equation, with A. Soffer. JCP Vol. 225, Issue 2, p.p.

1218-1232. arXiv:math/0609183

  • [2] Multiscale Resolution of Shortwave-Longwave Interaction, with A. Soffer. CPAM (accepted). arXiv:0705.3501
  • [3] Stable Open Boundaries for Anisotropic Waves, with A. Soffer (submitted). arXiv:0805.2929
  • All papers available from my webpage: http://cims.nyu.edu/~stucchio/

81 Monday, July 7, 2008