Cours : Dynamique Non-Lin eaire Laurette TUCKERMAN - - PowerPoint PPT Presentation

Cours : Dynamique Non-Lin´ eaire Laurette TUCKERMAN laurette@pmmh.espci.fr

• VII. Reaction-Diffusion Equations:
• 1. Excitability
• 2. Turing patterns
• 3. Lyapunov functionals
• 4. Spatial analysis and fronts

Reaction-Diffusion Systems

∂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

Excitability

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

Turing patterns

Instability of homogeneous solutions (¯ u, ¯ v) to reaction-diffusion systems 0 = f(¯ u, ¯ v) 0 = g(¯ u, ¯ v)

• =

⇒ 0 = f(¯ u, ¯ v) + Du∆¯ u 0 = g(¯ u, ¯ v) + Dv∆¯ v

• What about stability? Does diffusion damp spatial variations?

Linear stability analysis: u(x, t) = ¯ u + ˜ ueσt+ik·x v(x, t) = ¯ v + ˜ veσt+ik·x

• =

⇒ σ˜ u = fu˜ u + fv˜ v − Duk2˜ u σ˜ v = gu˜ u + gv˜ v − Dvk2˜ v

• Mk ≡

fu − Duk2 fv gu gv − Dvk2

• =

fu fv gu gv

• − k2

Du 0 Dv

• If Du = Dv ≡ D, then

σk± = σ0± − k2D ≤ σ0± (¯ u, ¯ v) stable to homogeneous perturbations = ⇒ (¯ u, ¯ v) stable to inhomogeneous perturbations. Diffusion is stabilizing.

Alan Turing (famous WW II UK cryptologist, founder of computer science) 1952: homogeneous state can be unstable if Du = Dv For instability, need Trk > 0 or Detk < 0

For instability, need Trk > 0 or Detk < 0 Homogeneous stability ⇐ ⇒ Tr0 = fu + gv < 0 and Det0 = fugv − fvgu > 0

• Trk = fu + gv − (Du + Dv)k2 = Tr0 − (Du + Dv)k2 < Tr0 < 0

So for instability, need Detk < 0 Detk = fugv − fvgu + DuDvk4 − (Dvfu + Dugv)k2 = Det0

>0

+ DuDvk4

• >0, dominates for k≫1

−(Dvfu + Dugv)k2 Find negative minimum for intermediate k2: 0 = d Detk dk2

• k∗

= 2DuDvk2

∗ − (Dvfu + Dugv)

k2

∗ = Dvfu + Dugv

2DuDv = ⇒ need Dvfu + Dugv > 0

Need Detk < 0 at k2

∗ = (Dvfu + Dugv)/(2DuDv):

0 > Detk|k∗ = Det0 + DuDvk4

∗ − (Dvfu + Dugv)k2 ∗

= Det0 + (Dvfu + Dugv)2 4DuDv − 2(Dvfu + Dugv)2 4DuDv = Det0 − (Dvfu + Dugv)2 4DuDv 0 > 4DuDv(fugv − fvgu) − (Dvfu + Dugv)2 Collecting the four conditions: Tr0 = fu + gv < 0 Det0 = fugv − fvgu > 0 2DuDvk2

∗ = Dvfu + Dugv > 0

4DuDv Detk|k∗ = 4DuDv(fugv − fvgu) − (Dvfu + Dugv)2 < 0

Turing patterns were first produced experimentally: –in 1990 by de Kepper et al. at Univ. of Bordeaux –in 1992 by Swinney et al. at Univ. of Texas at Austin Turing pattern in a chlorite- iodide-malonic acid chemical laboratory experiment. From R.D. Vigil, Q. Ouyang & H.L. Swinney, Turing patterns in a simple gel reactor, Physica A 188, 17 (1992) Might be mechanism for: –differentiation within embryos –formation of patterns on animal coats, e.g. zebras and leopards

Lyapunov functionals

1D systems: no limit cycles, usually just convergence to fixed point Generalize to multidimensional variational, potential, or gradient flows: du dt = −∇Φ ⇐ ⇒ dui dt = − ∂Φ ∂ui For gradient flow, Jacobian is Hessian matrix: H = −    ∂2Φ/(∂u1∂u1) ∂2Φ/(∂u1∂u2) . . . ∂2Φ/(∂u2∂u1) ∂2Φ/(∂u2∂u2) . . . . . . . . . . . .    H symmetric = ⇒ no complex eigenvalues = ⇒ no Hopf bifurcations dΦ dt =

• i

∂Φ ∂ui dui dt = −

• i

∂Φ ∂ui ∂Φ ∂ui = −|∇Φ|2 Φ decreases monotonically, either to −∞ or to point where du/dt = −∇Φ = 0 = ⇒ no limit cycles

Generalize to reaction-diffusion systems involving potential Φ(u): ∂u ∂t = −∇Φ + ∂2u ∂x2

• n xlo ≤ x ≤ xhi

Boundary conditions: Dirichlet u(xlo) = ulo u(xhi) = uhi

