Numerical Methods for Differential Equations 1.- Numerical Methods - PowerPoint PPT Presentation

Numerical Methods for Differential Equations 1.- Numerical Methods for ODEs Luis M. Abia, J. C. López Marcos, O. Angulo abia@mac.uva.es University of Valladolid Valladolid, (Spain) Euro Summer School Lipari (Sicilia, Italy) 13-26 September 2009– p. 1/64

Contents Numerical Methods for ODEs and DDEs (2 units, L. M. Abia) Runge-Kutta and Multistep Methods. Absolute Stability. Numerical Methods for DDEs. Numerical Methods for Structured Population Models (2 units, O. Angulo) Euro Summer School Lipari (Sicilia, Italy) 13-26 September 2009– p. 2/64

The Problem y ′ ( x ) = f ( x, y ( x )) , a ≤ x ≤ b, y ( a ) = A , y ( x ) = [ y 1 ( x ) , . . . , y d ( x )] T ∈ R d , R d → I f = [ f 1 , . . . , f d ] T ∈ R d . R d , f : [ a, b ] × I SUCH THAT: R d → I R d continuous f : [ a, b ] × I f satisfies a Lipschitz condition with respect to y � f ( x, y 1 ) − f ( x, y 2 ) � ≤ L � y 1 − y 2 � for all x ∈ [ a, b ] , y 1 , y 2 ∈ R d = ⇒ (Picard) Existence and uniqueness of a solution y ( x ) in [ a, b ] and � y 1 ( x ) − y 2 ( x ) � ≤ e L ( x − a ) � y 1 ( a ) − y 2 ( a ) � , x ∈ [ a, b ] for any two solutions y 1 ( x ) , y 2 ( x ) . Euro Summer School Lipari (Sicilia, Italy) 13-26 September 2009– p. 3/64

2 120 1.8 L = 4 , L = 4 L = 4 , L = − 4 100 1.6 1.4 80 1.2 1 60 y y y ′ = 4 y y ′ = − 4 y 0.8 40 0.6 0.4 20 0.2 0 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 x x A more accurate estimate is obtained if f satisfies a one side Lipschitz condition with respect to y < f ( x, y 1 ) − f ( x, y 2 ) , y 1 − y 2 > ≤ L � y 1 − y 2 � 2 x ∈ [ a, b ] , y 1 , y 2 ∈ R d . for all Then � y 1 ( x ) − y 2 ( x ) � ≤ e L ( x − a ) � y 1 ( a ) − y 2 ( a ) � , x ∈ [ a, b ] for any two solutions y 1 ( x ) , y 2 ( x ) . Euro Summer School Lipari (Sicilia, Italy) 13-26 September 2009– p. 4/64

Discrete Variable Methods a = x 0 < x 1 < x 2 < . . . < x N = b . x n +1 = x n + h n +1 , n = 0 , . . . , N − 1 , h = m´ ax 0 ≤ n ≤ N − 1 h n +1 . { y n } N n =0 , Numerical Solution, y n ≈ y ( x n ) 1 y 0 y 1 0.9 y’=−2 x y 2 y(x) 0.8 0.7 y 2 0.6 y 0.5 h 5 0.4 y 3 0.3 y 4 0.2 y 5 h 0 h 1 h 2 h 3 h 4 0.1 x_0 x_1 x_2 x_3 x_4 x_5 x y 0 = A , y n +1 = y n + h n +1 f ( x n , y n ) , n = 0 , . . . , N − 1 , (Euler Method) Euro Summer School Lipari (Sicilia, Italy) 13-26 September 2009– p. 5/64

Discrete Variable Methods Methods can be explicit or implicit : y n +1 = y n + h n +1 f ( x n , y n ) , n = 0 , . . . , N − 1 , y 0 ≈ A (Euler Method) y n +1 = y n + h n +1 f ( x n , y n +1 ) , n = 0 , . . . , N − 1 , y 0 ≈ A (Backward Euler) One-step methods y n +1 = y n + h n +1 Φ ( x n , y n , y n +1 ; f , h n +1 ) , y 0 dado . Multistep methods k � y n +1 = α j y n +1 − j + h Φ ( x n , y n +1 , . . . , y n + k − 1 ; f , h ) , j =1 or k � y n +1 = α n,j y n +1 − j + h n +1 Φ ( x n , y n +1 , . . . , y n + k − 1 ; f , ∆ n ) , j =1 with ∆ n = { x n +1 , . . . , x n − k +1 } , and y 0 , . . . , y k − 1 dados . Others Euro Summer School Lipari (Sicilia, Italy) 13-26 September 2009– p. 6/64

