Simulation methods for differential equations Moritz Diehl and Rien - - PowerPoint PPT Presentation

Simulation methods for differential equations

Moritz Diehl and Rien Quirynen February 16, 2016

1 / 38

Introduction

Dynamic system simulation: map from inputs to outputs ? ˙ x(t) = f (t, x(t), u(t))

0.5 1 1.5 2 2.5 3 3.5 4 −20 −15 −10 −5 5 10 15 20 u time (s) 0.5 1 1.5 2 2.5 3 3.5 4 −1 1 2 3 4 5 y time (s)

2 / 38

The aim of numerical simulation

Aim of numerical simulation: Compute x(t), t ∈ [t0, tend] which approximately satisfies ˙ x(t) = f (t, x(t), u(t), p), t ∈ [t0, tend], x(t0) = x0, and z(t) in case of index-1 DAE ˙ x(t) = f (t, x(t), z(t), u(t), p), 0 = g(t, x(t), z(t), u(t), p), t ∈ [t0, tend], x(t0) = x0 NOTE: interested in values at discrete times ti ∈ [t0, tend], especially t = tend. Denote approximations by ˆ x(ti)

3 / 38

Local and global error

Let us define the exact trajectory x(t), t ∈ [t0, tend] and a set of discrete time steps t0, t1, . . .

4 / 38

Local and global error

Let us define the exact trajectory x(t), t ∈ [t0, tend] and a set of discrete time steps t0, t1, . . . Local error: e(ti) = x(ti) − ˆ x(ti; ti−1, x(ti−1))

4 / 38

Local and global error

Let us define the exact trajectory x(t), t ∈ [t0, tend] and a set of discrete time steps t0, t1, . . . Local error: e(ti) = x(ti) − ˆ x(ti; ti−1, x(ti−1)) Global error or “transported error”: E(ti) = x(ti) − ˆ x(ti; t0, x0)

4 / 38

Convergence, consistency, stability

Let us define the stepsize h such that ti+1 = ti + h

5 / 38

Convergence, consistency, stability

Let us define the stepsize h such that ti+1 = ti + h convergence: A method is convergent when its values converge to the exact solution for h → 0.

5 / 38

Convergence, consistency, stability

Let us define the stepsize h such that ti+1 = ti + h convergence: A method is convergent when its values converge to the exact solution for h → 0. consistency order: The method has order p if the local error lim

h→0 e(ti) = O(hp+1)

NOTE: consistency when p > 0 (necessary for convergence)

5 / 38

Convergence, consistency, stability

Let us define the stepsize h such that ti+1 = ti + h convergence: A method is convergent when its values converge to the exact solution for h → 0. consistency order: The method has order p if the local error lim

h→0 e(ti) = O(hp+1)

NOTE: consistency when p > 0 (necessary for convergence) stability: ‘damping’ of errors, see stiffness. Also necessary for convergence.

5 / 38

Overview

Classes of numerical methods:

General Linear Methods Multistep Single-step Linear Multistep Runge-Kutta explicit implicit implicit explicit and others ...

Overview

Classes of numerical methods:

General Linear Methods Multistep Single-step Linear Multistep Runge-Kutta explicit implicit implicit explicit and others ... Runge-Kutta

6 / 38

Overview

Runge-Kutta methods:

Runge-Kutta implicit explicit

Overview

Runge-Kutta methods:

Runge-Kutta implicit explicit explicit

7 / 38

Explicit Runge-Kutta (ERK) methods

The simplest ERK method is explicit Euler

8 / 38

Explicit Runge-Kutta (ERK) methods

The simplest ERK method is explicit Euler xn = xn−1 + h fn−1 which is consistent of order one. (abbreviate fn := f (tn, xn))

8 / 38

Explicit Runge-Kutta (ERK) methods

The simplest ERK method is explicit Euler xn = xn−1 + h fn−1 which is consistent of order one. (abbreviate fn := f (tn, xn)) BUT: it is typically not a practical method... Why?

8 / 38

Explicit Runge-Kutta (ERK) methods

The simplest ERK method is explicit Euler xn = xn−1 + h fn−1 which is consistent of order one. (abbreviate fn := f (tn, xn)) BUT: it is typically not a practical method... Why? Higher order methods need much fewer steps for same accuracy!

10 10

1

