ODE solvers

Bader G., Aseher U. (1987) A New Basis Implementation for a Mixed Order Boundary Value ODE Solver, SIAM J. Sci. Stat. Comput. 8 483-500. 

This book describes and gives Fortran subroutines for a wide variety of ODE solvers. More to the point, it gives numerical recipes for practically anything you will ever need to compute. Volume 2 is also available online. It discusses Fortran 90 in the context of parallel computing. C, Pascal, and Basic versions of Volume 1 can be purchased. 

The extension to multiple reactions is done by writing Equation (3.1) (or the more complicated versions of Equation (3.1) that will soon be developed) for each of the N components. The component reaction rates are found from Equation (2.7) in exactly the same ways as in a batch reactor. The result is an initial value problem consisting of N simultaneous, first-order ODEs that can be solved using your favorite ODE solver. The same kind of problem was solved in Chapter 2, but the independent variable is now z rather than t. 

The partial differential equations describing the catalyst particle are discretized with central finite difference formulae with respect to the spatial coordinate [50]. Typically, around 10-20 discretization points are enough for the particle. The ordinary differential equations (ODEs) created are solved with respect to time together with the ODEs of the bulk phase. Since the system is stiff, the computer code of Hindmarsh [51] is used as the ODE solver. In general, the simulations progressed without numerical problems. The final values of the rate constants, along with their temperature dependencies, can be obtained with nonlinear regression analysis. The differential equations were solved in situ with the backward... 

Hindmarsh, A.C. (1983) ODEPACK-A systematized collection of ODE solvers, in Scientific Computing (eds R. Stepleman et al), IMACS/North Holland Publishing Company, pp. 55-64. 

The classic methods use an ODE solver in combination with an optimization algorithm and solve the problem sequentially. This solution strategy is referred to as a sequential solution and optimization approach, since for each iteration the optimization variables are set and then the differential equation constraints are integrated. Though straightforward, this approach is generally inefficient because it requires the accurate solution of the model equations at each iteration within the optimization, even when iterates are far from the final optimal solution. 

Thus, due to limitations on the available computer memory, DNS of homogeneous turbulent reacting flows has been limited to Sc 1 (i.e., gas-phase reactions). Moreover, because explicit ODE solvers (e.g., Runge-Kutta) are usually employed for time stepping, numerical stability puts an upper limit on reaction rate k. Although more complex... 

ODE solver. Relative to non-stiff ODE solvers, stiff ODE solvers typically use implicit methods, which require the numerical inversion of an Ns x Ns Jacobian matrix, and thus are considerably more expensive. In a transported PDF simulation lasting T time units, the composition variables must be updated /Vsm, = T/At 106 times for each notional particle. Since the number of notional particles will be of the order of A p 106, the total number of times that (6.245) must be solved during a transported PDF simulation can be as high as A p x A sim 1012. Thus, the computational cost associated with treating the chemical source term becomes the critical issue when dealing with detailed chemistry. 

Fractional time stepping is widely used in reacting-flow simulations (Boris and Oran 2000) in order to isolate terms in the transport equations so that they can be treated with the most efficient numerical methods. For non-premixed reactions, the fractional-time-stepping approach will yield acceptable accuracy if A t r . Note that since the exact solution to the mixing step is known (see (6.248)), the stiff ODE solver is only needed for (6.249), which, because it can be solved independently for each notional particle, is uncoupled. This fact can be exploited to treat the chemical source term efficiently using chemical lookup tables. 

However, in practice, the integral must be evaluated using a stiff ODE solver or chemical lookup tables (see Section 6.9). Because transported PDF simulations are typically used for reacting flows with complex chemistry, the chemical-reaction step will often dominate the overall computational cost. It is thus important to consider carefully the computational efficiency of the chemical-reaction step when implementing a transported PDF simulation. 

We demonstrate the use of Matlab s numerical integration routines (ODE-solvers) and apply them to a representative collection of interesting mechanisms of increasing complexity, such as an autocatalytic reaction, predator-prey kinetics, oscillating reactions and chaotic systems. This section demonstrates the educational usefulness of data modelling. 

