Systems of odes

This works rather well with odes, as also seen in Fig. 4.7. For use with pdes, however, it may be considered too much trouble to program, especially as there are easier options, for example extrapolation, which produce results that are just as good. Also, if BDF is nonetheless chosen, it was found in Sect. 4.8.1 and proved mathematically in Appendix B, that the simple start with a simple time correction produces about equally good results for much less effort. 

All the techniques described above can also be applied to the numerical solution of systems of odes, and here we are getting closer to what happens when we solve pdes, because in effect, one reduces them to ode systems when discretising them. 

Instead of a single variable y, there is now a number n of variables, 

In principle, all the methods described above for single odes can be used for the solution of such a system, when extended suitably. In the case of explicit methods such as Euler or RK, this is very simple to implement, whereas with implicit methods such as BI or the trapezium method, there are some choices to be made. 

For brevity, the Euler method will be treated as a special case of RK, considered as RK1. The method is then to start by calculating a vector of k values, one for each y element. Discretising directly from (4.51), this is 

Systems of coupled ODEs might arise from reformulating a higher-order differential equation to a system of first-order differential equations, or as a description of a system that consists of coupled variables. Systems of differential equations can be solved as an extention of the methodology for a single differential equation. The principle is shown in Example 6.4, which considers a stirred tank reactor. 

