Introduction
Strange attractor, deterministic chaos, self organization—interpretation of these concepts, begun relatively recently, is still under way [17]. Lorenz's discovery of the extraordinary behavior of the numerical solution of a system of ordinary differential equations [8] marks the dawn of the era of chaotic dynamics, a field which is the product of classical physics and computeraided mathematics. Chaotic dynamics has very substantially broadened our conceptions of space and time, determinism and randomness.
By chaotic behavior is meant a substantial, arbitrarily large change in the character of the behavior of a system as a result of an arbitrarily small change in the external parameters or initial conditions. Such a definition makes it problematic to use the apparatus of differential calculus to describe chaotic motion mathematically. This is because the requirement that the solutions of differential equations must be continuous and robust leads to the condition that only infinitesimal dispersal of trajectories can result from a small change in the initial conditions of any parameter. The appearance of computerized calculations on the scene, while significantly extending the possible range of numerical integration, cast doubt on the robustness of iterative methods and raised the question of the finite accuracy of computer calculations. Whether Lorenz's attractor is a product of the finite accuracy of calculations on a computer, or whether it represents a new class of solutions of differential equations, is a question which, in our view, calls for analysis. We are faced here with a special situation—the traditional chain linking a system of differential equations with an analytic solution cannot be automatically applied to cases where there is no analytic solution. The finite accuracy of calculations, which from the definition of chaotic behavior cannot fail to influence the solution, has severed this link and made of the computer an essential element of the chain: mathematical model—computer (numerical method, accuracy)—solution and analysis of solution. This problem leads us to the paradoxes of the ancient Greeks: is it mathematically and logically possible to formulate a contradictionfree description of the process of approach to an object if the distance to this object contains an infinite number of segments? Newton's answer, the calculus of infinitesimals—the theory of differential equations, to date the only mathematical apparatus for describing dynamics—will seemingly have to be modified before it can be used to describe the nonlinear dynamics of systems with chaotic behavior. The problem of the existence of a lim Dx/Dt for Dt ® 0 is particularly crucial in the study of chaos. Is the language of differential calculus which appears in the mathematical formulation of all dynamical laws truly so universal? Or is it possible to construct a mathematical tool or even a new calculus for describing dynamics without using derivatives and free of the contradiction between the continuous time and space of differential equations and the discrete process of calculation? Equally urgently, how is one to define the role and place of chaotic behavior in the deterministic world of Newton and the probabilisticstatistical world of Boltzmann?
Below we present a new theory and concordant mathematical model for describing the dynamics of complex chemical transformations. As it is algebraic and iterative, it is free of the contradictions mentioned above. It describes all known types of dynamic behavior in chemical reactions, including chaotic behavior.
The p theorem of the theory of similarity and dimensionality and the stoichiometry of chemical reactions
From the point of view of the theory of dimensionality, all physical quantities fall into two caregories, primary and secondary. Primary quantities are those that are determined by direct measurement, yielding a numerical value: u = V/c_{o} , where V is the primary quantity and c_{o} is the standard unit of measure. Take the average height of a man, V = 175 cm; the numerical value is obviously u = 175 (c_{o} = 1 cm) and depends on the choice of unit. On the other hand, the ratio of similar quantities, i.e., those having the same dimensionality, is independent of the unit. Thus the ratio of the numerical value of a man's height to the length of his arm does not depend on the choice of unit of measurement. This concept of direct measurement is called the "absoluteness of relation".
Secondary physical quantities are expressed in terms of primary quantities. To properly construct a secondary quantity, one has to determine what primary quantities it depends on and how (form of the dependence). The equations with whose help secondary quantities are expressed in terms of primary ones are called determination equations.
It is natural to demand that the concept of absolute "absoluteness of relation" ratio also hold for secondary quantities. Hence the determination equations should have the wholly concrete form:
where u1, u2, .... uM are primary quantities, X is a secondary quantity and f is a homogeneous function. The property of homogeneity is defined by the condition that any transformation of the variables will entrain a similar transformation of the function. In other words, if we change the system of units (multiply u_{j} by some arbitrary coefficients 1/t_{i}), then X will also change, as a result of the homogeneity of ƒ:
The condition of homogeneity of function f leads to the function f(u) :
where C is a constant. Thus to obtain a new physical characteristic (secondary quantity), one need only multiply the primary quantities with each other and raise them to some power.
Having defined the concepts of primary and secondary quantities characterizing the process under investigation, we will now look for relations between these quantities which would ensure invariance under transformations of primary quantities of the type:
where t_{j} are arbitrary coefficients.
The solution of this problem is given by the p theorem of the theory of similarity and dimensionality, and in Brandt's formulation [9] it is as follows. Consider the system of secondary quantities:
where  a_{ij}  is a deterministic matrix of dimension N x M of rank r. Let the function F(X) be homogeneous with respect to the quantities u_{j} and X_{i}. Then the equation
is equivalent to the following
where the first r arguments are units and p1, p2,...pL (L = NM) are L independent products formed from X1, X2,...XN. The theorem has the following consequence: let the quantities X1, X2,...XN have a determinant matrix  a_{ij}  of rank M:
where R is a nondegenerate matrix of dimension MxM , and Q is a matrix of dimension (NM)xM. Then if F is a homogeneous function, Eq. 4 is equivalent to Eq. 5, and
are NM independent quantities with the deterministic matrix  u_{li} :
where I is the unit matrix (NM)x(NM). From Eq. 6 and Eq. 8,
which expresses the condition of independence of pl from uj.
Thus the p theorem enables us to describe systems using the quantities pl , which are invariant under transformations of the primary quantities uj, while formula (8) makes it possible to construct invariant complexes of p_{l} from the quantities X_{i} on the given deterministic matrix  a_{ij} .
The invariance of the complexes of p_{1} reflects the following fact about modeling the behavior of systems: since the choice of primary quantities is arbitrary, for any giset of the quantities u_{j} the system of secondary quantities X_{i} will depend on the values of the quantities u_{j}; this is not a reflection of the essence of the phenomena taking place in the system but rather a consequence of the definite dependence of the primary and secondary quantities expressed by Eq. 3. To construct quantities that do not depend on the choice of the quantities u_{j}, one must construct quantities characterizing the internal properties of the system under investigation. These quantities p_{1}, being criteria of similarity, will be the same for similar systems.
Let us illustrate this
with a simple example. In describing the kinematics of motion, let us specify
a set of four quantities: the distance S, time T, velocity V and acceleration
W. The determinant matrix 
aij 
of the dimensionalities of these quantities is:
s[m]  t[c]  

