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

Part V Stability analysis 2 / 26

Part 5. Section 1 Introduction to stability of numerical methods 1 Introduction to stability of numerical methods 2 Linear stability analysis for one-step methods 3 Linear stability analysis for multi-step methods 4 Stiffness 3 / 26

Stability properties: a graphical view Note: there are several kinds of stability. Impose C p Consequence x ( t ) Problem Impose C m Consequence x n Method From a generic point of view we have: Impose a certain conditions C p on IVP which force the exact solution x ( t ) to exhibit a certain stability Apply a numerical method on IVP Question: what conditions must be imposed on the method such that the approximate solution ( x n ) n ∈ N has the same stability property? 4 / 26

Total stability of IVP Consider, a perturbed IVP y = f ( t , y ) + δ ( t ) ˙ with y (0) = y 0 + δ 0 and t ∈ [0 , b ] ( δ ( t ) , δ 0 ) denotes the perturbations Definition: totally stable IVP From ( δ ( t ) , δ 0 ) and ( δ ∗ ( t ) , δ ∗ 0 ) two perturbations y ( t ) and y ∗ ( t ) the associated solutions if ∀ t ∈ [0 , b ] , ∀ ε > 0 , ∃ K > 0 , � δ ( t ) − δ ∗ ( t ) �≤ ε ∧ � δ 0 − δ ∗ 0 �≤ ε ⇒� y ( t ) − y ∗ ( t ) �≤ K ε then IVP is totally stable . 5 / 26

Zero stability of numerical methods We consider the application of numerical method on a perturbed IVP so we have a perturbed numerical scheme Definition: zero-stability From δ n and δ ∗ n two discrete-time perturbation y n and y ∗ n the associated numerical solution if ∀ n ∈ [0 , N ] , ∀ ε > 0 , ∃ K > 0 , ∀ h ∈ (0 , h 0 ] � δ n − δ ∗ n �≤ ε ⇒� y n − y ∗ n �≤ K ε then the method is zero-stable In a different point of view, we want to solve ˙ y = 0 with y (0) = y 0 and so numerical method should produce as a solution y ( t ) = y 0 . (It is obvious for RK methods) 6 / 26

Zero stability for multi-step methods First and second characteristic polynomials for linear multi-step methods are k k � � α i z i β i z i ρ ( z ) = and σ ( z ) = i =0 i =0 Root condition A linear multi-step method satisfies the root condition is the roots of the first characteristic polynomial ρ have modulus less than or equal to one and those of modulus one are simple. Theorem A multi-step method is zero stable is it satisfies the root condition. Theorem No zero-stable linear k -step method can have order exceeding k + 1 7 / 26

Consistency of numerical methods We denote by Φ f ( t n , y n ; h ) a Runge-Kutta method such that y n +1 = y n + h Φ f ( t n , y n ; h ) If Φ f is such that h → 0 Φ f ( t n , y n ; h ) = f ( t n , y n ) . lim then the Runge-Kutta method is consistent to the IVP. As a consequence, the truncation error is such that: h → 0 y ( t n +1 ) − y n − h Φ f ( t n , y n ; h ) = 0 lim Consistency for s -stage RK methods A necessary and sufficient condition is that s � b i = 1 i =1 8 / 26

Convergence of numerical methods A Runge-Kutta method is said convergent if h → 0 y n = y ( t n ) lim Problem Lipschitz Problem totally stable condition y n con- consistency and Method verges to y ( t ) zero-stability 9 / 26

Part 5. Section 2 Linear stability analysis for one-step methods 1 Introduction to stability of numerical methods 2 Linear stability analysis for one-step methods 3 Linear stability analysis for multi-step methods 4 Stiffness 10 / 26

