Numerical methods for dynamical systems Alexandre Chapoutot ENSTA - PowerPoint PPT Presentation

Numerical methods for dynamical systems Alexandre Chapoutot ENSTA Paris master CPS IP Paris 2020-2021

Differential equations Many classes Ordinary Differential Equations (ODE) y ( t ) = f ( t , y ( t )) ˙ Differential-Algebraic equations (DAE) y ( t ) = f ( t , y ( t ) , x ( t )) ˙ 0 = g ( t , y ( t ) , x ( t )) Delay Differential Equations (DDE) ˙ y ( t ) = f ( t , y ( t ) , y ( t − τ )) and others: partial differential equations, etc. Remark This talk focuses on ODE 2 / 42

High order vs first order and non-autonomous vs autonomous High order vs first order � � � � ˙ y 1 y 2 ¨ y = f ( y , ˙ y ) ⇔ = with y 1 = y and y 2 = ˙ y . ˙ f ( y 1 , y 2 ) y 2 Non-autonomous vs autonomous � � � � ˙ 1 t y = f ( t , y ) ⇔ ˙ ˙ z = = = g ( z ) . y ˙ f ( t , y ) 3 / 42

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 4 / 42

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 . 5 / 42

One-step methods: Runge-Kutta family 1 One-step methods: Runge-Kutta family Building Runge-Kutta methods 2 3 Variable step-size methods 4 Solving algebraic equations in IRK 5 Implementation in Python 6 Special cases : symplectic integrator 6 / 42

Examples of Runge-Kutta methods Single-step fixed step-size explicit Runge-Kutta method e.g. explicit Trapezoidal method (or Heun’s method) 1 is defined by: 0 k 1 = f ( t ℓ , y ℓ ) , k 2 = f ( t ℓ + 1 h , y ℓ + h 1 k 1 ) 1 1 � 1 2 k 1 + 1 � y i +1 = y ℓ + h 2 k 2 1 1 2 2 y expl. trap. rule Intuition y 1 y = t 2 + y 2 ˙ y 0 = 0 . 46 1 h = 1 . 0 k 2 dotted line is the exact solution. k 1 y 0 1 t 1 2 1 example coming from “Geometric Numerical Integration”, Hairer, Lubich and Wanner, 2006. 7 / 42

Examples of Runge-Kutta methods Single-step variable step-size explicit Runge-Kutta method e.g. Bogacki-Shampine ( ode23 ) is defined by: k 1 = f ( t n , y n ) k 2 = f ( t n + 1 2 h n , y n + 1 0 2 h k 1 ) 1 1 2 2 k 3 = f ( t n + 3 4 h n , y n + 3 3 3 4 h k 2 ) 0 4 4 2 1 4 1 � 2 9 k 1 + 1 3 k 2 + 4 � 9 3 9 y n +1 = y n + h 9 k 3 2 1 4 9 3 9 k 4 = f ( t n + 1 h n , y n +1 ) 7 1 1 1 � 7 24 4 3 8 24 k 1 + 1 4 k 2 + 1 3 k 3 + 1 � z n +1 = y n + h 8 k 4 Remark: the step-size h is adapted following � y n +1 − z n +1 � 1 example coming from “Geometric Numerical Integration”, Hairer, Lubich and Wanner, 2006. 7 / 42

Examples of Runge-Kutta methods Single-step fixed step-size implicit Runge-Kutta method e.g. Runge-Kutta Gauss method (order 4) is defined by: √ √ � � � � � � �� 1 3 1 1 3 k 1 = f t n + 2 − h n , y n + h 4 k 1 + 4 − k 2 (1a) 6 6 √ √ � � � �� � �� 1 3 1 3 k 1 + 1 k 2 = f t n + 2 + h n , y n + h 4 + 4 k 2 (1b) 6 6 � 1 2 k 1 + 1 � y n +1 = y n + h 2 k 2 (1c) Remark: A non-linear system of equations must be solved at each step. 1 example coming from “Geometric Numerical Integration”, Hairer, Lubich and Wanner, 2006. 7 / 42

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 · · · b 1 b 2 b s i b ′ b ′ b ′ · · · (optional) 1 2 s Which induces the following recurrence formula: � � s s � � k i = f t n + c i h n , y n + h a ij k j y n +1 = y n + h b i k i (2) j =1 i =1 Explicit method (ERK) if a ij = 0 is i ≤ j Diagonal Implicit method (DIRK) if a ij = 0 is i ≤ j and at least one a ii � = 0 Singly Diagonal implicit method (SDIRK) if a ij = 0 is i ≤ j and all a ii = γ are identical. Implicit method (IRK) otherwise 8 / 42

Building Runge-Kutta methods 1 One-step methods: Runge-Kutta family Building Runge-Kutta methods 2 3 Variable step-size methods 4 Solving algebraic equations in IRK 5 Implementation in Python 6 Special cases : symplectic integrator 9 / 42

Building RK methods: Order condition Every numerical method member of the Runge-Kutta family follows the condition order . Order condition This condition states that a method of this family is of order p if and only if the p + 1 first coefficients of the Taylor expansion of the true solution and the Taylor expansion of the numerical methods are equal. In other terms, a RK method has order p if y ( t n ; y n − 1 ) − y n = h p +1 Ψ f ( y n ) + O ( h p +2 ) 10 / 42