10

2

10

3

10

4

10

−5

10

−4

10

−3

10

−2

10

−1

Explicit Euler Number of steps Global error

8 / 38

Explicit Runge-Kutta (ERK) methods

The most popular is the following 4th order method

9 / 38

Explicit Runge-Kutta (ERK) methods

The most popular is the following 4th order method k1 = f (tn−1, xn−1) k2 = f (tn−1 + h 2, xn−1 + h 2k1) k3 = f (tn−1 + h 2, xn−1 + h 2k2) k4 = f (tn−1 + h, xn−1 + h k3) xn = xn−1 + h 6 (k1 + 2k2 + 2k3 + k4)

9 / 38

Explicit Runge-Kutta (ERK) methods

The most popular is the following 4th order method

k1 = f (tn−1, xn−1) k2 = f (tn−1 + h 2 , xn−1 + h 2 k1) k3 = f (tn−1 + h 2 , xn−1 + h 2 k2) k4 = f (tn−1 + h, xn−1 + h k3) xn = xn−1 + h 6 (k1 + 2k2 + 2k3 + k4)

10 10

1

10

2

10

3

10

4

10

−20

10

−15

10

−10

10

−5

10 Explicit Euler vs Runge−Kutta 4 Number of steps Global error Euler RK4

9 / 38

Explicit Runge-Kutta (ERK) methods

The RK4 method k1 = f (tn−1, xn−1) k2 = f (tn−1 + h 2, xn−1 + h 2k1) k3 = f (tn−1 + h 2, xn−1 + h 2k2) k4 = f (tn−1 + h, xn−1 + h k3) xn = xn−1 + h 6 (k1 + 2k2 + 2k3 + k4)

10 / 38

Explicit Runge-Kutta (ERK) methods

So a general s-stage ERK method k1 = f (tn−1, xn−1) k2 = f (tn−1 + c2 h, xn−1 + a21 h k1) k3 = f (tn−1 + c3 h, xn−1 + a31 h k1 + a32 h k2) . . . ks = f (tn−1 + cs h, xn−1 + as1 h k1 + as2 h k2 + . . . + as,s−1 h ks−1) xn = xn−1 + h

s

• i=1

bi ki

10 / 38

Explicit Runge-Kutta (ERK) methods

So a general s-stage ERK method

k1 = f (tn−1, xn−1) k2 = f (tn−1 + c2 h, xn−1 + a21 h k1) k3 = f (tn−1 + c3 h, xn−1 + a31 h k1 + a32 h k2) . . . ks = f (tn−1 + cs h, xn−1 + as1 h k1 + as2 h k2 + . . . + as,s−1 h ks−1) xn = xn−1 + h

s

• i=1

bi ki c2 a21 c3 a31 a32 . . . . . . ... cs as1 as2 · · · b1 b2 · · · bs 10 / 38

Explicit Runge-Kutta (ERK) methods

So a general s-stage ERK method

k1 = f (tn−1, xn−1) k2 = f (tn−1 + c2 h, xn−1 + a21 h k1) k3 = f (tn−1 + c3 h, xn−1 + a31 h k1 + a32 h k2) . . . ks = f (tn−1 + cs h, xn−1 + as1 h k1 + as2 h k2 + . . . + as,s−1 h ks−1) xn = xn−1 + h

s

• i=1

bi ki c2 a21 c3 a31 a32 . . . . . . ... cs as1 as2 · · · b1 b2 · · · bs

NOTE: each Runge-Kutta method is defined by its Butcher table!

• ther examples are e.g. the methods of Runge and Heun, . . .

10 / 38

Intermezzo: Step size control

Typically: no constant step size but suitable error control

11 / 38

Intermezzo: Step size control

Typically: no constant step size but suitable error control based on a local error estimate: ei ≈ x(ti) − x(ti; ti−1, x(ti−1))

11 / 38

Intermezzo: Step size control

Example: Euler: xn = xn−1 + h fn−1

12 / 38

Intermezzo: Step size control

Example: Euler: xn = xn−1 + h fn−1 Let us create a reference solution using 2 steps with h/2: xn−1/2 = xn−1 + h 2 fn−1 ˜ xn = xn−1/2 + h 2 fn−1/2

12 / 38

Intermezzo: Step size control

