Ordinary Differential Equations (ODEs)
Download
Report
Transcript 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]
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
The Jacobian matrix reads
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,
their position) instead of 10’000
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
XP
k
X Y 2 P
k2
A and B are held constant
3
A X
2X 2Z
k
2X
4
AP
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