One-Step Methods for ODEs Local solution at x = x n , u ′ = f ( x, u ) , x > x n , u ( x n ) = y n . Local error at the step x n → x n +1 le n +1 := u ( x n +1 ) − y n +1 Example: Taylor Series Method (d=1 for simplicity): u ( x n +1 ) ≈ y n +1 = y n + hu ′ ( x n ) + h 2 2 u ′′ ( x n ) + h 3 6 u ′′′ ( x n ) , where u ′ ( x n ) = f, u ′′ ( x n ) = f x + f y f, u ′′′ ( x n ) = f xx + 2 f xy f + f yy f 2 + f y f x + f 2 y f, with all the functions evaluated at ( x n , y n ) . Euro Summer School Lipari (Sicilia, Italy) 13-26 September 2009– p. 7/64

General Explicit Runge-Kutta Methods A general explicit s -stages Runge-Kutta method is given by the formulae s � = y n + h b i f ( x n + c i h, Y i ) , y n +1 i =1 i − 1 � = y n + h a ij f ( x n + c j h, Y j ) , i = 1 , . . . , s, Y i j =1 The quantities Y i , i = 1 , . . . , s , are called inner stages of the method. The method is completely defined by the parameters c 1 0 c = [ c 1 , . . . , c s ] T , ( abscisae ) c 2 a 21 0 . . . ... b = [ b 1 , . . . , b s ] T , ( weights ) . . . , . . . � s A = ( a ij ) , j =1 a ij = c i , i = 1 , . . . , s c s a s 1 a s 2 · · · 0 b 1 b 2 · · · b s Euro Summer School Lipari (Sicilia, Italy) 13-26 September 2009– p. 8/64

