Transcript Ordinary Differential Equations (ODEs)

Ordinary Differential Equations (ODEs)
dy (t )
 f (t , y )
dt
Daniel Baur
ETH Zurich, Institut für Chemie- und Bioingenieurwissenschaften
ETH Hönggerberg / HCI F128 – Zürich
E-Mail: [email protected]
http://www.morbidelli-group.ethz.ch/education/index
Daniel Baur / Numerical Methods for Chemical Engineers / Explicit ODE Solvers
1
Problem Definition
 We are going to look at initial value problems of explicit
ODE systems, which means that they can be cast in the
following form
dy (t )
 f (t , y )
dt
y (t0 )  y0
Daniel Baur / Numerical Methods for Chemical Engineers / Explicit ODE Solvers
2
Higher Order ODEs
 Higher order explicit ODEs can be cast into the first order
form by using the following «trick»
d 2 y1
 f (t , y1 )
2
dt
d 2 y1 d  dy1 
 

2
dt
dt  dt 
dy1
y2 :
dt
 dy1
 dt  y2

2
d
y
d
 2  y1  f (t , y )
1
2
 dt
dt
 Therefore, a system of ODEs of order n can be reduced to
first order by adding n-1 variables and equations
Daniel Baur / Numerical Methods for Chemical Engineers / Explicit ODE Solvers
3
Autonomous Systems
 A system is called autonomous if the independent variable
(t) does not explicitly appear in the system
dy (t )
dt
 f ( y)
 It is always possible to transform a non-autonomous
system into an autonomous one
dy (t )
 f (t , y )
dt
 y1  y

 y2  t
 dy1
 dt  f ( y2 , y1 )

 dy 2  1
 dt
Daniel Baur / Numerical Methods for Chemical Engineers / Explicit ODE Solvers
4
Example: A 1-D First Order ODE
 Consider the following initial value problem
 dy (t )
 2y t

 dt
 y (0)  1
y (t ) 
1
3e 2t  2t  1

4
7
 This ODE and its solution follow a general form
(Johann Bernoulli, Basel, 17th – 18th century)
6
5
y
dy (t )
 r (t )  p (t ) y (t )
dt
  p ( s )ds 
p ( s )ds


y (t )  e
k

r
(
t
)
e
d
t





4
3
2
1
0
0.1
0.2
0.3
0.4
Daniel Baur / Numerical Methods for Chemical Engineers / Explicit ODE Solvers
0.5
t
0.6
0.7
0.8
0.9
1
5
Numerical Solution of a 1-D First Order ODE
 dy (t )
 2y t

 dt
 y (0)  1
1
y (t )   3e 2t  2t  1
4
 Since we know the first derivative, Taylor expansion
suggests itself as a first approach
7
yn 1
dy
yn 
dt

6
 (tn 1  tn )  O  tn 1  tn 
5
yn
2

x
 yn  hyn  O  h 2   yn  h  f (tn , yn )  O  h 2 
4
3
 This is known as the Euler method, named after Leonhard
Euler (Basel, 18th century)
2
Euler Method
Exact Solution
1
0
0.1
0.2
0.3
Daniel Baur / Numerical Methods for Chemical Engineers / Explicit ODE Solvers
0.4
0.5
t
0.6
0.7
0.8
0.9
1
6
Some Nomenclature
 On the next slides, the following notation will be used
yn
yn
Value of y at point tn
h
Integration step
Value of the derivative of y at point tn
Daniel Baur / Numerical Methods for Chemical Engineers / Explicit ODE Solvers
7
Categorization of Algorithms
 One step integration algorithms:
yn 1  yn  h  f (tn , yn )
 Multi-step integration algorithms:
yn 1  yn 1  2h  f (tn , yn )
 Explicit integration algorithms:
k1  h  f (tn , yn )

k2  h  f (tn  h, yn  k1 )
 y  y  (k  k ) / 2
n
1
2
 n 1
 Implicit integration algorithms:
yn 1  yn  h  f (tn , yn 1 )
Daniel Baur / Numerical Methods for Chemical Engineers / Explicit ODE Solvers
8
Accuracy of Integration Algorithms
 The local error is always proportional to a known power of
the integration step
yn 2
yn 1  yn  h  f (tn , yn )    h
2
 An algorithm is of order p if it can correctly integrate a
polynomial of order p
 An algorithm of order p always has a local error in the order
of hp+1
 An algorithm is called convergent if for h  0 the global
error decreases (assuming there are no rounding errors)
Daniel Baur / Numerical Methods for Chemical Engineers / Explicit ODE Solvers
9
Examples
 Forward Euler Algorithm
 Explicit Runge-Kutta (p = 2)
yn 1  yn  h  f (tn , yn )
k1  h  f (tn , yn )

k2  h  f (tn  h, yn  k1 )
 y  y  (k  k ) / 2
