J. Phys. Chem. 1981, 85, 2861-2862
2861
COMMENTS Comments on a Paper by Edelson and Thomas on Sensltivlty Analysis of Oscillating Reactlons
Sir: The paper in question1is concerned with determining the dependence of the period of an oscillating chemical reaction on parameters such as the rate constants f f k in the kinetic model. The formula presented there is
-=
(4)
(
:)t
where T is the period, [n] is the concentration of any one of the species, and the subscripts on the right indicate the times at which the various derivatives in that formula are to be evaluated. The authors include initial conditions as well as rate constants as parameters in their problem. The derivation of (4), however, does not allow for the fact that, to obtain a strictly periodic solution of the given mass action kinetic system, the initial conditions, or at least a subset of them, cannot be chosen arbitrarily; rather they must lie on the orbit of the periodic solution. It is convenient to distinguish these restricted “initial condition” parameters from the others; we denote them by PI, p2, ... and the others, including the rate constants, by al,a2,... Then for strict periodicity, /3 is a function of a,and the statement that the solution being considered is periodic can be written as [nl(t, @(a))= [nl(t d a ) , a, @(a)) This is analogous to the authors’ (2). Differentiating this with respect to yields
ar _ aak -
(z)t(z)t+v
based on a Floquet-type analysis applied to the linear equations with periodic coefficients satisfied by the various functions appearing in (*). The rate at which the desired constant is approached is essentially the rate at which the orbit of the given periodic solution is approached by solutions close to it. Stiff systems are characterized by a rapid approach, and consequently small delay times before (4)becomes valid. It might be argued that (4)is always valid, on the following basis. For fixed t > 0, the solutions of (1)are little affected by the initial conditions /3 as long as they lie in a small neighborhood of the orbit. Therefore we can suppose that /3 is not constrained by a,and disregard this dependence in the above derivation of (*). This would lead to (4). This reasoning is erroneous; saying that one quantity y depends little on another quantity x because x is restricted to a small range does not mean that dy/dx is small. And (*) shows that derivatives are of primary concern here. To summarize, the validity of the procedure outlined in the cited paper depends on an unstated assumption which should be made clear: the time t in (4) must be large enough. How large it should be depends on the degree of stiffness of the original system (1). Stiffer systems allow smaller values of t. In particular cases, a practical check on the validity can be made while computing the sensitivity functions a[n]/aak;the right-hand side of (4) will be constant when t is in the valid region. Department of Mathematics University of Arizona Tucson, Arizona 8572 1
Paul C. Fife
Received: May 14, 198 1; In Final Form: June 8, 198 1
-
Sensltlvlty Analysis of Oscillating Reactlons. Note on the Transient Effect
-I-
Sir: In our previous paper1 we derived an expression for
( The last term here is not included in the published expression (4). Since it cannot be easily calculated due to the unknown dependence of /3 on a,it is important to know when it can be neglected. One approach which suggests itself immediately is to compute the right side of (4)for t values up until such a time that it becomes independent oft. This constant value will then be the correct value for aT/aak. This approach, in fact, is borne out by theoretical considerations, which indicate that when the given periodic solution is stable, then the last term in (*) approaches 0 as t m, and the first approaches a constant. These considerations are
-
(1) Edelson, D.;Thomas, V. M. J. Phys. Chem. 1981, 85, 1555-8. 0022-3654/a~/20a5-2a81$01.25/0
the sensitivity coefficients of the period of an oscillating reaction with respect to the rate parameters and initial conditions of the model mechanism. These coefficients are given in terms of the sensitivities of one of the oscillating species concentrations at corresponding points on adjacent cycles (ref 1, eq 4). Fife2 has commented that a rigorous treatment of our proposed technique should include a correction term which is required to account for the initial transient induction period of the reactions. His remarks are based on the following considerations: (1) The delay equation (eq 2) in our paper is valid for all t > 0 only if the starting conditions lie exactly on the periodic orbit. (2) If the initial conditions do not lie exactly on this orbit, there is a transient in the oscillating components which corresponds to the path taken to reach the orbit. (1) Edelson, D.; Thomas,V. M. J. Phys. Chem. 1981,85, 1555. (2) Fife, P. C. J. Phys. Chem. Preceding article in this issue.
0 1981 American Chemical Soclety
2882
Additions and Corrections
The Journal of Physical Chemistry, Vol. 85, No. 19, 1981
OREGONATOR
Figure 1. Computer-produced projections of the three-dimensionai limit cycle for the first 75 s (approximately 8 cycles) of the Oregonator of ref 1. The symbols mark the following times (in seconds) on the plot: n, 1 x 0,1 x 10-4; A, 1 x io-? +, 1 x IO-*;X, 1 x io-’; 0 , 1; v, 10.
The duration of this transient is dependent on the initial conditions and its presence could lead to an error in our calculations. In principle we agree with Fife’s analysis. It is necessary to assure that the transient has died out (to within a sufficient numerical accuracy) by the smallest time used
in the sensitivity calculation. Fife presents an expression for the transient correction term, It would be possible to evaluate this, or, as Fife suggests, to make additional calculations of & / d a as a function of time until a constant value is obtained. The former would involve the additional computation of the sensitivity of the oscillating component used for the analysis ([Z] for the Oregonator) to the initial [Y], concentrations of all the oscillating components ([XI, [Z]). Extending the computation of & / d a to longer times would require a deeper probe into the sources of error involved in integrating the Green’s function kernel, as mentioned in our paper.l A more practical approach is to examine the numerical values of the solution vector to ascertain that they recycle accurately at the times used. Figure 1is a plot of the limit cycle of the Oregonator for the rate constants and initial conditions used in ref 1. Despite the apparent large distance between the initial conditions and the limit cycle (greatly distorted on the expanded logarithmic scale) the limit cycle is reached within 1ms of real time, and thereafter is retraced to much better precision (