Building RK methods: Order condition Taylor expansion of the exact and the numerical solutions At a time instant t n the Taylor expansion of the true solution with the Lagrange remainder states that there exists ξ ∈ ] t n , t n +1 [ such that: p h i � i ! y ( i ) ( t n ; y 0 ) + O ( h p +1 ) n y ( t n +1 ; y 0 ) = y ( t n ; y 0 ) + i =1 p h i � i ! f ( i − 1) ( t n , y ( t n ; y 0 )) + O ( h p +1 ) n = y ( t n ; y 0 ) + i =1 The Taylor expansion (very complex expression) of the numerical solution is given by expanding, around ( t n , y n ), the expression: s � y n +1 = y n + h b i k i i =1 Consequence of the condition order The definition of RK methods (Butcher table coefficients) is based on the solution of under-determined system of algebraic equations. 11 / 42

Example: 3-stages explicit RK method (scalar IVP) One considers a scalar ODE ˙ y = f ( t , y ) with f : R × R → R One tries to determine the coefficients b i ( i = 1 , 2 , 3), c 2 , c 3 , a 32 such that y n +1 = y n + h ( b 1 k 1 + b 2 k 2 + b 3 k 3 ) k 1 = f ( t n , y n ) k 2 = f ( t n + c 2 h , y n + hc 2 k 1 ) k 3 = f ( t n + c 3 h , y n + h ( c 3 − a 32 ) k 1 + ha 32 k 2 Some notations (evaluation at point ( t n , y ( t n )): f tt = ∂ 2 f ( t , y ) f t = ∂ f ( t , y ) f ty = ∂ f ( t , y ) f = f ( t , y ) · · · ∂ t 2 ∂ t ∂ t ∂ y Note: in Butcher tableau we always have the row-sum condition s � c i = a ij , i = 1 , 2 , . . . , s . j =1 12 / 42

Example: 3-stages explicit RK method (scalar IVP) Taylor expansion of y ( t n +1 ), the exact solution, around t n : y ( t n +1 ) = y ( t n ) + hy (1) ( t n ) + h 2 2 y (2) ( t n ) + h 3 6 y (3) ( t n ) + O ( h 4 ) Moreover, y (1) ( t n ) = f y (2) ( t n ) = f t + f y ˙ y = f t + ff y y (3) ( t n ) = f tt + f ty f + f ( f ty + f yy f ) + f y ( f y + ff y ) = f tt + 2 ff ty + f 2 f yy + f y ( f t + ff y ) With F = f t + ff y and G = f tt + 2 ff ty + f 2 f ty , one has: y ( t n +1 ) = y ( t n ) + hf + h 2 2 F + h 3 6 ( Ff y + G ) + O ( h 4 ) 13 / 42

Example: 3-stages explicit RK method (scalar IVP) Taylor expansion k i around t n k 2 = f + hc 2 ( f t + k 1 f y ) + h 2 2 c 2 f tt + 2 k 1 f ty + k 2 + O ( h 3 ) � � 1 f yy 2 = f + hc 2 F + h 2 2 c 2 2 G + O ( h 3 ) k 3 = f + h { c 3 f t + [( c 3 − a 32 ) k 1 + a 32 k 2 ] f y } + h 2 � c 2 3 f tt + 2 c 3 [( c 3 − a 32 ) k 1 + a 32 k 2 ] f ty 2 + [( c 3 − a 32 ) k 1 + a 32 k 2 ] 2 f yy � + O ( h 3 ) = f + hc 3 F + h 2 ( c 2 a 32 Ff y + 1 2 c 2 3 G + O ( h 3 ) (substituting k 1 = f and k 2 ) Taylor expansion of y n +1 (localizing assumption y n = y ( t n )) y n +1 = y ( t n ) + h ( b 1 + b 2 + b 3 ) f + h 2 ( b 2 c 2 + b 3 c 3 ) F + h 3 � 2 b 3 c 2 a 32 Ff y + ( b 2 c 2 2 + b 3 c 2 3 ) G � + O ( h 4 ) 2 14 / 42

Example: 3-stages explicit RK method (scalar IVP) Building one stage method We fix b 2 = b 3 = 0, so one gets y n +1 = y ( t n ) + hb 1 f + O ( h 2 ) In consequence b 1 = 1 (by identification) so one gets Euler’s method (order 1) 15 / 42

Example: 3-stages explicit RK method (scalar IVP) Building two stages method We fix b 3 = 0, so one gets y n +1 = y ( t n ) + h ( b 1 + b 2 ) f + h 2 b 2 c 2 F + 1 2 h 3 b 2 c 2 2 G + O ( h 3 ) In consequence to get order 2 methods, we need to solve b 2 c 2 = 1 b 1 + b 2 = 1 2 Remark: there is a (singly) infinite number of solutions. Two particular solutions of order 2: 0 0 1 1 1 1 2 1 1 0 1 2 2 16 / 42

Example: 3-stages explicit RK method (scalar IVP) Building three stages method In consequence to get order 3 methods, we need to solve b 2 c 2 + b 3 c 3 = 1 b 1 + b 2 + b 3 = 1 2 3 = 1 b 3 c 2 a 32 = 1 b 2 c 2 2 + b 3 c 2 3 6 Remark: there is a (doubly) infinite number of solutions. Two particular solutions of order 3: 0 0 1 1 1 1 3 3 2 2 2 2 0 1 − 1 2 3 3 1 3 1 2 1 0 4 4 6 3 6 17 / 42