One-step Methods for ODEs Classical low order Runge-Kutta Methods are easily obtained from � x n + h u ( x n + h ) = u ( x n ) + f ( x, u ( x )) dx, x n using quadrature rules. For example, Modified Euler method y n +1 = y n + h f ( x n + 1 2 h, y n + 1 2 h f ( x n , y n )) , Improved Euler Method y n +1 = y n + h 2 ( f ( x n , y n ) + f ( x n + h, y n + h f ( x n , y n )) Euro Summer School Lipari (Sicilia, Italy) 13-26 September 2009– p. 9/64

Other RK formulae The general approach is to match up to a given order the terms in the Taylor serie expansions for u ( x n + h ) and y n +1 ( h ) . For example, for the scalar equation y ′ = f ( x, y ) c 1 0 c 2 a 21 0 and for an explicit 3-stages RK-method , we have c 3 a 31 a 32 0 b 1 b 2 b 3 u ( x n + h ) = y n + hf + 1 2 h 2 ( f x + ff y ) + 1 6 h 3 ( f y ( f x + ff y ) + ( f xx + 2 ff xy + f 2 f yy )) + O ( h 4 ) , and y n +1 = y n + h ( b 1 + b 2 + b 3 ) f + h 2 ( b 2 c 2 + b 3 c 3 )( f x + ff y ) + 1 2 h 3 [2 b 3 c 2 a 32 ( f x + ff y ) f y + ( b 2 c 2 2 + b 3 c 2 3 )( f xx + 2 ff xy + f 2 f yy )] + O ( h 4 ) . where f x , f y , f xx , . . . , are evaluated at ( x n , y n ) . Euro Summer School Lipari (Sicilia, Italy) 13-26 September 2009– p. 10/64

Explicit 2-stage RK methods ( b 3 = c 3 = a 31 = a 32 = 0 , b 1 , b 2 , a 21 = c 2 ): 0 0 0 b 1 + b 2 = 1 , c/ 2 c/ 2 0 c = 2 , Modified Euler , , b 2 c 2 = 1 c = 1 , Improved Euler 2 1 − 1 /c 1 /c Explicit 3-stage RK methods ( b 1 , b 2 , b 3 , c 3 , a 31 , a 32 , a 21 ) 0 0 b 1 + b 2 + b 3 = 1 , b 2 c 2 + b 3 c 3 = 1 1 / 3 1 / 3 1 / 2 1 / 2 2 , , b 2 c 2 2 + b 3 c 2 3 = 1 2 / 3 0 2 / 3 1 − 1 2 3 b 3 c 2 a 32 = 1 1 / 4 0 3 / 4 1 / 6 2 / 3 1 / 6 6 Euro Summer School Lipari (Sicilia, Italy) 13-26 September 2009– p. 11/64

4-Stages Classical RK-methods Order conditions b 1 + b 2 + b 3 + b 4 = 1 , b 2 c 2 + b 3 c 3 + b 4 c 4 = 1 / 2 , b 2 c 2 2 + b 3 c 2 3 + b 4 c 2 4 = 1 / 3 , b 3 a 32 c 2 + b 4 ( a 42 c 2 + a 43 c 3 ) = 1 / 6 , b 2 c 3 2 + b 3 c 3 3 + b 4 c 3 4 = 1 / 4 , b 3 c 3 a 32 c 2 + b 4 c 4 ( a 42 c 2 + a 43 c 3 ) = 1 / 8 , b 3 a 32 c 2 2 + b 4 ( a 42 c 2 2 + a 43 c 2 3 ) = 1 / 2 , b 4 a 43 a 32 c 2 = 1 / 24 0 0 1 / 2 1 / 2 1 / 3 1 / 3 1 / 2 0 1 / 2 2 / 3 − 1 / 3 1 . 1 0 0 1 1 1 − 1 1 1 / 6 2 / 6 2 / 6 1 / 6 1 / 8 3 / 8 3 / 8 1 / 8 Euro Summer School Lipari (Sicilia, Italy) 13-26 September 2009– p. 12/64

The algebraic theory of order for RK methods y ′ = f ( y ) , f = [ f 1 , . . . , f d ] . a ≤ x ≤ b, We introduce the notation (from tensorial calculus) ∂f i ∂ 2 f i ∂y j = f i ∂y j ∂y k = f i ,j , ,j,k , 1 ≤ i, j, k ≤ d. and The convention that when an index is repeated in an expression then we assume that the expression is to be summed over all values of the index. For example, ( u ′ ) i ( x n ) = f i ( y n ) , i = 1 , . . . , d, d � ( u ′′ ) i ( x n ) = f i ,j ( y n ) f j ( y n ) = f i ,j ( y n ) f j ( y n ) , i = 1 , . . . , d. j =1 ,j,k f j f k + f i ( u ′′′ ) i ( x n ) = f i ,j f j ,k f k , ,j,k,l f j f k f l + f i ,l f l f k + f i ( u (4) ) i ( x n ) = f i ,j,k f j ,j,k f j f k ,l f l ,k f k + f i ,k,l f k f l + f i + f i ,j,l f l f j ,j f j ,j f j ,k f k ,l f l , Euro Summer School Lipari (Sicilia, Italy) 13-26 September 2009– p. 13/64

Order Conditions for RK Methods t q t ❏ ❏ ❏ ❏ t n t p t t ❏ ✡ ✡ ❏ ✡ ✡ ❏ ❏ t m t l t t ❏ ✡ ❏ ✡ ❏ ✡ ❏ ✡ t j t k t t ❏ ✡ ❏ ✡ ❏ ✡ ❏ ✡ t i t root root We associate an elementary differential of f to each monotone labelling of a rooted tree, as for example, in d � ( F ) i ( λτ ) = f i ,j,k f k f j ,l,m f m f l ,p,n f n f p ,q f q . j,k,l,m,n,p,q =1 Then � u ( q ) ( x n + h ) = α ( τ ) F ( τ )( y n ) τ ∈ T q Here T q is the set of all the rooted trees with order (number of nodes) q , and α ( τ ) represents the number of different monotone labelling for the rooted tree τ . Euro Summer School Lipari (Sicilia, Italy) 13-26 September 2009– p. 14/64

Order Conditions for RK Methods t q t ❏ ❏ ❏ ❏ t n t p t t ❏ ✡ ✡ ❏ ✡ ✡ ❏ ❏ t m t l t t ❏ ✡ ❏ ✡ ❏ ✡ ❏ ✡ t j t k t t ❏ ✡ ❏ ✡ ❏ ✡ ❏ ✡ t i t root root For each monotone labelling λτ of a rooted tree τ of order q with labels j 1 < j 2 < . . . < j q we introduce the coefficients s � � Ψ i ( τ ) = a ij a ik a jm a jl a ln a lp a pq = c i a ij c j a jl c l a lp c p , j,k,l,m,n,p,q =1 j,l,p and we put Ψ ( τ ) = [Ψ 1 ( λτ ) , . . . , Ψ s ( λτ )] . We denote τ = [ τ 1 , τ 2 , . . . , τ M ] if τ is obtained by connecting the roots of the rooted trees τ 1 , . . . , τ M to a new node that is going to become the root of the rooted tree τ . Then we define the density γ ( τ ) of τ , by recurrence, γ ( τ ) := ρ ( τ ) γ ( τ 1 ) γ ( τ 2 ) . . . γ ( τ M ) , ρ ( τ ) = number of nodes Euro Summer School Lipari (Sicilia, Italy) 13-26 September 2009– p. 15/64