n
1
2
 n 1
7
7
h = 0.01
h = 0.01
6
6
h = 0.1
5
h = 0.2
4
x
x
5
h = 0.2
4
Error at t = 1 (h = 0.01)
3
3
• Euler Forward: 1.72%
2
2
• RK (p=2):
1
0
0.1
0.2
0.3
0.4
0.5
t
0.6
0.7
0.8
0.9
1
0.01%
1
0
0.1
0.2
0.3
0.4
Daniel Baur / Numerical Methods for Chemical Engineers / Explicit ODE Solvers
0.5
t
0.6
0.7
0.8
0.9
1
10
Conditioning of the Problem
 A system of ODEs is called well conditioned if a small
error in the function (e.g. parameters) and/or in the initial
conditions generate a small error in the solution (without
rounding errors!)
 Example of perturbations / errors:
y  f (t , y )   (t )
y0  y0   (t )
The conditioning of a problem has nothing
to do with the stability of an algorithm!
Daniel Baur / Numerical Methods for Chemical Engineers / Explicit ODE Solvers
11
Example of Bad Conditioning
 Consider the ODE
y  9 y  10 e  t
4
 The general solution reads
10
y (t )  c e9t  e  t
3
10
 If the initial condition y(0) = 1 applies, the constant is 0
10
x
y (t )  e
2
t
 On the other hand, if y(0)10 = 1.0001, the solution reads
1
y (t )  0.0001e9t  e  t
x(0) = 1.0001
x(0) = 1
0
10
-1
10
0
0.2
0.4
0.6
0.8
1
t
Daniel Baur / Numerical Methods for Chemical Engineers / Explicit ODE Solvers
1.2
1.4
1.6
1.8
2
12
Examples for Conditioning
 Consider these initial value problems
y  9 x  10 e  t
y (0)  1
y   e t
y   x
y (0)  1
y (0)  1
 They all have the same solution
y (t )  e  t
 But if we perturb the initial condition (y(0) = 1.0001)
y (t )  e  t  0.0001e9t

Diverging Error
y (t )  e  t  0.0001

Constant Error
y (t )  e  t  0.0001e  t

Converging Error
Daniel Baur / Numerical Methods for Chemical Engineers / Explicit ODE Solvers
13
Measuring the Conditioning
 We can tell a priori if a system is ill conditioned:
If all the eigenvalues of the Jacobian of a system of ODEs
have a positive real part, then the system is ill
conditioned.


Re   J   0
fi
with J ik 
 Ill conditioned system
y j
y  9 y  10 e  t
J 9
y   e t
y   x
J  0  Neutral
J  1  Well Conditioned
 Ill Conditioned
Daniel Baur / Numerical Methods for Chemical Engineers / Explicit ODE Solvers
14
Stability of the Algorithms
 The Forward Euler algorithm is stable only if
  1  hf y  1
where κ is the amplification factor
 Therefore
1. If fy > 0 (the ODE is ill conditioned)
 |κ| > 1 for h > 0
2. If fy < 0 (the ODE is well conditioned)  |κ| < 1 but only if
2
h
fy
Daniel Baur / Numerical Methods for Chemical Engineers / Explicit ODE Solvers
15
Stability and Precision
 When dealing with ODE systems, it is common practice to
reduce the problem of stability to the ODE
y(t )  max  y (t )
where λmax is the largest eigenvalue of the Jacobian of the
system. Note that eigenvalues are typically complex.
 For example, take a look at the Forward Euler method
 The method is stable if
 The local error of the method is
 With the maximum step size
  1  hf y  1
 h
2
fy
y 2 et
2
  h  h
2
2
2
L
hmax  
  max
 2 emax t
max
L
Daniel Baur / Numerical Methods for Chemical Engineers / Explicit ODE Solvers
16
Explicit Runge-Kutta Methods
 A very commonly used method is RK with p = 4
k1  hf (tn , yn )
k2  hf (tn  0.5h, yn  0.5k1 )
k3  hf (tn  0.5h, yn  0.5k2 )
k4  hf (tn  h, yn  k3 )
yn 1  yn 
1
 k1  2k2  2k3  k4 
6
 Advantages
 RK is typically efficient and requires small memory
 It is easy to change the integration step
 Disadvantages
 Efficiency is lower than multi-step methods
 It is difficult to evaluate the local error
Daniel Baur / Numerical Methods for Chemical Engineers / Explicit ODE Solvers
17
Multi-Step Algorithms
 Multi-step methods use not only the current function value
yn to compute the next value yn+1, but also function values
at previous times (yn-1, ...)
 One example are the Adams-Bashforth methods