• r Neumann (homogeneous)

∂u ∂x(xlo) = 0 ∂u ∂x(xhi) = 0

Define free energy or Lyapunov functional: F(u) ≡ xhi

xlo

dx

• Φ(u(x, t))
• potential energy

+ 1 2

• ∂u(x, t))

∂x

• 2
• kinetic energy
• Seek quantity analogous to gradient:

F (x + dx) = F(x) + ∇F(x) · dx + O(|dx|)2 for all dx The functional derivative δF/δu is defined to be such that F(u + δu) = F(u) + xhi

xlo

dx δF δu · δu + O(δu)2 for every δu

Expand: F(u + δu) = xhi

xlo

dx

• Φ(u + δu) + 1

2

• ∂(u + δu)

∂x

• 2

= xhi

xlo

dx

• Φ(u) + ∇Φ(u) · δu + . . . + 1

2

• ∂u

∂x + ∂δu ∂x + . . .

• 2

= xhi

xlo

dx

• Φ(u) + 1

2

• ∂u

∂x

• 2

+ xhi

xlo

dx

• ∇Φ(u) · δu + ∂u

∂x · ∂δu ∂x

• + O(δu)2

Integrate by parts: xhi

xlo

dx ∂u ∂x · ∂δu ∂x = ∂u ∂x · δu xhi

xlo

− xhi

xlo

dx∂2u ∂x2 · δu Surface term vanishes since ∂u

∂x(xlo) = ∂u ∂x(xhi) = 0

for Neumann BCs δu(xlo) = δu(xhi) = 0 for Dirichlet BCs

F(u+δu)= xhi

xlo

dx

• Φ(u) +
• ∂u

∂x

• 2

+ xhi

xlo

dx

• ∇Φ(u) − ∂2u

∂x2

• · δu+O(δu)2

The functional derivative δF/δu is defined to be such that F(u + δu) = F(u) + xhi

xlo

dx δF δu · δu + O(δu)2 for every δu = ⇒ xhi

xlo

dx δF δu · δu = xhi

xlo

dx

• ∇Φ(u) − ∂2u

∂x2

• · δu

Choosing δu to be delta function centered on any x and pointing in any vector direction leads to pointwise equality: δF δu = ∇Φ(u) − ∂2u ∂x2 = −∂u ∂t

dF dt = lim

∆t→0

1 ∆t [F(t + δt) − F(t)] = lim

∆t→0

1 ∆t [F(u(t + ∆t)) − F(u(t))] = lim

∆t→0

1 ∆t

• F
• u(t) + ∂u

∂t ∆t + . . .

• − F(u(t))
• =

lim

∆t→0

1 ∆t

• F(u(t)) +

xhi

xlo

dxδF δu · ∂u ∂t ∆t + . . . − F(u(t))

• =

lim

∆t→0

1 ∆t xhi

xlo

dxδF δu · ∂u ∂t ∆t + . . .

• =

xhi

xlo

dxδF δu · ∂u ∂t = xhi

xlo

dx

• −∂u

∂t

• · ∂u

∂t = − xhi

xlo

dx

• ∂u

∂t

• 2

≤ 0 F decreases so limit cycles cannot occur. Can be applied in higher spatial dimensions via volume integration and Gauss’s Divergence Theorem.

Spatial Analysis and Fronts

∂u ∂t = −dΦ du + ∂2u ∂x2 Travelling wave solutions: u(x, t) = U(x − ct) with c = 0 for steady states ξ ≡ x − ct ∂u ∂t (x, t) = −c dU dξ (ξ) ∂2u ∂x2(x, t) = d2U dξ2 (ξ) Equation obeyed by steady states and travelling waves becomes −c du dξ = −dΦ du + d2u dξ2 = ⇒ d2u dξ2 = dΦ du−c du dξ Analogy between space and time = ⇒ x must be 1D

Spatial analysis or Mechanical analogy

d2u dξ2

• “acceleration”

= − d(−Φ) du “potential gradient” −c du dξ “friction” u position ξ time

du dξ

velocity −Φ potential E(ξ) ≡ −Φ + 1

2

• du

dξ

2 energy

d2u dξ2

acceleration −cdu

dξ

friction ˙ E = dE dξ = d dξ

• −Φ + 1

2 du dξ 2 = −dΦ du du dξ + du dξ d2u dξ2 =

• −dΦ

du + d2u dξ2 du dξ = −c du dξ 2    < 0 if c > 0 = 0 if c = 0 > 0 if c < 0 c < 0 ⇐ ⇒ “Increase in energy” “Negative friction”

• ⇐

⇒ just leftwards motion

If c = 0, then E constant with E = −Φ(u(ξ)) + 1 2 du dξ 2 E + Φ(u(ξ)) = 1 2 du dξ 2

• 2(E + Φ(u(ξ))) = du

dξ

• dξ =
• du
• 2(E + Φ(u))

[ξ] ξ

ξlo =

u(ξ)

ulo

du

• 2(E + Φ(u))