S  1  0 

= 4 
T  0  1  M  = 2 
V  1  1  L  = 2 
W  1  2  L  = 2 
(11)
The meaning of these expressions is obvious and is a consequence of the very general postulates of the theory.
Now that we have formulated the principal results of the theory of similarity and dimensionality, let us attempt to extend them to analysis of the chemical transformations. Take M substances in the known initial amounts u^{0}_{1} , u^{0}_{2} ,...u^{0}_{M} . In the terminology of the theory of dimensionality, the concentrations corresponding to these chemical substances or elements will be "primary" quantities. Having mixed these substances and assuming their chemical transformation—accompanied by the formation of a series of new substances—we must now introduce the concept of "secondary" quantities formed as a result of the reaction with the concentrations X_{i}. The determination matrix for the substances formed as a result of the reaction will look like this:
u_{1}  u_{2}  . . .  u_{M}  

X_{1}  a_{11}  a_{12}  . . .  a_{1M} 
X_{2}  a_{21}  a_{22}  . . .  a_{2M} 
:  :  :  :  
X_{N}  a_{N1}  a_{N2}  . . .  a_{NM} 
H  I  

X_{1}  H  1  0  N=5  
X_{2}  I  0  1  M=2  
X_{3}  HI  1  1  L=3  
X_{4}  H_{4}  2  0  
X_{5}  I_{4}  0  2 
1 1 1 0 0   
uli   =   2 0 0 1 0   
0 2 0 0 1  
(15)
The requirement that the quantities p_{1} be invariant (be independent of the initial concentrations of the quantities uj ) is fulfilled, since the matrix product Eq. 9 is equal to zero.
Now it is an easy matter to recognize in the quantities pl the equilibrium constants of the corresponding chemical reactions:
H + I  =  HI 
2I  =  I_{2} 
2H  =  H_{2} 
N  
S  u_{li} A_{i}= 0 
i=1 
N  
S  a_{ij} X_{i}= b_{j} 
i=1 
In deriving the relations within the bounds of the theory of similitude and dimensionality, we used notions such as "absoluteness of relation ratio" of quantities and p theorem of the theory of similarity and dimensionality. The fact that Eq. 7 is identical in form with the mass action law, which is obtained in chemical thermodynamics from the minimum free energy condition, supports the analogy of this law with the relations obtained in the theory of similarity and dimensionality. Unlike dimensionless invariants Eq. 12, chemical invariants possess dimensionality, hence in this case we are speaking of a purely formal analogy. However, if we assume that invariance—in the sense of independence from the concentrations of the components in the chemical systems—is just as fundamental a property of physicochemical systems as is nondimensionality, we arrive at the following ideas:
1. The form of Eq. 7 must be conserved whether one is describing reactions in the ideal gas phase or in any condensed phase (liquid, solid, plasma). In order to adequately describe equilibrium states in these phases, it is necessary and sufficient to write down the entire set of physicochemical reactions (constituents). These descriptions must include defects, clusters , structural characteristics etc., in the form of a system of reactions (Eq. 16). Instances of the use of this approach are known and have been successful [11, 12].
2. The form of Eq. 7 must be conserved not merely in the equilibrium state but also at any arbitrary instant. It would in fact be strange if the requirement we used regarding independence from the initial concentrations had to be obeyed at some definite instant t ® ¥ in the equilibrium state. This condition must be met independently of time; that is, according to the logic of the theory of similarity and dimensionality, it is necessary that our "invariants" be independent of time, just as they are independent of the initial concentrations of the components. For this condition to be realized, it is necessary to introduce time—or a function of time—among the "elements" and obtain a system of equations of the type of Eq. 7 in which time appears—that is, use equation Eq. 7 to describe the kinetics of chemical transformations.
A new model for the kinetics of chemical transformations
Eqs. 7
and 17 form a system of N nonlinear algebraic
equations having a unique solution for X_{i}
> 0 and defining a complex chemical equilibrium. This system has been thoroughly
studied and admits of reduction of dimensionality. Eqs. 7
and 17 reduce to the system of M equations [13]:
(18)
Since usually M << N, the computational advantages of (18) are obvious.
As concluded in the preceding section, when formulating p _{l}, when X_{i} has time dependence, one must include an additional element—the time function q (t). Let us take Eq. 13 and introduce q (t) among the components I and H of the chemical system. Then the matrix  aij  becomes
H  I  q (t)  

