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

Temporal Discretization

du dt = f(u) u, f in RN Goal: turn differential equation into difference equation First order methods: Forwards (Explicit) Euler: uFE(t + ∆t) = uFE(t) + ∆tf(uFE(t)) explicit: “=” is an assignment statement Backwards (Implicit) Euler: uBE(t + ∆t) = uBE(t) + ∆tf(uBE(t + ∆t)) implicit: “=” is an equation to be solved for uBE(t + ∆t)

Taylor series:

u(t + ∆t) = u(t) + ∆tdu dt + ∆t2 2 d2u dt2 + . . . = u(t) + ∆tf(u(t)) + ∆t2 2 f ′(u(t))f(u(t)) + . . . Forwards (Explicit) Euler: uFE(t + ∆t) = uFE(t) + ∆tf(uFE(t)) Backwards (Implicit) Euler: uBE(t + ∆t) = uFE(t) + ∆tf(uBE(t + ∆t)) = uBE(t) + ∆t

• f(uBE(t)) + ∆tf ′(uBE(t))f(uBE(t)) + . . .
• = uBE(t) + ∆tf(uBE(t)) + ∆t2f ′(uBE(t))f(uBE(t)) + . . .

First order: ∆t terms match Taylor series but not ∆t2 terms

Second order methods: Adams-Bashforth (explicit) uAB(t+∆t) = uAB(t)+∆t 3 2f(uAB(t)) − 1 2f(uAB(t − ∆t))

• Crank-Nicolson (implicit) also called trapezoidal

uCN(t+∆t) = uCN(t)+∆t 1 2f(uCN(t)) + 1 2f(uCN(t + ∆t))

• Backwards Differentiation (implicit)

uBD(t+∆t) = 4 3uBD(t)−1 3uBD(t−∆t)+2 3∆tf(uBD(t+∆t)) Second order: ∆t, ∆t2 terms match Taylor series but not ∆t3

For fixed T , take T/∆t steps If one-step error is ∆tp, total error at time T is (T/∆t)∆tp = T ∆tp−1 First-order methods have one-step error ∆t2 and fixed-time error ∆t Second-order methods have one-step error ∆t3 and fixed-time error ∆t2 Have assumed f(u(t)), but can also have f(u(t), t) Second or higher order differential equations: Write u0 = u, u1 = u′

0, u2 = u′′ 0, etc.

du dt = f(u) u = (u0, u1, . . .), f = (f0, f1, . . .)

Stability: example of heat equation with periodic boundary conditions

∂tu = ∂xxu u(x, t) =

kmax

• k=1

ˆ uk(t) sin kx ∂tˆ uk = −k2ˆ uk

EXACT SOLUTION

ˆ uk(t + ∆t) = e−k2∆tˆ uk(t)

EXPLICIT EULER

ˆ uFE

k (t + ∆t) = ˆ

uFE

k (t) − k2∆tˆ

uFE

k (t)

= (1 − k2∆t)ˆ uFE

k (t)

As kmax → ∞, ∆tmax =

2 k2

max → 0

IMPLICIT EULER

ˆ uBE

k (t + ∆t) = ˆ

uBE

k (t) − k2∆tˆ

uBE

k (t + ∆t)

(1 + k2∆t)ˆ uBE

k (t + ∆t) = ˆ

uBE

k (t)

ˆ uBE

k (t + ∆t) = (1 + k2∆t)−1ˆ

uBE

k (t)

Matrix version: uBE(t + ∆t) = (I − ∆tL)−1uBE(t)

Crank-Nicolson: uCN(t + ∆t) = uCN(t) + ∆t 2

• f(uCN(t)) + f(uCN(t + ∆t))
• uCN(t + ∆t) − ∆t

2 f(uCN(t + ∆t)) = uCN(t) + ∆t 2 f(uCN(t))

• 1 + k2∆t

2

• ˆ

uCN

k

(t + ∆t) =

• 1 − k2∆t

2

• ˆ

uCN

k

