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

∂ t U = LU + N ( U ) Time stepping: 0 = LU + N ( U ) Steady state solving: 0 = F ( U ) Newton’s method 0 = F ( U − u ) ≈ F ( U ) − DF ( U ) u � DF ( U ) u = F ( U ) U ← U − u

Newton’s Method converges quadratically U n +1 = U n − F ( U n ) F ′ ( U n ) 2 F ′′ ( U n )( U − U n ) 2 + . . . F ( U ) = 0 = F ( U n ) + F ′ ( U n )( U − U n ) + 1 F ′ ( U n ) + F ′ ( U n ) F ′′ ( U n ) 0 = F ( U n ) F ′ ( U n ) ( U − U n ) 2 + . . . F ′ ( U n )( U − U n ) + 1 2 F ′′ ( U n ) 0 = F ( U n ) F ′ ( U n ) ( U − U n ) 2 + . . . F ′ ( U n ) + ( U − U n ) + 1 2 F ′′ ( U n ) F ′ ( U n ) ( U − U n ) 2 + . . . 0 = − U n +1 + U + 1 2 F ′′ ( U n ) F ′ ( U n ) ( U − U n ) 2 + . . . U n +1 − U = 1 2 F ′′ ( U n ) F ′ ( U n ) ǫ 2 + . . . ǫ n +1 = 1 2 Typical sequence: ǫ = 10 − 1 , 10 − 2 , 10 − 4 , 10 − 8 , 10 − 16

Much faster than timestepping: U ( t ) = U + ce − λt U ( t n ) − U = ce − λt n U ( t n +1 ) − U = ce − λ ( t n +∆ t ) = e − λ ∆ t ( U ( t n ) − U )) Linear convergence: ǫ n +1 ∼ cǫ n with | c | � 1 since ∆ t ≪ 1 In addition to converging faster than timestepping, Newton’s method can converge to unstable states.