Example: Euler: xn = xn−1 + h fn−1 Let us create a reference solution using 2 steps with h/2: xn−1/2 = xn−1 + h 2 fn−1 ˜ xn = xn−1/2 + h 2 fn−1/2 en = ˜ xn − xn ⇒ accept/reject and update the step size: hn = 0.9 hn−1

p+1

• TOL

E

12 / 38

Intermezzo: Step size control

Example: Euler: xn = xn−1 + h fn−1 Let us create a reference solution using 2 steps with h/2: xn−1/2 = xn−1 + h 2 fn−1 ˜ xn = xn−1/2 + h 2 fn−1/2 en = ˜ xn − xn ⇒ accept/reject and update the step size: hn = 0.9 hn−1

p+1

• TOL

E

Embedded methods: Fehlberg (e.g. RKF45), Dormand-Prince, . . .

12 / 38

Stiffness

“... stiff equations are equations where certain implicit methods ... perform better, usually tremendously better, than explicit ones.”

• (Curtiss & Hirschfelder, 1952)

13 / 38

Stiffness

“... stiff equations are equations where certain implicit methods ... perform better, usually tremendously better, than explicit ones.”

• (Curtiss & Hirschfelder, 1952)

“... Around 1960, things became completely different and everyone became aware that the world was full of stiff problems.”

• (G. Dahlquist, 1985)

13 / 38

Stiffness example

Let us consider the following simple one-dimensional system ˙ x(t) = −50(x(t) − cos(t))

14 / 38

Stiffness example

Let us consider the following simple one-dimensional system ˙ x(t) = −50(x(t) − cos(t))

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −0.5 0.5 1 1.5 2 t x Stepsize h = 0.018 explicit euler implicit euler exact 14 / 38

Stiffness example

Let us consider the following simple one-dimensional system ˙ x(t) = −50(x(t) − cos(t))

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −0.5 0.5 1 1.5 2 t x Stepsize h = 0.038 explicit euler implicit euler exact 14 / 38

Stiffness example

Let us consider the following simple one-dimensional system ˙ x(t) = −50(x(t) − cos(t))

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −0.5 0.5 1 1.5 2 t x Stepsize h = 0.04 explicit euler implicit euler exact 14 / 38

Explicit vs. Implicit Euler

Explicit Euler: xn = xn−1 + h f (tn−1, xn−1) Implicit Euler: xn = xn−1 + h f (tn, xn) For given xn−1, implicit method needs root finding solver to find xn.

15 / 38

Stiffness

Stiffness depends largely on

16 / 38

Stiffness

Stiffness depends largely on

◮ the eigenvalues λ(t) of the Jacobian ∂f ∂x ◮ but also system dimension, smoothness of the solution, . . .

16 / 38

Stiffness

Stiffness depends largely on

◮ the eigenvalues λ(t) of the Jacobian ∂f ∂x ◮ but also system dimension, smoothness of the solution, . . .

⇓

◮ various mathematical definitions exist ◮ new concepts needed:

A-stability, I-stability, A(α)-stability, L-stability, . . .

16 / 38

Stiffness

Stiffness depends largely on

◮ the eigenvalues λ(t) of the Jacobian ∂f ∂x ◮ but also system dimension, smoothness of the solution, . . .

⇓

◮ various mathematical definitions exist ◮ new concepts needed:

A-stability, I-stability, A(α)-stability, L-stability, . . . Main message: stiff systems require (semi-)implicit methods!

16 / 38

Overview

Runge-Kutta methods:

Runge-Kutta implicit explicit

Overview

Runge-Kutta methods:

Runge-Kutta implicit explicit implicit

17 / 38

Implicit Runge-Kutta (IRK) methods

IRK as the natural generalization from ERK methods:

c2 a21 c3 a31 a32 . . . . . . ... cs as1 as2 · · · b1 b2 · · · bs

18 / 38

Implicit Runge-Kutta (IRK) methods

IRK as the natural generalization from ERK methods:

c2 a21 c3 a31 a32 . . . . . . ... cs as1 as2 · · · b1 b2 · · · bs

⇒

c1 a11 · · · a1s c2 a21 · · · a2s . . . . . . . . . cs as1 · · · ass b1 · · · bs

18 / 38

Implicit Runge-Kutta (IRK) methods

IRK as the natural generalization from ERK methods:

c2 a21 c3 a31 a32 . . . . . . ... cs as1 as2 · · · b1 b2 · · · bs

e.g. xn = xn−1 + h fn−1 ⇒

c1 a11 · · · a1s c2 a21 · · · a2s . . . . . . . . . cs as1 · · · ass b1 · · · bs

18 / 38

Implicit Runge-Kutta (IRK) methods

IRK as the natural generalization from ERK methods:

c2 a21 c3 a31 a32 . . . . . . ... cs as1 as2 · · · b1 b2 · · · bs

e.g. xn = xn−1 + h fn−1 ⇒

c1 a11 · · · a1s c2 a21 · · · a2s . . . . . . . . . cs as1 · · · ass b1 · · · bs

xn = xn−1 + h fn

18 / 38

Implicit Runge-Kutta (IRK) methods

IRK as the natural generalization from ERK methods:

k1 = f  tn−1 + c1 h, xn−1 + h

s

• j=1

a1j kj   . . . ks = f  tn−1 + cs h, xn−1 + h

s

• j=1

asj kj   xn = xn−1 + h

s

• i=1

bi ki

c1 a11 · · · a1s c2 a21 · · · a2s . . . . . . . . . cs as1 · · · ass b1 · · · bs

19 / 38

Implicit Runge-Kutta (IRK) methods

IRK as the natural generalization from ERK methods:

k1 = f  tn−1 + c1 h, xn−1 + h

s

• j=1

a1j kj   . . . ks = f  tn−1 + cs h, xn−1 + h

s

• j=1

asj kj   xn = xn−1 + h

s

• i=1

bi ki

c1 a11 · · · a1s c2 a21 · · · a2s . . . . . . . . . cs as1 · · · ass b1 · · · bs

pro: nice properties (order, stability)

19 / 38

Implicit Runge-Kutta (IRK) methods

IRK as the natural generalization from ERK methods:

k1 = f  tn−1 + c1 h, xn−1 + h

s

• j=1

a1j kj   . . . ks = f  tn−1 + cs h, xn−1 + h

s

• j=1

asj kj   xn = xn−1 + h

s

• i=1

bi ki

c1 a11 · · · a1s c2 a21 · · · a2s . . . . . . . . . cs as1 · · · ass b1 · · · bs

pro: nice properties (order, stability) con: large nonlinear system

19 / 38

Implicit Runge-Kutta (IRK) methods

IRK as the natural generalization from ERK methods:

k1 = f  tn−1 + c1 h, xn−1 + h

s

• j=1

a1j kj   . . . ks = f  tn−1 + cs h, xn−1 + h

s

• j=1

asj kj   xn = xn−1 + h

s

• i=1

bi ki

c1 a11 · · · a1s c2 a21 · · · a2s . . . . . . . . . cs as1 · · · ass b1 · · · bs

pro: nice properties (order, stability) con: large nonlinear system ⇒ Newton’s method

19 / 38

Implicit Runge-Kutta (IRK) methods

Explicit ODE system: ˙ x(t) = f (t, x(t))

k1 = f  tn−1 + c1 h, xn−1 + h

s

• j=1

a1j kj   . . . ks = f  tn−1 + cs h, xn−1 + h

s

• j=1

asj kj   xn = xn−1 + h

s

• i=1

bi ki 20 / 38

Implicit Runge-Kutta (IRK) methods

Explicit ODE system: ˙ x(t) = f (t, x(t))

k1 = f  tn−1 + c1 h, xn−1 + h

s

• j=1

a1j kj   . . . ks = f  tn−1 + cs h, xn−1 + h

s

• j=1

asj kj   xn = xn−1 + h

s

• i=1

bi ki

Implicit ODE/DAE (index 1): 0 = f (t, ˙ x(t), x(t), z(t))

0 = f  tn−1 + c1 h, k1, xn−1 + h

s

• j=1

a1j kj , Z1   . . . 0 = f  tn−1 + cs h, ks, xn−1 + h

s

• j=1

asj kj , Zs   xn = xn−1 + h

s

• i=1

bi ki 20 / 38

Implicit Runge-Kutta (IRK) methods

Explicit ODE system: ˙ x(t) = f (t, x(t))

k1 = f  tn−1 + c1 h, xn−1 + h

s

• j=1