X_{1}  H  1  0  0 
X _{2}  I  0  1  0 
X _{3}  q (t)  0  0  1 
X _{4}  HI  1  1  g _{1} 
X _{5}  H _{2}  2  0  g _{2} 
X _{6}  I _{2}  0  2  g _{3} 
  1 1 g_{1} 1 0 0    
 u_{li}  =    2 0 g_{2} 0 1 0   
  0 2 g_{3} 0 0 1   
Like q
(t), the concentrations X_{i}
are functions of the time t. Rewriting Eq.
20 in the form
(21)
several properties of the function q (t) for closed chemical systems become apparent. Assuming that the initial values of the concentrations are X_{1} (0) = b_{1}, X_{2} (0) = b_{2}, X_{3} (0) = 0, X_{4} (0) = 0 and X_{5} (0) =0, function q (t) should satisfy the conditions qg (0) = 0; qg (¥) = 1. Based on these restrictions, the following function was proposed for q (t) [14]:
(22)
the being empirical parameters characterizing the rate of the lth reaction.
Thus assuming that the
parameters p _{l},
a_{l}
and b_{j} are given, to obtain the kinetic
curves x_{i} (t) of the system of a complex
reaction specified by the stoichiometric matrix 
u_{li}
, it is
necessary to solve the following system of equations (or the equivalent
system (Eqs. 18)):
(23)
for all values of t_{s},
s = 0, 1, 2... for the specified initial conditions
(24)
The procedure for determining the parameters p _{l} and a_{l} from a statistical analysis of kinetic experiments is given in [15].
It is interesting to compare the solutions obtained for diverse chemical reaction mechanisms (u_{li} ) with the corresponding solutions obtained by integrating the system of differential equations of the kinetic mass action law. These calculations have been carried out and indicate good qualitative correspondence between the solutions for all types of reaction in closed systems [16].
Additional independent
support for this theory is provided by the simple time dependence proposed
by M. Garfinkle in 1982 for the chemical affinity A. After analyzing experimental
data on the kinetics of over 100 reactions in closed systems, he demonstrated
that all can be described satisfactorily by introducing a dependence for
A [17, 18]:
(25)
A_{r}, t_{k} being empirical parameters. It is easy to show that the function q (t) is related to the chemical affinity A by a simple expression, ln q (t) = A (t). The results of analysis of real experimental kinetic data, as well as the ideas we have formulated regarding "chemical" invariance, support the validity of the ideas put forward here for closed chemical systems.
The dynamics of chemical reactions
In the preceding section we constructed a new mathematical model for describing the kinetics of chemical reactions and used it successfully to describe kinetics in closed systems.
However, problems arise
when we seek to use the new kinetic model (NKM) to describe chemical oscillation
reactions. As is well known, the kinetic mass action law (KMAL) equations
leading to oscillations are obtained when the reaction mechanism is made
to include reciprocal relations (autocatalysis, or mutual influence of
the concentrations of the component parts of the system on the acceleration
or deceleration of a given reaction). Take for example:
(26)
This means that [C] accelerates the process of its own formation.
How should the equations of the new kinetic model be changed to allow description of processes of the type illustrated by Eq. 26?
The answer lies in the
method of construction of chemical invariants p_{l}.
Every time we turn our attention to a new factor (new chemical element
in the system, time t). we must change the determinant matrix 
a_{ij}
so as to assure that the invariants p,
subsequently constructed in accordance with Eq. 8,
will not depend on this factor, which appears here as a new variable. As
we did for t, for which we introduced the function q
(t_{s}) in 
a_{ij},
we should introduce a new function of the concentrations Xi,
Y(X_{i}),
which exert a definite influence on the course of the given reaction ([C]
in Eq. 26), as functions of the times t_{sp}
(henceforth written as (sp)):
A  B  q(t_{s})  y(C (t_{sp}))  

