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 u , f in R N dt = f ( u ) Goal: turn differential equation into difference equation First order methods: Forwards (Explicit) Euler: u FE ( t + ∆ t ) = u FE ( t ) + ∆ tf ( u FE ( t )) explicit: “=” is an assignment statement Backwards (Implicit) Euler: u BE ( t + ∆ t ) = u BE ( t ) + ∆ tf ( u BE ( t + ∆ t )) implicit: “=” is an equation to be solved for u BE ( t + ∆ t )

Taylor series: dt + ∆ t 2 d 2 u u ( t + ∆ t ) = u ( t ) + ∆ tdu dt 2 + . . . 2 = u ( t ) + ∆ tf ( u ( t )) + ∆ t 2 2 f ′ ( u ( t )) f ( u ( t )) + . . . Forwards (Explicit) Euler: u FE ( t + ∆ t ) = u FE ( t ) + ∆ tf ( u FE ( t )) Backwards (Implicit) Euler: u BE ( t + ∆ t ) = u FE ( t ) + ∆ tf ( u BE ( t + ∆ t )) = u BE ( t ) + ∆ t f ( u BE ( t )) + ∆ tf ′ ( u BE ( t )) f ( u BE ( t )) + . . . � � = u BE ( t ) + ∆ tf ( u BE ( t )) + ∆ t 2 f ′ ( u BE ( t )) f ( u BE ( t )) + . . . First order: ∆ t terms match Taylor series but not ∆ t 2 terms

Second order methods: Adams-Bashforth (explicit) � 3 2 f ( u AB ( t )) − 1 � u AB ( t +∆ t ) = u AB ( t )+∆ t 2 f ( u AB ( t − ∆ t )) Crank-Nicolson (implicit) also called trapezoidal � 1 2 f ( u CN ( t )) + 1 � u CN ( t +∆ t ) = u CN ( t )+∆ t 2 f ( u CN ( t + ∆ t )) Backwards Differentiation (implicit) u BD ( t +∆ t ) = 4 3 u BD ( t ) − 1 3 u BD ( t − ∆ t )+2 3∆ tf ( u BD ( t +∆ t )) Second order: ∆ t , ∆ t 2 terms match Taylor series but not ∆ t 3

For fixed T , take T/ ∆ t steps If one-step error is ∆ t p , total error at time T is ( T/ ∆ t )∆ t p = T ∆ t p − 1 First-order methods have one-step error ∆ t 2 and fixed-time error ∆ t Second-order methods have one-step error ∆ t 3 and fixed-time error ∆ t 2 Have assumed f ( u ( t )) , but can also have f ( u ( t ) , t ) Second or higher order differential equations: Write u 0 = u , u 1 = u ′ 0 , u 2 = u ′′ 0 , etc. du dt = f ( u ) u = ( u 0 , u 1 , . . . ) , f = ( f 0 , f 1 , . . . )

Stability: example of heat equation with periodic boundary conditions ∂ t u = ∂ xx u k max � u ( x, t ) = u k ( t ) sin kx ˆ k =1 u k = − k 2 ˆ ∂ t ˆ u k E XACT SOLUTION u k ( t + ∆ t ) = e − k 2 ∆ t ˆ ˆ u k ( t )

E XPLICIT E ULER u FE u FE k ( t ) − k 2 ∆ t ˆ u FE ˆ k ( t + ∆ t ) = ˆ k ( t ) = (1 − k 2 ∆ t )ˆ u FE k ( t ) 2 As k max → ∞ , ∆ t max = max → 0 k 2 I MPLICIT E ULER u BE u BE k ( t ) − k 2 ∆ t ˆ u BE ˆ k ( t + ∆ t ) = ˆ k ( t + ∆ t ) (1 + k 2 ∆ t )ˆ u BE u BE k ( t + ∆ t ) = ˆ k ( t ) u BE k ( t + ∆ t ) = (1 + k 2 ∆ t ) − 1 ˆ u BE ˆ k ( t ) Matrix version: u BE ( t + ∆ t ) = ( I − ∆ tL ) − 1 u BE ( t )

Crank-Nicolson: u CN ( t + ∆ t ) = u CN ( t ) + ∆ t f ( u CN ( t )) + f ( u CN ( t + ∆ t )) � � 2 u CN ( t + ∆ t ) − ∆ t 2 f ( u CN ( t + ∆ t )) = u CN ( t ) + ∆ t 2 f ( u CN ( t )) 1 + k 2 ∆ t 1 − k 2 ∆ t � � � � u CN u CN ˆ ( t + ∆ t ) = ˆ ( t ) k k 2 2 ( t + ∆ t ) = 1 − k 2 ∆ t u CN 2 u CN ˆ ˆ ( t ) k k 1 + k 2 ∆ t 2 | Amp. factor | < 1 but spurious large- k behavior: slow oscillatory decay

A-stable methods: du u exact ( t ) = e − k 2 t dt = − k 2 u = ⇒ = ⇒ t →∞ u exact ( t ) = 0 lim A-stable method = ⇒ t →∞ u numerical ( t ) = 0 lim Define q ≡ − k 2 ∆ t and generalize to q complex Define amplification factor Φ( q ) Consider behavior for R e ( q ) < 0 A-stable: | Φ( q ) | < 1 L-stable: | Φ( q ) | → 0 as q → ∞