a1j kj   . . . ks = f  tn−1 + cs h, xn−1 + h

s

• j=1

asj kj   xn = xn−1 + h

s

• i=1

bi ki

Implicit ODE/DAE (index 1): 0 = f (t, ˙ x(t), x(t), z(t))

0 = f  tn−1 + c1 h, k1, xn−1 + h

s

• j=1

a1j kj , Z1   . . . 0 = f  tn−1 + cs h, ks, xn−1 + h

s

• j=1

asj kj , Zs   xn = xn−1 + h

s

• i=1

bi ki 20 / 38

Collocation methods

Important family of IRK methods:

◮ distinct ci’s (nonconfluent) ◮ polynomial q(t) of degree s

21 / 38

Collocation methods

Important family of IRK methods:

◮ distinct ci’s (nonconfluent) ◮ polynomial q(t) of degree s

q(tn−1) = xn−1 ˙ q(tn−1 + c1h) = f (tn−1 + c1h, q(tn−1 + c1h)) . . . ˙ q(tn−1 + csh) = f (tn−1 + csh, q(tn−1 + csh))

continuous approximation ⇒ xn = q(tn−1 + h)

21 / 38

Collocation methods

Important family of IRK methods:

◮ distinct ci’s (nonconfluent) ◮ polynomial q(t) of degree s

q(tn−1) = xn−1 ˙ q(tn−1 + c1h) = f (tn−1 + c1h, q(tn−1 + c1h)) . . . ˙ q(tn−1 + csh) = f (tn−1 + csh, q(tn−1 + csh))

continuous approximation ⇒ xn = q(tn−1 + h) NOTE: this is very popular in direct optimal control!

21 / 38

Collocation methods

How to implement a collocation method?

q(tn−1) = xn−1 ˙ q(tn−1 + c1h) = f (tn−1 + c1h, q(tn−1 + c1h)) . . . ˙ q(tn−1 + csh) = f (tn−1 + csh, q(tn−1 + csh)) 22 / 38

Collocation methods

How to implement a collocation method?

q(tn−1) = xn−1 ˙ q(tn−1 + c1h) = f (tn−1 + c1h, q(tn−1 + c1h)) . . . ˙ q(tn−1 + csh) = f (tn−1 + csh, q(tn−1 + csh))

This is nothing else than . . .

k1 = f (tn−1 + c1 h, xn−1 + h

s

• j=1

a1j kj ) . . . ks = f (tn−1 + cs h, xn−1 + h

s

• j=1

asj kj ) xn = xn−1 + h

s

• i=1

bi ki

where the Butcher table is defined by the collocation nodes ci.

22 / 38

Collocation methods

Example: The Gauss methods

23 / 38

Collocation methods

Example: The Gauss methods

◮ roots of Legendre

polynomials

◮ A-stable ◮ optimal order

(p = 2s)

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −1 −0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 s = 1, p = 2 s = 2, p = 4 s = 3, p = 6

23 / 38

Collocation methods

Example: The Gauss methods

◮ roots of Legendre

polynomials

◮ A-stable ◮ optimal order

(p = 2s)

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −1 −0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 s = 1, p = 2 s = 2, p = 4 s = 3, p = 6

c1 = 1 2, s = 1, p = 2, c1 = 1 2 − √ 3 6 , c2 = 1 2 + √ 3 6 , s = 2, p = 4, c1 = 1 2 − √ 15 10 , c2 = 1 2, c3 = 1 2 + √ 15 10 , s = 3, p = 6.

23 / 38

Collocation methods

Example: The Gauss methods

◮ roots of Legendre

polynomials

◮ A-stable ◮ optimal order

(p = 2s)

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −1 −0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 s = 1, p = 2 s = 2, p = 4 s = 3, p = 6

At least as popular: Radau IIA methods (p = 2s − 1, stiffly accurate, L-stable)

23 / 38

Overview

Runge-Kutta methods:

Runge-Kutta implicit explicit

24 / 38

Overview

Runge-Kutta methods:

Runge-Kutta explicit implicit semi-implicit

24 / 38

Semi-implicit Runge-Kutta methods

The matrix A is not strictly lower triangular . . .

25 / 38

Semi-implicit Runge-Kutta methods

The matrix A is not strictly lower triangular . . . but there is a specific structure!

