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)

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

Un+1 = Un − F (Un) F ′(Un) F (U) = 0 = F (Un) + F ′(Un)(U − Un) + 1

2F ′′(Un)(U − Un)2 + . . .

0 = F (Un) F ′(Un) + F ′(Un) F ′(Un)(U − Un) + 1

2

F ′′(Un) F ′(Un) (U − Un)2 + . . . 0 = F (Un) F ′(Un) + (U−Un) + 1

2

F ′′(Un) F ′(Un) (U − Un)2 + . . . 0 = −Un+1 + U + 1

2

F ′′(Un) F ′(Un) (U − Un)2 + . . . Un+1 − U = 1

2

F ′′(Un) F ′(Un) (U − Un)2 + . . . ǫn+1 = 1

2

F ′′(Un) F ′(Un) ǫ2 + . . . Typical sequence: ǫ = 10−1, 10−2, 10−4, 10−8, 10−16

Much faster than timestepping: U(t) = U + ce−λt U(tn) − U = ce−λtn U(tn+1) − U = ce−λ(tn+∆t) = e−λ∆t(U(tn) − 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 ˙ ǫ = f(¯ x) + f ′(¯ x)ǫ + 1 2f ′′(¯ x)ǫ2 + · · · ≈ f ′(¯ x)ǫ ǫ(t) = etf ′(¯

x)ǫ(0)

increases if f ′(¯ x) > 0 decreases if f ′(¯ x) < 0

Saddle-node Bifurcations

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

f(x, µ) = c00 + c10x + c01µ + c20x2 + . . . general quadratic polynomial = ±˜ µ ± ˜ x2 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

∂tu = µu −

• q2

c + ∆

2 u − u3 Derived by J. Swift and P.C. Hohenberg (Phys. Rev. A 15, 319 (1977)) to describe pattern formation in convection For u ≪ 1, ∂tu = µu −

• q2

c + ∆

2 u u ∼ exp(ikx + σt) σu = µu −

• q2

c − k22 u =

⇒ σ = µ, k = qc Add quadratic term to obtain hexagons ∂tu = µu −

• q2

c + ∆

2 u + g1u2 − u3 Include qc and q′

c = 1 to obtain quasipatterns

∂tu = µu −

• q2

c + ∆

2 (1 + ∆)2 u + g1u2 − u3

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: ∂tU = F (U) = µU −

• q2

c + ∆

2 U − U 3 Equation for steady state: 0 = F (U) = µU −

• q2

c + ∆

2 U − U 3 Loop: calculate and compare with ǫ: ||F (U)|| ≡ ||µU −

• q2

c + ∆

2 U − U 3|| < ǫ ? If ||F (U)||✚

✚

<ǫ, then U not solution, so try U − u: 0 = µ(U − u) −

• q2

c + ∆

2 (U − u) − (U − u)3 = µU −

• q2

c + ∆

2 U − U 3 −

• µu −
• q2

c + ∆

2 u − 3U 2u − 3Uu2 − u3 Newton step: truncate at first order in u and solve for u(x):

• µ −
• q2

c + ∆

2 − 3U 2 u = µU −

• q2

c + ∆

2 U − U 3 Then replace and try again: U ← U − u

Continuation: going around saddle-nodes

Goal:

0 = RN(U) + LU 0 = p(U, R) − ¯ p where

• Ui some component

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 − RNUu − rN(U) − Lu + O(r, u)2 0 = p(U − u, R − r) − ¯ p = Ui − ¯ p − ui R − ¯ p − r

• RNU + L

N(U) 0 0 . . . 0 1 0 . . . 0 1

• u

r

• =

  RN(U) + LU Ui − ¯ p R − ¯ p

• 



• r

If p(U, R) = R (i.e. set Reynolds number), then set R = ¯ p, r = 0 and get previous case: [RNU + L] [u] = [RN(U) + LU] If p(U, R) = Ui, then must solve extended system for (u, r).

• (RNu + L)

N(U) 0 0 . . . 0 1 0 . . . 0 u r

• =

(RN(U) + LU) Ui − ¯ p

• Set ui = Ui − ¯

p Calculate (RNU + 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 Ui and con- sider Uj, j = i and R to be quadratic functions of Ui. Set new value of Ui, and evaluate new estimate of Uj and R Now, ∆R can change sign and can go around saddle-node.

Reaction-Diffusion Equations ∂tui = fi(u1, u2, . . .)

• reaction

+ Di∆ui diffusion Reactions fi couple different species ui at same location Diffusivity Di couples same species ui 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 ∂tu = f(u, v) + Du∆u ∂tu = f(u, v) ∂tv = g(u, v) + Dv∆v ∂tv = g(u, v) FitzHugh-Nagumo model Barkley model f(u, v) = u − u3/3 − v + I f(u, v) = 1

ǫ u(1 − u)

• u − v+b

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 stable if eigenvalues of fu fv gu gv

• have negative real parts

f(u, v) = 1

ǫ u(1 − u)

• u − v+b

a

• g(u, v) = u − v

∂tu = f = 0 separates ← − and − → O(ǫ−1) ∂tv = g = 0 separates ↑ and ↓ O(1) 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

TRAVELING WAVES: U(x − Ct, y, z)

Goal : 0 = C∂xU + N(U) + LU 0 = ∂xU(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∂xU + N(U) + LU − C∂xu − c∂xU − NUu − Lu 0 = ∂xU(x∗) − ∂xu(x∗) C∂x + NU + L ∂xU ∂x|x=x∗ u c

• =

C∂xU + N(U) + LU ∂x|x=x∗U

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 Lu + NUu = LU + N(U)