Linear stability We consider the IVP: y = λ y ˙ with λ ∈ C , ℜ ( λ ) < 0 Applying a RK method, we get y n +1 = R (ˆ ˆ h ) y n with h = λ h R (ˆ h ) is called the stability function of the method. Stability function of RK methods h ) = det ( I − ˆ hA + ˆ l b t ) h 1 R (ˆ det ( I − ˆ hA ) So, lim n →∞ x n = 0 when | R (ˆ h ) | < 1 11 / 26

Linear stability of ERK – 1 The stability function for s -stage ( s = 1 , 2 , 3 , 4 ⇒ p = s ) ERK is reduced to a polynomial function: h + 1 h 2 + · · · + 1 R (ˆ h ) = 1 + ˆ ˆ ˆ h s 2! s ! p=1 p=2 3 3 2 2 1 1 0 0 -1 -1 -2 -2 -3 -3 -3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 3 p=3 p=4 3 3 2 2 1 1 0 0 -1 -1 -2 -2 -3 -3 -3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 3 12 / 26

Linear stability of ERK – 1 The stability function for s -stage ( s > 4 ⇒ p < s ) ERK is reduced to a polynomial function: s h + 1 h 2 + · · · + 1 h p + R (ˆ h ) = 1 + ˆ ˆ ˆ � γ q ˆ h q 2! p ! q = p +1 with γ q depending only on the coefficients of the ERK methods. For example, for RKF45 ( s = 5 and p = 4) h + 1 h 2 + 1 h 3 + 1 1 h 4 + R (ˆ h ) = 1 + ˆ ˆ ˆ ˆ ˆ h 5 2! 6 24 104 DOPIR54 ( s = 6 and p = 5) h + 1 h 2 + 1 h 3 + 1 1 1 h 4 + h 5 + R (ˆ h ) = 1 + ˆ ˆ ˆ ˆ ˆ ˆ h 6 2! 6 24 120 600 13 / 26

Part 5. Section 3 Linear stability analysis for multi-step methods 1 Introduction to stability of numerical methods 2 Linear stability analysis for one-step methods 3 Linear stability analysis for multi-step methods 4 Stiffness 14 / 26

Linear stability of Adams-Bashworth methods We consider the scalar linear IVP y = λ y ˙ with λ ∈ C , ℜ ( λ ) < 0 For linear problem, the stability polynomial of a multi-step method is π ( r , ˆ h ) = ρ ( r ) − ˆ ˆ h σ ( r ) with h = λ h Stability Domains of AB 1 1 0.8 2 0.6 0.4 4 3 I m { λ · h } 0.2 5 6 0 −0.2 −0.4 −0.6 −0.8 −1 −2.5 −2 −1.5 −1 −0.5 0 0.5 R e { λ · h } 15 / 26

Linear stability of Adams-Moulton methods We consider the scalar linear IVP y = λ y ˙ with λ ∈ C , ℜ ( λ ) < 0 For linear problem, the stability polynomial of a multi-step method is π ( r , ˆ h ) = ρ ( r ) − ˆ ˆ h σ ( r ) with h = λ h Stability Domains of AM 4 2 3 3 2 4 1 1 I m { λ · h } 5 6 0 −1 −2 −3 −4 −7 −6 −5 −4 −3 −2 −1 0 1 2 3 R e { λ · h } 16 / 26

Linear stability of Adams-Bashworth-Moulton methods We consider the IVP: x = λ x ˙ with λ ∈ C , ℜ ( λ ) < 0 Stability Domains of ABM 1.5 1 4 5 3 0.5 6 I m { λ · h } 0 −0.5 −1 −1.5 −2.5 −2 −1.5 −1 −0.5 0 0.5 1 R e { λ · h } 17 / 26

Linear stability of BDF We consider the scalar linear IVP y = λ y ˙ with λ ∈ C , ℜ ( λ ) < 0 For linear problem, the stability polynomial of a multi-step method is π ( r , ˆ h ) = ρ ( r ) − ˆ ˆ h σ ( r ) with h = λ h Stability Domains of BDF 10 5 5 4 I m { λ · h } 3 2 1 0 −5 −10 −5 0 5 10 15 20 R e { λ · h } 18 / 26

Part 5. Section 4 Stiffness 1 Introduction to stability of numerical methods 2 Linear stability analysis for one-step methods 3 Linear stability analysis for multi-step methods 4 Stiffness 19 / 26

Stiff versus non-stiff problems Problem 1 � � � � � � � � y 1 ˙ − 2 1 y 1 2 sin( t ) = + y 2 ˙ 1 − 2 y 2 2(cos( t ) − sin ( t )) Problem 2 � � � � � � � � y 1 ˙ − 2 1 y 1 2 sin( t ) = + y 2 ˙ 998 − 999 y 2 999(cos( t ) − sin ( t )) Both have the same exact solution: � � � � � � � � � � y 1 ( t ) 1 sin( t ) y 1 (0) 2 = 2 exp( − t ) + with initial values = y 2 ( t ) 1 cos( t ) y 2 (0) 3 20 / 26

Simulation results 21 / 26

Stiff linear ODE: a definition We consider linear constant coefficients IVP of the form: y = A y + φ ( t ) ˙ assuming that all eigenvalues λ are such that ℜ ( λ ) < 0 We denote by | ℜ ( λ ) | = max 1 ≤ i ≤ n | ℜ ( λ i ) | | ℜ ( λ ) | = min 1 ≤ i ≤ n | ℜ ( λ i ) | the stiffness ratio is defined by | ℜ ( λ ) | / | ℜ ( λ ) | Stiffness definition - 1 (Lambert) A linear constant coefficients system is stiff iff all eigenvalues are such that ℜ ( λ ) < 0 and the stiffness ratio is large. 22 / 26

Others stiffness definitions Definition 2 (Lambert) Stiffness occurs when stability requirements, rather than those of accuracy, constrain the step size. Definition 3 (Lambert) Stiffness occurs when some components of the solution decay much more quickly than others. Global definition (Lambert) If a numerical method with a finite region of absolute stability, applied to a system with any initial values, if forced to use in a certain interval of integration a step size which is excessively small in relation to the smoothness of the exact solution in that interval, then the system is said to be stiff in that interval. 23 / 26

Linear stability definition for stiff systems - 1 A -stability A method is A-stable if R s ⊇ { ˆ h : ℜ (ˆ h ) < 0 } A ( α )-stability A method is A ( α ) -stable , α ∈ ]0 , π/ 2[, if R s ⊇ { ˆ h : − α < π − arg(ˆ h ) < α } Stiffly stability A method is stiffly stable if R S ⊇ R 1 ∪ R 2 such that R 1 = { ˆ h : ℜ (ˆ h ) < − a } and R 2 = { ˆ h : − a ≤ ℜ (ˆ h ) ≤ 0 , − c ≤ ℑ (ˆ h ) ≤ c } with a and c two positive real numbers. 24 / 26

Linear stability definition for stiff systems - 2 L -stability A one step method is L -stable if it is A -stable and when applied to stable scalar test equations ˙ y = λ y it yields y n +1 = ℜ ( h λ ) x n where | ℜ ( h λ ) |→ 0 as ℜ ( h λ ) → −∞ Relation between the stability definitions L -stability ⇒ A -stability ⇒ stiffly stability ⇒ A ( α )-stability 25 / 26