Stiff equations model problem

The standard discretization for the equations (9) in molecular dynamics is the (explicit) Verlet method. Stability considerations imply that the Verlet method must be applied with a step-size restriction k < e = j2jK,. Various methods have been suggested to avoid this step-size barrier. The most popular is to replace the stiff spring by a holonomic constraint, as in (4). For our first model problem, this leads to the equations d... 

Stiffness occures in a problem if there are two or more very different time scales on which the dependent variables are changing. Since at least one component of the solution is "fast", a small step size must be selected. There is, however, also a "slow" variable, and the time interval of interest is large, requiring to perform a large number of small steps. Such models are common in many areas, e.g., in chemical reaction kinetics, and solving stiff equations is a challenging problem of scientific computing. 

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 combination of time marching and Newton s method can be illustrated via a very simple model problem [277]. Consider two reactions, R + A B + P and R + B 2P, where in the first a reactant R reacts with a compound A to produce a compound B and a product P. Then, in the second reaction, R further reacts with B to produce two moles of P. If the reaction rates are significantly different, this will lead to a stiff system. For the sake of our example, presume that the mole fraction of R is fixed at a value of 0.1, and that the rate constants for the reactions are k = 1011 and ki = 1012, respectively. Furthermore take the equilibrium constants for the two reactions to be AT] =5 and K.2 = 15. With these parameters set, the mole fractions of A and B (A and B) are governed by the following system of equations. (The value of P is determined from the fact that the mole fractions must sum to unity.)... 

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.)... 

So-called stiff differential equation models are particularly challenging to solve. Stiff models have dynamic behavior that encompasses a wide range of time scales. An example would be fast kinetics combined with long fluid-residence times in a chemical reactor. Gear s method is perhaps the most commonly used technique for solving these types of problems. 

An expanded formulation of the steady-state permeation model has been presented. Two numerical problems - stiffness and an ill-conditioned boundary value problem - are encountered in solving the system equations. These problems can be circumvented by matching forward and reverse integrations at a point near the inlet (n = 0) but outside the combustion zone. The model predicts a... 

On the other hand, these disparate scales often allow us to approximate the complete mechanistic description with simpler rate expressions that retain the essential features of the full problem on the time scale or in the concentration range of interest. Although these approximations were often used in earlier days to allow easier model solution, that is not their primary purpose today. Most models, even stiff differential equation models with fairly disparate time scales, can be solved efficiently with modern ODE solvers. On the other hand, the physical insight provided by these approximations remains valuable. Moreover the reduction of complex mechanisms removes from consideration many parameters that would be difficult to estimate from available data. The next two sections describe, two of the most widely... 

Consider the problem of characterizing the movement of a simple oscillator A small mass m placed at the extreme of a stiff rod of length L oscillates under ideal conditions (friction neither in the mechanical components nor in the movement itself) due to the gravity. If we indicate with F the tension on the rod and with mg the force of gravity, the equations modeling the motion, in Cartesian coordinates, are given by... 

The file ex43.m solves this set of differential equations using ode23s. The model differential equations are defined in m-file mode 143.m. The program ex43a. m solves the problem using the stiff equation solver ode 15. s. 

The principle of the perfectly-mixed stirred tank has been discussed previously in Sec. 1.2.2, and this provides essential building block for modelling applications. In this section, the concept is applied to tank type reactor systems and stagewise mass transfer applications, such that the resulting model equations often appear in the form of linked sets of first-order difference differential equations. Solution by digital simulation works well for small problems, in which the number of equations are relatively small and where the problem is not compounded by stiffness or by the need for iterative procedures. For these reasons, the dynamic modelling of the continuous distillation columns in this section is intended only as a demonstration of method, rather than as a realistic attempt at solution. For the solution of complex distillation problems, the reader is referred to commercial dynamic simulation packages. 

The modeling of complex solids has greatly advanced since the advent, around 1960, of the finite element method [196], Here the material is divided into a number of subdomains, termed elements, with associated nodes. The elements are considered to consist of materials, the constitutive equations of which are well known, and, upon change of the system, the nodes suffer nodal displacements and concomitant generalized nodal forces. The method involves construction of a global stiffness matrix that comprises the contributions from all elements, the relevant boundary conditions and body and thermal forces a typical problem is then to compute the nodal displacements (i. e., the local strains) by solving the system K u = F, where K is the stiffness matrix, u the... 

Dynamic simulations are also possible, and these require solving differential equations, sometimes with algebraic constraints. If some parts of the process change extremely quickly when there is a disturbance, that part of the process may be modeled in the steady state for the disturbance at any instant. Such situations are called stiff, and the methods for them are discussed in Numerical Solution of Ordinary Differential Equations as Initial-Value Problems. It must be realized, though, that a dynamic calculation can also be time-consuming and sometimes the allowable units are lumped-parameter models that are simplifications of the equations used for the steady-state analysis. Thus, as always, the assumptions need to be examined critically before accepting the computer results. 

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 first two sections of Chapter 5 give a practical introduction to dynamic models and their numerical solution. In addition to some classical methods, an efficient procedure is presented for solving systems of stiff differential equations frequently encountered in chemistry and biology. Sensitivity analysis of dynamic models and their reduction based on quasy-steady-state approximation are discussed. The second central problem of this chapter is estimating parameters in ordinary differential equations. An efficient short-cut method designed specifically for PC s is presented and applied to parameter estimation, numerical deconvolution and input determination. Application examples concern enzyme kinetics and pharmacokinetic compartmental modelling. 

Thus due to potential stiffness, integrating a BVP such as the one in equation (5.37) in the positive uj direction from 0 to 1 may not be wise in all cases. In fact, backward integration is much more stable for our simplified model since in backward integration, the eigenvalues switch signs and then the problem is no longer stiff according to the definition. 

The use of the compressibility term can be described as follows. The greater the stiffness a system model has, the more quickly the flow reacts to a change in pressure, and vice versa. For instance, if all fluids in the system are incompressible, and quasi-steady assumptions are used, then a step change to a valve should result in an instantaneous equilibrium of flows and pressures throughout the entire system. This makes for a stiff numerical solution, and is thus computationally intense. This pressure-flow solution technique allows for some compressibility to relax the problem. The equilibrium time of a quasi-steady model can be modified by changing this parameter, for instance this term could be set such that equilibrium occurs after 2 to 3 seconds for the entire model. However, quantitative results less than this timescale would then potentially not be captured accurately. As a final note, this technique can also incorporate flow elements that use the momentum equation (non-quasi-steady), but its strength is more suited by quasi-steady flow assumptions. 

Mathematical models that contain ordinary differential equations face an inherent computational difficulty associated with the stiffness of the equations. Stiffness of ordinary differential equations depends on the relative magnitudes of the response modes or the characteristic time constants of the system being modeled. In solid fuel conversion problems where particles of varying sizes are considered the differential equations for the thermal transients of the particles are usually stiff. Estab-... 

Warner [176] has given a comprehensive discussion of the principal approaches to the solution of stiff differential equations, including a hundred references among the most pertinent books, papers and application packages directed at simulating kinetic models. Emphasis has been put not only on numerical and software problems such as robustness, improving the linear equation solvers, using sparse matrix techniques, etc., but also on the availability of a chemical compiler, i.e. a powerful interface between kineticist and computer. 

---