Laurette TUCKERMAN laurette@pmmh.espci.fr Numerical Methods for - - PowerPoint PPT Presentation

Laurette TUCKERMAN laurette@pmmh.espci.fr

Numerical Methods for Differential Equations in Physics

Time stepping: ∂tU = LU + N(U) Steady state solving: 0 = LU + N(U) Linear Stability Analysis: λu = Lu + NUu Navier-Stokes Equations ∂tU = −(U · ∇)U − ∇P + ν∆U = −(I − ∇∇−2∇·)(U · ∇)U + ν∆U = N(U) + L U NUu ≡ −(U · ∇)u − (u · ∇)U AUu = NUu + Lu Must solve λu = Lu + NUu

Linear Stability Analysis λu = Lu + NUu How to calculate eigenpairs (λ, u)?

1) Direct: Diagonalisation = QR decomposition Storage: M 2 Time: M 3 3D case with Mx = My = Mz = 102 = ⇒ M = 106 M 2 = 1012 M 3 = 1018 2) Iterative: Calculate a few desired eigenpairs. Use only matrix-vector products u → Au To diagonalise an arbitrary matrix, Each product u → Au requires M 2 operations Generating M eigenpairs requires M iterations

• M 3

Can gain: If A is structured or sparse, then u → Au takes ∼ M ops Aim method at desired eigenvalues.

Power method

Aψj = λjψj where |λ1| > |λ2| > · · · u =

• j

ujψj Au =

• j

ujλjψj . . . ANu =

• j

ujλN

j ψj = λN 1 u1

• ψ1 + ψ2

u2 u1 λ2 λ1 N + ψ3 u3 u1 λ3 λ1 N + . . .

• ↑

↑ converges to first evec first error term Rayleigh quotient: AN−1u, ANu AN−1u, AN−1u ≈ λN−1

1

u1ψ1, λN

1 u1ψ1

λN−1

1

u1ψ1, λN−1

1

u1ψ1 = λ1 What if we want more than one eigenvalue? What if we want an eigenvalue which is not the dominant one (largest |λ|)?

MATRIX TRANSFORMATIONS

If A u = λ u then f(A) u = f(λ) u f(A) = eA∆t f(A) = A−1

f(A) =

j fjAj

fj chosen dynam- ically to extract desired eigenvalues: principle of ARPACK (Sorensen et al.)

EXPONENTIAL POWER METHOD

un+1 = (I − ∆tL)−1(I + ∆tNU)un ≈ e∆t(L+NU)un Approximation valid for ∆t ≪ 1 Time-stepping linearized evolution equation Enhancement factor at each iteration is

• e∆tλ1

e∆tλ2

• 1

where λ1 > λ2 > · · ·

INVERSE POWER METHOD

un+1 = (L + NU)−1un Solve matrix equation using iterative method if necessary Enhancement factor at each iteration is

• λ2

λ1

• ≫ 1 for λ1 ≈ 0

Can shift to find eigenvalues closest to s

• λ2 − s

λ1 − s

• ≫ 1

for λ1 ≈ s

AXISYMMETRIC SPHERICAL COUETTE FLOW

AXISYMMETRIC SPHERICAL COUETTE FLOW

AXISYMMETRIC SPHERICAL COUETTE FLOW

Basic flow at Re = 650 Leading eigenvector

Eigenvalues of Spherical Couette Flow Obtained from inverse power power method with shift

BOSE-EINSTEIN CONDENSATION

Ultra-cold coherent state of matter Predicted by Bose (1924) and Einstein (1925) Realized experimentally by Cornell, Ketterle, Wieman (1995) Nobel prize (2001) Gross-Pitaevskii / Nonlinear Schr¨

• dinger Equation

∂tΨ = i [ 1 2∆

• L

+ µ − V (r) − a|Ψ|2

• N

]Ψ V (x) = 1 2|ω · x|2 = 1 2(ωrr2 + ωzz2)

(cylindrical trap)

Spatial discretisation up to M = 102 × 102 × 102 = 106 Eigenvalues, energies determine decay rates of condensate