Fixed points and linear stability. ˙ x = f ( x ) unstable stable 0 = f (¯ x ) Fixed point ¯ x d dt (¯ x + ǫ ( t )) = f (¯ x + ǫ ) Linear stability of ¯ x x ) ǫ + 1 x ) ǫ 2 + · · · ≈ f ′ (¯ x ) + f ′ (¯ 2 f ′′ (¯ ǫ = f (¯ ˙ x ) ǫ � increases if f ′ (¯ x ) > 0 ǫ ( t ) = e tf ′ (¯ x ) ǫ (0) decreases if f ′ (¯ x ) < 0

Saddle-node Bifurcations x = f ( x ) = µ − x 2 ˙ Fixed points: x ± = ±√ µ ¯ for µ > 0 Stability: x ± = − 2( ±√ µ ) = ∓ 2 √ µ f ′ (¯ x ± ) = − 2¯ x + ) = f ′ ( √ µ ) = − 2 √ µ < 0 = f ′ (¯ ⇒ ¯ x + stable x − ) = f ′ ( −√ µ ) = 2 √ µ > 0 = f ′ (¯ ⇒ ¯ x − unstable

f ( x, µ ) = c 00 + c 10 x + c 01 µ + c 20 x 2 + . . . general quadratic polynomial x 2 = ± ˜ µ ± ˜ four cases, depending on signs of c ’s Newton’s method finds steady states independently of their stability Where might saddle-node bifurcations occur?

Swift-Hohenberg equation � 2 u − u 3 � q 2 ∂ t u = µu − c + ∆ Derived by J. Swift and P.C. Hohenberg (Phys. Rev. A 15, 319 (1977)) to describe pattern formation in convection For u ≪ 1 , � 2 u � q 2 ∂ t u = µu − c + ∆ u ∼ exp( ikx + σt ) c − k 2 � 2 u = � q 2 σu = µu − ⇒ σ = µ, k = q c Add quadratic term to obtain hexagons � 2 u + g 1 u 2 − u 3 � q 2 ∂ t u = µu − c + ∆ Include q c and q ′ c = 1 to obtain quasipatterns � 2 (1 + ∆) 2 u + g 1 u 2 − u 3 � q 2 ∂ t u = µu − c + ∆

2D Patterns produced by Swift-Hohenberg equation Stripes Hexagons Zigzag instability Quasicrystals

Snaking in 1D Swift-Hohenberg Equation

Thermosolutal Convection: Patterns with 1D snaking

Thermosolutal Convection: Patterns with 2D snaking

Newton’s method: example Swift-Hohenberg equation: � 2 U − U 3 � q 2 ∂ t U = F ( U ) = µU − c + ∆ Equation for steady state: � 2 U − U 3 � q 2 0 = F ( U ) = µU − c + ∆ Loop: calculate and compare with ǫ : � 2 U − U 3 || < ǫ ? � q 2 || F ( U ) || ≡ || µU − c + ∆ If || F ( U ) || ✚ <ǫ , then U not solution, so try U − u : ✚ � 2 ( U − u ) − ( U − u ) 3 � q 2 0 = µ ( U − u ) − c + ∆ � � 2 u − 3 U 2 u − 3 Uu 2 − u 3 � � 2 U − U 3 − � � q 2 q 2 = µU − c + ∆ µu − c + ∆ Newton step: truncate at first order in u and solve for u ( x ) : � � 2 − 3 U 2 � � 2 U − U 3 � � q 2 q 2 µ − c + ∆ u = µU − c + ∆ Then replace and try again: U ← U − u

Continuation: going around saddle-nodes

Goal: 0 = RN ( U ) + LU � � U i some component 0 = p ( U, R ) − ¯ p where R Newton step: ( U, R ) not solution, so try ( U − u, R − r ) 0 = ( R − r ) N ( U − u ) + L ( U − u ) = RN ( U ) + LU − RN U u − rN ( U ) − Lu + O ( r, u ) 2 � U i − ¯ � p − u i 0 = p ( U − u, R − r ) − ¯ p = R − ¯ p − r

  � u RN ( U ) + LU � � � � U i − ¯ RN U + L N ( U ) � = p   r 0 0 . . . 0 1 0 . . . 0 1 R − ¯ p � �� � or If p ( U, R ) = R (i.e. set Reynolds number), then set R = ¯ p , r = 0 and get previous case: [ RN U + L ] [ u ] = [ RN ( U ) + LU ] If p ( U, R ) = U i , then must solve extended system for ( u, r ) . � � u � ( RN ( U ) + LU ) � � � ( RN u + L ) N ( U ) = 0 0 . . . 0 1 0 . . . 0 0 r U i − ¯ p  Set u i = U i − ¯ p  Calculate ( RN U + L ) u  Add N ( U ) r

It may be sufficient to extrapolate quadratically. Far from saddle-node bifurcation, U is considered to be a function of R . To get an initial guess for U at a new R , extrapolate. Zeroth order extrapolation: U ( R (2) ) initial guess = U ( R (1) ) Linear extrapolation: U ( R (2) ) initial guess = U ( R (1) ) + ( U ( R (1) ) − U ( R (0) )) R (2) − R (1) R (1) − R (0) Quadratic extrapolation: Fit quadratic polynomials through U ( R (0) ) , U ( R (1) ) , U ( R (2) ) as func- tions of R and evaluate polynomial at new value R (3) . Close to saddle-node bifurcation, choose distinguished value U i and con- sider U j , j � = i and R to be quadratic functions of U i . Set new value of U i , and evaluate new estimate of U j and R Now, ∆ R can change sign and can go around saddle-node.

Reaction-Diffusion Equations ∂ t u i = f i ( u 1 , u 2 , . . . ) + D i ∆ u i � �� � � �� � diffusion reaction Reactions f i couple different species u i at same location Diffusivity D i couples same species u i at different locations Describe oscillating chemical reactions, such as famous Belousov-Zhabotinskii reaction, discovered by two Soviet scientists in 1950s-1960s. Also describe phenomena in –biology (population biology, epidemiology, neurosciences) –social sciences (economics, demography) –physics

Two species Spatially homogeneous ∂ t u = f ( u, v ) + D u ∆ u ∂ t u = f ( u, v ) ∂ t v = g ( u, v ) + D v ∆ v ∂ t v = g ( u, v ) FitzHugh-Nagumo model Barkley model � � f ( u, v ) = 1 u − v + b f ( u, v ) = u − u 3 / 3 − v + I ǫ u (1 − u ) a g ( u, v ) = 0 . 08 ( u + 0 . 7 − 0 . 8 v ) g ( u, v ) = u − v u -nullclines f ( u, v ) = 0 , v -nullclines g ( u, v ) = 0 , • steady states � f u f v � stable if eigenvalues of have negative real parts g u g v

� � f ( u, v ) = 1 u − v + b ǫ u (1 − u ) g ( u, v ) = u − v a O ( ǫ − 1 ) ∂ t u = f = 0 ← − and − → separates ∂ t v = g = 0 ↑ ↓ O (1) separates and u = 1 excited phase u = 0 v ∼ 1 refractory phase u = 0 v ≪ 1 excitable phase u = ( v + b ) /a excitation threshold

Waves in Excitable Medium Spatial variation + diffusion + excitability = ⇒ propagating waves Excitable media in physiology: –neurons –cardiac tissue (the heart) Pacemaker periodically emits electrical signals, propagated to rest of heart

Simulations from Barkley model , Scholarpedia Spiral waves in 2D Spiral waves in 3D

T RAVELING WAVES : U ( x − Ct, y, z ) � 0 = C∂ x U + N ( U ) + LU Goal : 0 = ∂ x U ( x ∗ ) example of phase condition Newton step: ( U, C ) not solution, so try ( U − u, C − c ) 0 = ( C − c ) ∂ x ( U − u ) + N ( U − u ) + L ( U − u ) = C∂ x U + N ( U ) + LU − C∂ x u − c∂ x U − N U u − Lu 0 = ∂ x U ( x ∗ ) − ∂ x u ( x ∗ ) � C∂ x + N U + L ∂ x U � � u � C∂ x U + N ( U ) + LU � � = ∂ x | x = x ∗ 0 c ∂ x | x = x ∗ 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 Lu + N U u = LU + N ( U )