- Email: [email protected]
Iterative solution of the master equation for thermal unimolecular reactions by the method of variable successive overrelaxation
Volume 131. number 6
CHEMICAL PHYSICS LETTERS
21 November 1986
ITERATIVE SOLUTION OF THE MASTER EQUATION FOR THERMAL UNIMOLECULAR REACTIONS BY THE METHOD OF VARIABLE SUCCESSIVE OVERRELAXATION Sung Hoon KANG
and Kyung-Hoon
JUNG
Department of Chemistry, Korea Advanced Institute of Scienceand Technology, P.O. Box 150 Chongyang, Seoul 131, Korea
Received 4 August 1986
A new algorithm which we have &led the variable successive overrelaxation (VSOR) method is suggested to solve iteratively the master equation for thermal unimolecular reactions. The rate of convergence using this method is at least six times faster than using the well-known Tardy-Rabinovitch algorithm. A comparative study of various iterative techniques has been carried out on a well-known system (the thermal isomerization of cyclopropane-1 ,I-&).
1. Introduction (1)
The study of thermal unimolecular reactions with various collision partners has given useful information on collisional energy transfer from highly excited vibrational energy levels. These studies require the solution of a master equation describing the time evolution of the population of the various energy levels in the molecule. While analytical solutions of the master equation have been obtained in special cases by Troe [ 11, Forst [2], and Pritchard [3], it is more usual to fmd numerical solutions for any assumed transition probabilities by an iterative method developed by Tardy and Rabinovitch [4] or by many-shot expansion techniques [S] . These widely used iterative methods suffer, however, from very slow convergence and often require more than 100 iterations, especially for we&-collision systems. In this work, we present a new iterative algorithm called the variable successive overrelaxation (VSOR) method, which shows a marked increase in the rate of convergence and is computationally very efficient. Considerable emphasis has been placed on the convergence analysis of the VSOR method applied to the numerical solution of the master equation for thermal reaction. A master equation for thermal unimolecular reactions can be written as 496
where ni is the population of energy level i in the molecule, fi the rate of external input,pij the probability of a collisional transition from energy levelj to i, ki the microscopic rate constant, and w the energy-independent collision frequency. The pii must satisfy the completeness condition, i.e. “ipii = 1, and detailed balance, i.e. pi&q = pj,$‘, where $9 denotes the equilibrium population of energy level i. Thermal reaction systems can be characterized under steady-state conditions by dni/dt = 0, and fi = 0, i.e. no external input. Furthermore, the steady&ate populations nf” between the zero-point energy and . some energy level Ek below the critical energy EO for reaction are good approximations to the equilibrium state populations nlg [ 1,461, i.e. nf” = nlg for i < k. Hence, eq. (1) for thermal systems under steady-state conditions reduces to
d”iss= wig piinyq + wlIZkpijn,?---nfS-kiniSS=O* dt
Applying the detailed balance condition eq. (2) then reads
(2) for the pii,
0 009-2614/86/$03.50 0 Elsevier Science Publishers (North-Holland Physics Publishing Division)
B.V.
Volume I3 1, number 6
bi =
j$kpiin,,,
=c j&k
= n;’ t kin;Slw -
p..neg = nlq c irk
I1 ’
c Pijnj”
j>k
2.1. Convergence analysis
p.. ”
for i> k.
(3)
This can be represented in matrix notation [4,5] as B= [(I-P)tdK]NSS=JNSS,
(4)
where I is the unit matrix and J the reduced transport matrix. All column vectors (B,NsS) and square matrices (P, K , J) have as many elements as there are energy levels under consideration.
2. Method For simplicity, eq. (4) is transformed into an equivalent equation by a prescaling procedure. Introducing a diagonal matrix F whose elements are the equilibrium reactant populations, the master equation (4) is then given by B = J(FF-l)NSS=(JF)(F-lNSS) =(JF)T
X=
FJTjY=AX,
A
i-l (bi -
c aij#+’ j=l
-
Xk+l = HXktC,
k=O,l,...,
(7)
where H is an iteration matrix. It can easily be seen thatifamatrixAisgivenbyA=M -N,thenfor H =M-lNandC=M-lB,X=HXtCifandonly if AX = B. Let the spectral radius p(H) of a matrix H be the maximum absolute value of all eigenvalues of H [7] : the basic convergence criterion for eq. (7) is that the iterative method is convergent for an arbitrary starting vector X0 if and only ifp(H) < 1 [7,8]. In order to obtain the properties of the iteration L-U,whereDis matrix H,, weassumeA=Ddefined as a diagonal matrix obtained from A by replacing all its off-diagonal elements by zero, and - L and - U represent the strictly lower and upper triangular parts of A , respectively. Then, if r is a diagonal matrix with elements rl, y2, .... y,, , eq. (6) can be rearranged into the matrix form DXk+l
= F J T , X = F-lNs,
Xf+l= (1 - +yi)-%$ t- ;l
It is necessary to prove that a given iterative method meets the requirements of sequential convergence for an arbitrary initial vector. To simplify the argument, consider the general iterative formula
(5)
and the superscript T denotes the transpose matrix. From detailed balance, it follows that both F-l J and A = FJ T are real symmetric matrices and X represents a column vector with elements xi = ny/J;:i = nfs/nfq. In order to solve the linear equation given above, we use the following iterative procedure,
where
2 I November 1986
CHEMICAL PHYSICS LETTERS
5 j=i+l
@) aijxJ
3
(6)
where the superscript k denotes the kth iteration and ri the relaxation parameter. Eq. (6) reduces to the Gauss-Siedel (GS) method if ri = 1 and to the successive overrelaxation @OR) method if ri = 7 + 1 [7,8]. In the present study, unlike in the SOR method, the relaxation parameter, ri, is varied from equation to equation in the range 0 < ri < 2: hence the name “variable overrelaxation method”.
-~LXk+l=(I-I’)DXktrUXktI-B.
This can
(8)
be rearranged as
Xk+l =H 7Xk + (D - rL)-lrB
,
(9)
where the iteration matrix H,=(D-rL)-‘[(I-r)DtrU].
Furthermore, it can be seen that M = P-1( D -lY’L) andN =r-‘[(I-T)Dt W] WhenAissplitinto A = M - N . We now, have to prove that the spectral radius of the iteration matrix H, given in eq. (9) is less than unity, i.e.p(H?)< 1. For the purpose of the proof for the convergence of the VSOR method, a well-known mathematical theorem [9] was utilized with no further proof, which states that if a Hermitian matrix A is given by A = M - N and M* t N is positive definite, then p(H) < 1 if and only if A is positive definite, where M* denotes the conjugate transpose matrix of M. Investigating the properties of the matrix A In the corresponding linear equation (S), A is obviously Hermitian, since it is a real symmetric matrix. From the completeness condition for the transition proba497
VolumeI31,number6
bility matrix P, it is readily seen that A has a weak diagonal dominance such that Iuiil > Zj+ laii I for i = 1) 2, .... n and for at least one i, laiil> Zjizilajil. Also; all diagonal elements of A are positive and there exists a positive diagonal matrix D such that A D is lower semistrictly diagonallydominant,i.e.aiidi > Zji+ilaiildj for i = 1,2, .... n, mdaiidi> Z&‘lajildi for i = 2,3, .... n. This property of Ais a necessary and sufficient condition that A is a non-singular matrix [8, p, 1341381. Making use of an important characteristic property of a Hermitian matrix A with non-negative diagonal elements and weak diagonal dominance (ref. [7], p. 41, theorem 5.5),‘we know that the matrix A is positive definite. Hence from the result mentioned above, it remains to prove that M* t N is positive definite in order for p(H,) to be less than unity. In the VSOR method, M = I’-1( D - r L) and N =I+[(&r)DWU]. Hence, M*tN=~‘-lD)*-L*tr-1D-D
=r-l(2br)D.
tIJ
(10)
Since 0 < ri < 2 for all i in the VSOR method, M * + N is positive definite. On developing the present method, the requirements of both completeness and detailed balance for collisional transition probabilities were essential to get the prescaled master equation (7) to obtain the various important properties of the matrix A and to make the rigorous convergence analysis. The SOR method using a constant relaxation parameter can be considered as a special case of the VSOR method. A is a Stieltjes matrix since A is positive definite, Uii> 0, and aji G 0 for i fj [7]. It is known (ref. [7] , ch. 12) that if A is a Stieltjes matrix, then a satisfactory value for the optimum relaxation parameter rb in the SOR method is rb = 2 {l - [ 1 P(C)~]L/?}-~, where C = I - (diag A)-l A. Since evaluating the eigenvalues of a large matrix C is not practical, the monotonicity theorem was used to obtain a good estimate of p(C). This theorem says that if A is positive definite and AL is obtained from A by deleting certain rows and corresponding columns of A,thenp(CL)
49&
21 November1986
CHEMICALPHYSICSLETTERS
below the critical energy, since C L is independent of the collision frequency for thermal reaction systems. when the practical choice of 7b is made by using the procedure mentioned above, rapid convergence occurs in the region of lower energy levels, but the reduced population xf in the higher energy states tends to be oscillatory during the iteration process due to a rather large VdUe of rb. Thus, one can attempt to obtain an optimum ri for each level by reducing yb monotonically, as given by ri = 1 t (Tb l)w(w + ki)-1. Any element ri lies in the range 1 < ri < 2 since 0 < O(W + kr)-’ < 1. The optimal value for 7b in the present work was found to be 1.48 by the trial and error method [7, p, 2091.
3. Application As an illustrative example of a thermal unimolecular reaction,‘Ge present the twochannel thermal isomerization of cyclopropane-1 J-d2 [lo]. Note that reliable and significant information on collisional energy transfer can be obtained from multichannel unimolecular reactions which produce kinetic data with less experimental error and show a more sensitive dependence on the collisional energy transfer parameters. The microscopic rate constant ki was caIculated from the RRKM expression. Vibrational frequencies, moments of inertia of the reactant and activated complexes, and critical energies for the cyclopropane-l ,I-d2 system were obtained from ref. [lo]. The Jacobi (J), Gauss-Siedal (GS), and SOR Table 1 Number of iterations required to solve the steady-state mapter equation (5) by various iteration methods a) k/k,
1.84x1O-2 8.14x1O-2 2.48x 10-l 5.27x 10-l 8.03x 10-r
Number of iterations
VSOR
TR
J
GS
SOR
24 26 28 29 29
144 156 167 176 196
122 133 141 153 164
66 12 II 80 83
26 30 31 53 96
a) For two-channel cyclopropane-1 ,l -d2 isomerization at 773 K with an exponential collisional transfer model for which the average energy transferred downward (AE)down is 300 cm-‘.
Volume 13 1, number 6
CHEMICAL PHYSICS LETTERS
schemes can be applied to the iterative solution of the prescaled steady-state master equation (5) with absolute confidence in sequential convergence for each iterative method, since the matrix A in eq. (5) is a symmetric positive-definite matrix with positive diagonal and negative off-diagonal elements [7,8]. The number of iterations required to solve the steadystate master equation is compared for these methods together with the present technique in table 1. The initial values for the relative populations were $ = w(w + k&l and convergence was considered to have been achieved when the value of lx:+1 - ~$1 lxik+’1-l became less than 10e5. The numerical results indicate that the convergence rate of the VSOR method is about six times as fast as that of the TR method and even faster than the other iteration methods. Using the VSOR procedure for the stepladder, Gaussian, and Poisson function models of the collisional transition probability, also shows a substantial improvement in convergence rates. The rate of convergence for collision systems with moderate and large average step sizes was found to be greater than for weak collision systems. It can therefore be seen that the VSOR method gives a considerable reduction in the CPU time required to solve a steady-state master equation for thermal unimolecular reactions.
21 November 1986
Acknowledgement The financial support of the Korea Science and Engineering Foundation is gratefully acknowledged. We would also like to thank Professor B.D. Choi, Department of Applied Mathematics, KAIST, for helpful discussions.
References [l] J. Troe, Ber. Bunsenges. Physik. Chem. 77 (1973) 665; 78 (1974) 478;J. Chem. Phys. 66 (1977) 4745. [2] AP. Penner and W. Forst, Chem. Phys. ll(1975) 243, J. Chem. Phys. 67 (1977) 5296; 72 (1980) 1435. [ 31 H.O. Pritchard, Accounts Chem. Res. 9 (1976) 99. [ 41 D.C. Tardy and B.S. Rabinovitch, J. Chem. Phys. 45 (1966) 3720;48 (1968) 1282. [5] W.G. Valance and E.W. S&lag, J. Chem. Phys. 45 (1966) 216,428O. [6] D.C. Tardy and B.S. Rabinovitch, Chem. Rev. 77 (1977) 369. [ 71 D.M. Young, Iterative solution of large linear systems (Academic Press, New York, 1971). [ 81 A. Berman and R.J. Plemmons, Nonnegative matrices in the mathematical sciences (Academic Press, New York, 1979). [ 91 J.M. Oreta and R.J. Plemmons, Linear Algebra Appl. 28 (1979) 177. [lo] I.E. Klein, B.S. Rabinovitch and K.-H. Jung, J. Chem. Phys. 67 (1977) 3833; I.E. Klein and B.S. Rabinovitch, Chem. Phys. 35 (1978) 439.
499
CHEMICAL PHYSICS LETTERS
21 November 1986
ITERATIVE SOLUTION OF THE MASTER EQUATION FOR THERMAL UNIMOLECULAR REACTIONS BY THE METHOD OF VARIABLE SUCCESSIVE OVERRELAXATION Sung Hoon KANG
and Kyung-Hoon
JUNG
Department of Chemistry, Korea Advanced Institute of Scienceand Technology, P.O. Box 150 Chongyang, Seoul 131, Korea
Received 4 August 1986
A new algorithm which we have &led the variable successive overrelaxation (VSOR) method is suggested to solve iteratively the master equation for thermal unimolecular reactions. The rate of convergence using this method is at least six times faster than using the well-known Tardy-Rabinovitch algorithm. A comparative study of various iterative techniques has been carried out on a well-known system (the thermal isomerization of cyclopropane-1 ,I-&).
1. Introduction (1)
The study of thermal unimolecular reactions with various collision partners has given useful information on collisional energy transfer from highly excited vibrational energy levels. These studies require the solution of a master equation describing the time evolution of the population of the various energy levels in the molecule. While analytical solutions of the master equation have been obtained in special cases by Troe [ 11, Forst [2], and Pritchard [3], it is more usual to fmd numerical solutions for any assumed transition probabilities by an iterative method developed by Tardy and Rabinovitch [4] or by many-shot expansion techniques [S] . These widely used iterative methods suffer, however, from very slow convergence and often require more than 100 iterations, especially for we&-collision systems. In this work, we present a new iterative algorithm called the variable successive overrelaxation (VSOR) method, which shows a marked increase in the rate of convergence and is computationally very efficient. Considerable emphasis has been placed on the convergence analysis of the VSOR method applied to the numerical solution of the master equation for thermal reaction. A master equation for thermal unimolecular reactions can be written as 496
where ni is the population of energy level i in the molecule, fi the rate of external input,pij the probability of a collisional transition from energy levelj to i, ki the microscopic rate constant, and w the energy-independent collision frequency. The pii must satisfy the completeness condition, i.e. “ipii = 1, and detailed balance, i.e. pi&q = pj,$‘, where $9 denotes the equilibrium population of energy level i. Thermal reaction systems can be characterized under steady-state conditions by dni/dt = 0, and fi = 0, i.e. no external input. Furthermore, the steady&ate populations nf” between the zero-point energy and . some energy level Ek below the critical energy EO for reaction are good approximations to the equilibrium state populations nlg [ 1,461, i.e. nf” = nlg for i < k. Hence, eq. (1) for thermal systems under steady-state conditions reduces to
d”iss= wig piinyq + wlIZkpijn,?---nfS-kiniSS=O* dt
Applying the detailed balance condition eq. (2) then reads
(2) for the pii,
0 009-2614/86/$03.50 0 Elsevier Science Publishers (North-Holland Physics Publishing Division)
B.V.
Volume I3 1, number 6
bi =
j$kpiin,,,
=c j&k
= n;’ t kin;Slw -
p..neg = nlq c irk
I1 ’
c Pijnj”
j>k
2.1. Convergence analysis
p.. ”
for i> k.
(3)
This can be represented in matrix notation [4,5] as B= [(I-P)tdK]NSS=JNSS,
(4)
where I is the unit matrix and J the reduced transport matrix. All column vectors (B,NsS) and square matrices (P, K , J) have as many elements as there are energy levels under consideration.
2. Method For simplicity, eq. (4) is transformed into an equivalent equation by a prescaling procedure. Introducing a diagonal matrix F whose elements are the equilibrium reactant populations, the master equation (4) is then given by B = J(FF-l)NSS=(JF)(F-lNSS) =(JF)T
X=
FJTjY=AX,
A
i-l (bi -
c aij#+’ j=l
-
Xk+l = HXktC,
k=O,l,...,
(7)
where H is an iteration matrix. It can easily be seen thatifamatrixAisgivenbyA=M -N,thenfor H =M-lNandC=M-lB,X=HXtCifandonly if AX = B. Let the spectral radius p(H) of a matrix H be the maximum absolute value of all eigenvalues of H [7] : the basic convergence criterion for eq. (7) is that the iterative method is convergent for an arbitrary starting vector X0 if and only ifp(H) < 1 [7,8]. In order to obtain the properties of the iteration L-U,whereDis matrix H,, weassumeA=Ddefined as a diagonal matrix obtained from A by replacing all its off-diagonal elements by zero, and - L and - U represent the strictly lower and upper triangular parts of A , respectively. Then, if r is a diagonal matrix with elements rl, y2, .... y,, , eq. (6) can be rearranged into the matrix form DXk+l
= F J T , X = F-lNs,
Xf+l= (1 - +yi)-%$ t- ;l
It is necessary to prove that a given iterative method meets the requirements of sequential convergence for an arbitrary initial vector. To simplify the argument, consider the general iterative formula
(5)
and the superscript T denotes the transpose matrix. From detailed balance, it follows that both F-l J and A = FJ T are real symmetric matrices and X represents a column vector with elements xi = ny/J;:i = nfs/nfq. In order to solve the linear equation given above, we use the following iterative procedure,
where
2 I November 1986
CHEMICAL PHYSICS LETTERS
5 j=i+l
@) aijxJ
3
(6)
where the superscript k denotes the kth iteration and ri the relaxation parameter. Eq. (6) reduces to the Gauss-Siedel (GS) method if ri = 1 and to the successive overrelaxation @OR) method if ri = 7 + 1 [7,8]. In the present study, unlike in the SOR method, the relaxation parameter, ri, is varied from equation to equation in the range 0 < ri < 2: hence the name “variable overrelaxation method”.
-~LXk+l=(I-I’)DXktrUXktI-B.
This can
(8)
be rearranged as
Xk+l =H 7Xk + (D - rL)-lrB
,
(9)
where the iteration matrix H,=(D-rL)-‘[(I-r)DtrU].
Furthermore, it can be seen that M = P-1( D -lY’L) andN =r-‘[(I-T)Dt W] WhenAissplitinto A = M - N . We now, have to prove that the spectral radius of the iteration matrix H, given in eq. (9) is less than unity, i.e.p(H?)< 1. For the purpose of the proof for the convergence of the VSOR method, a well-known mathematical theorem [9] was utilized with no further proof, which states that if a Hermitian matrix A is given by A = M - N and M* t N is positive definite, then p(H) < 1 if and only if A is positive definite, where M* denotes the conjugate transpose matrix of M. Investigating the properties of the matrix A In the corresponding linear equation (S), A is obviously Hermitian, since it is a real symmetric matrix. From the completeness condition for the transition proba497
VolumeI31,number6
bility matrix P, it is readily seen that A has a weak diagonal dominance such that Iuiil > Zj+ laii I for i = 1) 2, .... n and for at least one i, laiil> Zjizilajil. Also; all diagonal elements of A are positive and there exists a positive diagonal matrix D such that A D is lower semistrictly diagonallydominant,i.e.aiidi > Zji+ilaiildj for i = 1,2, .... n, mdaiidi> Z&‘lajildi for i = 2,3, .... n. This property of Ais a necessary and sufficient condition that A is a non-singular matrix [8, p, 1341381. Making use of an important characteristic property of a Hermitian matrix A with non-negative diagonal elements and weak diagonal dominance (ref. [7], p. 41, theorem 5.5),‘we know that the matrix A is positive definite. Hence from the result mentioned above, it remains to prove that M* t N is positive definite in order for p(H,) to be less than unity. In the VSOR method, M = I’-1( D - r L) and N =I+[(&r)DWU]. Hence, M*tN=~‘-lD)*-L*tr-1D-D
=r-l(2br)D.
tIJ
(10)
Since 0 < ri < 2 for all i in the VSOR method, M * + N is positive definite. On developing the present method, the requirements of both completeness and detailed balance for collisional transition probabilities were essential to get the prescaled master equation (7) to obtain the various important properties of the matrix A and to make the rigorous convergence analysis. The SOR method using a constant relaxation parameter can be considered as a special case of the VSOR method. A is a Stieltjes matrix since A is positive definite, Uii> 0, and aji G 0 for i fj [7]. It is known (ref. [7] , ch. 12) that if A is a Stieltjes matrix, then a satisfactory value for the optimum relaxation parameter rb in the SOR method is rb = 2 {l - [ 1 P(C)~]L/?}-~, where C = I - (diag A)-l A. Since evaluating the eigenvalues of a large matrix C is not practical, the monotonicity theorem was used to obtain a good estimate of p(C). This theorem says that if A is positive definite and AL is obtained from A by deleting certain rows and corresponding columns of A,thenp(CL)
49&
21 November1986
CHEMICALPHYSICSLETTERS
below the critical energy, since C L is independent of the collision frequency for thermal reaction systems. when the practical choice of 7b is made by using the procedure mentioned above, rapid convergence occurs in the region of lower energy levels, but the reduced population xf in the higher energy states tends to be oscillatory during the iteration process due to a rather large VdUe of rb. Thus, one can attempt to obtain an optimum ri for each level by reducing yb monotonically, as given by ri = 1 t (Tb l)w(w + ki)-1. Any element ri lies in the range 1 < ri < 2 since 0 < O(W + kr)-’ < 1. The optimal value for 7b in the present work was found to be 1.48 by the trial and error method [7, p, 2091.
3. Application As an illustrative example of a thermal unimolecular reaction,‘Ge present the twochannel thermal isomerization of cyclopropane-1 J-d2 [lo]. Note that reliable and significant information on collisional energy transfer can be obtained from multichannel unimolecular reactions which produce kinetic data with less experimental error and show a more sensitive dependence on the collisional energy transfer parameters. The microscopic rate constant ki was caIculated from the RRKM expression. Vibrational frequencies, moments of inertia of the reactant and activated complexes, and critical energies for the cyclopropane-l ,I-d2 system were obtained from ref. [lo]. The Jacobi (J), Gauss-Siedal (GS), and SOR Table 1 Number of iterations required to solve the steady-state mapter equation (5) by various iteration methods a) k/k,
1.84x1O-2 8.14x1O-2 2.48x 10-l 5.27x 10-l 8.03x 10-r
Number of iterations
VSOR
TR
J
GS
SOR
24 26 28 29 29
144 156 167 176 196
122 133 141 153 164
66 12 II 80 83
26 30 31 53 96
a) For two-channel cyclopropane-1 ,l -d2 isomerization at 773 K with an exponential collisional transfer model for which the average energy transferred downward (AE)down is 300 cm-‘.
Volume 13 1, number 6
CHEMICAL PHYSICS LETTERS
schemes can be applied to the iterative solution of the prescaled steady-state master equation (5) with absolute confidence in sequential convergence for each iterative method, since the matrix A in eq. (5) is a symmetric positive-definite matrix with positive diagonal and negative off-diagonal elements [7,8]. The number of iterations required to solve the steadystate master equation is compared for these methods together with the present technique in table 1. The initial values for the relative populations were $ = w(w + k&l and convergence was considered to have been achieved when the value of lx:+1 - ~$1 lxik+’1-l became less than 10e5. The numerical results indicate that the convergence rate of the VSOR method is about six times as fast as that of the TR method and even faster than the other iteration methods. Using the VSOR procedure for the stepladder, Gaussian, and Poisson function models of the collisional transition probability, also shows a substantial improvement in convergence rates. The rate of convergence for collision systems with moderate and large average step sizes was found to be greater than for weak collision systems. It can therefore be seen that the VSOR method gives a considerable reduction in the CPU time required to solve a steady-state master equation for thermal unimolecular reactions.
21 November 1986
Acknowledgement The financial support of the Korea Science and Engineering Foundation is gratefully acknowledged. We would also like to thank Professor B.D. Choi, Department of Applied Mathematics, KAIST, for helpful discussions.
References [l] J. Troe, Ber. Bunsenges. Physik. Chem. 77 (1973) 665; 78 (1974) 478;J. Chem. Phys. 66 (1977) 4745. [2] AP. Penner and W. Forst, Chem. Phys. ll(1975) 243, J. Chem. Phys. 67 (1977) 5296; 72 (1980) 1435. [ 31 H.O. Pritchard, Accounts Chem. Res. 9 (1976) 99. [ 41 D.C. Tardy and B.S. Rabinovitch, J. Chem. Phys. 45 (1966) 3720;48 (1968) 1282. [5] W.G. Valance and E.W. S&lag, J. Chem. Phys. 45 (1966) 216,428O. [6] D.C. Tardy and B.S. Rabinovitch, Chem. Rev. 77 (1977) 369. [ 71 D.M. Young, Iterative solution of large linear systems (Academic Press, New York, 1971). [ 81 A. Berman and R.J. Plemmons, Nonnegative matrices in the mathematical sciences (Academic Press, New York, 1979). [ 91 J.M. Oreta and R.J. Plemmons, Linear Algebra Appl. 28 (1979) 177. [lo] I.E. Klein, B.S. Rabinovitch and K.-H. Jung, J. Chem. Phys. 67 (1977) 3833; I.E. Klein and B.S. Rabinovitch, Chem. Phys. 35 (1978) 439.
499