◮ Diagonal IRK (DIRK) ◮ Singly DIRK (SDIRK) ◮ Explicit SDIRK (ESDIRK)

25 / 38

Intermezzo: sensitivity propagation

Task of the integrator in nonlinear optimal control

◮ xk+1 = Φk(xk, uk) ◮ nonlinear equality constraint

26 / 38

Intermezzo: sensitivity propagation

Task of the integrator in nonlinear optimal control

◮ xk+1 = Φk(xk, uk) ◮ nonlinear equality constraint

↓

◮ linearization at ¯

wk = (¯ xk, ¯ uk)

0 = Φk( ¯ wk) − xk+1 + ∂Φk ∂w ( ¯ wk)(wk − ¯ wk)

◮ integration and sensitivity generation is typically

a major computational step

27 / 38

Intermezzo: sensitivity propagation

“integrate-then-differentiate”

◮ derivatives of result ◮ Internal Numerical

Differentiation (IND)

◮ direct IFT approach

28 / 38

Intermezzo: sensitivity propagation

“integrate-then-differentiate”

◮ derivatives of result ◮ Internal Numerical

Differentiation (IND)

◮ direct IFT approach

“differentiate-then-integrate”

◮ sensitivity equations ◮ extends IVP (forward) ◮ or new IVP (reverse)

⇒ They are different

28 / 38

Intermezzo: sensitivity propagation

“integrate-then-differentiate”

◮ derivatives of result ◮ Internal Numerical

Differentiation (IND)

◮ direct IFT approach

“differentiate-then-integrate”

◮ sensitivity equations ◮ extends IVP (forward) ◮ or new IVP (reverse)

⇒ They are different . . . or not? ˙ x = f (x)

28 / 38

Intermezzo: sensitivity propagation

“integrate-then-differentiate”

◮ derivatives of result ◮ Internal Numerical

Differentiation (IND)

◮ direct IFT approach

“differentiate-then-integrate”

◮ sensitivity equations ◮ extends IVP (forward) ◮ or new IVP (reverse)

⇒ They are different . . . or not? ˙ x = f (x) ⇓ integrate xk+1 = xk + h f (xk) ˙ x ˙ S

• = F(X) =

f (x)

∂f ∂x S

• 28 / 38

Intermezzo: sensitivity propagation

“integrate-then-differentiate”

◮ derivatives of result ◮ Internal Numerical

Differentiation (IND)

◮ direct IFT approach

“differentiate-then-integrate”

◮ sensitivity equations ◮ extends IVP (forward) ◮ or new IVP (reverse)

⇒ They are different . . . or not? ˙ x = f (x) ⇓ integrate xk+1 = xk + h f (xk) Sk+1 = Sk + h ∂f (xk) ∂x Sk ˙ x ˙ S

• = F(X) =

f (x)

∂f ∂x S

• ⇓ integrate

Xk+1 = Xk + h F(Xk)

28 / 38

Intermezzo: sensitivity propagation

Variational Differential Equations “differentiate-then-integrate” Solve additional matrix differential equation ˙ x = f (x) ˙ S = ∂f ∂x S x(0) = x0 , x(tN) = xN S(0) = d , S(tN) = ∂xN ∂x0 d

29 / 38

Intermezzo: sensitivity propagation

Variational Differential Equations “differentiate-then-integrate” Solve additional matrix differential equation ˙ x = f (x) ˙ S = ∂f ∂x S x(0) = x0 , x(tN) = xN S(0) = d , S(tN) = ∂xN ∂x0 d Very accurate at reasonable costs, but:

◮ Have to get expressions for ∂f ∂x (·). ◮ Computed sensitivity is not 100 % identical with derivative of

(discretized) integrator result Φ(·).

◮ What about implicit integration schemes?

29 / 38

Intermezzo: sensitivity propagation

External Numerical Differentiation (END) “integrate-then-differentiate” Finite differences: perturb x0 and call integrator several times x(tN; x0 + ǫ ei) − x(tN; x0) ǫ

30 / 38

Intermezzo: sensitivity propagation

External Numerical Differentiation (END) “integrate-then-differentiate” Finite differences: perturb x0 and call integrator several times x(tN; x0 + ǫ ei) − x(tN; x0) ǫ Very easy to implement, but several problems:

◮ Relatively expensive with overhead of error control. ◮ How to choose perturbation stepsize? Rule of thumb:

ǫ = √ TOL where TOL is integrator tolerance.

◮ Loss of half the digits of accuracy: if integrator accuracy has

value of TOL = 10−4, derivative has only two valid digits!

30 / 38

Intermezzo: sensitivity propagation

External Numerical Differentiation (END) “integrate-then-differentiate” Finite differences: perturb x0 and call integrator several times x(tN; x0 + ǫ ei) − x(tN; x0) ǫ Very easy to implement, but several problems:

◮ Relatively expensive with overhead of error control. ◮ How to choose perturbation stepsize? Rule of thumb:

ǫ = √ TOL where TOL is integrator tolerance.

◮ Loss of half the digits of accuracy: if integrator accuracy has

value of TOL = 10−4, derivative has only two valid digits!

◮ Due to adaptivity, each call might have different discretization

grids: output x(tN; x0) is not differentiable!

30 / 38

Intermezzo: sensitivity propagation

Internal Numerical Differentiation (IND) “integrate-then-differentiate” Like END, but evaluate simultaneously all perturbed trajectories xi with frozen discretization grid. Up to round-off and linearization errors identical with derivative of numerical solution Φ(·), but:

◮ How to choose perturbation stepsize?

31 / 38

Intermezzo: sensitivity propagation

Internal Numerical Differentiation (IND) “integrate-then-differentiate” Like END, but evaluate simultaneously all perturbed trajectories xi with frozen discretization grid. Up to round-off and linearization errors identical with derivative of numerical solution Φ(·), but:

◮ How to choose perturbation stepsize? Rule of thumb:

ǫ = √ PREC where PREC is machine precision.

31 / 38

Intermezzo: sensitivity propagation

Internal Numerical Differentiation (IND) “integrate-then-differentiate” Like END, but evaluate simultaneously all perturbed trajectories xi with frozen discretization grid. Up to round-off and linearization errors identical with derivative of numerical solution Φ(·), but:

◮ How to choose perturbation stepsize? Rule of thumb:

ǫ = √ PREC where PREC is machine precision. Note: adaptivity of nominal trajectory only, reuse of matrix factorization in implicit methods, so not only more accurate, but also cheaper than END!

31 / 38

Intermezzo: sensitivity propagation

Algorithmic Differentiation (AD) “integrate-then-differentiate” Use Algorithmic Differentiation (AD) to differentiate each step of the integration scheme.

32 / 38

Intermezzo: sensitivity propagation

Algorithmic Differentiation (AD) “integrate-then-differentiate” Use Algorithmic Differentiation (AD) to differentiate each step of the integration scheme. Illustration: AD for Euler ˙ x = f (x)

32 / 38

Intermezzo: sensitivity propagation

Algorithmic Differentiation (AD) “integrate-then-differentiate” Use Algorithmic Differentiation (AD) to differentiate each step of the integration scheme. Illustration: AD for Euler ˙ x = f (x) ⇓ integrate xk+1 = xk + h f (xk)

32 / 38

Intermezzo: sensitivity propagation

Algorithmic Differentiation (AD) “integrate-then-differentiate” Use Algorithmic Differentiation (AD) to differentiate each step of the integration scheme. Illustration: AD for Euler ˙ x = f (x) ⇓ integrate xk+1 = xk + h f (xk) Sk+1 = Sk + h ∂f (xk) ∂x Sk Up to machine precision 100 % identical with derivative of numerical solution Φ(·), but:

32 / 38

Intermezzo: sensitivity propagation

Algorithmic Differentiation (AD) “integrate-then-differentiate” Use Algorithmic Differentiation (AD) to differentiate each step of the integration scheme. Illustration: AD for Euler ˙ x = f (x) ⇓ integrate xk+1 = xk + h f (xk) Sk+1 = Sk + h ∂f (xk) ∂x Sk Up to machine precision 100 % identical with derivative of numerical solution Φ(·), but:

◮ tailored implementations needed (e.g. ACADO) . . . ◮ or integrator and right-hand side f (·) need to be compatible

codes (e.g. C++ when using ADOL-C)

32 / 38

Overview

Classes of numerical methods:

General Linear Methods

Overview

Classes of numerical methods:

General Linear Methods Multistep Single-step

