Report

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