Crank-Nicolson: seek stability region in complex plane, where � � 1 + q/ 2 � ≡ � Φ CN ( q ) � � � � � ≤ 1 � � 1 − q/ 2 � | 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 R e ( q ) < 0 , have | Φ( q ) | < 1 � YES L-stable ⇐ ⇒ For R e ( 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: u AB ( t + ∆ t ) = u AB ( t ) + ∆ t 3 f ( u AB ( t )) − f ( u AB ( t − ∆ t )) � � 2 � 3 2 u AB ( t ) − 1 � = u AB ( t ) + q 2 u AB ( t − ∆ t ) � u AB ( t + ∆ t ) � 1 + 3 q − q � � � u AB ( t ) � = 2 2 u AB ( t ) u AB ( t − ∆ t ) 1 0 Require that both eigenvalues obey | λ | < 1 Eigenvalues λ obey characteristic polynomial � 1 + 3 q � ( − λ ) + q 0 = 2 − λ 2 1 + 3 q + q � � = λ 2 − λ 2 2 2 = λ 2 − λ q 3 λ − 1

Set λ = e iθ : 2 = λ 2 − λ 3 λ − 1 = e 2 iθ − e iθ q 3 e iθ − 1 = cos 2 θ − cos θ + i (sin 2 θ − sin θ ) (3 cos θ − 1) + 3 i 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: u BD ( t + ∆ t ) = 4 3 u BD ( t ) − 1 3 u BD ( t − ∆ t ) + 2 3∆ tf ( u BD ( t + ∆ t )) = 4 3 u BD ( t ) − 1 3 u BD ( t − ∆ t ) + 2 3 qu BD ( t + ∆ t ) � 1 − 2 � u BD ( t + ∆ t ) = 4 3 u BD ( t ) − 1 3 u BD ( t − ∆ t ) 3 q 1 � 4 3 u BD ( t ) − 1 � u BD ( t + ∆ t ) = 3 u BD ( t − ∆ t ) 1 − 2 q 3 � u BD ( t + ∆ t ) − 1 4 � � � � u BD ( t ) � 3 − 2 q 3 − 2 q = u BD ( t ) u BD ( t − ∆ t ) 1 0 Eigenvalues λ obey characteristic polynomial 4 1 � � 0 = 3 − 2 q − λ ( − λ ) + 3 − 2 q (2 q − 3) = − 4 λ + 1 λ 2

Require that both eigenvalues obey | λ | < 1 so set λ = e iθ (2 q − 3) = − 4 1 e iθ + e 2 iθ = − 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 s α j u n +1 − j = ∆ t � � β j f ( u n +1 − j ) j =0 j =0 Degrees of freedom { α j } , { β j } allow many routes to order- p accuracy. Scale by setting α 0 = 1 . Explicit ⇐ ⇒ β 0 = 0 : s u n +1 = − α j u n +1 − j + ∆ tβ j f ( u n +1 − j ) � � � j =1 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 order 1 order 2 order 3 order 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 e Lt = I + tL + t 2 2 L 2 + . . . = V V − 1 + tV Λ V − 1 + t 2 2 V Λ V − 1 V Λ V − 1 + I + t Λ + t 2 � � 2 Λ 2 + . . . V − 1 = V e Λ t V − 1 = V � λ 1 � � λ 1 � λ 2 � � 0 0 0 Λ 2 = 1 = λ 2 0 λ 2 0 λ 2 0 2 1 + tλ 1 + ( tλ 1 ) 2 � � + . . . 0 e t Λ = 2 1 + tλ 2 + ( tλ 2 ) 2 0 + . . . 2 � e tλ 1 � 0 = e tλ 2 0

Imaginary Eigenvalues � 0 − ω � L = λ ± = ± iω ω 0 � 0 − ω � � 0 − ω � − ω 2 � � 0 L 2 = = − ω 2 ω 0 ω 0 0 � 0 − ω � � − ω 2 ω 3 � � � 0 0 L 3 = = − ω 2 − ω 3 ω 0 0 0 e Lt = I + tL + t 2 2 L 2 + t 3 6 L 3 + . . . � 1 0 � 0 − ω � − ω 2 + t 2 + t 3 � � � � ω 3 � 0 0 = + t − ω 2 − ω 3 0 1 ω 0 0 0 2 6 � 1 − 1 2 ( tω ) 2 + . . . 6 ( tω ) 3 + . . . − tω + 1 � = 6 ( tω ) 3 + . . . 2 ( tω ) 2 + . . . tω − 1 1 − 1 � cos( ωt ) − sin( ωt ) � = sin( ωt ) cos( ωt )

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

Holds for any analytic function f ( L ) 1 � j ! f ( j ) (0) L j f ( L ) = j 1 V Λ V − 1 � j � j ! f ( j ) (0) � = j   1  V − 1 = V f (Λ) V − 1 � j ! f ( j ) (0)Λ j = V j where     λ 1 f ( λ 1 ) λ 2 ... f ( λ 2 )     f (Λ) = f  =     ...        λ N f ( λ N )

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