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: daniel.baur@chem.ethz.ch 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 Ay 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 Ay 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