Scientific Computing I Module 4: Numerical Methods for ODE Miriam - PowerPoint PPT Presentation

Lehrstuhl Informatik V Scientific Computing I Module 4: Numerical Methods for ODE Miriam Mehl based on Slides by Michael Bader (Winter 09/10) Winter 2011/2012 Miriam Mehl based on Slides by Michael Bader (Winter 09/10): Scientific Computing I Module 4: Numerical Methods for ODE, Winter 2011/2012 1

Lehrstuhl Informatik V Part I: Basic Numerical Methods Motivation: Direction Fields Euler’s Method Discretized Model vs. Discrete Model Implicit Euler Analysis of Numerical Schemes for ODE Local Discretisation Error Global Discretisation Error Order of Consistency/Convergence Miriam Mehl based on Slides by Michael Bader (Winter 09/10): Scientific Computing I Module 4: Numerical Methods for ODE, Winter 2011/2012 2

Lehrstuhl Informatik V Motivation: Direction Fields • given: initial value problem: ˙ y ( t ) = f ( t , y ( t )) , y ( t 0 ) = y 0 • easily computable: direction field 5 4 p(t) 3 2 1 0 0 2 4 6 8 10 t • idea: “follow the arrows” Miriam Mehl based on Slides by Michael Bader (Winter 09/10): Scientific Computing I Module 4: Numerical Methods for ODE, Winter 2011/2012 3

Lehrstuhl Informatik V “Following the Arrows” • direction field illustrates slope for given time t n and value y n : ˙ y n = f ( t n , y n ) • “follow arrows” = make a small step in the given direction: y n + 1 := y n + τ ˙ y n = y n + τ f ( t n , y n ) • motivates numerical scheme: y 0 := y 0 y n + 1 := y n + τ f ( t n , y n ) for n = 0 , 1 , 2 , . . . Miriam Mehl based on Slides by Michael Bader (Winter 09/10): Scientific Computing I Module 4: Numerical Methods for ODE, Winter 2011/2012 4

Lehrstuhl Informatik V Euler’s Method • numerical scheme is called Euler’s method : y n + 1 := y n + τ f ( t n , y n ) • results from finite difference approximation: y n + 1 − y n ≈ ˙ y n = f ( t n , y n ) τ (difference quotient instead of derivative) • or from truncation of Taylor expansion: y ( t n ) + O ( τ 2 ) y ( t n + 1 ) = y ( t n ) + τ ˙ Miriam Mehl based on Slides by Michael Bader (Winter 09/10): Scientific Computing I Module 4: Numerical Methods for ODE, Winter 2011/2012 5

Lehrstuhl Informatik V Euler’s Method and Direction Fields 60 50 p(t) 40 30 20 10 0 0 50 100 150 200 250 t use direction at the beginning of the timestep Miriam Mehl based on Slides by Michael Bader (Winter 09/10): Scientific Computing I Module 4: Numerical Methods for ODE, Winter 2011/2012 6

Lehrstuhl Informatik V Euler’s Method – 1D examples • model of Malthus, ˙ p ( t ) = α p ( t ) : p n + 1 := p n + τα p n • Logistic Growth, ˙ p ( t ) = α ( 1 − p ( t ) /β ) p ( t ) : � � 1 − p n p n + 1 := p n + τα p n β • Logistic growth with threshold: � � � 1 − p n 1 − p n � p n + 1 := p n + τα p n β δ Miriam Mehl based on Slides by Michael Bader (Winter 09/10): Scientific Computing I Module 4: Numerical Methods for ODE, Winter 2011/2012 7

Lehrstuhl Informatik V Euler’s Method in 2D • Euler’s method is easily extended to systems of ODE → use vector notation: y n + 1 := y n + τ f ( t n , y n ) • example: nonlinear extinction model � 71 � 8 − 23 12 p ( t ) − 25 ˙ p ( t ) = 12 q ( t ) p ( t ) � 73 � ˙ 8 − 25 12 p ( t ) − 23 q ( t ) = 12 q ( t ) q ( t ) • Euler’s method: � 71 � 8 − 23 12 p n − 25 p n + 1 = p n + τ 12 q n p n � 73 � 8 − 25 12 p n − 23 q n + 1 = q n + τ 12 q n q n Miriam Mehl based on Slides by Michael Bader (Winter 09/10): Scientific Computing I Module 4: Numerical Methods for ODE, Winter 2011/2012 8

Lehrstuhl Informatik V Discretized Model vs. Discrete Model • simplest example: model of Malthus p n + 1 := p n − τα p n , α > 0 • compare to discrete model: p n + 1 := p n − δ p n , δ > 0 with decay rate δ (“percentage”) • obvious restriction in the discrete model: δ < 1 • obvious restriction for τ in the discretized model? τα < 1 ⇒ τ < α − 1 • not that simple in non-linear models or systems of ODE! Miriam Mehl based on Slides by Michael Bader (Winter 09/10): Scientific Computing I Module 4: Numerical Methods for ODE, Winter 2011/2012 9

Lehrstuhl Informatik V Implicit Euler • Euler’s method (“explicit Euler”): y n + 1 := y n + τ f ( t n , y n ) • implicit Euler: y n + 1 := y n + τ f ( t n + 1 , y n + 1 ) • results from finite difference approximation: y n + 1 − y n ≈ ˙ y n + 1 = f ( t n + 1 , y n + 1 ) τ • explicit formula for y n + 1 not immediately available • to do: solve equation for y n + 1 Miriam Mehl based on Slides by Michael Bader (Winter 09/10): Scientific Computing I Module 4: Numerical Methods for ODE, Winter 2011/2012 10

