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: ∂ t U = LU + N ( U ) 0 = LU + N ( U ) Steady state solving: Linear Stability Analysis: λu = Lu + N U u Navier-Stokes Equations ∂ t U = − ( U · ∇ ) U − ∇ P + ν ∆ U = − ( I − ∇∇ − 2 ∇· )( U · ∇ ) U + ν ∆ U = N ( U ) + L U N U u ≡ − ( U · ∇ ) u − ( u · ∇ ) U A U u = N U u + Lu Must solve λu = Lu + N U u

Linear Stability Analysis λu = Lu + N U u How to calculate eigenpairs ( λ, u ) ? 1) Direct: Diagonalisation = QR decomposition Storage: M 2 Time: M 3 3D case with M x = M y = M z = 10 2 = ⇒ M = 10 6 M 2 = 10 12 M 3 = 10 18 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 � M 3 Generating M eigenpairs requires M iterations 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 = u j ψ j j � Au = u j λ j ψ j j . . . � � � λ 2 � N � λ 3 � N u 2 u 3 � A N u = u j λ N j ψ j = λ N 1 u 1 ψ 1 + ψ 2 + ψ 3 + . . . u 1 λ 1 u 1 λ 1 j ↑ ↑ converges to first evec first error term Rayleigh quotient: � λ N − 1 � A N − 1 u, A N u � u 1 ψ 1 , λ N 1 u 1 ψ 1 � 1 � A N − 1 u, A N − 1 u � ≈ = λ 1 � λ N − 1 u 1 ψ 1 , λ N − 1 u 1 ψ 1 � 1 1 What if we want more than one eigenvalue? What if we want an eigenvalue which is not the dominant one (largest | λ | )?

M ATRIX T RANSFORMATIONS If A u = λ u then f ( A ) u = f ( λ ) u f ( A ) = � j f j A j f j chosen dynam- ically to extract desired eigenvalues: principle of A RPACK f ( A ) = e A ∆ t f ( A ) = A − 1 (Sorensen et al.)

E XPONENTIAL P OWER M ETHOD u n +1 = ( I − ∆ tL ) − 1 ( I + ∆ tN U ) u n ≈ e ∆ t ( L + N U ) u n Approximation valid for ∆ t ≪ 1 Time-stepping linearized evolution equation Enhancement factor at each iteration is � � e ∆ tλ 1 � � � � � � 1 where λ 1 > λ 2 > · · · � � e ∆ tλ 2 �

I NVERSE P OWER M ETHOD u n +1 = ( L + N U ) − 1 u n Solve matrix equation using iterative method if necessary Enhancement factor at each iteration is � � λ 2 � � � ≫ 1 for λ 1 ≈ 0 � � λ 1 � Can shift to find eigenvalues closest to s � � λ 2 − s � � � ≫ 1 for λ 1 ≈ s � � � λ 1 − s

A XISYMMETRIC S PHERICAL C OUETTE F LOW

A XISYMMETRIC S PHERICAL C OUETTE F LOW

A XISYMMETRIC S PHERICAL C OUETTE F LOW Leading Basic flow at eigenvector Re = 650

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

B OSE -E INSTEIN C ONDENSATION 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¨ odinger Equation ∂ t Ψ = i [ 1 + µ − V ( r ) − a | Ψ | 2 2∆ ]Ψ � �� � ���� N L V (x) = 1 2 | ω · x | 2 = 1 2( ω r r 2 + ω z z 2 ) (cylindrical trap) Spatial discretisation up to M = 10 2 × 10 2 × 10 2 = 10 6 Eigenvalues, energies determine decay rates of condensate

2000 1300 1200 1600 E + E E 1 + E 1100 1 1200 1000 E E � � 900 800 8.0 2.0 2 2 � � + + 4.0 2 1.0 1 � 2 1 � 0.0 0.0 2 2 � � E G E G � � N N −1.0 N N −4.0 1400 1600 1800 2000 2200 800 1000 1200 1400 1600 1800 N N ω z = ω r / 5 (cigar) ω r = ω z / 5 (pancake) Hamiltonian saddle-node bifurcation of hyperbolic and elliptic fixed points 1 1 1 1 1 1 1 1 1 1 1 1 1 1

A B 0.075 0.075 0.05 0.05 0.025 0.025 p p 0 0 Q- Q Q - Q + + -0.025 -0.025 -0.05 -0.05 -0.075 -0.075 -0.2 -0.1 0 0.1 0.2 -0.2 -0.1 0 0.1 0.2 q q C D 0.075 ø 0.05 0.002 0.025 0.001 Q p + 0 q -0.4 -0.2 0.2 0.4 -0.025 -0.001 Q - -0.05 -0.002 -0.075 A -0.2 -0.1 0 0.1 0.2 C q B