(t) ˆ uCN

k

(t + ∆t) = 1 − k2∆t

2

1 + k2∆t

2

ˆ uCN

k

(t) |Amp. factor| < 1 but spurious large-k behavior: slow oscillatory decay

A-stable methods:

du dt = −k2u = ⇒ uexact(t) = e−k2t = ⇒ lim

t→∞ uexact(t) = 0

A-stable method = ⇒ lim

t→∞ unumerical(t) = 0

Define q ≡ −k2∆t and generalize to q complex Define amplification factor Φ(q) Consider behavior for Re(q) < 0 A-stable: |Φ(q)| < 1 L-stable: |Φ(q)| → 0 as q → ∞

Crank-Nicolson: seek stability region in complex plane, where

• ΦCN(q)
• ≡
• 1 + q/2

1 − q/2

• ≤ 1

|1 + q/2| ≤ |1 − q/2| q is closer to −2 than to 2: q is in left half plane A-stable ⇐ ⇒ stability region contains negative real axis ⇐ ⇒ For Re(q) < 0, have |Φ(q)| < 1 YES L-stable ⇐ ⇒ For Re(q) < 0, have |Φ(q)| → 0 as q → ∞ × NO

Forwards Euler Backwards Euler |1 + q| < 1 1/|1 − q| < 1 |1 − q| > 1 For q real, require −2 < q < 0 All q < 0 in stability region

Adams-Bashforth:

uAB(t + ∆t) = uAB(t) + ∆t 2

• 3f(uAB(t)) − f(uAB(t − ∆t))
• = uAB(t) + q

3 2uAB(t) − 1 2uAB(t − ∆t)

• uAB(t + ∆t)

uAB(t)

• =

1 + 3q

2

−q

2

1 uAB(t) uAB(t − ∆t)

• Require that both eigenvalues obey |λ| < 1

Eigenvalues λ obey characteristic polynomial 0 =

• 1 + 3q

2 − λ

• (−λ) + q

2 = λ2 − λ

• 1 + 3q

2

• + q

2 q 2 = λ2 − λ 3λ − 1

Set λ = eiθ: q 2 = λ2 − λ 3λ − 1 = e2iθ − eiθ 3eiθ − 1 = cos 2θ − cos θ + i(sin 2θ − sin θ) (3 cos θ − 1) + 3i sin θ = (cos 2θ − cos θ) (3 cos θ − 1) + (sin 2θ − sin θ)3 sin θ (3 cos θ − 1)2 + (3 sin θ)2 + i (cos 2θ − cos θ) (−3 sin θ) + (sin 2θ − sin θ)(3 cos θ − 1) (3 cos θ − 1)2 + (3 sin θ)2 not A-stable

Backwards Differentiation:

uBD(t + ∆t) = 4 3uBD(t) − 1 3uBD(t − ∆t) + 2 3∆tf(uBD(t + ∆t)) = 4 3uBD(t) − 1 3uBD(t − ∆t) + 2 3quBD(t + ∆t)

• 1 − 2

3q

• uBD(t + ∆t) = 4

3uBD(t) − 1 3uBD(t − ∆t) uBD(t + ∆t) = 1 1 − 2q

3

4 3uBD(t) − 1 3uBD(t − ∆t)

• uBD(t + ∆t)

uBD(t)

• =
• 4

3−2q −1 3−2q

1 uBD(t) uBD(t − ∆t)

• Eigenvalues λ obey characteristic polynomial

0 =

• 4

3 − 2q − λ

• (−λ) +

1 3 − 2q (2q − 3) = −4 λ + 1 λ2

Require that both eigenvalues obey |λ| < 1 so set λ = eiθ (2q − 3) = − 4 eiθ + 1 e2iθ = −4(cos θ − i sin θ) + (cos 2θ − i sin 2θ) q = 1 2 (3 − 4 cos θ + cos 2θ + i[4 sin θ − sin 2θ]) A-stable

General formalism:

s

• j=0

αjun+1−j = ∆t

