Ordinary and partial differential equations Victor Eijkhout 335/394 - PowerPoint PPT Presentation

Ordinary and partial differential equations Victor Eijkhout 335/394 fall 2011

ODEs and PDEs Time-evolving phenomena: IVP (Initial Value Problem), usually Ordinary Differential Equations Space-constraint phenomena: BVP (Boundary Value Problem), usually Partial Differential Equations 2

Ordinary Differential Equations 3

Simple Initial Value Problem Newtons’ Second Law: F = ma Acceleration depends linearly on the force exerted on the mass: dtu = d 2 a = d dt 2 x = F m Write velocity as vector → first order derivative in time: u ( t ) = ( x ( t ) , x ′ ( t )) t u implicitly depends on time, easier to solve: u ′ ( t ) = f ( u ( t )) , u (0) = u 0 4

Numerical treatment of differential equations Initial value problem: u ′ ( t ) = f ( t , u ( t )) , u (0) = u 0 , t > 0 Boundary value problem: u ′′ ( x ) = f ( x ) , x (0) = x 0 , x (1) = x 1 , x ∈ [0 , 1] General assumption: f has higher derivatives. IVP stability: solutions corresponding to different u 0 values converge as t → ∞ . Criterium:  > 0  unstable  ∂ ∂ u f ( t , u ) = = 0 neutrallystable   < 0 stable Simple example: f ( t , u ) = − λ u , then u ( t ) = u 0 e − λ t ; stable if λ > 0 5

Finite difference approximation We turn the continuous problem into a discrete one, by looking at finite time/space steps. Assume all functions are sufficiently smooth, and use Taylor series: u ( t + ∆ t ) = u ( t ) + u ′ ( t )∆ t + u ′′ ( t )∆ t 2 2! + u ′′′ ( t )∆ t 3 3! + · · · This gives for u ′ : u ′ ( t ) = u ( t + ∆ t ) − u ( t ) + O (∆ t 2 ) ∆ t So we approximate u ′ ( t ) ≈ u ( t + ∆ t ) − u ( t ) ∆ t and the “truncation error” is O (∆ t 2 ). 6

Finite differences 2 How does this help? In u ′ = f ( t , u ) substitute u ′ ( t ) → u ( t + ∆ t ) − u ( t ) ∆ t giving u ( t + ∆ t ) − u ( t ) = f ( t , u ( t )) ∆ t or u ( t + ∆ t ) = u ( t ) + ∆ t f ( t , u ( t )) Let t 0 = 0, t k +1 = t k + ∆ t = · · · = ( k + 1)∆ t , u ( t k ) = u k : u k +1 = u k + ∆ t f ( t k , u k ) Discretization ‘Explicit Euler’ or ‘Euler forward’. Does this compute something close to the true solution? ‘Discretization error’ 7

Some error analysis Local Truncation Error: assume computed solution is exact at step k , how wrong will we be at step k + 1? u ( t k ) + u ′ ( t k )∆ t + u ′′ ( t k )∆ t 2 u ( t k +1 ) = 2! + · · · u ( t k ) + f ( t k , u ( t k ))∆ t + u ′′ ( t k )∆ t 2 = 2! + · · · = u k + f ( t k , u k )∆ t u k +1 So u k +1 − u ( t k +1 ) L k +1 = u k − u ( t k ) + f ( t k , u k ) − f ( t k , u ( t k )) − u ′′ ( t k )∆ t 2 = 2! + · · · − u ′′ ( t k )∆ t 2 = 2! + · · · Global error is first order method: � � � � ∆ t 2 / 2! ∆ t 2 / 2! E k ≈ Σ k L k ≈ k ≈ ( T − t 0 / ∆ t ) ≈ O (∆ t ) 8

An Euler forward example Consider f ( t , u ) = − λ u , exact solution u ( t ) = u 0 e − λ t ; stable if λ > 0 Explicit Euler scheme u k +1 = u k − ∆ t λ u k = (1 − λ ∆ t ) u k = (1 − λ ∆ t ) k u 0 Then u k → 0 as k → ∞ ⇔ | 1 − λ ∆ t | < 1 ⇔ − 1 < 1 − λ ∆ t < 1 ⇔ − 2 < − λ ∆ t < 0 ⇔ 0 < λ ∆ t < 2 ⇔ ∆ t < 2 /λ Conditionally stable 9

