Numerical methods for dynamical systems Alexandre Chapoutot ENSTA Paris master CPS IP Paris 2020-2021
Part I Multi-step methods for IVP-ODE 2 / 43
I nitial V alue P roblem of O rdinary D ifferential E quations Consider an IVP for ODE, over the time interval [0 , t end ] y = f ( t , y ) ˙ with y (0) = y 0 IVP has a unique solution y ( t ; y 0 ) if f : R n → R n is Lipschitz in y ∀ t , ∀ y 1 , y 2 ∈ R n , ∃ L > 0 , � f ( t , y 1 ) − f ( t , y 2 ) �≤ L � y 1 − y 2 � . Goal of numerical integration Compute a sequence of time instants: t 0 = 0 < t 1 < · · · < t n = t end Compute a sequence of values: y 0 , y 1 , . . . , y n such that ∀ ℓ ∈ [0 , n ] , y ℓ ≈ y ( t ℓ ; y 0 ) . s.t. y ℓ +1 ≈ y ( t ℓ + h ; y ℓ ) with an error O ( h p +1 ) where h is the integration step-size p is the order of the method 3 / 43
Simulation algorithm Data: f the flow, y 0 initial condition, t 0 starting time, t end end time, h integration step-size t ← t 0 ; y ← y 0 ; while t < t end do Print( t , y ); y ← Euler( f , t , y , h ); t ← t + h ; end with, the Euler’s method defined by y n +1 = y n + hf ( t n , y n ) and t n +1 = t n + h . 4 / 43
Multi-step methods Recall: single-step methods solve IVP using one value y n and some values of f . A multi-step method approximate solution y n +1 of IVP using k previous values of the solution y n , y n − 1 , . . . , y n − k − 1 . Different methods implement this approach Adams-Bashforth method (explicit) Adams-Moulton method (implicit) Backward difference method (implicit) The general form of such method is k k � � α j y n + j = h β j f ( t n + j , y n + j ) . j =0 j =0 with α j and β j some constants and α k = 1 and | α 0 | + | β 0 | � = 0 5 / 43
Polynomial interpolation 1 Polynomial interpolation 2 Multi-step methods: Adams family Building Adams-Bashforth’s methods Building Adams-Moulton’s method Predictor-Corrector methods Implementation in Python 3 Multi-step methods: BDF 4 Order condition 5 Variable order and variable step-size multi-step methods 6 / 43
A quick remainder on polynomial interpolation Starting point: a function f ( t ) a sequence of n time instants t 1 , t 2 , . . . , t n . a sequence of points f 1 = f ( t 1 ), f 2 = f ( t 2 ), . . . , f n = f ( t n ) Goal Find a polynomial p of order n approximating f and passes through the ( n + 1) function values p ( t i ) = f i Theorem (Uniqueness of the Interpolating Polynomial) Given n unequal points x 1 , x 2 , . . . , x n and arbitrary values f 1 , f 2 , . . . , f n there is at most one polynomial p of degree less or equal to n − 1 such that p ( x i ) = f i , i = 1 , . . . , n . Note: different algorithms in function of the monomial basis 7 / 43
Polynomial interpolation: basis Standard basis We consider p ( x ) = a 0 + a 1 x + a 2 x 2 + a 3 x 3 + · · · + a n x n we have to find a i such that p ( x i ) = f ( x i ) so the Vandermond matrix x 2 x n 1 x 0 . . . a 0 f ( x 0 ) 0 0 x 2 x n 1 x 1 . . . a 1 f ( x 1 ) 1 1 = . . . . . . . . . . . . . . . . . . x 2 x n 1 x n . . . a n f ( x n ) n n Lagrange basis We consider p ( x ) = f ( x 0 ) ℓ 0 ( x ) + f ( x 1 ) ℓ 1 ( x ) + · · · + f ( x n ) ℓ n ( x ) such that n x − x i � ℓ i ( x ) = x i − x j j =0 , j � = i 8 / 43
Polynomial interpolation: error Interpolation error If f is n + 1 continuously differentiable on [ a , b ] then E n ( x ) = ( x − x 0 )( x − x 1 ) . . . ( x − x n ) f ( n +1) ( ξ ) ( n + 1)! with ξ ∈ ] a , b [ Comments : Vandermond matrix is not use as it is ill-conditioned Lagrange interpolation is useful when f change but not x i Remark For our purpose to define multi-step methods, equidistant time instants will be considered! 9 / 43
Multi-step methods: Adams family 1 Polynomial interpolation 2 Multi-step methods: Adams family Building Adams-Bashforth’s methods Building Adams-Moulton’s method Predictor-Corrector methods Implementation in Python 3 Multi-step methods: BDF 4 Order condition 5 Variable order and variable step-size multi-step methods 10 / 43
Adams-Bashforth method – 1 Integral form of IVP y = f ( t , y ) ˙ y ( t 0 ) = y 0 ⇔ � t � t n +1 y ( t ) = y 0 + f ( s , y ( s )) ds ⇒ y n +1 = y n + f ( s , y ( s )) ds t 0 t n Ingredients: We denote by t i = t n + ih the grid of points in time We assume given numerical approximations: y n , y n − 1 , . . . , y n − k +1 of the exact solution. we can use y i , i = n − k + 1 , . . . , n , to approximate f ( t , y ( t )) using f ( t i , y i ) ≡ f i . We can use polynomial interpolation with points: { ( t i , f i ) : i = n − k + 1 , . . . , n } to approximate integral. 11 / 43
Adams-Bashforth method – 2 We have n + 1 distinct (equidistant) points ( t 0 , f 0 ) , ( t 1 , f 1 ) , · · · , ( t n , f n ) with f i = f ( t i , y i ) Adams-Bashforth method is defined by � t n +1 � t n +1 n n � � y n +1 = y n + f i ℓ i ( s ) ds = y n + f i ℓ i ( s ) ds t n t n i =0 i =0 Example of first Adams-Bashforth methods of order k : k = 1: y n +1 = y n + hf n (explicit Euler method) k = 2: y n +1 = y n + h � 3 2 f n − 1 � 2 f n − 1 k = 3: y n +1 = y n + h � 23 12 f n − 16 5 � 12 f n − 1 + 12 f n − 2 k = 4: y n +1 = y n + h � 55 24 f n − 59 24 f n − 1 + 37 9 � 24 f n − 2 − 24 f n − 3 12 / 43
Adams-Bashforth’s method – 3 from sympy import ∗ t = Symbol ( ’ t ’ , r e a l=True , p o s i t i v e=True ) h = Symbol ( ’ h ’ , r e a l=True , p o s i t i v e=True ) tn = Symbol ( ’ t n ’ , r e a l=True , p o s i t i v e=True ) tnm3 = tn − 3 ∗ h tnm2 = tn − 2 ∗ h tnm1 = tn − h tnp1 = tn + h fnm3 = Symbol ( ’ f { n − 3 } ’ , r e a l=True ) fnm2 = Symbol ( ’ f { n − 2 } ’ , r e a l=True ) fnm1 = Symbol ( ’ f { n − 1 } ’ , r e a l=True ) fn = Symbol ( ’ f n ’ , r e a l=True ) yn = Symbol ( ’ y n ’ , r e a l=True ) ynp1 = Symbol ( ’ y { n+1 } ’ , r e a l=True ) p o i n t s o r d e r 1 = [ ( tn , fn ) ] p o i n t s o r d e r 2 = [ ( tnm1 , fnm1 ) , ( tn , fn ) ] p o i n t s o r d e r 3 = [ ( tnm2 , fnm2 ) , ( tnm1 , fnm1 ) , ( tn , fn ) ] p o i n t s o r d e r 4 = [ ( tnm3 , fnm3 ) , ( tnm2 , fnm2 ) , ( tnm1 , fnm1 ) , ( tn , fn ) ] 13 / 43
Adams-Bashforth’s method – 4 def l a g r a n g e b a s i s ( time , p o i n t s ) : acc = 1 f o r p o i n t i n p o i n t s : i f ( time != p o i n t [ 0 ] ) : acc = acc ∗ ( t − p o i n t [ 0 ] ) / ( time − p o i n t [ 0 ] ) e l s e : acc = p o i n t [ 1 ] ∗ acc return acc def l a g r a n g e ( p o i n t s ) : acc = 0 f o r p o i n t i n p o i n t s : acc = acc + l a g r a n g e b a s i s ( p o i n t [ 0 ] , p o i n t s ) return acc def build adams ( p o i n t s ) : p l = l a g r a n g e ( p o i n t s ) return s i m p l i f y ( i n t e g r a t e ( pl , ( t , tn , tnp1 ) ) ) p r i n t ( ”# # Order 1” ) formula1 = build adams ( p o i n t s o r d e r 1 ) p r i n t ( l a t e x (Eq( ynp1 , yn + formula1 ) ) ) 14 / 43
Explicit Adams-Bashforth formulae This is an explicit ODE solver Each integration step involves only one evaluation of f Using past values of f for order n we use n − 1 past values Adams-Bashforth algorithm of order n can only be used after n − 1 previous steps ( not self starting method ) 15 / 43
Adams-Moulton method – 1 We have n + 2 distinct (equidistant) points ( t 0 , f 0 ) , ( t 1 , f 1 ) , · · · , ( t n , f n ) , ( t n +1 , f n +1 ) with f i = f ( t i , y i ) Adams-Moulton method is defined by � t n +1 � t n +1 n + 1 n + 1 � � y n +1 = y n + f i ℓ i ( s ) ds = y n + f i ℓ i ( s ) ds t n t n i =0 i =0 Example of first Adams-Moulton methods of order k : k = 1: y n +1 = y n + h f n +1 (implicit Euler method) k = 2: y n +1 = h 2 ( f n + f n +1 ) + y n h k = 3: y n +1 = 12 (8 f n + 5 f n +1 − f n − 1 ) + y n h k = 4: y n +1 = 24 (19 f n + 9 f n +1 − 5 f n − 1 + f n − 2 ) + y n 16 / 43
Adams-Moulton method – 3 from sympy import ∗ t = Symbol ( ’ t ’ , r e a l=True , p o s i t i v e=True ) h = Symbol ( ’ h ’ , r e a l=True , p o s i t i v e=True ) tn = Symbol ( ’ t n ’ , r e a l=True , p o s i t i v e=True ) tnm2 = tn − 2 ∗ h tnm1 = tn − h tnp1 = tn + h fnm2 = Symbol ( ’ f { n − 2 } ’ , r e a l=True ) fnm1 = Symbol ( ’ f { n − 1 } ’ , r e a l=True ) fn = Symbol ( ’ f n ’ , r e a l=True ) fnp1 = Symbol ( ’ f { n+1 } ’ , r e a l=True ) yn = Symbol ( ’ y n ’ , r e a l=True ) ynp1 = Symbol ( ’ y { n+1 } ’ , r e a l=True ) p o i n t s o r d e r 1 = [ ( tnp1 , fnp1 ) ] p o i n t s o r d e r 2 = [ ( tn , fn ) , ( tnp1 , fnp1 ) ] p o i n t s o r d e r 3 = [ ( tnm1 , fnm1 ) , ( tn , fn ) , ( tnp1 , fnp1 ) ] p o i n t s o r d e r 4 = [ ( tnm2 , fnm2 ) , ( tnm1 , fnm1 ) , ( tn , fn ) , ( tnp1 , fnp1 ) ] 17 / 43
Recommend
More recommend