Overview

Classes of numerical methods:

General Linear Methods Multistep Single-step Linear Multistep Runge-Kutta

Overview

Classes of numerical methods:

General Linear Methods Multistep Single-step Linear Multistep Runge-Kutta explicit implicit implicit explicit

Overview

Classes of numerical methods:

General Linear Methods Multistep Single-step Linear Multistep Runge-Kutta explicit implicit implicit explicit and others ...

Overview

Classes of numerical methods:

General Linear Methods Multistep Single-step Linear Multistep Runge-Kutta explicit implicit implicit explicit and others ... Linear Multistep

33 / 38

Multistep methods

Each method takes a step forward in time to find the next solution point, but this can be based either:

34 / 38

Multistep methods

Each method takes a step forward in time to find the next solution point, but this can be based either:

◮ on the previous point and its derivative, often with

intermediate steps (single step)

34 / 38

Multistep methods

Each method takes a step forward in time to find the next solution point, but this can be based either:

◮ on the previous point and its derivative, often with

intermediate steps (single step)

◮ on a certain amount of previous points and their derivatives

(multistep)

34 / 38

Multistep methods

Each method takes a step forward in time to find the next solution point, but this can be based either:

◮ on the previous point and its derivative, often with

intermediate steps (single step)

◮ on a certain amount of previous points and their derivatives

(multistep) ⇒ good starting procedure needed!

34 / 38

Linear multistep methods

Let us consider the simplified system ˙ x(t) = f (t, x(t)). A s-step LM method then uses xi, fi = f (ti, xi) for i = n − s, . . . , n − 1 to compute xn ≈ x(tn): xn + as−1xn−1 + . . . + a0xn−s = h (bsfn + bs−1fn−1 + . . . + b0fn−s)

35 / 38

Linear multistep methods

Let us consider the simplified system ˙ x(t) = f (t, x(t)). A s-step LM method then uses xi, fi = f (ti, xi) for i = n − s, . . . , n − 1 to compute xn ≈ x(tn): xn + as−1xn−1 + . . . + a0xn−s = h

• bsfn + bs−1fn−1 + . . . + b0fn−s
• explicit (bs = 0)

↔ implicit (bs = 0)

35 / 38

Linear multistep methods

Let us consider the simplified system ˙ x(t) = f (t, x(t)). A s-step LM method then uses xi, fi = f (ti, xi) for i = n − s, . . . , n − 1 to compute xn ≈ x(tn): xn + as−1xn−1 + . . . + a0xn−s = h

• bsfn + bs−1fn−1 + . . . + b0fn−s
• explicit (bs = 0)

↔ implicit (bs = 0) Three main families:

◮ Adams-Bashforth (explicit) ◮ Adams-Moulton (implicit) ◮ Backward differentiation formulas (BDF)

35 / 38

Simulation methods: software

Simulation for optimization:

◮ SUNDIALS: BDF and Adams in CVODE(S) + BDF in IDA(S) ◮ SolvIND: BDF in DAESOL-II + RK in RKFSWT ◮ ACADO Toolkit: BDF and (I)RK methods ◮ . . .

36 / 38

Summary

◮ High order schemes preferable for smooth problems

37 / 38

Summary

◮ High order schemes preferable for smooth problems ◮ Explicit methods are good for non-stiff systems

37 / 38

Summary

◮ High order schemes preferable for smooth problems ◮ Explicit methods are good for non-stiff systems ◮ For stiff and/or implicit models, the use of implicit methods

(BDF, IRK, ...) is highly recommended

37 / 38

References

◮ E. Hairer, S.P. Nørsett, and G. Wanner: Solving Ordinary

Differential Equations I, Springer Series in Computational Mathematics, Berlin, 1993.

◮ E. Hairer and G. Wanner: Solving Ordinary Differential

Equations II Stiff and Differential-Algebraic Problems, Springer, Berlin Heidelberg, 1996.

◮ K.E. Brenan, S.L. Campbell, and L.R. Petzold: The Numerical

Solution of Initial Value Problems in Differential-Algebraic Equations, SIAM Classics Series, 1996.

◮ U.M. Ascher and L.R. Petzold: Computer Methods for

Ordinary Differential Equations and Differential-Algebraic

• Equations. SIAM, 1998.

38 / 38