KINETICS OF SMALL SYSTEMS Thesis for fine Degree 9‘ M. S. fiiCREEfiz‘AN SMTE “N""ERSZT Lowe ii E. Jacobs 1963 ‘JHESIS p Q f.f.'".'f.' a rt! v “ I EAST LANSING, MICHIGAN my} . u . " n._ \,'5. v _.A\‘.:.:‘ ABSTRACT KINETICS OF SMALL SYSTEMS by Lowell E. Jacobs Chemical kinetics is usually approached in a deterministic way, i. e. , once the state of the system is known it can also be determined at any other time. However this approach does not allow for fluctuations in the system. The present work approaches chemi- cal kinetics from a stochastic viewpoint which will allow fluctuations in the system to be recognized. KINETICS OF SMALL SYSTEMS By Lowell E. Jacobs A THESIS Submitted to -Michigan State University in partial fulfillment of the requirements for the degree of MASTER OF SCIENCE Department of Chemistry 1963 AC KNOWLED GMENT The author wishes to express his appreciation to Dr. D. A. McQuarrie for his guidance and assistance during this investigation. Appreciation is also extended to the National Institutes of Health for partial financial assistance. :4: 3:: >:: ::< a}: 1:: >;< >',< 2:: * >1: ii TABLE OF CONTENTS Page INTRODUCTION 0 O O O I O O O O O O O O I O O O O O O O O O O O 1 Examination of the System: A—->- B. . . . . . . . . . 3 Examination of the System 2A—>- B. . . . . . . . . . 5 MonteCarloA—-—>-B................. 7 MonteCarloZA———>—B................ ll DiscussionofGraphs.................. 12 BIBLIOGRAPHY O 0 O I O O O O O O O O O O O O O I O O O O O O O 13 APPENDIX 0 O O O O O O O O O O O O O O O O O O O O 0 O O O O 0 O 15 iii INTR ODU CT ION There has recently appeared a treatment of the classical thermodynamics of small systems, in which the well-known relations of classical thermodynamics have been extended to be applicable to small systems [I]. A similar extension is also possible in chemical kinetics and this thesis is part of a study in this direction. The extension of kinetics to small systems (e. g. , systems with one hundred particles or less) is effected by an application of stochastic processes to chemical kinetics. The classical approach to chemical kinetics can be called deterministic since once the state of the system is known at some time t1, the state at any other time is automatically fixed, and no fluctuations or deviations about this value are recognized. A However, in a stochastic model of kinetics the state of the system, for instance the concentration at time t, is represented by an integer- valued random variable and the fundamental solution to the problem is the probability density function of the random variable, from which various moments may be calculated. Fluctuations are then recognized through the second central moment, the variance. Delbriick [2] seems to have the earliest representation of a chemical reaction by a stochastic model, in whichrhe discussed fluctu- ations in an autocatalytic reaction. Singer [3] used this approach next in 1953 when he discussed apparent irreproducibilities in some chemi- cal reactions and gave several examples of systems in which fluctuations, or irreproducibilities, have been reported. The reactions treated by Singer were a chain reaction without branching and one with branching. Renyi [4] next used this idea about a year later when he discussed the reaction A f B —-> C. Bartholomay [5, 6] has treated the irreversible unimolecular reaction A —->- B and shown the connection between the two previous approaches and stationary Markov chains. Ishida [7] then used a multidimensional stochastic treatment to discuss the unimolecular-decomposition theory of Kassel. Bartholomay [8] also has discussed a stochastic model for the Michaelis-Menten formu- lation of an enzyme-substrate reaction. The contribution that an idea such as the stochastic representation of chemical kinetics could make to the field of biology seems very useful. Several papers have already been published which use the stochastic approach to discuss reactions in cells and transitions of one type of cell to another. Berger [9], for example, discussed the transition of reticulocytes to red blood cells using a stochastic interpretation. Shea [10] has discussed histology from a statistical standpoint, and Fraser [11] has used this approach to study genetic systems._ The most recent work in applying stochastic models to chemical kinetics is that of McQuarrie [12], in which he discussed the systems, A—>— B, A_<_—>-_ B, and A—L B, Alia—>— C. McQuarrie also has another paper [13] in which the stochastic model for the bimolecular reactions, zA—->- B, A + B ——>- c, and A-—->— B autocatalytically are discussed. The purpose of the present work was to apply a Monte Carlo process to a pure "death" reaction, similar to Kendall's procedure [14] for a simple "birth and death" reaction, and to calculate the individual fluctuations in the system in order that they might be compared graphically 'to the ensemble average fluctuations obtained from McQuarrie's formula. Kendall's procedure and the modifications applied to it will be discussed in detail later. . The general idea was simply to derive certain fundamental equa- tions of probability concerning the system and then to set up a mathe- matically equivalent sampling process in order to solve the equations. This is the general technique used in the Monte Carlo process. In this problem the sampling process and the calculations were carried out on the 160A computer at Michigan State University. The sampling process was simply a series of random digits obtained from The Rand Corporation [15]. The complete program appears on page 16. This amounts to a simulated kinetics experiment. The reactions under con- sideration are A ———>— B and 2A—9— B. Let us look at an example in order to see just what procedure was followed in this problem. Suppose we take the reaction A ——9- B and look at a system of five particles for five different runs. That means we have five particles in our system initially and we let them decay according to our probability equations. When this process is repeated five times we have a system of five particles for five runs. This is equivalent to an ensemble of five systems, each of which contains five particles initially. The systems that were constructed are listed on page 12 . From the results of computations on these systems the graphs of the mean and variance were drawn. Examination of the System: A ——>— B Let us now consider the reaction A—->- B [12] using~ a stochastic model. Let the integer-valued random variable, X(t), be . the number of A molecules in the system at time t. Then the stochastic model is. completely defined by the following three assumptions: (1) The probability of a transition (x) —-—>- (x-l) in the interval (t,t + At) is kxAt + 0(At), where k is a constant and 0(At) means that [0(At)/At—>- 0 as At ——>- 0. (2) The probability of a transition (x) —>- (x - j), where j > 1, in the interval (t, t, + At) is 0(At). (3) The reverse reaction occurs with probability zero. Now let Px(t) =tProb (X(t) = x} . That is, Px(t) is the prob- ability that the number of A molecules is x at time t. 'Let Px(t +f At) = ‘Prob {X(t + At) = x}. That is, Px(t + At) is the probability that the number of A molecules is x at time t + At. The state X(t + At) = x can be obtained in two ways. First, X(t) may equal x and no transition occurs in time At; and second, X(t) may equal x + l and a transition does occur in time At. Px(t + At) is thus a sum of two probabilities as follows: Px(t + At) = [k(x+1)At] Px+1(t) + [1 - k x At] Px(t) + 0(At) (1) where k(x + 1)At is the probability of a transition from the state X(t) = x + l to the state X(t + At) = x; Px+1(t) is the probability that X(t) = x + l; [1 - kx At] is the probability of no transition from the state X(t) = x; and Px(t) is the probability that X(t) = x. The term, 0(At), is added on-in order to compensate for any higher order terms. Transpose Px(t) and divide through by At to get: 13,,“ + At) - pxm At = “X + 1)Px+l(t) - kxPxM + 0(At)/At (2) Th lim of the left hand side of this e uation is b defi 't' d? e At->- 0 C! Y n1 ion Jdt , Therefore: dP dtx 3 k(x + 1) Px+1(t) ' kXPx“) (3) whichis the differential-difference equation of the process. Now, by using the generating function of Px(t) [16], that is, co _ F (s,t) = xgo Px(t) sx , (4) where I s I < 1 equation (3) may be transformed into a partial dif- ferential equation by multiplying through by 3x and summing over x from O to co. m dP x a) X m X E —2‘-dt s = k E (x + l) Px+l(t) s - k xPx(t)s (5) x=o x=o x=o and so by noting equation (4), one obtains 9F _ cDF at - k(1-S) ()s (6) The solution of this equation uSing the initial condition F(s, o) = s o, 15 kt xo F(s,t): [1+(s-1)e' ] (7) By using the following relationships: . 00 E {X(t)] = (%{)s=l E XPX(t) and x=o 2 Dz ( X(t)} = ( 8a sF2 )s=1 + (3:)s=1 - (g: 25:1 where E {X(t)} is the expectation value or mean of X(t), and D2 {X(t)} is the variance, the following equations result: E [xm] O '1“ (a) D’- {x003 Note that the mean value of the stochastic representation is ll X (D er-kt (1 - e-kt) (9) identical to the deterministic result, showing that the two representations are consistent in the mean. The stochastic model, however, also gives higher moments and so fluctuations can now be included in chemical kinetic s . Examination of the System 2A —->- B Now consider the reaction 2A -—->- B [13]. Here the procedure is similar to that for A —->- B, i. e. , the bimolecular reaction is considered to be a pure "death" process with a continuous time para- meter and transition probabilities. Fir st consider a system of - X0 .molecules of reactant A at time t: 0. - Let X(t) equal the random variable, i. e. , the number of A molecules is the system at time t. Then the stochastic model is completely defined by the following assumptions: (1) The probability of the transition x -—>- x - 2 in the time interval (t, t + At) is %—k'x(x - 1)At + 0(At) where k' is a constant and 0(At) is defined so that 0(At)/At ->- O as At—->- 0. The factor of i— eliminates counting collisions between the same two molecules twice. We shall let fk' = k in subsequent equations. (2) The probability of the transition x—->-— x - j where j > 2, in the time interval (t , t + At) is 0(At). (3) The probability of the transition x —+— x + j, where j > 0, in the time interval (t , t + At) is zero. Thus the following probability equation can be written, - Px(t + At) = [k(x+2)(x+l)At] Px+2(t) + [1 —kx(xa=1)At] Px(t) + 0(At) (10) where the notation is the same as before. - Now, by the standard procedure of transposing *Px(t), dividing through by At, and taking the limit, the following differential-difference equation results; .%?& -.- k(x+2)(x+l) Px+z(t) - kx(x = 1) Px(t)- (11) This is transformed into a partial differential equation for the generating function, F(s,t) , viz: 3F _ 53F Tt " (1- $2) "—z'a s (12) whose solution for the initial condition F(s, o) = s 0 is: x0 i F(s,t) = l. + Z AjCj T (3) exp {-i—jU-Ukt} (l3) j=2 where r ,xo - j+ 1) 25 - 1 f (x +1) ‘ 2, A ‘ [ l l 9 l (14) 2 — +1 + +1 J f‘ (X0 1 )f' (x 2.1 _p____) _i and Cj T (s) are Gegenbauer polynomials [17], and j= 2,4,6,-—-,xo (x0 even). It can be easily shown that: xo = 2: A. exp [-é— j(j - mitt] (15) \ j: 2 J X0 . . = .74 2 Aj [Lil—1)] exp [- i- j (j - 1m] (16) J: A study of second order reactions by the stochastic approach pro- vides an excellent method of determining to what extent the deterministic rate equation is applicable to systems where the inherent statistical fluctuations are relatively large. McQuarrie [12] has stated "that when dA a deterministic rate equation such as - 21?- = kAZ is considered there is a tacit assumption that A" = (AV, which is true only for a delta function type density function, i. e. , all central moments equal zero. " Dole and Inokuti [18] have studied the first and second order kinetics for systems in which the number of reacting species was of the order of 10 or less. They considered the conditions under which the ratio 1 3/ (A)2 approaches unity. Their calculations seemed to indicate that the ratio approaches unity much more rapidly than one might expect. Monte Carlo A —-->- B The procedure used in this problem was similar to that employed by Kendall [14]; the equations being modified to fit a pure "death" process. Kendall derives three differential-difference equations based on the probability equation for—the birth and death) process [16]: Px(t+ At) = ,1 (x + 1)At Px+1(t) + X(x-1)At Px_l(t) + [1-(i+(1)xAt]Px(t) (17) where X is the expected birth rate and p. is the expected death rate. The standard procedure of transposing Px(t), dividing through by At, and taking the limit, leads to: (18) )L (x+1) Px+1(t) + MX-l) Px-1(t) - (Hp) X Px(t) dPx(t) = dt IInagine now the portion of the system represented by the following diagram: J l 1 f I I -l x x+1 x we can write three By forbidding transitions out of the states xi 1 equations which give the probability of the system being in any one of the above three states at time t + At. Px+1(t+At)= ()txAt)Px(t) + Px+1(t) (19) Px(t + At) = [1 - (,1 H.) x At]Px(t) Px_1(t+ At) = p x Ath(t) + Px_1(t) These equations lead to: dP‘S'tfit) = x x M) (20) dde(tt) = _ (H + x) x PX“) J—L—dp - (t) = p. x Px(t) dt Introduce a new variable u where u = NT; N is the number of particles present and ’C’ is the length of time until the next "incident. " Then, from the second equation in (20) it follows that the distribution of ’Z,’ is: e'()‘ + “NT (x + u)NdT (o < T< oo ). Kendall [14] has shown that this "birth-and-death" process is mathematically equivalent to a scheme of sequential sampling using random numbers. This process of deriving equations to explain a particular phenomenon and then employing a mathematically equivalent scheme to obtain a solution is called the "Monte Carlo" method. _This was exactly the procedure used in this problem, a pure "death" process (i. e. ,' x 3' 0). Here we are concerned with the second equation in (20). The solution of this equation is mathematically equivalent to the follow- ing scheme of sequential sampling: random variables are drawn independently from a series of random numbers. After each step in the sequential sampling process the following quantities are calculated: Nr =- Nr-l ' 1 (21) + ur ( > 1) r r-l Nr-l _ tr=t The construction of the sample of u-variables is most easily effected by employing the transformation, x = e"u so that the x-variables must be uniformly distributed in- (0, 1). I Suppose we look at a typical set of calculations for a run with 10 particles. 10 i 33 Er IZir-l Zr _t_r-l 1 .32025 1.139 10 0.114 0.000 2 .27574 1.288 9 0.143 0.114 3 .79141 0.234 8 0.029 0.257 4 .33940 1.081 7 0.154 0.286 5 .53963 0.617 6 0.103 0.440 6 .39754 0.922 5 0.184 0.543 7 . 18502 1.687 4 0.422 0.727 8 .65081 0.430 3 0.143 1.149 9 .41346 0.883 2 0.441 1.292 10 .18978 1.662 1 1.662 1.733 11 .66021 0.415 0 00 3.395 Each interval is first assigned a random number X which is chosen independently. The u-variable is then obtained from the relation x = e’u. The values of N are listed from 10 to 0 since the reaction is unimolecular and hence decays one particle at a time. Next '7," is calculated from the equation 11 = NT. And finally t, the cumulative time, is obtained by adding the subsequent values of T . After we repeat this process for a series of runs we can tally the number of particles, N, in each run at various time. intervals and thus obtain a mean for the system at each time interval. . We can also obtain 2 the variance at each time interval from the equation 82 = 71— - N“ where Ni is the number of particles present at time ti' Thus we can now graph the mean and the variance for each system and compare it with the theoretical curve obtained from McQuarrie's formula [12]. 11 Monte Carlo 2A—+ B - For the reaction 2A —>— B we can again write three equations which give the probability that the system is in one of the states, x-2, x, orx+2 attime t+At. Px+2(t + At) = [1 - k (x+2)(x+1)At] Px+2(t) + k(x+4)(x+3) Px+4(t)At Px(t + At) = k(x+2)(x+1)Ath+z(t) + [1 — kx(x=1)At]Px(t) (22) Px_2(t + At) = k x (x-1)At Px(t) + [1-k(x-2)(x-3)At]Px_2(t) By restricting transitions from x :1: 2 we get: dPxiZ(t) _ dt - k (x+4)(x+3)Px+4(t) 513311- = -k x (x-l) Px(t) (23) dP t “"2329 = kx(x-1)Px(t) Now it is clear that the distribution of (E, is: e-kNm'lW kN(N-l)d’E’ (24) This necessitates a change in equations (21) to: Nr 3 Nr_1 '2 + “r (25) tr = tr-l N.,;(Nr-1) Now the sampling process and the calculations can be carried out in an exactly analogous manner to that for the reaction A ——->- B. . 12 Discus sion of Graphs ‘The graphs which resulted from this study are shown on the following pages. For the reaction A—->- B the graphs shown are (1) mean, XBAR, vs time and (2) Variance vs time and for the reaction 2A—>— B the graphs shown are (l) 1 . . XBAR vs time and (2) Variance vs time. The graphs correspond to various systems which consisted of differing numbers of runs and particles as outlined in the table , ,, below. . Inspection of the graphs shows how the individual fluctuations tend to follow the theoretical curves obtained from McQuarrie's formula. Reaction: A -—>- B - Reaction: 2A —->- B Systems Systems Particles Runs Particles Runs 5 5 6 5 5 25 6 10 5 50 6 20 5 100 6 50 10 5 20 5 10 10 20 10 10 25 20 20 10 50 20 50 20 5 100 4 20 10 100 10 20 25 20 50 100 5 100 10 10. 11. 12. l3. 14. 15. 16, BIBLIOGRAPHY Hill, T. L., J. Chem. Phys. 29.: 3182 (1962). _ "Thermodynamics of Small Systems" (W. A. Benjamin, Inc. , Delbriick,‘ M., J. Chem.-Phys. _§_, 120 (1940). Singer, K. 1,8 J. Roy. Statist. Soc. B15, 92 (1953). Renyi, A. ,1 Magyar Tudomanyos Ahad. Kem. Tudomanyok OsztalyAnak,‘.Kazlemenyei _71, 93 (1954). Bartholomay, A., Bull. Math. Biophys. _2_0_, 175 (1958). Bartholomay, A., Bull. Math. Biophys. El, 363 (1959). Ishida, K., Bull. Chem. Soc. Japan, 23, 1030 (1960). Bartholomay, A., Ann,_N. YJcad. Sci. 29., 897 (1962). Berger, P., J. Theor. Bio. _Z_,‘No. 3, 295 (1962). Shea, S. M., J. Theor. Bio. _3_, No. l, 111 (1962). Fraser, A. S.,- J. Theor. Bio. _2_, No. 3, 329 (1962). McQuarrie, D. A., J. Chem. Phys. _3_8, No. 2,433 (1963). McQuarrie, D. A., Jackimowski, C. J., Russell, M. E., J. Chem. Phys., to be published. Kendall, D. (3., J. Roy. Statist. Soc., B12, 116 (1950). The Rand Corporation, ("A Million Random Digits!“(Glencoe, Illinois: The Free Press, Publishers). Bharucha-Reid, "Elements of the Theory of Markov Processes and Their Applications" (New York: McGraw-Hill Book Company, Inc. , 1960). 13 17, 18. 19. 20. 21. 22. 23. 24. 14 Erdelyi, Magnus, ‘ Oberhettinger, and Tricomi, "Higher Transcen-L ._ dental. Functions" (New York: ’ McGraw-Hill Book Company, Inc. , 1953), p. 174. Dole,-'M.,. and Inokuti, M. ,8 J. Chem. Phys”. to be published. Kendall, D. c.,' J. Roy. Statist. Soc. B11, 230 (1950). Kendall, D. G., Proc. Third Berkeley Symposium on‘Math. 7 Statistics and Probability, i, 149 (1956). Meyer, H. A. (ed.),' "Symposium on Monte Carlo Methods" (New’York: John Wiley and Sons, Inc. , 1956). Beckenbach, E. F. , "Modern Mathematics for the Engineer" (New York: McGraw-Hill Book Company, 1956), p. 279. Fluendy and Smith, Quarterly Reviews 1_6_, 241 (1962). U. S. Department of Commerce-National Bureau of Standards. "Monte Carlo Method." (Washington 25, D. C.; U. S. Government Printing Office, 1951). APPENDIX 15 16 L. E. JACOBS 21-1016T A-B No. 4 8-1-63 100 101 230 15 102 200 502 503 504 505 506 507 508 DIMENSION X(l 1), U(11), N(11, 20), A(1 1), T(l l), NMOL(42, 20) READ 100, L, K FORMAT (I 3, 3X, 13) PRINT 1 FORMAT (lHO, 23H10 PARTICLES FOR 5 RUNS) DO 230 I=1, L READ 101, (N(I, IRUN), IRUN?I l, K) FORMAT (2613) CONTINUE DO 220 IRUN=1, K PRINT 15, IRUN FORMAT (1H0, 5H1RUN=I3 PRINT 2 FORMAT (1HO, 49H X U N A T) T(1)=0. 000 READ 102, (X(I), I=1, L) FORMAT (6F8. 8) DO 200 I=1, L U(I)=-LOGF(X(I)) A(I)=U(I)/N(I, IRUN) PRINT 3. X(I), U(I), N(I, IRUN), A(I), T(I) FORMAT (1H2. F11. 8, 3X, F8.4, 3X, 14, 3X, F8. 4, 3X, F8.4) T(I+1)=T(I) + A(I) TI=0. 000 NMOL(l, IRUN)=N(I, IRUN) 1:2 IY=1 DO 210 IFORTY=1,40 TI=TI + 0.1 IY=IY + 1 IF (TI-T(II) 502, 503, 504 NMOL (IY, IRUN)=N(I- 1, IRUN) GO TO 210 NMOL(IY, IRUN)=N(I, IRUN) GO TO 209 IF (TI-T(I+1)) 505, 506, 507 NMOL(IY, IRUN)=N(I, IRUN) GO TO 209 NMOL(IY, IRUN)=N(I+1, IRUN) GO TO 210 IF (TI-T(l+ 2)) 508, 509, 510 NMOL(IY, IRUN)=N(I+1, IRUN) GO TO 210 17 509 NMOL(IY, IRUN)=N(I+2, IRUN) GO TO 209 ‘ 510 IF (TI-T(I+3) 511,512,513 511 NMOL(IY, IRUN)=N(I+2, IRUN) GO TO 209 512 NMOL(IY, IRUN)=N(I+3, IRUN) GO TO 210 513 IF (TI-T(I+4)) 514, 515, 516 514 NMOL(IY, IRUN)=N(I+3, IRUN) GO TO 210 515 NMOL(IY, IRUN)=N(I+4, IRUN) GO TO 209 516 IF (TI-T(I+5)) 517, 518, 519 517 NMOL(IY, IRUN)=N(I+4, IRUN) GO TO 209 518 NMOL(IY, IRUN)=N(I+5; IRUN) GO TO 210 519 IF (TI-T(I+6)) 520, 521, 522 520 NMOL(IY, IRUN)=N(I+5, IRUN) GO TO 210 521 NMOL(IY, IRUN)=N(I+6, IRUN) GO TO 209 ,522 IF (TI-T(I+7)) 523, 524, 525 523 NMOL(IY, IRUN)=N(I+6, IRUN) GO YO 209 524 NMOL(IY, IRUN)=N(I+7, IRUN) GO TO 210 525 IF (TI-T(I+8)) 526, 527, 528 526 NMOL(IY, IRUN)=N(I+7,IRUN) GO TO 210 527 NMOL(IY, IRUN)=N(I+8, IRUN) GO TO 209 528 IF (TI-T(I+9)) 529, 530, 531 529 NMOL(IY, IRUN)=N(I+8, IRUN) GO TO 209 530 NMOL(IY, IRUN)=N(I+9, IRUN) GO TO 210 531 IF (TI-T(I+10)) 532, 533, 534 532 NMOL(IY, IRUN)=N(I+9, IRUN) GO TO 210 533 NMOL(IY, IRUN)=N(I+10, IRUN) GO TO 209 534 IF (TI-T(I+ll)) 535, 536, 537 535 NMOL(IY, IRUN)=N(I+10, IRUN) GO TO 209 18 536 NMOL(IY, IRUN)=N(I+11, IRUN) GO TO 210 537 - IF (TI-T(I+12)) 538,539,540 538 NMOL(IY, IRUN)=N(I+11,IRUN) GO To 210 _ 539 NMOL(IY, IRUN)=N(I+12, IRUN) GO TO 209 540 IF (TI=-T(I+13)) 541, 542, 543 541 NMOL(IY, IRUN)=N(I+12, IRUN) GO TO 209 542 NMOL(IY, IRUN)=N(I+13, IRUN) GO TO 210 543 IF (TI—T(I+14)) 544, 545, 545, 544 NMOL(IY, IRUN)=N(I+13, IRUN) GO TO 210 545 NMOL(IY, IRUN)=N(I+14, IRUN) 209 I=I+l 210 CONTINUE 220 CONTINUE PRINT 8 8 FORMAT(1H1, 50HL. E. JACOBS A-B No.4 10 PARTICLES FOR 5 RUNS 1-8) 5 . . PRINT 6 6 FORMAT (1H0, 65H TI SUMXSQ EXSQ XBAR XBARSQ 1 VAR STDEV) 1 TI=0. 000 D0 203 IY=1,41 SUMXSQ=OOO. 00 DO 201 IRUNzl. K 201 SUMXSQ=SUMXSQ + (NMOL(IY, IRUN)**2 EXSQ=SUMXSQ/K SUMeooo. 00 DO 202 IRUN=1, K 202 SUM=SUM + NMOL(IY, IRUN) XBAR=SUM/K XBARSQ=XBAR**2 VAR=EXSQ-XBARSQ STDEV=SQRTF(VAR) PRINT 7, TI, SUMXSQ, EXSQ, XBAR, XBARSQ, VAR, STDEV 7 FORMAT (1H2, F7. 3, 2X, F7. 0, 2X, F9. 2, 2X, F7. 2, 2X,- F9. 2, 2X, F7. 3, 2X,F7.3) 203 TI=TI + 0. 1 STOP END 19 * L.E. JACOBS 21-1016T 2A-B_No. 7 87-15-63. 100 101 230 15 102 200 502 503 504 505 _ 506 507 508 DIMENSION X_(ll),U(11),N(11,10),A(11),T(11),NMOL(46, 10) READ 100,- L,K, FORMAT (13, 3X, 13) PRINT 1 _ FORMAT (1H0, 24H20 PARTICLES FOR ~110 RUNS) DO 230 I=1, L READ 101, (N(I, IRUN), IRUN-:1, K) FORMAT (2613) CONTINUE Do 220 IRUN=1, K PRINT 15,IRUN FORMAT (1H0, 5HIRUN=I3) PRINT 2 FORMAT (1112,4914 x U N A T) T(1)=0. 000 READ 102, (X(I),I=1,L) FORMAT (6F8.8) DO 200 I=1, L U(I)=-LOGF(X(I)) A(I)=U(I)/(N(I, IRUN)*(N(I, IRUN)-1)) PRINT 3, X(I), U(1), N(I, IRUN), A(I), T(I) FORMAT (1H2, F11. 8, 3x, 118.4, 3x, 14, 3x, F8. 4, 3x, F8. 4) T(I+1)=T(I) + A(I) TI=0.000 NMOL(1,IRUN)=N(1, IRUN) I=2 IY=1 DO 210 IFORT=1,45 TI=TI + 0. 005 IY=IY + 1 IF (TI-T(I)) 502,503,504 NMOL (IY, IRUN)=N(I- 1, IRUN) GO To 210 NMOL(IY, IRUN)=N(I, IRUN) GO TO 209 (IF (TI-T(I+1)) 505,506,507 NMOL(IY, IRUN)=N(I, IRUN) GO TO 209 NMOL(IY, IRUN)=N(I+1, IRUN) GO TO 210 IF (TI-T(I+2)) 508,509,510 NMOL(IY, IRUN)=N(I+1, IRUN) GO TO 210 20 509 NMOL(IY, IRUN)=N(I+2, IRUN) GO TO 209 _ 510 IF (TI-T(I+3)) 511,512,513 511 ‘ NMOL(IY, IRUN)=N(I+2, IRUN) GO TO 209 . 1512 NMOL(IY, IRUN)=N(I+3, IRUN) GO TO 210 513 IF (TI-T(I+4)) 514, 515, 516 514 NMOL(IY, IRUN)+N(I+3, IRUN) GO TO 210 515 NMOL(IY, IRUN)=N(I+4, IRUN) GO TO 209 , 516 IF (TI-T(I+5) 517, 518, 519 517 NMOL(IY, IRUN)=N(I+4, IRUN) GO TO 209 518 NMOL(IY, IRUN)=N(I+5, IRUN) GO TO 210 519 IF (TI-T(I+6)) 520, 521, 522 520 NMOL(IY, IRUN)=N(I+5, IRUN) GO TO 210 ‘ 521 NMOL(IY, IRUN)=N(I+6, IRUN) GO TO 209 522 IF (TI-T(I+7)) 523, 524, 525 523 NMOL(IY, IRUN)=N(I+6, IRUN) GO TO 209 524 NMOL(IY, IRUN)=N(I+7, IRUN) GO TO 210 525 IF (TI-T(I+8)) 526, 527, 528 526 NMOL(IY, IRUN)=N(I+7, IRUN) GO TO 210 527 NMOL(IY, IRUN)=N(I+8, IRUN) GO TO 209 528 IF (TI-T'(I+9)) 529, 530, 531 529 NMOL(IY, IRUN)=N(I+8, IRUN) GO TO 209 530 NMOL(IY, IRUN)=N(I+9, IRUN) GO TO 210 531 ~ IF (TI-T(I+10)) 532, 533, 534 532 NMOL(IY, IRUN)=N(I+9, IRUN) GO TO 210 533 NMOL(IY, IRUN)=N(I+10, IRUN) GO TO 209 534 IF (TI-T(I+11)) 535,536,537 535 NMOL(IY, IRUN)=N(I+10, IRUN) GO TO 209 536 5.37 538 539 .540 541 542 543 544 545 209 210 220 201 202 203 10 21 NMOL(IY,- IRUN)=N(I+11, IRUN) GO TO 210 IF (TI-T(I+12)) 538, 539,540 NMOL(IY, IRUN)=N(I+11, IRUN) GO To 210 ' NMOL(IY‘,IRUN)=N(I+ 1 2 , IRUN) GO TO 209 , IF (TI-T(I+13)) 541,542,543 NMOL(IY, IRUN)=N(I+12, IRUN) GO To 209 NMOL(IY, IRUN)=N(I+13, IRUN) GO TO 210 IF (TI-T(I+14)) 544,545,545 NMOL(IY, IRUN)=N(I+13, IRUN) GO To 210 NMOL(IY, IRUN)=N(I+14, IRUN) I=I+1 CONTINUE CONTINUE PRINT 8 FORMAT (1H1, 50HL. E. JACOBS 2A-B No. 7 20 PARTICLES FOR RUNS 8-15) PRINT 6 FORMAT (1H0,65H TI SUMXSQ Exso XBAR XBARSQ VAR STDEV) TI=0. 000 Do 203 IY=1,46 SUMXSO=000. 00 DO 201IRUN=1,K SUMXSQ=SUMXSQ + (NMOL(IY, IRUN)**2 EXSQ=SUMXSQ/K SUM=000. 00 DO 202 IRUN=1-, K SUM=SUM + NMOL(IY, IRUN) XBAR=SUM/K . XBARSQ=XBAR* >1< 2 VAR=EXSO-XBARSQ STDEV=SORTF(VAR) PRINT 7, TI, SUMXSQ, EXSQ, XBAR, XBARSQ, VAR, STDEV FORMAT (1H2, F7. 3, 2x, F7. 0, 2x, F9. 2, 2x, F7. 2,2x, F9. 2, 2x F7.3,2X,F7.3) TI=TI+0. 005 STOP END LOO 2.00 Ti me O 5 Particles for 5 runs I 5 Particles for 25 runs A Theoretical curve ~ '.‘ “- 3.00 4.00 23 2.00 MB 0 5 Particles for 5 runs . 5 Particles for 25 runs A Theoretical curve .50 Variance .00 h .50 .00 ’ w .00 3.90 _ I.80 2.70 3.60 Time 24 AOB o 5 Particles for 50 runs I 5 Particles for 100 runs A Theoretical curve LOO 2.00 Time 3.00 4.00 25 O 5 Particles for 50 runs . 5 Particles for 100 runs A Theoretical curve 1.50 .00 Variance . 0| 0 '0900 .90 I.80 ' 2.70 3.60 Time 26 A-B O 10 Particles for 5 runs O 10 Particles for 10 runs [0.00 _ . I 10 Particles for 50 runs A Theoretical curve 5.00 O 8- U .0 X C IOOO O \ O .50 ‘ O O O O .00 100 . 2.00 3.00 4.00 _ Time 27 o 10 Particles for 5 runs . 10 Particles for 10 runs A Theoretical curve 300 Vaflance LOO .00 .00 .90 I80 2.70 3.60 Thus - 28 y 10 Particles for 25 runs ' 10 Particles for 50 runs A Theor etic a1 curve ' V p I 3.00 i" O Q Variance .00 ' .00 .90 |.80 2.70 3.60 Time __._ ,' _4_ "AN." 29 7 t % AoB 8.00 i _ . O 20 Particles for 5 runs O 20 Particles for 10 runs A Theoretical curve .03 O O Variance 4.00s .00 .90 l.80 2.70 3.60 .- Time . ~ 20.00 LOO OO LOO 30 “ 2.00 Time O 1' 20 Particles for 5 runs . 20 Particles for 25 runs I 20 Particles for 50 runs A Theoretical curve 3.00 . 4.00 31 6.00 Variance .00 .90 l.80 Time 7 20 Particles for '25 runs I 20 Particles for 50 runs A Theoretical curve 2.70 3.60 32 AeB a 100 Particles for 5 runs , . 100 Particles for 10 runs '00 A Theoretical curve I I. O .0 X .\q 5 O I , I O I . I O O .00 ' 1.00 ‘ 2.00 3.00 4.00 60' .5 (fl a: Variance Q 33 i.'80 Time 0 100 Particles for. 5 runs . 100 Particles for 10 runs A Theoretical curve A\‘ \ g A A ‘ A _. a... mam 2.70 3.60 34 . 10.00 , A... \ . ' o 10 Particles for 1 run I 10 Particles for 50 runs 750 L. O .0 X 5.00 2.50 .00 IIIIIII‘ .00 .90 1.80 ' 2.70 3.60 ' Time 35 20.00 A-B . 1 O 20 Particles for 1 run ) ' I 20 Particles for 50 runs 15.00 ‘ Xbar 5.00 .00 .90 l.80 2.70 3.50 ' Time 35 20.00 A-B _ i O 20 Particles for 1 run ) ' . 20 Particles for 50 runs 15.00 ' Xbar I0.00 = 5.00 .00 .90 * 1.80 2.70 3.60 ' Time 36 \ o 100 Particles for. 1 run I 100 Particles for 10 runs Xbar 50 25 o .x... .00 .90 l.80 2.70 3.60 Thne ’ 37 2A-B O 6 Particles for 5 runs . 6 Particles for 10‘ runs A Theoretical [curve 3.00 .100 .200 , Time .550 .450 Yiia‘r .550 . .250 .I5 0 000 D 38 2A-B .l00 Time 0 6 Particles for 5|runs V 6 Particles for 20 runs I 6 Particles for 50 runs A Theoretical curve .2 OO 39 2AoB V ‘6 Particles for 20 runs . 6 Particles for 50 nuns A Theoretical curve ‘ 3.00 2.00 '5 Varlance 0 .000 100 ' .200. Time 40 2A-B Q. 20 Particles for 5 runs I2,0o - ‘ . 20 Particles for 10 runs A Theoretical curve 8£KD Vaflance :5 O O '100 .000 J00 .200- Thne 41 ZA-B .450? .350? 6 C. O .250 u 'l. O . V .x. 150s .050* .000 100 Time 8 20 Particles for 5- runs V 20 Particles for 20 runs I 20 Particles for 50 runs A Theoretical curve .200 42 2A~B 9.00 v 20 Particles for 20 runs I 20 Particles for 50 runs A Theoretical curve Variance 3.00;; .000 .i00 .200 Time 43 ZA-B ' O 100 Particles for 4 runs I 100 Particles for 10 runs .400 .300 Xbar .200 .IOO ° .0000 .0500 .IOOO 44 ZA—B o 100 Particles for 4 runs . 100 Particles for 10 runs 30.00 20.00 . ,/ Variance l0.00 ' Q.%%OO .0500 .IOOO m 9 .650 .550 .450 his .350! .250 ' .000 45 2A-B J00 Thne o 6 Particles for 1 run I 6 Particles for 50 runs .200 .450 .35Ci I YBar .25Ci J5C) . O 05.000 46 2A-B J00 Time 0 20 Particles for 1 run I 20 Particles for 50 runs JEOO 47 .500 7 _ ”'3 o 100 Particles for 1 {run .400 .300 I Xbar .200” J00r “.0000 .0500 JOOO gums-m uafim 'l I l l I II I l l I I l I I ll llillll lllflllmllllmlllll IIHIIHHIHH 93 03083 2343 31