Validated simulation of ODEs Julien Alexandre dit Sandretto Department U2IS ENSTA Paris SSC310-2020
Contents Initial Value Problem of Ordinary Differential Equations Problem of integral computation Picard-Lindel¨ of Integration scheme Problem of wrapping effect Affine arithmetic Stepsize controller Validated integration in a contractor formalism Additive constraints Temporal constraints Bibliography Do it yourself Julien Alexandre dit Sandretto - Validated simulation of ODEs October 22, 2020- 2
Initial Value Problem of Ordinary Differential Equations Initial Value Problem of Ordinary Differential Equations Classical problem Consider an IVP for ODE, over the time interval [ 0 , T ] y = f ( y ) ˙ with y ( 0 ) = y 0 This IVP has a unique solution y ( t ; y 0 ) if f : R n → R n is Lipschitz. Interval IVP y = f ( y , p ) ˙ with y ( 0 ) ∈ [ y 0 ] and p ∈ [ p ] Julien Alexandre dit Sandretto - Validated simulation of ODEs October 22, 2020- 3
Initial Value Problem of Ordinary Differential Equations Numerical Integration � t How compute y ( t ) = y 0 + 0 f ( y ( s )) ds ? Goal of numerical integration ◮ Compute a sequence of time instants: t 0 = 0 < t 1 < · · · < t n = T ◮ Compute a sequence of values: y 0 , y 1 , . . . , y n such that ∀ i ∈ [ 0 , n ] , y i ≈ y ( t i ; y 0 ) . Goal of validated numerical integration ◮ Compute a sequence of time instants: t 0 = 0 < t 1 < · · · < t n = T ◮ Compute a sequence of values: [ y 0 ] , [ y 1 ] , . . . , [ y n ] such that ∀ i ∈ [ 0 , n ] , [ y i ] ∋ y ( t i ; y 0 ) . Julien Alexandre dit Sandretto - Validated simulation of ODEs October 22, 2020- 4
Problem of integral computation Problem of integral computation � h Discrete system given by y n + 1 = y n + 0 f ( y ( s )) ds � h Bounding of 0 f ( y ( s )) ds If y ( s ) is bounded s.t. y ( s ) ∈ [ x ] , ∀ s ∈ [ 0 , h ] , then � h f ([ x ]) ds ⊂ [ 0 , h ] · [ f ]([ x ]) 0 How bound y ( s ) ? Complex, it is what we are trying to compute ! We note by [˜ y n ] ⊃ { y ( s ) , s ∈ [ t n , t n + 1 ] } Julien Alexandre dit Sandretto - Validated simulation of ODEs October 22, 2020- 5
Picard-Lindel¨ of Picard-Lindel¨ of (or Cauchy-Lipschitz) Theorem (Banach fixed-point theorem) Let ( K , d ) a complete metric space and let g : K → K a contraction that is for all x, y in K there exists c ∈ ] 0 , 1 [ such that d ( g ( x ) , g ( y )) � c · d ( x , y ) , then g has a unique fixed-point in K. ———— We consider the space of continuously differentiable functions C 0 ([ t j , t j + 1 ] , R n ) and the Picard-Lindel¨ of operator � t p f ( y ) = t �→ y j + f ( y ( s )) ds , with y j = y ( t j ) (1) t j If this operator is a contraction then its solution is unique and its solution is the solution of IVP. Julien Alexandre dit Sandretto - Validated simulation of ODEs October 22, 2020- 6
Picard-Lindel¨ of Interval counterpart of Picard-Lindel¨ of With a first order integration scheme that is for f : R n → R n a continuous function and [ a ] ⊂ IR n , we have � a f ( s ) ds ∈ ( a − a ) f ([ a ]) = w ([ a ]) f ([ a ]) , (2) a we can define a simple enclosure function of Picard-Lindel¨ of such that [ p f ]([ r ]) def = [ y j ] + [ 0 , h ] · f ([ r ]) , (3) with h = t j + 1 − t j the step-size. In consequence, if one can find [ r ] such that [ p f ]([ r ]) ⊆ [ r ] then [˜ y j ] ⊆ [ r ] by the Banach fixed-point theorem. Julien Alexandre dit Sandretto - Validated simulation of ODEs October 22, 2020- 7
Picard-Lindel¨ of Interval counterpart of Picard-Lindel¨ of ~ [y j ] We can then build the Lohner 2-steps method: [y j ] 1. Find [˜ y j ] and h j with Picard-Lindel¨ of [y j+1 ] operator and Banach’s theorem 2. Compute [ y j + 1 ] with a validated integration t scheme: Taylor or Runge-Kutta t j t j+1 h j It is important to obtain [˜ y j ] and [ y j + 1 ] as tight as possible Integration scheme at order higher than one: Taylor for example Julien Alexandre dit Sandretto - Validated simulation of ODEs October 22, 2020- 8
Integration scheme Integration scheme Two main approaches: ◮ Taylor series (Vnode, CAPD, etc.): 1 h i f [ i ] ( y j ) + O ( h p + 1 ) with f [ i ] the i th term of y j + 1 = y j + � p serie expansion of f . O ( h p + 1 ) can be easily bounded by the Lagrange remainder of serie s.t. O ( h p + 1 ) = f [ p + 1 ] ( ξ ) , with ξ ∈ [˜ y j ] , and then O ( h p + 1 ) ∈ f [ p + 1 ] (˜ y j ) ◮ Runge-Kutta methods (DynIBEX): y j + 1 = Φ( y j , f , p ) + LTE , with Φ any RK method and LTE the local truncation error. Julien Alexandre dit Sandretto - Validated simulation of ODEs October 22, 2020- 9
Integration scheme Runge-Kutta methods s -stage Runge-Kutta methods are described by a Butcher tableau · · · c 1 a 11 a 12 a 1 s j . . . . . . . . . . . . c s a s 1 a s 2 · · · a ss i b 1 b 2 · · · b s Which induces the following recurrence: � s � s � � k i = f t j + c i h j , y j + h a il k l y j + 1 = y j + h b i k i l = 1 i = 1 ◮ Explicit method (ERK) if a il = 0 is i � l ◮ Diagonal Implicit method (DIRK) if a il = 0 is i � l and at least one a ii � = 0 ◮ Implicit method (IRK) otherwise Julien Alexandre dit Sandretto - Validated simulation of ODEs October 22, 2020- 10
Integration scheme Explicit methods Interval extensions 1. Computation of k 1 = f ( y j ) , k 2 = f ( y j + h · a 21 · k 1 ) , . . . , k i = f ( y j + h � i − 1 ℓ = 1 a i ℓ k ℓ ) , . . . , k s = f ( y j + h � s − 1 ℓ = 1 a s ℓ k ℓ ) 2. Computation of y j + 1 = y j + h � s i = 1 b i k i + LTE ⇒ with interval arithmetic (natural extension) Example of HEUN 0 0 0 1 1 0 1 / 2 1 / 2 [ k 1 ] = [ f ]([ y j ]) , [ k 2 ] = [ f ]([ y j ] + h [ k 1 ]) , [ y j + 1 ] = [ y j ] + h ([ k 1 ] + [ k 2 ]) / 2 Julien Alexandre dit Sandretto - Validated simulation of ODEs October 22, 2020- 11
Integration scheme Implicit Schemes Example of Radau IIA − 1 / 12 1 / 3 5 / 12 1 3 / 4 1 / 4 3 / 4 1 / 4 [ k 1 ] = [ f ]([ y j ] + h ( 5 [ k 1 ] / 12 − [ k 2 ] / 12 )) , [ k 2 ] = [ f ]([ y j ] + h ( 3 [ k 1 ] / 4 + [ k 2 ] / 4 ) We need to solve this system of implicit equations ! Julien Alexandre dit Sandretto - Validated simulation of ODEs October 22, 2020- 12
Integration scheme Solve implicit scheme with a contractor point of view k 1 is the approximate of f ( y ( t j + h / 3 )) , but by construction y ( t j + h / 3 ) ∈ [˜ y j ] , then [ k 1 ] ⊂ f ([˜ y j ]) (same for k 2 ) Algorithm based on contraction Require: f , [˜ y j ] , [ y j ] , LTE [ k 1 ] = [ f ]([˜ y j ]) and [ k 2 ] = [ f ]([˜ y j ]) while [ k 1 ] or [ k 2 ] improved do [ k 1 ] = [ k 1 ] ∩ [ f ]([ y j ] + h ( 5 [ k 1 ] / 12 − [ k 2 ] / 12 )) [ k 2 ] = [ k 2 ] ∩ [ f ]([ y j ] + h ( 3 [ k 1 ] / 4 + [ k 2 ] / 4 ) end while [ y j + 1 ] = [ y j ] + h ( 3 [ k 1 ] + [ k 2 ]) / 4 + LTE return [ y j + 1 ] Julien Alexandre dit Sandretto - Validated simulation of ODEs October 22, 2020- 13
Integration scheme How to compute the LTE ? h p + 1 � y ( t n ; y n − 1 ) − y n = C · � C ∈ R . with Order condition This condition states that a method of Runge-Kutta family is of order p iff ◮ the Taylor expansion of the exact solution ◮ and the Taylor expansion of the numerical methods have the same p + 1 first coefficients. Consequence The LTE is the difference of Lagrange remainders of two Taylor expansions Julien Alexandre dit Sandretto - Validated simulation of ODEs October 22, 2020- 14
Integration scheme A quick view of Runge-Kutta order condition theory Starting from y ( q ) = ( f ( y )) ( q − 1 ) and with the Chain rule, we have High order derivatives of exact solution y y = f ( y ) ˙ ¨ y = f ′ ( y )˙ f ′ ( y ) is a linear map y y ( 3 ) = f ′′ ( y )(˙ y ) + f ′ ( y )¨ f ′′ ( y ) is a bi-linear map y , ˙ y y ( 4 ) = f ′′′ ( y )(˙ y ) + f ′ ( y ) y ( 3 ) y ) + 3 f ′′ ( y )(¨ f ′′′ ( y ) is a tri-linear map y , ˙ y , ˙ y , ˙ . y ( 5 ) = f ( 4 ) ( y )(˙ . y , ˙ y , ˙ y , ˙ y ) + 6 f ′′′ ( y )(¨ y , ˙ y , ˙ y ) . + 4 f ′′ ( y )( y ( 3 ) , ˙ y ) + f ′ ( y ) y ( 4 ) y ) + 3 f ′′ ( y )(¨ y , ¨ . . . Julien Alexandre dit Sandretto - Validated simulation of ODEs October 22, 2020- 15
Integration scheme A quick view of Runge-Kutta order condition theory Inserting the value of ˙ y , ¨ y , . . . , we have: High order derivatives of exact solution y y = f ˙ ¨ y = f ′ ( f ) y ( 3 ) = f ′′ ( f , f ) + f ′ ( f ′ ( f )) y ( 4 ) = f ′′′ ( f , f , f ) + 3 f ′′ ( f ′ f , f ) + f ′ ( f ′′ ( f , f )) + f ′ ( f ′ ( f ′ ( f ))) . . . ◮ Elementary differentials , such as f ′′ ( f , f ) , are denoted by F ( τ ) Remark a tree structure is made apparent in these computations Julien Alexandre dit Sandretto - Validated simulation of ODEs October 22, 2020- 16
Integration scheme A quick view of Runge-Kutta order condition theory Rooted trees ◮ f is a leaf ◮ f ′ is a tree with one branch, . . . , f ( k ) is a tree with k branches Example f is associated to f f ′ f ′′ ( f ′ f , f ) f ′′ Remark: this tree is not unique e.g., symmetry Julien Alexandre dit Sandretto - Validated simulation of ODEs October 22, 2020- 17
Recommend
More recommend