(here p = 3)
4
14
x 10
70
Exact Solution
Forward Euler
Adams-Bashforth
12
60
h
yn 1  yn   23 f (tn , yn )  16 f (tn 1 , yn 1 )  5 f (tn  2 , yn  2 ) 
12
50
Absolute Error [%]
10
y
8
6
40
Forward Euler
Adams-Bashforth
30
 At the start of the integration, only y0 is known, therefore
we have to use some other method to calculate y1 and y2 to
get the algorithm started (for example the Euler method)
4
20
2
10
0
0
1
2
3
t
4
5
6
0
0
1
2
Daniel Baur / Numerical Methods for Chemical Engineers / Explicit ODE Solvers
3
t
4
5
6
18
Multi-Step Algorithms (Continued)
 Advatanges







Typically very efficient
Integration step is easy to change
Error estimation is simple
Order of the method is easily changed
It is easy to compute y(t) outside the grid points
Usually, fewer function evaluations are required
ODEs with a non-identity mass matrix A are easily solved
A  y  f (t , y )
Daniel Baur / Numerical Methods for Chemical Engineers / Explicit ODE Solvers
19
Stiff Problems
 A problem is considered stiff if e.g. the following apply:
 A linear constant coefficient system is stiff if all of its eigenvalues
have a negative real part and the stiffness ratio is large
dy (t )
dt
 A  y (t )  b ,
 max
SR 
 min
 Stability requirements, rather than accuracy considerations,
constrain the step length
 Some components of the solution show much faster dynamics than
others
 Explicit solvers require an impracticably small h to insure
stability, therefore implicit solvers are preferred
 This requires solving large (non-)linear algebraic systems, which will
be addressed in a later excercise
Daniel Baur / Numerical Methods for Chemical Engineers / Explicit ODE Solvers
20
How does Matlab do it? Non-Stiff Solvers
 ode45: Most general solver, best to try first; Uses an
explicit Runge-Kutta pair (Dormand-Prince, p = 4 and 5)
 ode23: More efficient than ode45 at crude tolerances,
better with moderate stiffness; Uses an explicit RungeKutta pair (Bogacki-Shampine, p = 2 and 3).
 ode113: Better at stringent tolerances and when the ODE
function file is very expensive to evaluate; Uses a variable
order Adams-Bashforth-Moulton method
Daniel Baur / Numerical Methods for Chemical Engineers / Explicit ODE Solvers
21
How does Matlab do it? Stiff Solvers
 ode15s: Use this if ode45 fails or is very inefficient and
you suspect that the problem is stiff; Uses multi-step
numerical differentiation formulas (approximates the
derivative in the next point by extrapolation from the
previous points)
 ode23s may be more efficient than ode15s at crude
tolerances and for some problems (especially if there are
largely different time-scales / dynamics involved); Uses a
modified Rosenbrock formula of order 2 (one-step method)
Daniel Baur / Numerical Methods for Chemical Engineers / Explicit ODE Solvers
22
Matlab Syntax Hints
 All ODE-solvers use the syntax
[T, Y] = ode45(@ode_fun, tspan, y0, …);
 ode_fun is a function handle taking two inputs, a scalar t (current
time point) and a vector y (current function values), and returning as
output a vector of the derivatives dy/dt in these points
 tspan is either a two component vector [tstart, tend] or a vector of
time points where the solution is needed [tstart, t1, t2, ...]
 y0 are the values of y at tstart
 Use parametrizing functions to pass more arguments to
ode_fun if needed
f_new = @(t,y)ode_fun(t, y, k, p);
 This can be done directly in the solver call
 [T, Y] = ode45(@(t,y)ode_fun(t, y, k, p), tspan,
y0);
Daniel Baur / Numerical Methods for Chemical Engineers / Explicit ODE Solvers
23
Exercise
 Consider a batch reactor where these reactions take place
k1
A 
 2B
k2
B 
C
dA
 k1A
dt
dB
 2k1A  k2 B
dt
Daniel Baur / Numerical Methods for Chemical Engineers / Explicit ODE Solvers
24
Assignment
1. Using the following values
 k1 = 1; k2 = 10; A0 = 1; B0 = 0
 Find the Jacobian of the system. Which are the eigenvalues of the
Jacobian? Is the system well conditioned?
 Use ode45 to solve system. Find the time where 99% of A has been
consumed and use this as the final integration time. Plot the results.
2. Find online templates for the Forward Euler method and
the RK method of order 4.
 Add the necessary lines to calculate a step, using the formulas on
slides 10 (Euler) and 17 (RK4). Which step size h is needed to have
a maximum error of 1% compared to ode45 at the end of the
integration? What stepsize was used by ode45?
3. Now change A0 = 0.9; B0 = 0.2; k2 = 1e5 and try ode45
again. What happened? What can you do to resolve this
issue? Plot the resulting concentration profiles.
Daniel Baur / Numerical Methods for Chemical Engineers / Explicit ODE Solvers
25