Implicit Euler Or ‘Euler backward’: u ( t − ∆ t ) = u ( t ) − u ′ ( t )∆ t + u ′′ ( t )∆ t 2 2! + · · · so u ′ ( t ) = u ( t ) − u ( t − ∆ t ) + u ′′ ( t )∆ t / 2 + · · · ∆ t Compute u ′ ( t ) = f ( t , u ( t )) as u ( t ) − u ( t − ∆ t ) = f ( t , u ( t )) ∆ t ⇒ u ( t ) = u ( t − ∆ t ) + ∆ tf ( t , u ( t )) ⇒ u k +1 = u k + ∆ tf ( t k +1 , u k +1 ) Implicit equation for u k +1 ! Let f ( t , u ) = − u 3 , then u k +1 = u k − ∆ tu 3 k +1 needs nonlinear solver 10

Stability of Implicit Euler Again the f ( t , u ) = − λ u example: u k +1 = u k − λ ∆ tu k +1 (1 + λ ∆ t ) u k +1 = u k � � � � k 1 1 u k +1 = u k = u 0 1+ λ ∆ t 1+ λ ∆ t If λ > 0 (stable equation), then u k → 0 for all values of λ and ∆ t : unconditionally stable. Pro: larger time steps possible, no worries Con: implicit equation needs to be solved 11

Higher order methods Runge-Kutta Adams-Bashforth Crank-Nicholson 12