X_{1}  A  1  0  0  0 
X_{2}  B  0  1  0  0 
q  0  0  1  0  
y  0  0  0  1  
X_{3}  C  1  1  g  b 
Taking into account Eqs.
8 and 22, the determinant
equation for the invariant p
of Eq. 26 is:
(28)
Determining the explicit
form of the Y
function is the present task facing the theory. Here we will limit ourselves
to the simple empirical dependence which, for any arbitrary mechanism of
reaction, can be written in the form:
(29)
d_{i l} being empirical parameters characterizing the degree of influence of x_{i} on the lth reaction.
Generalizing
(30)
Let us confine ourselves
to the special case of the function x_{i}
( t_{sp}),
p = s1, i.e.,(to instants of time precedings. Clearly, for ts
®¥
Eqs. 30 also determine the chemical equilibrium
of a closed chemical system. An open system is characterized by a time
of establishment of the stationary state (T). Processes in open systems
can be modeled using Eqs. 30 by modifying 0 <
t_{s} < T = constant. Our introduction
of an explicit function of time, q,
discrete time Y(sp),
has opened up the possibility of using a system of algebraic equations
to describe the kinetics of chemical reactions. By using the special delayed
function Y,
it is possible to employ Eqs. 30 to simulate the
complex (and also chaotic) dynamics of chemical reactions. Since Eqs.
30 represents an inexplicit system of multidimensional iterations,
for a given chemical reaction mechanism (matrix u_{li}
) and taking into account the reciprocal influence of the component parts
of the system (Y
function), Eqs. 30 fully determines X_{i}(t_{s}).
A particular case of Eqs. 30 is the following
representation in the form of a system of multidimensional iterations:
s
(31)
The next step in the development
of the new dynamic model is to apply it to describe the spatial behavior
of chemical reactions. In order to be able to apply Eqs.
30 to the construction of the spatial distribution X_{i}
(r_{g}), one must add to the righthand side
of Eqs. 30 the function f
(X_{i} (r_{gf}))
of all the concentrations X_{i} (r_{gf}),
calculated at points in space r_{f} in the
vicinity of rg. By analogy with the Y
function, function f
must be introduced by means of an additional column in matrix  a_{ij}
 . In order to do this the entire spatial region in which the chemical
reaction takes place is divided into cells with linear dimensions D
r_{g} (g = 1,2...). Assuming that the chemical
reaction in a cell with coordinates r is influenced by the concentrations
X_{i} (r_{g1})
from the cells nearest to r and by analogy with Eqs.
29, we write
(32)
where e_{il} are empirical parameters characterizing the influence of the concentrations in the neighborhood of x_{i} (r_{g}).
To illustrate the use of
the formalism set forth above, let us consider the chemical reaction dynamics
represented by the following mechanism:
A
® B
® C
(33)
Element A transforms into
B and B into C, while C accelerates the rate of reaction A B and B slows
its own transformation into C. In accordance with the formalism outlined
above (29, 30), the corresponding
equations describing the dynamics of this complex reaction become
(34)
X_{A,B,C} being collective measure of the constituents A,B, C in the system (particle number, concentration). Let us confine ourselves to the influence of X_{B} and X_{C} on the process at the instants s1 preceding the time t_{s}. Figs. 14 show different kinds of oscillatory solutions of Eqs. 34. Of particular interest to us is the quasicontinuous solution (Fig. 2) we discovered for t_{s} = T (open system). As it is discrete (does not have a derivative at any point), this solution imitates oscillatory processes. As is well known, obtaining oscillatory regimes of this kind has hitherto been the prerogative of the differential equations and special functions associated with them. In our view, the construction of an oscillatory process using a system of iteration that lacks a continuous variable (time, coordinate) points to the possibility of creating a discrete wave mechanics with applications in various fields of science. The model given by Eqs. 34 should make it possible to obtain diverse types of chaotic behavior, including different regimes which could be compared as required with the experimental data [2023]. We will not address this task here and instead will confine ourselves to a single example: the experimental bifurcation diagram obtained in [19] and that calculated by us using (34) (see Fig. 6). While not precisely identical, they are similar enough to encourage us in our hopes of applying our approach to the quantitative description of the BelousovZhabotinsky reaction.
To obtain the spatial behavior of a chemical reaction, according to Eqs. 32 we must treat the concentrations of the component parts of the system as functions of the spatial coordinates r. The corresponding equations and solutions, in the form of ring and spiral waves characteristic of the BZ reaction, are described in [24] and shown here in Fig. 5.
Conclusions
Above we formulated a new mathematical model of the dynamics of physicochemical transformations. However, it is natural to question the necessity of such a model: why is it that the equations of the kinetic mass action law, which are classics, become problematic when it comes to describing complex chaotic systems? There are several reasons for this:
1. The fact that computers are actively used to generate chaotic regimes characterized by high sensitivity to numerical procedure and finite computational accuracy makes it impossible to carry out constructive analysis and interpretation of the solutions obtained. The contradiction between the continuous spacetime variables of differential equations and the discrete method of calculation with computers makes it necessary to perform special analysis, which does not overcome this contradiction. The new dynamic model may be regarded as a contradictionfree consequence of such analysis.
2. From their inception, modern theoretical physics and chemistry have naturally maintained that systems have small reactions to small perturbations. Perturbation theory, the distinction between rapid and slow variables, the various simplification procedures and similar problems—all traditionally use the concept of a "small" system parameter with negligible effect. These conditions cannot be fulfilled when one turns to the description of systems with chaotic behavior. There is no such thing as a "smaller" parameter in such systems. An infinitesimal change in the initial conditions or parameters can lead to the exponential dispersal of trajectories. What earlier appeared to be an instability or peculiarity of solution from the viewpoint of differential calculation, assumes an entirely natural, regular character in the study of chaotic systems. This leads to the necessity of constructing a new calculus free of the above contradictions for systems with chaotic behavior.
Without claiming universality, we see the foregoing presentation as an attempt to construct an alternative language for describing the dynamics of systems of chemical reactions with complex behavior—including chaotic behavior. To be accepted our theory will have to change our way of conceiving time. Our model has a minimum of two types of time: t_{s} and s  p (p = s1, s2). ts could compare with the astronomic time traditionally used in differential equations. In principle this time is not very different from continuous time. But we obtain new time scales in our model when we introduce "internal" clocks by using the s  p iteration. The oscillations we obtained in Fig. 2 actually have a period which depends only on the mechanism of the interactions and on some parameters, unrelated to traditional time, frequency, etc. The system can only be characterized by its mechanism—its temporal behavior must be calculated. We are led to the same conclusion when we look at the spatial behavior: the geometry of the waves depends only on the mechanism involved. We obtained ring waves even when we started from a square. This means that the spatiotemporal behavior is a result of the type of nonlinear interactions present, and can be expressed by an algebraic iterational model.
Of course we are only beginning to understand this new conception of space and time, but it looks very attractive, especially when we compare our results with the experimental data.
To conclude this article, it is necessary to remark that the results reported here should lead to fundamental changes in our conception of randomness in physics. In reality the chaotic behavior we obtained using iterations is nothing more than the "random numbers generator" well known to all programmers engaged in simulations of stochastic processes. This means that if we succeed in constructing theories like the above (namely, a dynamics that leads to difference rather than to differential equations), it will no longer be necessary to take into consideration probabilistic concepts. Recent successes in the application of modern chaos theory give us grounds to hope that our dream will become a reality in the near future.
I would like to express my gratitude to the Hebrew University of Jerusalem, the Lady Davis Trust, and Prof. R. Levine, as well as to BenGurion University of the Negev, its President, Dr. A. Braverman, and the Doron Foundation for Education and Welfare, for their support of my activities.
Pictures
REFERENCES
1. Prigogine, I. and Lefever, R. Symmetry breaking instabilities in dissipative systems. J. Chem. Phys. 16951700 (1968).