Stiff problems

The additional number of differential equations and increased complexities of the equilibrium relationships may also be compounded by computational problems caused by widely differing magnitudes in the equibbrium constants for the various components. As discussed in Section 3.3.2, it is shown that this can lead to widely differing values in the equation time constants and hence to stiffness problems for the numerical solution. 

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 Runge-Kutta algorithm cannot handle so-called stiff problems. Computation times are astronomical and thus the algorithm is useless, for that class of ordinary differential equations, specialised stiff solvers have been developed. In our context, a system of ODEs sometimes becomes stiff if it comprises very fast and also very slow steps and/or very high and very low concentrations. As a typical example we model an oscillating reaction in The Belousov-Zhabotinsky (BZ) Reaction (p.95). 

The quasi steady state approximation is a powerful method of transforming systems of very stiff differential equations into non-stiff problems. It is the most important, although somewhat contradictive technique in chemical kinetics. Before a general discussion we present an example where the approximation certainly applies. 

As a bit of an aside, one can think of the algebraic constraint as an infinitely stiff problem. Referring to the stiff model problem (Section 15.2), stiff problems are characterized by a fast transient and a slowly varying solution. Regardless of the initial condition, a stiff problem will always decay to the slowly varying solution, and the stiffer the problem, the faster will be the decay (e.g., Fig. 15.1). The situation in a problem like that in Fig. 7.5 is that there is no transient in the y2 component because it is a constraint, and not a differential equation. If, however, the y2 equation is modeled as y 2 = — X(y2 — 1), then y2 = (y2(0) — l)e Xl. As A. becomes larger, the differential equation becomes stiffer, and as X —> oo, the differential equation becomes an algebraic constraint. 

The successful numerical solution of differential equations requires attention to two issues associated with error control through time step selection. One is accuracy and the other is stability. Accuracy requires a time step that is sufficiently small so that the numerical solution is close to the true solution. Numerical methods usually measure the accuracy in terms of the local truncation error, which depends on the details of the method and the time step. Stability requires a time step that is sufficiently small that numerical errors are damped, and not amplified. A problem is called stiff when the time step required to maintain stability is much smaller than that that would be required to deliver accuracy, if stability were not an issue. Generally speaking, implicit methods have very much better stability properties than explicit methods, and thus are much better suited to solving stiff problems. Since most chemical kinetic problems are stiff, implicit methods are usually the method of choice. 

In regions where the solution varies slowly, accuracy considerations alone would permit a large time step. However, for a stiff problem where nearby solutions vary rapidly, stability demands a very small time step, even in regions where the solution is changing slowly and the local truncation error (accuracy) can be controlled easily with a large time step. 

Implicit methods, which have far better stability properties than explicit methods, provide the computational approach to solving stiff problems. The simplest implicit method is the backward (implicit) Euler method, which is stated as... 

It is clear from this expression that the method is stable for all A. and all time steps hn that is, the method is unconditionally stable (for linear problems). A consequence of the strong stability is that the time step can be chosen primarily to maintain accuracy. In the slowly varying regions of stiff problems, the time steps can be very large compared those required to maintain stability for an explicit algorithm. 

When f is nonlinear, as it nearly always is, then an iteration is required to determine y +i. For stiff problems, the iterative solution is usually accomplished with a modified Newton method. We seek the solution yn+i to a nonlinear system that may be stated in residual form as... 

There are many high-quality, well-documented, software packages available to solve stiff problems in this form. However, one often encounters chemically reacting flow problems that are not easily posed as standard-form ODEs. In these cases problems can often be posed easily in a more general form, called differential-algebraic equations (DAE),... 

As an example, if only quasi-steady flow elements are used with volume pressure elements, a model s smallest volume size (for equal flows) will define the timescale of interest. Thus, if the modeler inserts a volume pressure element that has a timescale of one second, the modeler is implying that events which happen on this timescale are important. A set of differential equations and their solution are considered stiff or rigid when the final approach to the steady-state solution is rapid, compared to the entire transient period. In part, numerical aspects of the model will determine this, but also the size of the perturbation will have a significant impact on the stiffness of the problem. It is well known that implicit numerical methods are better suited towards solving a stiff problem. (Note, however, that The Mathwork s software for real-time hardware applications, Real-Time Workshop , requires an explicit method presumably in order to better guarantee consistent solution times.)... 

However, in batch distillation, the system is frequently very stiff, owing either to wide ranges in relative volatilities or large differences in tray and reboiler holdups. Therefore, if methods for non-stiff problems are applied to stiff problems (ODE models but having column holdup and/or energy balances), a very small integration step must be used to ensure that the solution remains stable (Meadow, 1963 Distefano, 1968 Boston et al., 1980 Holland and Liapis 1983, etc.). 