Boundary value problems Consider u ′′ ( x ) = f ( x , u , u ′ ) for x ∈ [ a , b ] where u ( a ) = u a , u ( b ) = u b in 1D and x ) for x ∈ Ω = [0 , 1] 2 with u (¯ − u xx (¯ x ) − u yy (¯ x ) = f (¯ x ) = u 0 on δ Ω. (1) in 2D. 13

Approximation of 2nd order derivatives Taylor series (write h for δ x ): u ( x + h ) = u ( x ) + u ′ ( x ) h + u ′′ ( x ) h 2 2! + u ′′′ ( x ) h 3 3! + u (4) ( x ) h 4 4! + u (5) ( x ) h 5 5! + · · · and u ( x − h ) = u ( x ) − u ′ ( x ) h + u ′′ ( x ) h 2 2! − u ′′′ ( x ) h 3 3! + u (4) ( x ) h 4 4! − u (5) ( x ) h 5 5! + · · · Subtract: u ( x + h ) + u ( x − h ) = 2 u ( x ) + u ′′ ( x ) h 2 + u (4) ( x ) h 4 12 + · · · so − u (4) ( x ) h 4 u ′′ ( x ) = u ( x + h ) − 2 u ( x ) + u ( x − h ) 12 + · · · h 2 Numerical scheme: − u ( x + h ) − 2 u ( x ) + u ( x − h ) = f ( x , u ( x ) , u ′ ( x )) h 2 (2nd order PDEs are very common!) 14

This leads to linear algebra − u ( x + h ) − 2 u ( x ) + u ( x − h ) = f ( x , u ( x ) , u ′ ( x )) h 2 Equally spaced points on [0 , 1]: x k = kh where h = 1 / n , then − u k +1 + 2 u k − u k − 1 = h 2 f ( x k , u k , u ′ for k = 1 , . . . , n − 1 k ) Written as matrix equation:       h 2 f 1 + u 0 2 − 1 u 1       h 2 f 2 − 1 2 − 1 u 2  =      . . ... ... ... . . . . 15

Matrix properties • Very sparse, banded • Symmetric (only because 2nd order problem) • Sign pattern: positive diagonal, nonpositive off-diagonal (true for many second order methods) • Positive definite (just like the continuous problem) 16

Initial Boundary value problem Heat conduction in a rod T ( x , t ) for x ∈ [ a , b ], t > 0: ∂ t T ( x , t ) − α ∂ 2 ∂ ∂ x 2 T ( x , t ) = q ( x , t ) • Initial condition: T ( x , 0) = T 0 ( x ) • Boundary conditions: T ( a , t ) = T a ( t ), T ( b , t ) = T b ( t ) • Material property: α > 0 is thermal diffusivity • Forcing function: q ( x , t ) is externally applied heating. The equation u ′′ ( x ) = f above is the steady state. 17

Discretization Space discretization: x 0 = a , x n = b , x j +1 = x j + ∆ x Time discretiation: t 0 = 0, t k +1 = t k + ∆ t Let T k j approximate T ( x j , t k ) Space: ∂ t T ( x j , t ) − α T ( x j − 1 , t ) − 2 T ( x j , t ) + T ( x j +1 , t ) ∂ = q ( x j , t ) ∆ x 2 Explicit time stepping: T k +1 − T k T k j − 1 − 2 T k j + T k j j j +1 = q k − α j ∆ x 2 ∆ t Implicit time stepping: T k +1 T k +1 j − 1 − 2 T k +1 + T k +1 − T k j j +1 j j = q k +1 − α ∆ x 2 j ∆ t 18

Computational form: explicit j + α ∆ t T k +1 = T k ∆ x 2 ( T k j − 1 − 2 T k j + T k j +1 ) + ∆ tq k j j This has an explicit form: � � I + α ∆ t k +1 = k + ∆ tq k T T ∆ x 2 � � � 19

Computational form: implicit − α ∆ t T k +1 ∆ x 2 ( T k j − 1 − 2 T k j + T k j +1 ) = T k j + ∆ tq k j j This has an implicit form: � � I − α ∆ t k +1 = T k + ∆ tq k ∆ x 2 K T � � � Needs to solve a linear system in every time step 20

Stability of explicit scheme Let q ≡ 0; assume T k j = β k e i ℓ x j ; for stability we require | β | < 1: j + α ∆ t T k +1 T k ∆ x 2 ( T k j − 1 − 2 T k j + T k = j +1 ) j β k e i ℓ x j + α ∆ t ∆ x 2 ( β k e i ℓ x j − 1 − 2 β k e i ℓ x j + β k e i ℓ x j +1 ) ⇒ β k +1 e i ℓ x j = 1 + 2 α ∆ t ∆ x 2 [1 2( e i ℓ ∆ x + e − ℓ ∆ x ) − 1] ⇒ β = 1 + 2 α ∆ t ∆ x 2 (cos( ℓ ∆ x ) − 1) = 21

β k +1 = 1 + 2 α ∆ t ∆ x 2 (cos( ℓ ∆ x ) − 1) β k To get | β | < 1: • 2 α ∆ t ∆ x 2 (cos( ℓ ∆ x ) − 1) < 0: automatic • 2 α ∆ t ∆ x 2 (cos( ℓ ∆ x ) − 1) > − 2: needs 2 α ∆ t ∆ x 2 < 1, that is ∆ t < ∆ x 2 2 α big restriction on size of time steps 22

Stability of implicit scheme − α ∆ t T k +1 ∆ x 2 ( T k +1 − 2 T k +1 + T k +1 T k j +1 ) = j j 1 j j ⇒ β k +1 e i ℓ ∆ x − α ∆ t ∆ x 2 ( β k +1 e i ℓ x j − 1 − 2 β k +1 e i ℓ x j + β k +1 e i ℓ x j +1 ) β k e i ℓ x j = 1 + 2 α ∆ t ⇒ β − 1 = ∆ x 2 (1 − cos( ℓ ∆ x )) 1 β = 1 + 2 α ∆ t ∆ x 2 (1 − cos( ℓ ∆ x )) Noting that 1 − cos( ℓ ∆ x ) > 0, the condition | β | < 1 always satisfied: method always stable. 23

Sparse matrix in 2D case Sparse matrices so far were tridiagonal: only in 1D case. Two-dimensional: − u xx − u yy = f on unit square [0 , 1] 2 Difference equation: 4 u ( x , y ) − u ( x + h , y ) − u ( x − h , y ) − u ( x , y + h ) − u ( x , y − h ) = h 2 f ( x , y ) 24