V is the reactor volume, UA T — T ) is the external cooling for controlling the reaction temperature, and ri(—A//)Kis the heat caused by chemical reaction. As these equations are coupled, they must be advanced simultaneously. For simplicity, and for pedagogical reasons, we will use the explicit Euler method to solve this problem  

Solving Equations (6.31) and (6.32) obviously requires that the slopes /i(Ca, , T ) and /2(Ca, , Tn) are calculated before Ca and T are updated. Using a fourth-order Runge-Kutta algorithm means that eight slopes need to be calculated in each step. 

The methods described earlier in this chapter apply to higher-order differential equations, 

---

It is convenient to analyse tliese rate equations from a dynamical systems point of view similar to tliat used in classical mechanics where one follows tire trajectories of particles in phase space. For tire chemical rate law (C3.6.2) tire phase space , conventionally denoted by F, is -dimensional and tire chemical concentrations, CpC2,- are taken as ortliogonal coordinates of F, ratlier tlian tire particle positions and velocities used as tire coordinates in mechanics. In analogy to classical mechanical systems, as tire concentrations evolve in time tliey will trace out a trajectory in F. Since tire velocity functions in tire system of ODEs (C3.6.2) do not depend explicitly on time, a given initial condition in F will always produce tire same trajectory. The vector R of velocity functions in (C3.6.2) defines a phase-space (or trajectory) flow and in it is often convenient to tliink of tliese ODEs as describing tire motion of a fluid in F with velocity field/ (c p). 

This statement is not exactly true - the slightly different system of ODEs is defined by an asymptotic expansion in powers of At which is generally divergent. 

If for example we discretize the region over which the PDE is to be solved into M grid blocks, use of finite differences (or any other discretization scheme) to approximate the spatial derivatives in Equation 10.1 yields the following system of ODEs ... 

In summary, in the equilibrium-chemistry limit, the computational problem associated with turbulent reacting flows is greatly simplified by employing the presumed mixture-fraction PDF method. Indeed, because the chemical source term usually leads to a stiff system of ODEs (see (5.151)) that are solved off-line, the equilibrium-chemistry limit significantly reduces the computational load needed to solve a turbulent-reacting-flow problem. In a CFD code, a second-order transport model for inert scalars such as those discussed in Chapter 3 is utilized to find ( ) and and the equifibrium com-... 

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 total concentration of catalyst is much smaller than the concentrations of the reactants or products. Note, that in real systems, the reactions are reversible and usually there are more intermediates, but for the present purpose, this minimal reaction mechanism is sufficient. The system of ODEs ... 

Let me make my own personal preference clear from the outset. I have solved literally hundreds of systems of ODEs for chemical engineering systems over my 30 years of experience, and t have found only one or two situations where the plain old simple-minded first-order Euler algorithm was not the best choice for the problem. We will show some comparisons of different types of algorithms on different problems in this chapter and the next. 

Burrage, K., and Petzold, L. R., On order reduction for Runge-Kutta methods applied to differential/algebraic systems and to stiff systems of ODEs," Lawrence Livermore National Laboratory, UCR-98046 preprint (1988). 

Once the PDE has been semi-discretized (i.e., discretize the spatial derivatives but not the timelike derivatives) to form a system of ODEs, the ODEs can be solved by high-level software packages. In the standard form there are many such packages available, with relatively fewer for DAEs (see Section 15.3.3). In the method of lines, the spatial differencing must be done by the user, while time discretization and error control is handled by the ODE software. Overall, the effort to develop a new simulation is reduced, since a good deal of existing high-level software can be used. 

In many cases ordinary differential equations (ODEs) provide adequate models of chemical reactors. When partial differential equations become necessary, their discretization will again lead to large systems of ODEs. Numerical methods for the location, continuation and stability analysis of periodic and quasi-periodic trajectories of systems of coupled nonlinear ODEs (both autonomous and nonautonomous) are extensively used in this work. We are not concerned with the numerical description of deterministic chaotic trajectories where they occur, we have merely inferred them from bifurcation sequences known to lead to deterministic chaos. Extensive literature, as well as a wide choice of algorithms, is available for the numerical analysis of periodic trajectories (Keller, 1976,1977 Curry, 1979 Doedel, 1981 Seydel, 1981 Schwartz, 1983 Kubicek and Hlavacek, 1983 Aluko and Chang, 1984). 

A quasi-periodic solution to a system of ODEs is characterized by at least two frequencies that are incommensurate (their ratio is an irrational number) (Bohr, 1947 Besicovitch, 1954). Several such frequencies may be present on high-order tori, but for the two-dimensional forced systems we examine, we may have no more than two distinct frequencies (a two-torus, T2). A quasi-periodic solution is typically bora when a pair of complex conjugate FMs of a periodic trajectory leave the unit circle at some angle , where /2ir is irrational. Such a solution is also expected when we periodically perturb an autonomously oscillating system with a frequency incommensurate to its natural frequency. 

As this book is mainly intended for undergraduates, we only treat those chemical/biological processes and problems that can be modeled with ODEs. PDEs are significantly more complicated to understand and solve. However, we will often have to solve systems of ODEs, rather than one single ODE. Such ODE systems contain several DEs in one and the same independent variable, but they generally involve several functions or state variables in their formulation. In fact, systems of ODEs (and matrix DEs) occur quite naturally in chemical/biological engineering problems. 

The most simple isothermal case was studied first by finding the analytical solution to the BVP. Possible numerical difficulties of this case were analyzed via the eigenvalues of the corresponding constant coefficient linear system of ODEs in order to make our readers aware of the potential harm of stiffness when trying to solve DEs. The isothermal case with nonlinear kinetics was later solved numerically in a MAT-LAB program. 

This chapter has reviewed existing results in addressing the analysis and control of multiple-time-scale systems, modeled by singularly perturbed systems of ODEs. Several important concepts were introduced, amongst which the classification of perturbations to ODE systems into regular and singular, with the latter subdivided into standard and nonstandard forms. In each case, we discussed the derivation of reduced-order representations for the fast dynamics (in a newly defined stretched time scale, or boundary layer) and the corresponding equilibrium manifold, and for the slow dynamics. Illustrative examples were provided in each case. 

The process response is presented in Figure 4.6. Observe that all the state variables exhibit a fast transient, followed by a slow approach to steady state, which is indicative of the two-time-scale behavior of the system, and is consistent with our observation that processes with impurities and purge are modeled by systems of ODEs that are in a nonstandard singularly perturbed form. 

Equation (7.17) is a description of the fast dynamics of the high-purity distillation column. It involves only the stage temperatures and it can be easily verified that the system of ODEs describing the fast dynamics (as well as the quasi-steady-state conditions that result from setting the left-hand side of (7.17) to zero) are linearly independent. The constraints arising from the fast dynamics can therefore be solved (typically numerically) for the quasi-steady-state values of the stage temperatures, T = [7 (. 7. .. Tp Tg], which can then be substituted into the ODE system (7.8) in order to obtain a description of the dynamics after the fast temperature transient ... 

The initialisation of variables in a system of equations is important. While, in systems of ODEs all of the state variables must be initialised, in DAE systems only some of the variables need to be initialised, which is the same as the number of differential variables for index one system. The other variables can be determined using the algebraic equations. It is inconvenient for the user to be required to initialise all of the variables as this might require the solution of a set of nonlinear algebraic equations. Pantelides (1988) developed a procedure for consistent initialisation of DAE systems. Readers are directed to this reference for further details. 

We begin with single odes. At the end of this chapter, systems of odes are dealt with, as they are in fact one way of handling pdes, using the Method of Lines (MOL, see Chap. 9), which has a system of odes as an intermediate stage, or something close to it. 

This section is left to the last because it pertains to systems of odes (or DAEs), although Rosenbrock methods are a kind of Runge-Kutta method... 

Before moving on to real Rosenbrock methods, consider again (4.66). The left-hand side contains a term in fy if we are dealing with a system of odes, this is called the Jacobian of the system. It is often constant, evaluable in advance. It will be seen in Chap. 9 that unless the diffusion problem has nonlinear concentration terms (for example from higher-order homogeneous reactions), the Jacobian is constant. If not, it must be evaluated at every step. 

This will be the basis for what follows. The diffusion equation, discretised on the right-hand side as in (8.8), is now a system of odes in the concentration vector C, of the form... 

The two methods are BDF and extrapolation. Both methods are used for the numerical solution of odes and are described in Chap. 4. The extension to the solution of pdes is most easily understood if the pde is semidiscretised that is, if we only discretise the right-hand side of the diffusion equation, thus producing a set of odes. This is the Method of Lines or MOL. Once we have such a set, as seen in (8.9), the methods for systems of odes can be applied, after adding boundary conditions. 

Yet another, quite different, approach to solving a system of odes, such as one obtains as an intermediate step when using, for example, MOL or OC, is the eigenvalue-eigenvector method. Its use for electrochemical simulations was described in two papers in 1989 and 1990 [255,332]. The method has some drawbacks, and does not appear to have seen much use since these two papers. It does have one unique feature there is no discretisation of time. A solution is generated by the algorithm, at any chosen time. So, although the method may at times be fairly inefficient, if one wants a current or concentrations at only one or a few time points, this could be faster than a time march with the usually small time intervals. 

The method starts with a system of odes, represented as in (9.54),... 

Equation (14.6), when discretised for, say, a potential jump experiment (Cottrell), gives rise to a system of odes, depending on the method of discretisation used. For example, using the EX method, we have for the v fh equation... 

---