800 1200 1600 2000 800 1000 1200 1400 1600 1800 −4.0 0.0 4.0 8.0

• 2

• 1
• 2
• 1
• 2

900 1000 1100 1200 1300 1400 1600 1800 2000 2200 −1.0 0.0 1.0 2.0

• 1

• 2

• 2

• 2
• 1

ωz = ωr/5 (cigar) ωr = ωz/5 (pancake) Hamiltonian saddle-node bifurcation

• f hyperbolic and elliptic fixed points

• 0.2
• 0.1

0.1 0.2

• 0.075
• 0.05
• 0.025

0.025 0.05 0.075

• 0.2
• 0.1

0.1 0.2

• 0.075
• 0.05
• 0.025

0.025 0.05 0.075

• 0.2
• 0.1

0.1 0.2

• 0.075
• 0.05
• 0.025

0.025 0.05 0.075

• 0.4
• 0.2

0.2 0.4

• 0.002
• 0.001

0.001 0.002

q p q B p q p A C D

q ø

A B C

Q Q Q Q Q +

• Q-

+ +

Hamiltonian Systems

f(A) = A f(A) = eA∆t

f(A) = A−1 f(A) = A−2

STEADY STATE SOLVING VIA NEWTON’S METHOD:

0 = LΨ + N(Ψ)

LINEAR STABILITY OF STEADY STATE Ψ:

∂tψ = i [(1 2∆ + µ − V (r))ψ − aΨ2(2ψ + ψ∗)] A ψR ψI

• ≡
• −(L + DN I)

L + DN R ψR ψI

• DN R ≡ µ − V (x) − 3aΨ2

DN I ≡ µ − V (x) − aΨ2

A2 ψR ψI

• =

−(L + DN I)(L + DN R) −(L + DN R)(L + DN I) ψR ψI

• Square, Shift and Invert:

(A2 − s2I)ψn+1 = ψn Sometimes it is easier to solve a preconditioned version: L−2(A2 − s2I)ψn+1 = L−2ψn

Floquet theory

Linear equations with constant coefficients: a¨ x + b ˙ x + cx = 0 = ⇒ x(t) = α1eλ1t + α2eλ2t where aλ2 + bλ + c = 0 ˙ x = cx = ⇒ x(t) = ectx(0)

N

• n=0

cnx(n) = 0 = ⇒ x(t) =

N

• n=1

αneλnt where

N

• n=0

cnλn = 0 Generalize to linear equations with periodic coefficients: a(t)¨ x + b(t) ˙ x + c(t)x = 0 = ⇒ x(t) = α1(t)eλ1t + α2(t)eλ2t a(t), b(t), c(t) have period T = ⇒ α1(t), α2(t) have period T

Floquet theory continued

a(t)¨ x + b(t) ˙ x + c(t)x = 0 = ⇒ x(t) = α1(t)eλ1t + α2(t)eλ2t α1(t), α2(t) Floquet functions λ1, λ2 Floquet exponents eλ1T, eλ2T Floquet multipliers λ1, λ2 not roots of polynomial = ⇒ must calculate numerically or asymptotically ˙ x = c(t)x = ⇒ x(t) = eλtα(t)

N

• n=0

cn(t)x(n) = 0 = ⇒ x(t) =

N

• n=1

eλntαn(t)

Floquet theory and linear stability analysis

Dynamical system: ˙ x = f(x) Limit cycle solution: x(t + T ) = x(t) with ˙ x(t) = f(x(t)) Stability of x(t) : x(t) = x(t) + ǫ(t) ˙ x + ˙ ǫ = f(x(t)) + f ′(x(t))ǫ(t) + . . . ˙ ǫ = f ′(x(t))ǫ(t) Floquet form! ǫ(t) = eλtα(t) with α(t) of period T Re(λ) > 0 = ⇒ x(t) unstable Re(λ) < 0 = ⇒ x(t) stable λ complex = ⇒ most unstable or least stable pert has period different from x(t)

Region of stability

for exponent λ for multiplier eλT Imaginary part non-unique = ⇒ choose Im(λ) ∈ (−πi/T, πi/T ]

