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

▶
scientific computing i
SMART_READER_LITE
LIVE PREVIEW

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

Jan 06, 2023 334 likes •687 views

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

slide-1
SLIDE 1

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

slide-2
SLIDE 2

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

slide-3
SLIDE 3

Lehrstuhl Informatik V

Motivation: Direction Fields

  • given: initial value problem:

˙ y(t) = f(t, y(t)), y(t0) = y0

  • easily computable: direction field

p(t) t 5 10 4 3 8 2 1 6 4 2

  • 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

slide-4
SLIDE 4

Lehrstuhl Informatik V

“Following the Arrows”

  • direction field illustrates slope for given time tn and value yn:

˙ yn = f(tn, yn)

  • “follow arrows” = make a small step in the given direction:

yn+1 := yn + τ ˙ yn = yn + τ f(tn, yn)

  • motivates numerical scheme:

y0 := y0 yn+1 := yn + τ f(tn, yn) 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

slide-5
SLIDE 5

Lehrstuhl Informatik V

Euler’s Method

  • numerical scheme is called Euler’s method:

yn+1 := yn + τ f(tn, yn)

  • results from finite difference approximation:

yn+1 − yn τ ≈ ˙ yn = f(tn, yn) (difference quotient instead of derivative)

  • or from truncation of Taylor expansion:

y(tn+1) = y(tn) + τ ˙ y(tn) + O(τ 2)

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

slide-6
SLIDE 6

Lehrstuhl Informatik V

Euler’s Method and Direction Fields

p(t) 60 50 40 t 30 20 250 10 200 150 100 50

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

slide-7
SLIDE 7

Lehrstuhl Informatik V

Euler’s Method – 1D examples

  • model of Malthus, ˙

p(t) = αp(t): pn+1 := pn + τα pn

  • Logistic Growth, ˙

p(t) = α (1 − p(t)/β) p(t): pn+1 := pn + τα

  • 1 − pn

β

  • pn
  • Logistic growth with threshold:

pn+1 := pn + τα

  • 1 − pn

β 1 − pn δ

  • pn

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

slide-8
SLIDE 8

Lehrstuhl Informatik V

Euler’s Method in 2D

  • Euler’s method is easily extended to systems of ODE → use

vector notation: yn+1 := yn + τ f(tn, yn)

  • example: nonlinear extinction model

˙ p(t) = 71

8 − 23 12p(t) − 25 12q(t)

  • p(t)

˙ q(t) = 73

8 − 25 12p(t) − 23 12q(t)

  • q(t)
  • Euler’s method:

pn+1 = pn + τ 71

8 − 23 12pn − 25 12qn

  • pn

qn+1 = qn + τ 73

8 − 25 12pn − 23 12qn

  • qn

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

slide-9
SLIDE 9

Lehrstuhl Informatik V

Discretized Model vs. Discrete Model

  • simplest example: model of Malthus

pn+1 := pn − τα pn, α > 0

  • compare to discrete model:

pn+1 := pn − δpn, δ > 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

slide-10
SLIDE 10

Lehrstuhl Informatik V

Implicit Euler

  • Euler’s method (“explicit Euler”):

yn+1 := yn + τ f(tn, yn)

  • implicit Euler:

yn+1 := yn + τ f(tn+1, yn+1)

  • results from finite difference approximation:

yn+1 − yn τ ≈ ˙ yn+1 = f(tn+1, yn+1)

  • explicit formula for yn+1 not immediately available
  • to do: solve equation for yn+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

slide-11
SLIDE 11

Lehrstuhl Informatik V

Implicit Euler and Direction Fields

p(t) 60 50 40 t 30 20 250 10 200 150 100 50

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

slide-12
SLIDE 12

Lehrstuhl Informatik V

Implicit Euler – Examples

  • example: Model of Malthus

pn+1 := pn − τα pn+1 ⇒ pn+1 = 1 1 + τα pn

  • correct (discrete) model?

α > 0 : then 0 < (1 − τα)−1 < 1 for any τ α < 0 : then τ < α−1 required!

  • 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

slide-13
SLIDE 13

Lehrstuhl Informatik V

Implicit Euler – 2D Example

  • example: arms race

pn+1 = pn + τ (b1 + a11pn+1 + a12qn+1) qn+1 = qn + τ (b2 + a21pn+1 + a22qn+1)

  • solve linear system of equations:

(1 − τa11)pn+1 − τa12qn+1 = pn + τb1 −τa21pn+1 + (1 − τa22)qn+1 = qn + τb2 (for each time step n → n + 1)

  • in vector notation: (I − τA)yn+1 = yn + τ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

slide-14
SLIDE 14

Lehrstuhl Informatik V

Local Discretisation Error

  • local influence of using differences instead of derivatives
  • example: Euler’s method

l(τ) = max

[a,b]

  • y(t + τ) − y(t)

τ − f(t, y(t))

  • = τ