Lehrstuhl Informatik V Implicit Euler and Direction Fields 60 50 p(t) 40 30 20 10 0 0 50 100 150 200 250 t use direction at end of the timestep Miriam Mehl based on Slides by Michael Bader (Winter 09/10): Scientific Computing I Module 4: Numerical Methods for ODE, Winter 2011/2012 11

Lehrstuhl Informatik V Implicit Euler – Examples • example: Model of Malthus 1 p n + 1 := p n − τα p n + 1 ⇒ p n + 1 = 1 + τα p n • correct (discrete) model? 0 < ( 1 − τα ) − 1 < 1 for any τ α > 0 : then τ < α − 1 required! α < 0 : then • in physics α > 0 is more frequent! (damped systems, friction, . . . ) • implicit schemes preferred when explicit schemes require very small τ Miriam Mehl based on Slides by Michael Bader (Winter 09/10): Scientific Computing I Module 4: Numerical Methods for ODE, Winter 2011/2012 12

Lehrstuhl Informatik V Implicit Euler – 2D Example • example: arms race = p n + τ ( b 1 + a 11 p n + 1 + a 12 q n + 1 ) p n + 1 q n + 1 = q n + τ ( b 2 + a 21 p n + 1 + a 22 q n + 1 ) • solve linear system of equations: ( 1 − τ a 11 ) p n + 1 − τ a 12 q n + 1 = p n + τ b 1 − τ a 21 p n + 1 + ( 1 − τ a 22 ) q n + 1 = q n + τ b 2 (for each time step n → n + 1) • in vector notation: ( I − τ A ) y n + 1 = y n + τ b Miriam Mehl based on Slides by Michael Bader (Winter 09/10): Scientific Computing I Module 4: Numerical Methods for ODE, Winter 2011/2012 13

Lehrstuhl Informatik V Local Discretisation Error • local influence of using differences instead of derivatives • example: Euler’s method �� � � y ( t + τ ) − y ( t ) = τ � � y n || + O ( τ 2 ) 2 || ¨ l ( τ ) = max − f ( t , y ( t )) � � τ � � [ a , b ] • memory hook: insert exact solution y ( t ) into y n + 1 − y n − ˙ y n τ A numerical scheme is called consistent , if l ( τ ) → 0 for τ → 0 Miriam Mehl based on Slides by Michael Bader (Winter 09/10): Scientific Computing I Module 4: Numerical Methods for ODE, Winter 2011/2012 14

Lehrstuhl Informatik V Global Discretisation Error • compare numerical solution with exact solution • example: Euler’s method e ( τ ) = k = 0 ,..., k max {� y k − y ( t k ) �} max y ( t ) exact solution; y k solution of the discretized equations (depends on τ ) A numerical scheme is called convergent , if e ( τ ) → 0 for τ → 0 Miriam Mehl based on Slides by Michael Bader (Winter 09/10): Scientific Computing I Module 4: Numerical Methods for ODE, Winter 2011/2012 15

Lehrstuhl Informatik V Global Discretisation Error Euler e n = y n − y ( t n ) , y n + 1 = y n + τ f ( t n , y n ) , y ( t n ) + τ f ( t n , y ( t n )) + τ 2 ¨ y ( t n ) + O ( τ 3 ) , y ( t n + 1 ) = 2 ˙ | e n | + τ M | e n | + N τ 2 ⇒ | e n + 1 | ≤ ( 1 + M τ ) | e n | + N τ 2 = ( 1 + τ M ) 2 | e n − 1 | + ( 1 + τ M ) N τ 2 + N τ 2 ≤ + N τ 2 e n τ M − 1 . . . ≤ e n τ M | e 0 | ≤ τ M ���� = 0 N τ e n τ M − 1 = = O ( τ ) . M Miriam Mehl based on Slides by Michael Bader (Winter 09/10): Scientific Computing I Module 4: Numerical Methods for ODE, Winter 2011/2012 16

Lehrstuhl Informatik V Order of Consistency/Convergence A numerical scheme is called consistent of order p ( p -th order consistent), if l ( τ ) = O ( τ p ) A numerical scheme is called convergent of order p ( p -th order convergent), if e ( τ ) = O ( τ p ) We have shown that the explicit Euler method is consistent and convergent of order 1. Miriam Mehl based on Slides by Michael Bader (Winter 09/10): Scientific Computing I Module 4: Numerical Methods for ODE, Winter 2011/2012 17

Lehrstuhl Informatik V Part II: Advanced Numerical Methods Runge-Kutta-Methods 2nd-order Runge-Kutta 4th-order Runge-Kutta Multistep Methods Adams-Bashforth Adams-Moulton Problems for Numerical Methods for ODE Ill-Conditioned Problems Stability Stiff Equations Summary Miriam Mehl based on Slides by Michael Bader (Winter 09/10): Scientific Computing I Module 4: Numerical Methods for ODE, Winter 2011/2012 18

Lehrstuhl Informatik V Runge-Kutta-Methods • 1st idea: use additional evaluations of f : p � β i f ( t i n , y i n ) with t i y n + 1 = y n + τ n ∈ [ t n ; t n + 1 ] i = 1 open questions: Where to obtain y i n , i = 1 , . . . , p ? How to choose β i , i = 1 , . . . , p ? • 2nd idea: numerical approximations for missing values of y : p � α i , j f ( t j n , y j y i n := y n + τ n ) j = 1 explicit Runge-Kutta: α i , j = 0 if j ≥ i . • 3rd idea: choose β i and α i , j such that order of consistency is maximal (use qudrature rules) Miriam Mehl based on Slides by Michael Bader (Winter 09/10): Scientific Computing I Module 4: Numerical Methods for ODE, Winter 2011/2012 19