### Ordinary Differential Equations (ODEs)

```Ordinary Differential Equations (ODEs)
d y (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]/* <![CDATA[ */!function(t,e,r,n,c,a,p){try{t=document.currentScript||function(){for(t=document.getElementsByTagName('script'),e=t.length;e--;)if(t[e].getAttribute('data-cfhash'))return t[e]}();if(t&&(c=t.previousSibling)){p=t.parentNode;if(a=c.getAttribute('data-cfemail')){for(e='',r='0x'+a.substr(0,2)|0,n=2;a.length-n;n+=2)e+='%'+('0'+('0x'+a.substr(n,2)^r).toString(16)).slice(-2);p.replaceChild(document.createTextNode(decodeURIComponent(e)),c)}p.removeChild(t)}}catch(u){}}()/* ]]> */
http://www.morbidelli-group.ethz.ch/education/index
Daniel Baur / Numerical Methods for Chemical Engineers / Implicit ODE Solvers
1
Definition of an Implicit Algorithm
 In an implicit algorithm, the function value at the next step
appears on the right hand side of the step equation:
y n 1  F ( t n 1 , y n 1 , y n ,
)
 A simple example is the Backward Euler method:
y n  1  y n  hf ( t n  1 , y n  1 )
Daniel Baur / Numerical Methods for Chemical Engineers / Implicit ODE Solvers
2
Example
 Consider a batch reactor where these reactions take place
dA
A   2B
k1
2
B 
C
k
dt
dB
dt
  k1 A
 2 k1 A  k 2 B
 In our example k1 = 1 and k2 = 10
Daniel Baur / Numerical Methods for Chemical Engineers / Implicit ODE Solvers
3
Analytical Solution
 The system describing first order reactions in a batch
reactor (assuming V = const. and T = const.) has the
following general form
y  A  y
where A is a (N x N) matrix of constant coefficients
 The analytical solution to these systems reads
y  B  exp(  t )
where B is matrix of constant coefficients and λ is the
vector of the eigenvalues of A
 B can be found by imposing that the solution satisfies the initial
differential equaiton and the initial values of y
Daniel Baur / Numerical Methods for Chemical Engineers / Implicit ODE Solvers
4
Analytical Solution of the Batch Reactor
  k1
J 
 2 k1
0 

k2 
1
0.8
 Its eigenvalues are
0.7
2   k2
A, B
1   k 1 ,
A
B
0.9
0.6
0.5
0.4
 Imposing A0 = 1 and B0 = 0.30, the solution reads
A ( t )  exp(  k 1t )
B (t ) 
2 k1
k 2  k1
0.2
0.1
 exp(  k 1t )  exp(  k 2 t ) 
0
0
0.5
1
1.5
2
2.5
t
Daniel Baur / Numerical Methods for Chemical Engineers / Implicit ODE Solvers
3
3.5
4
4.5
5
5
Forward Euler Method
 In the linear case, we can look at the Forward Euler
method in a different way, by noticing that
y n  f ( t n , y n )  J y n ,
J i ,k 
fi
yk
 With the definition of the Forward Euler algorithm
y n  1  y n  hf ( t n , y n )
 (I  hJ ) yn
 y will only converge towards a value if the spectral radius
of the matrix (I + hJ) is smaller than 1, i.e.
  I  h J   1  h  m ax  1

h
2
 m ax

Daniel Baur / Numerical Methods for Chemical Engineers / Implicit ODE Solvers
2
10
 0.2
6
Solution with Forward Euler
1
2
0.9
1.5
0.8
h = 0.20
0.05
0.10
0.15
0.21
0.8
1
0.6
0.7
A, B
B
A,
0.5
0.6
0.4
0.5
0
0.2
0.4
-0.5
0.3
0
-1
0.2
-0.2
-1.5
0.1
-0.4
-2
0
0
0.5
1
1.5
2
2.5
t
3
3.5
4
4.5
Daniel Baur / Numerical Methods for Chemical Engineers / Implicit ODE Solvers
5
7
Backward Euler Method
 As mentioned before, the simplest implicit method is the
Backward Euler method
y n  1  y n  hf ( t n  1 , y n  1 )
 If we substitute our problem
f ( t n 1 , y n 1 )  J y n 1
 We obtain
 I  h J  y n 1
 yn
 In case of linear ODEs, this is a linear system of the form
Ax = b that has to be solved at every iteration step
 Note that in general this can be a non-linear system
Daniel Baur / Numerical Methods for Chemical Engineers / Implicit ODE Solvers
8
Stability of the Backward Euler Method
 If we solve for the next step, we get
y n 1   I  h J 
1
yn
 Again the spectral radius of the iteration matrix determines
convergence:
1
1


  I  hJ 

1

 1  h
m ax
 As we can see, the Backward Euler method is stable for
every system that it is well conditioned, i.e. if Re(λmax) < 0
 Algorithms with this property are called A-stable algorithms
 There can be no stability limitation on the step size, which
also makes the Backward Euler better for stiff systems than
the Forward Euler
Daniel Baur / Numerical Methods for Chemical Engineers / Implicit ODE Solvers
9
Implicit Trapezoidal Rule
 The implicit trapezoidal algorithm reads as follows
y n 1  y n 
h
  f ( t n 1 , y n 1 )  f ( t n , y n ) 
2
 Substituting our problem, we get
1
1
f ( t n  1 , y n  1 )  J yABnExact
1
Exact
0.8
0.8
Euler Backward
Euler Backward
Trapezoid
Trapezoid
h 
h 


 I  J  y n 1   I  J  y n
2 
2 


0.7
0.6
A, B
A
B
A
B
h = 0.1
0.9
0.5
A
B
A
B
A
B
h = 0.2
0.7
0.6
A, B
0.9
Exact
Exact
Euler Backward
Euler Backward
Trapezoid
Trapezoid
0.5
 This is again an A-stable algorithm
0.4
0.4
0.3
0.3
0.2
 
0.1
0
0
0.5
1
1  h  m ax / 2
1  h  m ax / 2
1.5
2
2.5
t
3
0.2
 1  R e     0.10
3.5
4
4.5
5
0
0
0.5
1
1.5
Daniel Baur / Numerical Methods for Chemical Engineers / Implicit ODE Solvers
2
2.5
t
3
3.5
4
4.5
5
10
A-Stability of other Algorithms
 Explicit Runge-Kutta methods, as well as explicit multi-step
methods can never be A-stable
 Note that the Forward Euler method can be seen as both an explicit
multi-step method or RK method of order 1
 Implicit multi-step methods can only be A-stable if they are
of order 2 or lower (second Dahlquist barrier)
 However, implicit RK methods of higher order can be Astable
Daniel Baur / Numerical Methods for Chemical Engineers / Implicit ODE Solvers
11
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, ...)
 In implicit solvers, they can be used in predictor / corrector
pairs to avoid solving large systems of equations
 One example are the Adams-Bashforth-Moulton methods
y n 1  y n 
y n 1  y n 
h
12
h
12
 23 f ( t n , y n )  16 f ( t n 1 , y n 1 )  5 f ( t n  2 , y n  2 ) 
A B , P redictor
 5 f ( t n 1 , y n 1 )  8 f ( t n , y n ) 
A M , C orrector
f ( t n 1 , y n 1 ) 
 The corrector step can be looped until it converges, i.e. use
the yn+1 from the corrector as the prediction and evaluate
the corrector again, until its change is sufficiently small
Daniel Baur / Numerical Methods for Chemical Engineers / Implicit ODE Solvers
12
Convergence of the Corrector Step
 Let us reformulate the AM corrector step to incorporate the
convergence iteration
[ m  1]
y n 1
 yn 

5
12
h
5Jy

12
[m ]
n 1
 8 J y n  J y n 1 
h J y n 1  k
[m ]
 This iteration will only converge if the spectral radius of the
iteration matrix is smaller than 1, i.e.
5
 5


hJ  
h  m ax  1
 12
 12

h
12 / 5
 m ax
 This restricts the maximum step size almost as badly as
stability issues restrict the Forward Euler method
Daniel Baur / Numerical Methods for Chemical Engineers / Implicit ODE Solvers
13
Example of the AB/AM Algorithm
 Startup was done with the Implicit Trapezoid rule; h = 0.2
With Convergence Loop
No Convergence Loop
2.5
1
0.9
2
0.8
1.5
0.7
0.6
A, B
A, B
1
0.5
0.5
0.4
0
0.3
-0.5
0.2
-1
-1.5
0.1
0
0.5
1
1.5
2
2.5
t
3
3.5
4
4.5
5
0
0
0.5
1
1.5
2
2.5
t
3
3.5
4
4.5
5
 The convergence loop more than triples the overall time
needed, but it increases stability
Daniel Baur / Numerical Methods for Chemical Engineers / Implicit ODE Solvers
14
Implicit Algorithms and Stiff Problems
 Since implicit algorithms are generally more stable than
explicit algorithms (some are even A-stable!), they fare
much better with stiff problems, where the step size is often
restricted by stability issues
 For non-A-stable implicit algorithms, the main goal is usually to get
the largest possible stability region, since this is the main advantage
of implicit algorithms
 However, this stability comes at the price of larger
computational demand per step, which is needed to solve
the arising algebraic equation systems
 There are however highly specialized algorithms to solve systems
arising from implicit solvers, which can take into account special
features of the problem like sparseness or bandedness
Daniel Baur / Numerical Methods for Chemical Engineers / Implicit ODE Solvers
15
Sparse and Banded Matrices
 It is computationally very advatageous if sparse or in the
best case even banded matrices can be used:
0
0
10
10
20
20
30
30
40
40
50
50
60
60
70
70
80
80
90
90
100
100
0
20
40
60
80
100
0
20
nz = 510
40
60
80
100
nz = 298
Sparse: Storing and operating on only
Banded: All non-zero elements are
510 non-zero elements (times two for
grouped in a band around the diagonal,
which further simplifies computations
Daniel Baur / Numerical Methods for Chemical Engineers / Implicit ODE Solvers
16
Variable Step Size and Order Algorithms
 Since the step size h is of such importance for stability and
accuracy, sophisticated algorithms adjust it during runtime
 This requires error estimation, which is usually done by
comparing the result to what the same algorithm produces
with a higher order
 Some algorithms even adjust their order p dynamically
Daniel Baur / Numerical Methods for Chemical Engineers / Implicit ODE Solvers
17
Assignment 1
1. Implement the Implicit Trapezoid method (see slide 10) for
the batch reactor given on slide 3.
 Simply use left division «\» to solve the arising linear systems.
 What is the maximum step size h that still insures stability?
 What step size h is needed to get a maximum error of 0.1% at all
time points, compared to the analytical solution?
 Plot the solutions.
Daniel Baur / Numerical Methods for Chemical Engineers / Implicit ODE Solvers
18
Exercise 2
 The following reaction scheme is known as the Oregonator
1
A  Y 
 XP
k
X  Y   2 P
k2
A and B are held constant
3
A  X 
 2X  2Z
k
2X
4

 AP
k
P and Q are constantly
withdrawn, therefore ≈ 0
5
B  Z 
 Y
k
Daniel Baur / Numerical Methods for Chemical Engineers / Implicit ODE Solvers
19
Exercise 2 (Continued)
 This leads to the following system of ODEs
dX
dt
dY
dt
dZ
dt
 k1 A Y  k 2 X Y  k 3 A X  2 k 4 X
2
  k1 A Y  k 2 X Y  k 5 B Z
 2k3BX  k5BZ
Daniel Baur / Numerical Methods for Chemical Engineers / Implicit ODE Solvers
20
Assignment 2
1. Solve the system using a suitable built-in Matlab solver
(which one? Note the dynamics!), using the following
parameters:
 A = B = 0.5 = const.;
k1 = 1.34; k2 = 1e9; k3 = 8e3; k4 = 4e7; k5 = 1;
X0 = Y0 = Z0 = 0.2;
tspan = [0, 300];
 Plot the solution. What is special, given the fact that this is a model
for a chemical reaction?
Daniel Baur / Numerical Methods for Chemical Engineers / Implicit ODE Solvers
21
```