2||¨ yn|| + O(τ 2)

  • memory hook: insert exact solution y(t) into

yn+1 − yn τ − ˙ yn 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

slide-15
SLIDE 15

Lehrstuhl Informatik V

Global Discretisation Error

  • compare numerical solution with exact solution
  • example: Euler’s method

e(τ) = max

k=0,...,kmax {yk − y(tk)}

y(t) exact solution; yk 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

slide-16
SLIDE 16

Lehrstuhl Informatik V

Global Discretisation Error Euler

en = yn − y(tn), yn+1 = yn + τf(tn, yn), y(tn+1) = y(tn) + τf(tn, y(tn)) + τ 2 ¨ y(tn) 2 + O(τ 3), ⇒ |en+1| ˙ ≤ |en| + τM|en| + Nτ 2 = (1 + Mτ)|en| + Nτ 2 ≤ (1 + τM)2|en−1| + (1 + τM)Nτ 2 + Nτ 2 ≤ . . . ≤ enτM |e0|

  • =0

+Nτ 2 enτM − 1 τM = Nτ enτM − 1 M = O(τ).

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

slide-17
SLIDE 17

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

slide-18
SLIDE 18

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

slide-19
SLIDE 19

Lehrstuhl Informatik V

Runge-Kutta-Methods

  • 1st idea: use additional evaluations of f:

yn+1 = yn + τ

p

  • i=1

βif(ti

n, yi n) with ti n ∈ [tn; tn+1]

  • pen questions: Where to obtain yi

n, i = 1, . . . , p? How to choose

βi, i = 1, . . . , p?

  • 2nd idea: numerical approximations for missing values of y:

yi

n := yn + τ p

  • j=1

αi,jf(tj

n, yj n)

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

slide-20
SLIDE 20

Lehrstuhl Informatik V

Runge-Kutta-Methods of 2nd Order

  • example: 2nd-order Runge-Kutta:

y1

n

= yn, y2

n

= yn + τf(tn, y1

n ),

yn+1 = yn + τ 1 2(f(tn, y1

n) + f(tn+1, y2 n)).

(“method of Heun”)

  • further example: modified Euler (also 2nd order)

y1

n

= yn, y2

n

= yn + τ 2f(tn, y1

n),

yn+1 = yn + τ f

  • tn + τ

2, y2

n

  • Miriam Mehl based on Slides by Michael Bader (Winter 09/10): Scientific Computing I

Module 4: Numerical Methods for ODE, Winter 2011/2012 20

slide-21
SLIDE 21

Lehrstuhl Informatik V

Runge-Kutta-Method of 4th order

classical 4th-order Runge-Kutta:

  • intermediate steps:

t1

n

= tn, y1

n = yn,

t2

n

= tn + τ 2, y2

n = yn + τ

2f(t1

n, y1 n ),

t3

n

= tn + τ 2, y3

n = yn + τ

2f(t2

n, y2 n ),

t4

n

= tn+1, y4

n = yn + τf(t3 n, y3 n ).

  • explicit scheme:

yn+1 = yn + τ 6

  • f(t1

n, y1 n ) + 2f(t2 n, y2 n ) + 2f(t3 n, y3 n ) + f(t4 n, y4 n )

  • How would you translate this method into efficient pseudo-code?

Remember that unneccessary evaluations of f can be very costly!!!

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

slide-22
SLIDE 22

Lehrstuhl Informatik V

Pseudo-Code Runge-Kutta-Method of 4th order

f[1] = f(t,y[n]); dt2 = dt/2; t2 = t + dt2; f[2] = f(t2, y[n] + dt2*f[1]); f[3] = f(t2, y[n] + dt2*f[2]); t = t + dt; f[4] = f(t, y[n] + dt*f[3]); y[n+1] = y[n] + dt * (f[1] + 2*f[2] + 2*f[3] + f[4]);

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

slide-23
SLIDE 23

Lehrstuhl Informatik V

Multistep Methods

  • 1st idea: use previous steps for computation:

yn+1 = g(yn, yn−1, . . . , yn−q+1)

  • 2nd idea: use integral form of ODE

˙ y(t) = f(t, y(t))

tn+1

  • tn

˙ y(t)dt =

tn+1

  • tn

f(t, y(t))dt y(tn+1) − y(tn) =

tn+1

  • tn

f(t, y(t))dt =?

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

slide-24
SLIDE 24

Lehrstuhl Informatik V

Multistep and Numerical Quadrature

  • 3rd idea: use numerical method for integration

→ interpolate f using a polynomial p: y(tn+1) − y(tn) =

tn+1

  • tn

f(t, y(t))dt ≈

tn+1

  • tn

p(t)dt where p(tj) = f(tj, y(tj)) for j = n − s + 1, . . . , n.

  • compute integral and obtain quadrature rule:

yn+1 = yn +

n

  • j=n−s+1

αj f(tj, yj)

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

slide-25
SLIDE 25