In contrast to our preferred standard mode in this book, we do not develop a Matlab function for the task of numerical integration of the differential equations pertinent to chemical kinetics. While it would be fairly easy to develop basic functions that work reliably and efficiently with most mechanisms, it was decided not to include such functions since Matlab, in its basic edition, supplies a good suite of fully fledged ODE solvers. ODE solvers play a very important role in many applications outside chemistry and thus high level routines are readily available. An important aspect for fast computation is the automatic adjustment of the step-size, depending on the required accuracy. Also, it is important to differentiate between stiff and non-stiff problems. Proper discussion of the difference between the two is clearly outside the scope of this book, however, we indicate the stiffness of problems in a series of examples discussed later. So, instead of developing our own ODE solver in Matlab, we will learn how to use the routines supplied by Matlab. This will be done in a quite extensive series of examples. 

The organisation of the Matlab ODE solvers requires some explanation. For this example, the core is a function, ode autocat.m, that returns the derivatives of the concentrations at any particular time or better, for any set of concentration of the reacting species. Essentially it is the Matlab code for equation (3.83). 

This function is called numerous times from the Matlab ODE solver. In the example it is the ode45 which is the standard Runge-Kutta algorithm. ode45 requires as parameters the file name of the inner function, ode autocat. m, the vector of initial concentrations, cO, the rate constants, k, and the total amount of time for which the reaction should be modelled (20 time units in the example). The solver returns the vector t at which the concentrations were calculated and the concentrations themselves, the matrix C. Note that due to the adaptive step size control, the concentrations are computed at times t which are not predefined. 

We are using the solvers here in their very basic version. Many additional parameters can be controlled, such as maximal step size or required accuracy. We refer to the original documentation for more information about these topics. In the above program, autocat.m, the 20 represents the total time. The ODE solver calculates the optimal step size automatically and returns the time vector t with the concentrations C. The ODE solver can also be forced to return concentrations at specific times by passing the complete vector of times instead of only the total time. 

With the availability of numerical ODE solvers, exercises of the kind just presented are superfluous. While the results are the same for large parts of the data, numerical integration delivers a complete analysis that covers the whole reaction from time 0 to the end. 

We return to spectrophotometric absorbance data and Beer s law Y=CA. The matrix C of concentration profiles, as computed by the relevant function, contains the profiles of all species that are part of the model. Remember, in kinetics this relevant function is an ODE solver, e.g. Runge-Kutta, while for equilibrium investigations, such as titrations, it is the Newton-Raphson routine that computes C. The following considerations hold for both kinetic and titrimetric absorbance data Y. 

Despite its success, the embedded model approach still requires repeated solution of the process model (and sensitivities). For large processes or for processes that require the solution of rigorous underlying procedures, this approach can become expensive. Moreover, for stiff or otherwise difficult systems, this approach is only as reliable as the ODE solver. The embedded model approach also offers only indirect ways of handling time-dependent constraints. Finally, the optimal solution of this approach is only as good as its control variable parameterization, which often can only be improved by a priori information about the specific problem. Consequently, we now consider the simultaneous approach to (16) as an alternative to solution methods for (17). 

For these time periods, the ODEs and active algebraic constraints influence the state and control variables. For these active sets, we therefore need to be able to analyze and implicitly solve the DAE system. To represent the control profiles at the same level of approximation as for the state profiles, approximation and stability properties for DAE (rather than ODE) solvers must be considered. Moreover, the variational conditions for problem (16), with different active constraint sets over time, lead to a multizone set of DAE systems. Consequently, the analogous Kuhn-Tucker conditions from (27) must have stability and approximation properties capable of handling ail of these DAE systems. 

Fig. 15. ODE solver for state differential equations using collocation on finite elements with information processed from element to element.

Even for simple reactions analytical solution is impossible. Instead we use an ODE solver to integrate the equations. The boundary values are the inlet flow and inlet concentrations. 

P. Valkd and S. Vajda, An extended ODE solver for sensitivity calculations, Computers and Chemistry, 8 (1984) 255-271. 

Develop two method-of-lines simulations to solve this problem. In the first, formulate the problem as standard-form ordinary differential equations, y7 = ff(f, y). In the second, formulate the problem in differential-algebraic (DAE) form, 0 = g(t, y, y ). Standard-form stiff, ordinary-differential-equation (ODE) solvers are readily avalaible. DAE solvers are less readily available, but Dassl is a good choice. The Fortran source code for Dassl is available at http //wwwjietlib.org. 

---