The preceding concept of stability is not sufficient when stiff problems are considered and it is necessary to introduce the concept of absolute stability. A numerical method is said to be absolutely stable if the global discretization error remains bounded for a given step size h when the number, N, of steps tends to infinity. 

Finite difference — Finite difference is an iterative numerical procedure that has been used to quantify current-voltage-time relationships for numerous electrochemical systems whose analyses have resisted analytic solution [i]. There are two generic classes of finite difference analysis 1. explicit finite difference (EFD), where a new set of parameters at t + At is computed based on the known values of the relevant parameters at t and 2. implicit finite difference (IFD), where a new set of parameters at t + At is computed based on the known values of the relevant parameters at t and on the yet-to-be-determined values at t + At. EFD is simple to encode and adequate for the solution of many problems of interest. IFD is somewhat more complicated to encode but the resulting codes are dramatically more efficient and more accurate - IFD is particularly applicable to the solution of stiff problems which involve a wide dynamic range of space scales and/or time scales. 

In this method, instead of working with 4> one works with a new set of variables that are combinations of the fundamental solutions 4> and 3. These new variables all vary with y at the identical rate, removing the stiffness problem. For the Orr-Sommerfeld equation the new variables are (for details about these new variables see Drazin Reid (1981), Sengupta (1992)),... 

The method of characteristics, the distance method of lines (continuous-time discrete-space), and the time method of lines (continuous-space discrete-time) were used to solve the solids stream partial differential equations. Numerical stiffness was not considered a problem for the method of characteristics and time method of lines calculations. For the distance method of lines, a possible numerical stiffness problem was solved by using a simple sifting procedure. A variable-step fifth-order Runge-Kutta-Fehlberg method was used to integrate the differential equations for both the solids and the gas streams. 

For the numerical study of the whole spectrum (for g R fixed), [79] uses a spectral tau-Chebychev discretization in y and the Arnold method (see [88]) to solve the generalized eigenvalue problem (see [89]). This numerical method is based on the orthonormalization of the Krylov space of the iterates of the inverse of the matrix A B. This method has been used more recently in [90]. It has been proven efficient in the stiff problems arising in the study of spectral stability of viscoelastic fluids. 

In ref. 144 the author presents the construction of a non-standard explicit algorithm for initial-value problems. The order of the developed method is two and also is A-stable. The new proposed method is proven to be suitable for solving different kind of initial-value problems such as non-singular problems, singular problems, stiff problems and singularly perturbed problems. Some numerical experiments are considered in order to check the behaviour of the method when applied to a variety of initial-value problems. 

In section 2.2.4, stop condition was used to predict the maximum yield in a chemical reaction. In section 2.2.5, a stiff problem was solved using Maple s default numerical solver. We concluded that the conventional numerical methods fail for this stiff problem. Maple s stiff solver was found to be superior for this stiff problem. Generally, one has to use a stiff solver only if the conventional methods fail, as stiff solvers take more time to solve ordinary IVPs than the conventional solvers. 

In section 5.2.4, a stiff nonlinear PDE was solved using numerical method of lines. This stiff problem was handled by calling Maple s stiff solver. The temperature explodes after a certain time. The numerical method of lines (NMOL) technique was then extended to coupled nonlinear parabolic PDEs in section 5.2.5. By comparing with the analytical solution, we observed that NMOL predicts the behavior accurately. 

Both analytical and numerical methods of lines are presented in this chapter for elliptic partial differential equations. Semianalytical method, presented in this chapter is very powerful technique, and is valid for elliptic Partial differential equations with mixed boundaries also (Subramanian and White, 1999). Numerical method of lines presented in this chapter should be used with precaution, as it may not work for stiff problems. A total of seven examples were presented in this chapter. 

The extension of Gillespie s algorithm to spatially distributed systems is straightforward. A lattice is used to represent binding sites of adsorbates, which correspond to local minima of the potential energy surface. The discrete nature of KMC coupled with possible separation of time scales of various processes could render KMC inefficient. The work of Bortz et al. on the n-fold or continuous time MC CTMC) method can lead to computational speedup of the KMC method, which, however, has been underutilized most probably because of its difficult implementation. This method classifies all atoms in a finite number of classes according to their transition probability. Probabilities are computed a priori and each event is successful, in contrast to the Metropolis method (and other null event algorithms) whose fraction of unsuccessful (null) events increases drastically at low temperatures and for stiff problems. In conjunction with efficient search within a class and dynamic variation of atom coordi-nates, " the CPU time can be practically independent of lattice size. After each event, the time is incremented by a continuous amount. 

---