Lehrstuhl Informatik V

Adams-Bashforth

  • s = 1 ⇒ use yn only (leads to Euler’s method):

p(t) = f(tn, yn), yn+1 = yn + τ f(tn, yn)

  • s = 2 ⇒ use yn−1 and yn:

p(t) = tn − t τ f(tn−1, yn−1) + t − tn−1 τ f(tn, yn), yn+1 = yn + τ 2

  • 3f(tn, yn) − f(tn−1, yn−1)
  • usually consistent of s-th order
  • modified at start (no previous values available)

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

slide-26
SLIDE 26

Lehrstuhl Informatik V

Adams-Moulton

  • use idea of Adams-Bashforth, but:

include value yn+1 ⇒ implicit scheme

  • first order: implicit Euler

p(t) = f(tn+1, yn+1), yn+1 = yn + τ f(tn+1, yn+1)

  • second order:

yn+1 = yn + τ 2 (f(tn, yn) + f(tn+1, yn+1))

  • how to obtain yn+1?
  • solve (nonlinear) equation ⇒ difficult!
  • easier and more common: predictor-corrector approach

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

slide-27
SLIDE 27

Lehrstuhl Informatik V

Problems for Numerical Methods for ODE

Possible problems:

  • Ill-Conditioned Problems:

small changes in the input ⇒ big changes in the exact solution of the ODE

  • Instability:

big errors in the numerical solution compared to the exact solution (for arbitrarily small time steps although the method is consistent)

  • Stiffness:

small time steps required for acceptable errors in the approximate solution (although the exact solution is smooth)

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

slide-28
SLIDE 28

Lehrstuhl Informatik V

Ill-Conditioned Problems

  • small changes in input entail completely different results
  • numerical treatment of such problems is always difficult!
  • discriminate:
  • only at critical points?
  • everywhere?
  • possible risks:
  • non-precise input
  • round-off errors,. . .
  • question: what are you interested in?
  • really the solution for specific initial condition?
  • statistical info on the solution?
  • general behaviour (patterns)?

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

slide-29
SLIDE 29

Lehrstuhl Informatik V

Stability

Example: ˙ y(t) = −2y(t) + 1, y(0) = 1

  • exact solution: y(t) = 1

2(e−2t + 1)

  • well-conditioned: yε(0) = 1 + ε ⇒ yε(t) − y(t) = εe−2t
  • use midpoint rule (multistep scheme):

yn+1 = yn−1 + 2τ · f(xn, yn)

  • leads to numerical scheme:

yn+1 = yn−1 + 2τ (1 − 2yn)

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

slide-30
SLIDE 30

Lehrstuhl Informatik V

Stability (2)

Observation:

  • 2-step rule:

yn+1 = yn−1 + 2τ (1 − 2yn) start with exact initial values: y0 = y(0) and y1 = y(τ)

  • numerical results for different sizes of τ:
  • τ = 1.0 ⇒ y9 = −4945.5, y10 = 20953.9
  • τ = 0.1 ⇒ y79 = −1725.3, y80 = 2105.7
  • τ = 0.01 ⇒ y999 = −154.6, y1000 = 158.7
  • midpoint rule is 2nd-order consistent, but does not converge

here: oscillations or instable behaviour

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

slide-31
SLIDE 31

Lehrstuhl Informatik V

Stability (3)

  • reason: difference equation generates spurious solutions
  • analysis: roots µi of characteristic polynomial

y2 = y0 + 4τ(1 − y); all |µi| < 1? Stability of ODE schemes:

  • single step schemes: always stable
  • multistep schemes: additional stability conditions
  • in general:

consistency + stability = convergence

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

slide-32
SLIDE 32

Lehrstuhl Informatik V

Stiff Equations

Example: ˙ y(t) = −1000y(t) + 1000, y(0) = y0 = 2

  • exact solution: y(t) = e−1000t + 1
  • explicit Euler (stable):

yk+1 = yk + τ(−1000yk + 1000) = (1 − 1000τ)yk + 1000τ = (1 − 1000τ)k+1 + 1

  • oscillations and divergence for δt > 0.002
  • Why that? Consistency and stability are asymptotic terms!

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

slide-33
SLIDE 33

Lehrstuhl Informatik V

Stiff Equations – Summary

Typical situation:

  • one term in the ODE demands very small time step
  • but does not contribute much to the solution

Remedy: use implicit (or semi-implicit) methods

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

slide-34
SLIDE 34

Lehrstuhl Informatik V

Summary

Runge-Kutta-methods:

  • multiple evaluations of f (expensive, if f is expensive to compute)
  • stable, well-behaved, easy to implement

Multistep methods:

  • higher order, but only evaluations of f (interesting, if f is

expensive to compute)

  • stability problems; behave “like wild horses”
  • in practice: do not use uniform τ and s

Implicit methods:

  • for stiff equations
  • most often used as corrector scheme

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