Stiff system of differential equations

A general numerical algorithm of the boundary layer type for stiff systems of differential equations has been proposed by Miranker [173] and applied to a few kinetic problems by Aiken and Lapidus [174,175]. The principle of the method will be briefly described in the case of the following system of differential equations, involving stiff variable x and non-stiff variable y. 

The Runge-Kutta-Fehlberg is a further modification of the Runge-Kutta fourth-order method. It uses a fifth function evaluation to determine the appropriate step size. This method appears to be very efficient for non-stiff systems of differential equations. Additional details regarding this method and a computer listing can be found in a report by Fehlberg and in the chapter by Watt and Shampine. ... 

Methods for Stiff Equations. Stiff systems of differential equations occur commonly with physiologically based models. That is, very fast processes and very... 

One problem with the availability of a large set of standard libraries is that at times it is possible to use the wrong module for a given task. A common situation involves the solution of differential equations. Some systems of differential equations contain derivatives that vary over wide scales, and these are known as stiff systems of differential equations. Therefore, a stiff differential equation solver should be used in these cases otherwise, substantial numerical errors or convergence problems will result. 

In porous reactors several scales appear in a natural way. Firstly, very often slow and fast chemical reactions take place simultaneously this is the temporal aspect. Therefore, stiff systems of differential equations have to be studied theoretically and treated numerically. Secondly, the complicated geometric micro-structure of the pores or capillary tubes and in certain cases also of macropores, cracks, and fractures make it necessary to study chemical processes on micro- and macro-scales simultaneously this is the spatial aspect. Both aspects mentioned may be intertwined when slow and fast processes take place in different domains in space. 

As most of the models in chemical engineering lead to highly nonlinear and stiff systems of differential equations, chemical engineering is a source of challenging numerical problems which cover the whole field of numerical mathematics. 

Thus, complex chemical processes are represented as a number of simple reactions that are very inhomogeneous on a time scale. Generally, it is impossible to separate the fast processes and the slow ones from each other, so that a continuous time monitoring of the total kinetic process is needed to understand the essence of the phenomenon. Mathematical models provide an adequate tool for the scanning of the kinetic curves. Fig. 1(a) shows a typical example of curves where two time scales are present. These time scales differ up to an order of 10 from each other. If one considers the process on the logarithmic scale, then just three different time scales may be identified, see Fig. 1(b). The presence of both fast and slow variables is explained by the occurrence of either large or small factors in the dynamical equations. For example, this is the case for so-called stiff systems of differential equations. 

Some systems may show stiff properties, especially those for oxidations. Here the system of differential equations to be integrated are not stiff . Even at calculated runaway temperature, ordinary integration methods can be used. The reason is that equilibrium seems to moderate the extent of the runaway temperature for the reversible reaction. 

The error in Runge-Kutta calculations depends on h, the step size. In systems of differential equations that are said to be stiff, the value of h must be quite small to attain acceptable accuracy. This slows the calculation intolerably. Stiffness in a set of differential equations arises in general when the time constants vary widely in magnitude for different steps. The complications of stiffness for problems in chemical kinetics were first recognized by Curtiss and Hirschfelder.27... 

We strongly suggest the use of the reduced sensitivity whenever we are dealing with differential equation models. Even if the system of differential equations is non-stiff at the optimum (when k=k ), when the parameters are far from their optimal values, the equations may become stiff temporarily for a few iterations of the Gauss-Newton method. Furthermore, since this transformation also results in better conditioning of the normal equations, we propose its use at all times. This transformation has been implemented in the program for ODE systems provided with this book. 

It is instructive to study a much simpler mathematical equation that exhibits the essential features of boundary-layer behavior. There is a certain analogy between stiffness in initial-value problems and boundary-layer behavior in steady boundary-value problems. Stiffness occurs when a system of differential equations represents coupled phenomena with vastly different characteristic time scales. In the case of boundary layers, the governing equations involve multiple physical phenomena that occur on vastly different length scales. Consider, for example, the following contrived second-order, linear, boundary-value problem ... 

Dassl, solves stiff systems of differential-algebraic equations (DAE) using backward differentiation techniques [13,46]. The solution of coupled parabolic partial differential equations (PDE) by techniques like the method of lines is often formulated as a system of DAEs. It automatically controls integration errors and stability by varying time steps and method order. 