s

• j=0

βjf(un+1−j)

Degrees of freedom {αj}, {βj} allow many routes to order-p accuracy. Scale by setting α0 = 1. Explicit ⇐ ⇒ β0 = 0: un+1 =

s

• j=1
• −αjun+1−j + ∆tβjf(un+1−j)
• Adams-Bashforth (explicit): α0 = 1, α1 = −1, αj≥2 = 0, β0 = 0

Select β1≤j≤p to achieve p-order accuracy. Adams-Moulton (implicit): α0 = 1, α1 = −1, αj≥2 = 0 Select β0≤j≤p−1 to achieve p-order accuracy. Backwards-Differentiation (implicit): α0 = 1, βj≥1 = 0 Select β0, α1≤j≤p to achieve p-order accuracy.

Stability regions of Adams-Bashforth formulas:

Forwards Euler

Stability regions of Adams-Moulton formulas:

Backwards Euler Crank-Nicolson

Stability regions of Backwards Differentiation formulas:

Backwards Euler

• rder 1
• rder 2
• rder 3
• rder 4

In general, increasing accuracy = ⇒ decreasing stability A-stable methods must be implicit and at most second order. Explicit methods approximate exponential by polynomial: Necessarily grow as q → ±∞ Implicit methods approximate the exponential by a rational. What about using the exponential itself? This is sometimes possible if the operator: –is linear –can be cheaply diagonalized

Exponential of a matrix with real eigenvalues

eLt = I + tL + t2 2 L2 + . . . = V V −1 + tV ΛV −1 + t2 2 V ΛV −1 V ΛV −1+ = V

• I + tΛ + t2

2 Λ2 + . . .

• V −1 = V eΛtV −1

Λ2 = λ1 λ2 λ1 λ2

• =

λ2

1

λ2

2

• etΛ =
• 1 + tλ1 + (tλ1)2

2

+ . . . 1 + tλ2 + (tλ2)2

2

+ . . .

• =

etλ1 etλ2

Imaginary Eigenvalues

L = 0 −ω ω

• λ± = ±iω

L2 = 0 −ω ω 0 −ω ω

• =

−ω2 −ω2

• L3 =

0 −ω ω −ω2 −ω2

• =
• ω3

−ω3

• eLt = I + tL + t2

2 L2 + t3 6 L3 + . . . = 1 0 0 1

• + t

0 −ω ω

• + t2

2 −ω2 −ω2

• + t3

6

• ω3

−ω3

• =

1 − 1

2(tω)2 + . . .

−tω + 1

6(tω)3 + . . .

tω − 1

6(tω)3 + . . .

1 − 1

2(tω)2 + . . .

• =

cos(ωt) − sin(ωt) sin(ωt) cos(ωt)

Complex Eigenvalues: exp

• t

µ −ω ω µ

• = eµt

cos(ωt) − sin(ωt) sin(ωt) cos(ωt)

• Mixed spectrum:

Λ =      λ1 0 µ −ω ω µ 0 λ4      Re(λ1) ≥ Re(λ2) = Re(λ3) ≥ Re(λ4) exp(Λt) =      eλ1t eµt cos(ωt) −eµt sin(ωt) eµt sin(ωt) eµt cos(ωt) 0 eλ4t     

Holds for any analytic function f(L) f(L) =

• j

1 j!f (j)(0)Lj =

• j

1 j!f (j)(0)

• V ΛV −1j

= V  

j

1 j!f (j)(0)Λj   V −1 = V f(Λ)V −1 where f(Λ) = f      λ1 λ2 ... λN      =      f(λ1) f(λ2) ... f(λN)     

du dt = Lu = ⇒ u(t) = eLtu(0) u(0) ↓ multiply by matrix V −1 V −1u(0) ↓ multiply by diagonal matrix eΛt etΛV −1u(0) ↓ multiply by matrix V V etΛV −1u(0) = eLtu(0)

Example: Heat Equation

∂tu = ∂2

xxu

u(x, t) =