Hamiltonian Systems f ( A ) = A f ( A ) = e A ∆ t

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

S TEADY STATE SOLVING VIA N EWTON ’ S M ETHOD : 0 = L Ψ + N (Ψ) L INEAR STABILITY OF STEADY STATE Ψ : ∂ t ψ = i [(1 2∆ + µ − V ( r )) ψ − a Ψ 2 (2 ψ + ψ ∗ )] � ψ R � � ψ R � � � − ( L + DN I ) 0 A ≡ ψ I L + DN R ψ I 0 DN R ≡ µ − V (x) − 3 a Ψ 2 DN I ≡ µ − V (x) − a Ψ 2

� ψ R � A 2 = ψ I � − ( L + DN I )( L + DN R ) � � ψ R � 0 − ( L + DN R )( L + DN I ) ψ I 0 Square, Shift and Invert: ( A 2 − s 2 I ) ψ n +1 = ψ n Sometimes it is easier to solve a preconditioned version: L − 2 ( A 2 − s 2 I ) ψ n +1 = L − 2 ψ n

Floquet theory Linear equations with constant coefficients: ⇒ x ( t ) = α 1 e λ 1 t + α 2 e λ 2 t a ¨ x + b ˙ x + cx = 0 = aλ 2 + bλ + c = 0 where ⇒ x ( t ) = e ct x (0) x = cx = ˙ N N � � c n x ( n ) = 0 = α n e λ n t ⇒ x ( t ) = n =0 n =1 N � c n λ n = 0 where n =0 Generalize to linear equations with periodic coefficients: ⇒ x ( t ) = α 1 ( t ) e λ 1 t + α 2 ( t ) e λ 2 t a ( t )¨ x + b ( t ) ˙ x + c ( t ) x = 0 = a ( t ) , b ( t ) , c ( t ) have period T = ⇒ α 1 ( t ) , α 2 ( t ) have period T

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

Floquet theory and linear stability analysis Dynamical system: x = f ( x ) ˙ x ( t + T ) = x ( t ) with ˙ x ( t ) = f ( x ( t )) Limit cycle solution: Stability of x ( t ) : x ( t ) = x ( t ) + ǫ ( t ) ˙ f ( x ( t )) + f ′ ( x ( t )) ǫ ( t ) + . . . x + ˙ ǫ = f ′ ( x ( t )) ǫ ( t ) ǫ ˙ = e λt α ( t ) Floquet form! ǫ ( 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 multiplier e λT for exponent λ Imaginary part non-unique = ⇒ choose Im ( λ ) ∈ ( − πi/T, πi/T ] In R N , 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 with downstream recirculation zone Ideal flow Von K´ arm´ an vortex street ( Re ≥ 46 ) Off Chilean coast Laboratory experiment past Juan Fernandez islands (Taneda, 1982)

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

Stability analysis of von K´ arm´ an vortex street 2D limit cycle U 2 D ( x, y, t mod T ) obeys: ∂ t U 2 D = − (U 2 D · ∇ )U 2 D − ∇ P 2 D + 1 Re ∆U 2 D Perturbation obeys: ∂ t u 3 D = − (U 2 D ( t ) ·∇ )u 3 D − (u 3 D ·∇ )U 2 D ( t ) −∇ p 3 D + 1 Re ∆u 3 D Equation homogeneous in z , periodic in t = ⇒ u 3 D ( x, y, z, t ) ∼ e iβz e λ β t f β ( x, y, t mod T ) Fix β , calculate largest µ = e λ β T via linearized Navier-Stokes

From Barkley & Henderson, J. Fluid Mech. (1996) mode A: Re c = 188 . 5 mode B: Re c = 259 β c = 1 . 585 = ⇒ λ c ≈ 4 β c = 7 . 64 = ⇒ λ c ≈ 1 Temporally: µ = 1 = ⇒ steady bifurcation to limit cycle with same periodicity as U 2 D 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( i k · x) � e i k · x ˆ Height ζ ( x, y, t ) = ζ k ( t ) k Oscillating gravity = ⇒ temporal Floquet problem, T = 2 π/ω � e λ j ˆ k t f j ζ k ( t ) = k ( t ) j

� e i k · x ˆ Height ζ ( x, y, t ) = ζ k ( t ) k Ideal fluids (no viscosity), sinusoidal forcing = ⇒ Mathieu eq. t ˆ 0 [1 − a cos( ωt )] ˆ ∂ 2 ζ k + ω 2 ζ k = 0 ω 2 0 combines g , densities, surface tension, wavenumber k 2 � e λ j ˆ k t f j ζ k ( t ) = k ( t ) j =1 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 0 1 harmonic same as forcing ω/ 2 − 1 subharmonic twice forcing period