In this sense, differential equations appear more tractable since they do not require particle tracking. Indeed, the solution of the coupled equations of mass, momentum and energy balance including the material equation, properly described on a suitable finite element mesh, theoretically provides the material lines. Nevertheless, the correct description of the basic experiments often requires the use of strong nonlinear terms. Such improvements may be unsatisfying from the numerical point of view since they can lead to stiff systems of nonlinear equations and to many convergence related problems. 

The basic K)ints of this study were 1) the use of the Gear s stiff type method for the direct integration of the system of differential equations without the recourse to the steady state assumption 2) the experimental determination of the rate constants for the most crucial steps (i.e., benzoyloxy radical addition and hydrogen abstraction from the model compound of (III) and hydrr n abstraction from the ethylen-propylene mcaety), when not derivable from the literature data. 

The establishment of a detailed kinetic model provides an opportunity for the numerical prediction of the behaviour of a chemical system under conditions that may not be accessible by experimental means. However, large-scale models with many variables may require considerable computer resource for their implementation, especially under non-isothermal conditions, for which stiffness of the system of differential equations for mass and energy to be integrated is a problem. Computation in a spatial domain, for which partial differential expressions are appropriate, becomes considerably more demanding. There are also many important fluid mechanical problems in reactive systems, the detailed kinetic representation of the chemistry for which would be highly desirable, but cannot yet be computed economically. In such circumstances there is a place for the use of reduced or simplified kinetic models, as discussed in Chapter 7. Thus,... 

The system of differential equations is integrated using CVODE numerical integration package. CVODE is a solver for stiff and nonstiff ordinary differential equation systems [60]. The fraction of dose absorbed is calculated as the sum of all drug amounts crossing the apical membrane as a function of time, divided by the dose, or by the sum of all doses if multiple dosing is used. 

The rigidity of the system of differential equations depends on the equations, their initial conditions, and the nnmerical method. Nonstiff methods can be employed to solve stiff problems, but these require much more computational time. One example is the propagation of flame fronts. The stiffness of the equations has to do with complex chemical process and differences in time scales. 

Luss and Amundson (1968) have studied the dynamics of catalytic fluidized beds. The system is a good example of a stiff set of differential equations. Catalytic fluidized beds are utilized for a variety of reactions such as oxidation of naphthalene and ethylene and the production of alkyl chlorides. A batch fluidization reactor is usually built as a cylindrical shell with a support for the catalyst bed. The reactants enter from the bottom through a cone and cause the catalyst particles to be fluidized in the reactor. The reactants leave through a cyclone in which the entrained solids are separated and returned to the bed. 

Once again, a pseudo-steady-state approximation may be adopted to reduce the stiffness of the system of differential equations for short polymerization times. 

The differential equations which arise in almost all chemical kinetic studies of complex reaction schemes are "stiff differential equations". In chemical kinetics the stiffness is caused by the huge differences in the reaction rate constants of the various elementary reactions. It is impossible to solve such a system of differential equations by the usual Rung-Kutta methods. Therefore we used a program, described by Gear (J ) as a special multi step predictor-corrector method with self adjusting optimum step size control. 

The differential eqns (10) and (11) were solved numerically by the method of Lee (J 4) which is well suited to solve non-linear boundary value problems with non-constant coefficients. However, as will be discussed later it was not possible to obtain convergence for certain parameter combinations, particularly at high values of )c. It is assumed that this has to be attributed to stiffness of the system of differential equations. Since... 

The system of differential equations, provided by Eqs. (9.27), is of the gyroscopic undamped type with time-invariant matrices, see Gasch and Knothe [78]. The mass matrix M" is s unmetric, the g3uoscopic matrix M (12) is antimetric and the stiffness matrix P (12) is S3mimetric as well as positive definite. The s unmetry properties are inherited from the continuous beam description and become obvious, for instance, in Eqs. (8.34) by means of the matrices I, I, I" from Eqs. (7.72). The solution to the problem at hand consists of three parts that are to be determined separately in the following. 

The errors induced within methods based on timescale separations will be discussed in more detail in Sect 7.8 below. On the other hand, since equilibrations will exist within the groups, the introduction of such families is likely to lead to the elimination of fast timescales, thus reducing the stiffness of the reduced system of differential equations with resultant increases in simulation speed. 0 D, for example, has an atmospheric lifetime of the order of 10 s (see Sect. 6.3), and therefore, its presence within a scheme can lead to large stiffness ratios when treated as an individual species. Within reactive flow models, further computational gains may also be made by advecting these families within the transport step rather than individual species, thereby reducing the number of transported variables. 

---