Transcript ODE_Implicit - The Morbidelli Group

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]
http://www.morbidelli-group.ethz.ch/education/index
Daniel Baur / Numerical Methods for Chemical Engineers / Implicit ODE Solvers
1
Runge-Kutta methods
 If we approximate the slope field more accurately between
two iteration steps, we can generate methods higher local
and global accuracy orders
 Approximating the slope field in m stages:
k1  f  t n , y n 
k 2  f  t n  hc 2 , y n  ha 21 k 1 
c
N odes
k 3  f  t n  hc 3 , y n  h  a 31 k 1  a 32 k 2  
km

 m 1
 f  t n  hc m , y n  h   a m i k i
 i 1

y n 1



i
bi
W eights
a ij
R unge-K utta M atrix
ki
S lopes
m
N um ber of stages
 m

 yn  h   b jk j 
 j 1

Daniel Baur / Numerical Methods for Chemical Engineers / Implicit ODE Solvers
2
Butcher tableau
 A convenient way of of writing down explicit RK methods is
the Butcher tableau
Example: Heun, order 2
0
c2
a 21
c3
a 31
0
1
a 32
1
1 2
cm
a m1
am 2
a m m 1
b1
b2
b m 1
1 2
k1  f  t n , y n 
bm
k 2  f  t n  h , y n  hk 1 
y n 1  y n  h 1 2 k1  1 2 k 2 
Daniel Baur / Numerical Methods for Chemical Engineers / Implicit ODE Solvers
3
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 )
 A system of equations has to be solved at every iteration
step; Why even use implicit algorithms?
Daniel Baur / Numerical Methods for Chemical Engineers / Implicit ODE Solvers
4
Example
 Consider a batch reactor where two reactions take place
dE
E   2I
k1
2
I 
P
k
dt
dI
dt
  k1 E
 2 k1 E  k 2 I
 This can be written in matrix form:
dy
dt
 Ay
  k1
A  
 2 k1
0 

k2 
E 
y  
I
Daniel Baur / Numerical Methods for Chemical Engineers / Implicit ODE Solvers
5
Stability of numerical ODE solvers
 Linear constant coefficient ODE systems such as this are
used as «test-equations» for solvers
dy
 Ay

dt
f  tn , yn   A  yn
 Example: Forward Euler


 I  h  A   y
y n 1  y n  h  f t n , y n  y n  h  A  y n
n
 This iteration will only converge to a value if the spectral
radius of the iteration matrix is smaller than 1:
  I  h A    I  hA 
m ax
 1  h A 
1
m ax
Daniel Baur / Numerical Methods for Chemical Engineers / Implicit ODE Solvers
6
Stability region of the forward Euler algorithm
 The forward Euler algorithm is stable if
1  h A 
1
m ax
 Since eigenvalues are complex in general, we need to use
the definition of the complex absolute value
a  Re  
  a  bi
 
a b
2
b  Im   
2
 Therefore the algorithm is stable if
1  ha 
2
  hb   1
2
Daniel Baur / Numerical Methods for Chemical Engineers / Implicit ODE Solvers
7
Visualizing the stability region
 We can plot the value of our stability equation as a function
of ha and hb
z  ha , hb  
1  ha 
2
  hb 
2
Daniel Baur / Numerical Methods for Chemical Engineers / Implicit ODE Solvers
8
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 )  A y n 1
 We obtain
 I  h A  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
9
Stability of the Backward Euler Method
 If we solve for the next step, we get
y n 1   I  h A 
1
yn
 Again the spectral radius of the iteration matrix determines
convergence:
1
1


  I  hA 

1

 1  h
m ax
 Therefore, the algorithm is stable if
1
1  ha 
2
  hb 
1
2
Daniel Baur / Numerical Methods for Chemical Engineers / Implicit ODE Solvers
10
Visualizing the stability region
1
1  ha 
2
  hb 
1
2
 We see that if a = Re(λ)max < 0, the stability is
independent on the step size
 This property is called A-stability
Daniel Baur / Numerical Methods for Chemical Engineers / Implicit ODE Solvers
11
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
12
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
d y (t )
dt
 A  y (t )  b ,
SR 


m ax
m in
 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
Daniel Baur / Numerical Methods for Chemical Engineers / Explicit ODE Solvers
13
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
14
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,
their position) instead of 10’000
which further simplifies computations
Daniel Baur / Numerical Methods for Chemical Engineers / Implicit ODE Solvers
15
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
16
Assignment
1. Consider the batch reactor given on slide 6. What is the
maximum step size h that would give a stable forward
Euler solution?
 Remember that the eigenvalues of a matrix can be calculated by
solving det(A – λI) = 0 for λ
 The determinant of a 2x2 matrix is
  a11
det  
  a 21
a12  
   a11 a 22  a 21 a12
a 22  
2. What is the maximum step size h that would give a stable
backward Euler solution?
Daniel Baur / Numerical Methods for Chemical Engineers / Implicit ODE Solvers
17
Assignment (continued)
3. From the Butcher tableaus, implement the 2nd order Heun
method and «the» 4th order RK method.
 Solve the batch reactor and Plot the solutions.
 The Butcher tableaus for the methods are:
0
0
Heun 2
1
1
1 2
1 2
1 2
1 2
1 2
0
1 2
1
0
0
1
1 6
1 3
1 3
RK 4
1 6
 Use k1 = 1; k2 = 2; E0 = 1; I0 = 0; tSpan = [0,10];
 The analytical solution reads
E ( t )  exp(  k 1t )
I( t )

2 k1
k 2  k1
 exp(  k 1t )  exp(  k 2 t ) 
 What step size h is needed to get a maximum error of 0.1% at all
time points, compared to the analytical solutions?
Daniel Baur / Numerical Methods for Chemical Engineers / Implicit ODE Solvers
18
Assignment (continued)
4. Change k2 to 10000.
 Do the explicit algorithms still work? Try to solve the system with
ode45. Does it perform well? What could be the problem?
5. Implement the backward Euler algorithm to solve the new
system with k2 = 10000 and plot the solutions.
 Note that you cannot just put the backward Euler formula into Matlab
as is! Use fsolve to solve for yn+1.
 Use ones(size(y0)) as initial guess for fsolve
 You might want to use these to avoid spamming you command line:
options = optimset('display','off');
y(i+1,:) = fsolve( ........, options);
 What can you say about the behavior of the chemical reaction when
you plot it?
Daniel Baur / Numerical Methods for Chemical Engineers / Implicit ODE Solvers
19