In RN, x(t) stable iff real parts of all λj are negative Monodromy matrix: ˙ M = Df(x(t))M with M(t = 0) = I Floquet multipliers/functions = eigenvalues/vectors of M(T )

Cylinder wake

Ideal flow with downstream recirculation zone Von K´ arm´ an vortex street (Re ≥ 46) Laboratory experiment (Taneda, 1982) Off Chilean coast past Juan Fernandez islands

von K´ arman vortex street: Re = U∞d/ν ≥ 46 spatially: temporally: two-dimensional (x, y) periodic, St = fd/U∞ (homogeneous in z) appears spontaneously U2D(x, y, t mod T )

Stability analysis of von K´ arm´ an vortex street

2D limit cycle U2D(x, y, tmod T ) obeys: ∂tU2D = −(U2D · ∇)U2D − ∇P2D + 1 Re∆U2D Perturbation obeys: ∂tu3D = −(U2D(t)·∇)u3D−(u3D·∇)U2D(t)−∇p3D+ 1 Re∆u3D Equation homogeneous in z, periodic in t = ⇒ u3D(x, y, z, t) ∼ eiβzeλβtfβ(x, y, t mod T ) Fix β, calculate largest µ = eλβT via linearized Navier-Stokes

From Barkley & Henderson, J. Fluid Mech. (1996)

mode A: Rec = 188.5 βc = 1.585 = ⇒ λc ≈ 4 mode B: Rec = 259 βc = 7.64 = ⇒ λc ≈ 1 Temporally: µ = 1 = ⇒ steady bifurcation to limit cycle with same periodicity as U2D Spatially: circle pitchfork (any phase in z)

mode A at Re = 210 mode B at Re = 250

From M.C. Thompson, Monash University, Australia (http://mec-mail.eng.monash.edu.au/∼mct/mct/docs/cylinder.html)

Faraday instability

Faraday (1831): Vertical vibration of fluid layer = ⇒ stripes, squares, hexagons In 1990s: first fluid-dynamical quasicrystals: Edwards & Fauve Kudrolli, Pier & Gollub

• J. Fluid Mech. (1994)

Physica D (1998)

Oscillating frame of reference = ⇒ “oscillating gravity” G(t) = g (1 − a cos(ωt)) G(t) = g (1 − a [cos(mωt) + δ cos(nωt + φ0)]) Flat surface becomes linearly unstable for sufficiently high a Consider domain to be horizontally infinite (homogeneous) = ⇒ solutions exponential/trigonometric in x = (x, y) Seek bounded solutions = ⇒ trigonometric: exp(ik · x) Height ζ(x, y, t) =

• k

eik·xˆ ζk(t) Oscillating gravity = ⇒ temporal Floquet problem, T = 2π/ω ˆ ζk(t) =

• j

eλj

ktf j

k(t)

Height ζ(x, y, t) =

• k

eik·xˆ ζk(t) Ideal fluids (no viscosity), sinusoidal forcing = ⇒ Mathieu eq. ∂2

t ˆ

ζk + ω2

0 [1 − a cos(ωt)] ˆ

ζk = 0 ω2

0 combines g, densities, surface tension, wavenumber k

ˆ ζk(t) =

2

• j=1

eλj

ktf j

k(t)

Re(λj

k) > 0 for some j, k =

⇒ ˆ ζk ր = ⇒ flat surface unstable = ⇒ Faraday waves with wavelength 2π/k Im(λj

k) eλT

waves period 1 harmonic same as forcing ω/2 −1 subharmonic twice forcing period

Floquet functions

X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X

t / T (ζ-< ζ >)/

1 2 3

• 0.003
• 0.002
• 0.001

0.001 0.002 0.003

λ

X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X

t / T (ζ-< ζ >)/

1 2 3

• 0.002

0.002 0.004 0.006

λ

X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X

t / T (ζ-< ζ >)/

1 2 3

• 0.03
• 0.02
• 0.01

0.01 0.02 0.03

λ

within tongue 1 /2 within tongue 2 /2 within tongue 3 /2 subharmonic harmonic subharmonic µ = −1 µ = +1 µ = −1

Hexagonal patterns in Faraday instability