kmax

• k=1

ˆ uk(t) sin kx ∂tˆ uk = −k2ˆ uk Begin with u(xj, t = 0) on a grid of values {x0, x1, . . . xn} Take Fourier transform: {u(xj, t = 0)} = ⇒ {ˆ uk(t = 0)} ˆ uk(t) = e−k2tˆ uk(t = 0) Take inverse Fourier transform: {ˆ uk(t)} = ⇒ {u(xj, t)} Complete answer valid for all time.

Not usually possible, since most problems are: –not linear –not easily diagonalized but exponential integration can be combined with other ideas Time-splitting for du

dt = Lu + N(u)

Backwards Euler / Forwards Euler: u(t + ∆t) = u(t) + ∆tLu(t + ∆t) + ∆tN(u(t)) (I − ∆tL)u(t + ∆t) = u(t) + ∆tN(u(t)) u(t + ∆t) = (I − ∆tL)−1 [u(t) + ∆tN(u(t)] Exponential / Forwards Euler: u(t + ∆t) = eL∆tu(t) + eL∆t/2∆tN(u(t)) Exponential / Adams-Bashforth: u(t+∆t) = eL∆tu(t)+eL∆t/2∆t 2 (3N(u(t)) − N(u(t − ∆t)))

du dt = Lu + N(u) Equivalent integral equation: u(t + ∆t) = eL∆t

• u(t) +

t+∆t

t

dτ e−L(τ−t)N(u(τ))

• Approximate e−L(τ−t) by a constant:

value at τ = t : t+∆t

t

dτ N(u(τ)) value at τ = t + ∆t : e−L∆t t+∆t

t

dτ N(u(τ)) value at τ = t + ∆t/2 : e−L∆t/2 t+∆t

t

dτ N(u(τ)) average value : 1 − e−L∆t − 1 L∆t t+∆t

t

dτ N(u(τ))

Integration scheme for N: duN dt = N(uN) ⇔uN(t + ∆t) − u(t) ≈ t+∆t

t

dτ N(uN(τ)) Resulting “slaved exponential” integration scheme u(t + ∆t) = eL∆tu(t) + eL∆t − 1 L∆t (uN(t + ∆t) − u(t)) If N(u) = N is constant, scheme gives exact soln. for all t: u(t) = eLtu(0) + eLt − 1 L N As t → ∞, if L has negative eigs, then u(t) → −L−1N If N is not constant but has a slower timescale than L, then numerical scheme is accurate for all ∆t on timescale of N.

Runge-Kutta Methods

RK2 (two evaluations of f per timestep, order 2) U1 = u(t) U2 = u(t) + ∆tf(U1) u(t + ∆t) = u(t) + ∆t 2 (f(U1) + f(U2)) RK4 (four evaluations of f per timestep, order 4) F1 = f(u(t)) F2 = f

• u(t) + ∆t

2 F1

• F3 = f
• u(t) + ∆t

2 F2

• F4 = f (u(t) + ∆tF3)

u(t + ∆t) = u(t) + ∆t 6 (F1 + 2F2 + 2F3 + F4) Many others

For conservative (Hamiltonian, area-preserving) systems, numerical method should also preserve area. Oscillator: d2u dt2 = −ω2u ⇐ ⇒ du

dt = −ωv dv dt = ωu

d dt u v

• =

0 −ω ω u v

• r

d dt(u + iv) = iω(u + iv) Set q = iω∆t and consider |Φ(q)|

Forwards Euler: ΦFE(iω∆t) = 1 + iω∆t |ΦFE(iω∆t)|2 = 1 + (ω∆t)2 > 1 Increases area Backwards Euler: ΦBE(iω∆t) = 1 1 − iω∆t |ΦBE(iω∆t)|2 = 1 1 + (ω∆t)2 < 1 Decreases area Crank-Nicolson: ΦCN(iω∆t) = 1 + iω∆t/2 1 − iω∆t/2 |ΦCN(iω∆t)|2 = 1 Preserves area