= elliptic integral if Φ(u) = u3 = ⇒ ξ(u) = ⇒ u(ξ) yields results but no intuition

Dynamical systems approach with ξ as time

v ≡ du dξ = ⇒ ˙ u = v ˙ v = dΦ

du − cv

If c = 0, then system is Hamiltonian: H = −Φ + 1 2v2 = ⇒ ˙ u = ∂H

∂v

˙ v = −∂H

∂u

Add diffusion to supercritical pitchfork = ⇒ Ginzburg-Landau equation: ∂u ∂t = µu − u3 + ∂2u ∂x2 Steady states 0 = µu − u3 + d2u dx2 Integrate to obtain the potential: −dΦ du = µu − u3 = ⇒ −Φ = µ 2 u2 − 1 4u4

Steady states: d2u dx2 = dΦ du = ⇒ ˙ u = v ˙ v = dΦ

du

Fixed points of new dynamical system: 0 = v 0 = dΦ du = −µ¯ u + ¯ u3 = ⇒ ¯ u = 0

• r

¯ u = ±√µ Same ¯ u as without diffusion, but stability under new dynamics is different: J = 1 Φ′′ 0

• =
• 1

3¯ u2 − µ 0

• =
• 0 1

−µ 0

• r
• 0 1

2µ 0

• Hamiltonian ⇐

⇒ Tr(J ) =

∂2H ∂u∂v − ∂2H ∂v∂u = 0 ⇐

⇒ eigs are ±λ λ(−λ) = −Φ′′ = ⇒ λ± = ± √ Φ′′ = ±√−µ for ¯ u = 0 ±√2µ for ¯ u = ±√µ λ = ±iω = ⇒ center = elliptic fixed point λ = ±σ = ⇒ saddle = hyperbolic fixed point

µ = −1 µ = +1

µ = +1

Types of Trajectories

µ = −1 µ = +1 unbounded, crossing between left to right

• unbounded, staying on left or on right
• periodic
• front (limiting case of periodic)
• Periodic:

Trajectories in the (u, ˙ u) phase plane are elliptical. Particle oscillates back and forth in potential well. Fronts: Trajectory leaves ¯ u = −√µ at zero velocity, arrives exactly at ¯ u = √µ with zero velocity, since there is no friction. Profile has u = −√µ on left, narrow transition region, u = √µ on right. Type of trajectory is determined by the initial conditions (temporal point of view) or the boundary conditions (spatial point of view). Periodic bound- ary conditions on a domain of fixed wavelength select the periodic profile. Boundary conditions u(±∞) = ±√µ lead to front solution.

Front solutions connect two maxima of −Φ, i.e. hyperbolic unstable fixed points of the transformed dynamical system. These correspond to stable spatially homogeneous solutions to the original reaction-diffusion system: du dt = −dΦ du Stability determined by −d2Φ du2 (¯ u) < 0 > 0

• =

⇒ ¯ u

• stable

unstable

• Thus, homogeneous stable steady states are maxima of −Φ.

Nonzero c

Front between u−∞ and u+∞, which are maxima of −Φ(u) and hence stable solutions to spatially homogeneous equations, Dirichlet BCs u(ξ = ±∞) = u±∞ = ⇒ Neumann BCs du dξ (ξ = ±∞) = 0 Travelling wave solutions: 0 = cdu dξ − dΦ du + d2u dξ2 Multiply by du/dξ: 0 = c du dξ 2 − dΦ(u(ξ)) dξ + 1 2 d dξ du dξ 2 Integrate over ξ interval: 0 = c +∞

−∞

dξ du dξ 2 − +∞

−∞

dξ dΦ(u(ξ)) dξ + +∞

−∞

dξ 1 2 d dξ du dξ 2 = c +∞

−∞

dξ du dξ 2 − [Φ]+∞

−∞ + 1

2 du dξ 2+∞

−∞

⇐ vanishes because of Neumann BCs

c = Φ+∞ − Φ−∞ +∞

−∞ dξ

• du

dξ

2 where Φ±∞ ≡ Φ(u±∞) Front velocity c > 0 if Φ−∞ < Φ+∞, i.e. if −Φ−∞ > −Φ+∞. Front moves from left to right = ⇒ u−∞, −Φ−∞ domain invades u+∞, −Φ+∞ domain Front motion increases size of domain with greater −Φ. Mechanical analogy: Trajectory goes from u−∞, −Φ−∞ to u+∞, with lower potential −Φ+∞. For “velocity” du/dξ and “kinetic energy” to vanish at both endpoints, energy must be lost via friction. Hence c is positive. “Negative friction” is possible since c < 0 just means that the front moves towards the left.

Trajectory Phase portrait from lower left hill to higher right hill Former center has become focus uses “negative friction” to increase its energy surrounded by spiralling trajectories

Perturbed Ginzburg-Landau equation 0 = cdu dξ + µu − u3 − 0.1 + d2u dξ2 Potential −Φ = 1 2µu2 − 1 4u4 − 0.1u has two maxima of different heights