THE UNIVERSITY OF MICHIGAN INDUSTRY PROGRAM OF THE COLLEGE OF ENGINEERING A METHOD FOR THE ANALYSIS OF MULTICOMPONENT EXPONENTIAL DECAY CURVES Donald G. 9ardner Jeanne C. Gardner George Laush W. Wayne Meinke fE U'NIVERSITY OF MICHIGAN ENGINEERING LIBRARY March, 1959 IP-358

N ^ ^ \l I Nf^ 4"\

ACKNOWLEDGMENT The authors are indebted to Professors Alan Perlis and John Carr III for many valuable discussions of the problem, and to Professor W. Kehl for making available the University of Pittsburgh's IBM 650 computer. ii

TABLE OF CONTENTS Page ACKNOWLEDGMENT.............................................. ii LIST OF TABLES................................................... iv LIST OF FIGURES................................................ v I. INTRODUCTION................................................ 1 II. PREVIOUS METHODS........................................... 2 III. SOLUTION BY FOURIER TRANSFORMS............................ 4 IV. A FEW DETAILS ON THE NUMERICAL SOLUTION................... 7 A. Cut-Off Error.......................................... 7 B. Numerical Integration................................. 9 V. RESULTS OF THE NUMERICAL EVALUATION................. 11 A. Effect of Cut-Off with Respect to.................... 14 B. Effect of Cut-Off with Respect to x.................... 26 C. Effect of Poor Extrapolation of the Initial Data....... 30 D. Effect of Scatter in the Initial Data................ 33 E. Effect of Decreased Accuracy in the Numerical Integration........................................ 33 VI. DISCUSSION................................................. 35 iii

LIST OF TABLES Table Page I Input Data for the Accurate Decay Curve f(t) 100e-0.02t Together with % Deviation Factors Used to Simulate a Decay Curve with Scatter.................. 12-13 II Comparison of the Results Graphically Determined from Figure 7 for the Four-Component Decay Curve with the Actual Values......................................... 24 iv

LIST OF FIGURES Figure Page 1 Effect of increasing po from 2.0 to 4.0 in the analysis of a single-component decay curve. X = 0.02, xo = 7.0, A x = 0.25, and A g = 0.1. Ordinate units are arbitrary........................... 15 2 Effect of increasing io from 6.0 to 8.0 in the analysis of a single-component decay curve. X = 0.02, xo = 7.0, A x = 0.25, and A j. = 0.1. Ordinate units are arbitrary.. 16 3 Detail of the peak shown in Fig. 2 for io = 8.0, and the results for the analysis of a single-component decay curve with no = 9.0. X = 0.02, xo = 7.0, A x = 0.25, and A p = 0.1. Ordinate units are arbitrary.. 17 4 Effect of increasing po from 3.0 to 6.0 in the analysis of a two-component decay curve. X1 = 0.1, l = 0.01, xo = 7.0, A x = 0.25, and A p = 0.1. Ordinate units are arbitrary............................................. 19 5 Effect of increasing lo from 6.0 to 8.0 in the analysis of a two-component decay curve. X1 = 0.02, X2 = 0.01, xo = 7.0, A x = 0.25, and A. = 0.1. Ordinate units are arbitrary............................. 20 6 Effect of increasing po from 6.0 to 8.0 in the analysis of a three-component decay curve. X1 = 0.1, X = 0.02, 3 = 0.01, xO = 7.0, A x = 0.25, and A p = 0.1. Ordinate units are arbitrary....................................... 22 7 Results of the analysis of a four-component decay curve with Po = 8.0. X1 = 0.5, X = 0.1, 3 = 0.02, 4 = 0.01, o 7.0, A x = 0.25, and A = 0.1. Ordinate units are arbitrary................................................. 23 8 Plot of the function eXf(eX)versus x for a single component decay curve with X = 0.02. Cut-offs at xo = 5.25 and xo = 6.25 shown............................. 27 9 Effect of cut-off at xo = 5.25 in the analysis of a single-component decay curve with o0 decreasing from 6.0 to 4.0. X = 0.02, A x = 0.25, and A p = O.1. Ordinate units are arbitrary.............................. 28 10 Effect of cut-off at xo = 6.25 in the analysis of a single-component decay curve with pO increasing from 6.0 to 8.0. X = 0.02, A x = 0.25, and A. = 0.1. Ordinate units are arbitrary............................ 29 v

LIST OF FIGURES (CONT'D) Figure Page 11 High and low extrapolations of the single-component decay curve for X = 0.02............................... 31 12 Effect of poor extrapolation of the initial data in the analysis of a single-component decay curve. X = 0.02, xo = 7.0, A x = 0.25, and A la = 0.1. Ordinate units are arbitrary.............................. 32 13 Effect of scatter in the initial data in the analysis of a single-component decay curve. X = 0.02, x = 7.0, A x = 0.25, and A p = 0.1. Ordinate units are arbitrary.. 34 vi

I. INTRODUCTION There are several types of problems in science in which experimental observations may best be represented by a linear combination of exponentials of the form: n -t f(t) = Ni et In these problems the parameters Ni and Xi have biological or physical significance. Therefore, in fitting a function of the above form to the data it is not sufficient that the function merely approximate the data closely, but it is also necessary that the parameters be accurately estimated. It should be noted that in Equation (1) the exponentials are all assumed to be separate and unrelated, i.e., none of the components are produced as the result of the decay of another component. The problem may be stated as follows: a function f(t) is approximated by experimentally determining an estimate of f(t) at a finite number of values of t. From this discontinuous set of data it is desired to obtain n(total number of components), and estimates of the Ni's and the X is. The essential difficulties in the solution of this problem are that we are dealing with a series of non-linear equations, that the data are only approximating the function f(t) over a finite range in t and that the exponential series possesses strongly non-orthogonal properties. With respect to this problem, the authors have principally been concerned with the analysis of multicomponent radioactive decay curves. However, similar problems arise in the study of a) first order chemical kinetics, b) certain diffusion problems, such as neutrons in a moderator,

-2c) some order-disorder transitions in solid state physics, d) dielectric properties of certain compounds, e) relaxation properties of organic polymers, f) pulses in electrical networks, g) survival and mortality experiments in the biological sciences and h) servoproblems of the guided missile type, to mention a few examples. Since lately in many cases it has become technically feasible and even convenient to obtain experimental data of reasonably good accuracy, the method of analysis of these data assumes greater importance. The purpose of this paper is to describe a mathematical method of analysis which appears to possess certain advantages over previous methods. The method has been evaluated on the IBM 650 computer located at the University of Pittsburgh. II. PREVIOUS METHODS By far the most common method used to resolve a decay curve into its components is the graphical approach. Here the data are plotted on semilog paper, and the curve resolved by a repeated subtraction of straight lines. The limitations of the method are apparent and need not be enumerated here. The method is, however, certainly the easiest to perform. The method may be considerably refined by employing a least-squares technique to fit the straight lines, and some error estimation becomes available. The difficulties inherent in the subtraction procedure, however, still remain. Mathematical approaches to this problem have been suggested by Proney(l), Hudson(2), (1) F. B. Hildebrand, Introduction to Numerical Analysis (McGraw-Hill Book Company, Inc., New York, 1956). (2) G. E. Hudson, Am. J. Phys., 21, 362 (1953).

-3Householder(3), Cornell(4), and Ziegler(5). To the authors' knowledge none of these methods, except perhaps that due to Zeigler, are in any sort of routine use. In none of these methods does the number of components "fall out" of the analysis, although in two cases(3'5) tests are included to determine the number needed to adequately fit the data. As Lanczos(6) has pointed out, there are a number of simple and straight-forward mathematical solutions to the problem of separating exponentials, but unfortunately enormous practical problems arise when they are applied to experimental data from physical experiments. The principal reason for this is the exceedingly non-orthogonal behavior of the exponential functions. The end result is that the initial data must be extremely accurate if more than two or three exponentials are to be separated. In most cases the accuracy required is far beyond that usually available. In Lanczos' opinion no amount of least-square or other statistical treatment can make up for the extreme sensitivity of the parameters to very small changes in the initial data. In the method to be described here a transformation of the initial function to the complex plane is made, wherein the new function exhibits entirely different properties. The hope that this method may succeed where other methods, which deal with purely real numbers, fail is based on the presence of periodic functions which arise during the analysis. (3) A. S. Householder, U.S. Atomic Energy Commission Report ORNL-455 (Feb. 1950). (4) R. G. Cornell, U.S. Atomic Energy Commission Report ORNL-2120 (Sept. 1956). (5) G. R. Keepin, T. S. Wimett, and R. K. Zeigler, J. Nuclear Energy, 6, 1 (1957). (6) C. Lanczos, Applied Analysis (Prentice-Hall, Inc., New York, 1956).

-4III. SOLUTION BY FOURIER TRANSFORMS The function f(t) in Equation (1) is in the form of a Dirichlet series which may be expressed as a Stieltjes integral. n f(t)= N e-itJ e-t dh(X) i=l (2) The function f(t) may also be expressed in the form of a Laplace integral equation. f(t) = e g(X) dX (3) Here h(X) is a step function and g(X) is a sum of delta functions. Due to the error inherent in the experimental estimate of f(t) and in the numerical computations necessary to obtain g(X), a plot of g(X) versus X will appear in the form of a frequency spectrum. The presence of a peak in the spectrum indicates a component, the abscissa value at the center of a peak is the decay constant Xi, while the height of the peak is proportional to the coefficient Ni. The problem then is to determine g(X) given the experimentally determined function f(t). The method advocated here is based on a well-known general approach for solving linear inte(7) gral equations.( We apply this approach to the specific case of the Laplace integral equation in a manner after Perlis.(8) A somewhat similar treatment has been described by Paley and Wiener.(9) (7) E. C. Titchmarsh, Introduction to the Theory of Fourier Integrals (Oxford University Press, London, 1937). (8) A. J. Perlis, U.S. Atomic Energy Commission Report NP-786 (Sept. 1948). (9) R. E. Paley and N. Wiener, Fourier Transforms in the Complex Domain (American Mathematical Society, 1934 )

-5We begin with Equation (3),,00 f(t)=J e-t g(X)d (3) and proceed to transform the variables X and t. Let X = e-Y and t = eX eTen f (ex) = e-e g(e-Y)e-Y dy Multiply both sides by eX 00 e f(eX)= e-e(x Y) e( )g(e'Y)dy -0 (5) Now F(P- e 0eX f(eX)e'dx 2 7 oo (6) where F(u) is the Fourier transform of exf(eX). Combining Equations (5) and (6) we obtain F()= /2ees r{ e g(e-y)dy}eip(S+Y)ds V27T -oO^o (7) Let s = x- y, or x = s + y. Then F(,) =4. g(eY)ei dy e e ds (8)

-6Rearranging we have F ( I e-e (-e(Y) e(x Y)g(e-y)dye ix dx vr2 ^oo ^ (9) Before proceeding further it is important to keep in mind that the term es = eX-Y was formed by combining in Equation (5) the term e-Y, obtained by differentiating X = e-Y, with eX. The kernel e-e(xY)e(-Y) in Equation (5) was kept intact and later separated from the function g(e-Y)dy. This last function is related to the original variables as follows: g(e-Y)dy =g(X)/X dA (10) Hence, when we eventually obtain g(e-Y) as a function of y from Equation (9), this will be equivalent to a plot of g(X)/X as a function of X. Returning to Equation (9) we see that the right hand side is the product of the Fourier transform G(p) of g(e-Y) and the Fourier transform K(i) of e-ees. Therefore, F(//) =/ G()K( ) K(11) and G(I F(u) G(.)=, - i/27r K(p) (12) Taking the inverse Fourier transform of G(k) we obtain g(e-Y)- e'YPM dp (13)

-7In this case K(p) can be evaluated analytically, and it turns out to be the Euler integral for the complex Gamma function: K(A)=-I r(l+i) (14) (14) IV. A FEW DETAILS ON THE NUMERICAL SOLUTION Briefly, the method of solution was shown to consist.of essentially only two integrations. First the Fourier transform F(p) is found using Equation (6). This is divided by the complex Gamma function given in Equation (14). Finally g(e-Y) as a function of y is found using the inverse Fourier transform shown in Equation (13). From the latter, a plot of g(X)/X versus X may be immediately obtained. Since it is convenient to use equidistant values of ui in determining F(p), K(4) can most easily be found from tabulations (10). The function g(e-Y) may be found using the same t values or an equidistant subset thereof. Two additional subjects warrant comment: these are the cut-off error and the set-up of the numerical integrations. A. Cut-Off Error It is clear that in Equations (6) and (13) one cannot numerically integrate from -a to a. Consider Equation (6). Here we must introduce the limits + xo which are the cut-off points of the integral. Instead of Equation (6) we now have xo,/ F(F )= exf(e )ei dx + E(xo,#) -Xo (15) (10) National Bureau of Standards, Applied Mathematics Series 34 (1954). "Tables of the Gamma Function for Complex Arguments."

-8The calculated 2 F((l) will be in error by at least the amount E(xou.). The major difficulty of the method is now apparent. We are trying to simulate a curve with an abrupt cut-off at xo by a sum of exponentlals that extend to x = oo. This has the effect of adding into F(p.) Fourier components which extend the range in 1i on which F(k) maintains appreciable value. Since K(pL) diminishes rapidly with increasing lt, for some value of 1 the quotient F(p)/K([i) in Equation (13) will begin to grow without bound. The'end result has been found to be error ripples in the plot of g(e-Y) versus y which tend to obscure the results. Hence, it is necessary that a finite x exist such that the value of E(xo,p) is sufficiently small that a good solution is possible. It is most unfortunate that once F(p) has been warped by the cut-off, the warping cannot be removed by any subsequent mathematical treatment. If it is not possible to experimentally follow the decay curve for a long enough time, then the data must be treated in some appropriate manner. For example, the longest-lived component might be extrapolated to give the necessary information. Even a somewhat inaccurate extrapolation will usually yield far better results than if no extrapolation is made. Alternatively, it might be possible to subtract off the longestlived component before the analysis. Finally, a drastic and unrecommended step would be to introduce a convergence factor into Equation (13). It was noted that the cut-off errors in Equation (6) tend to increase the height of the error ripples in the final results. A cutoff at A a + p. in Equation (13), on the other hand, has been shown to primarily affect the frequency of the error ripples and the breadth of the true peaks. The larger the value of p., the narrower and more well defined are the component peaks. Tne maximum useable value of pL will depend on

-9how good the initial data are and on the cut-off at x. If po is chosen too large, then the cut-off at xo will cause the amplitude of the error ripples to increase. If Lo is chosen too small, there will be an unnecessary loss in resolution of the peaks in the final result. The greater the value in p. at which F(i)/K(p) remains well behaved, the finer the resolution can be made in the final results. B. Numerical Integration It is convenient to adjust the units of t in the initial data such that the decay constants Xi fall in the range from 0 to 1. In other words, the half-life of the shortest-lived component should not be less than 0.693 units of t. This is no restriction since the range of t is infinite and the scale with which t is measured can be arbitrary. Next, each value of f(t) is multiplied by its value of t and the results plotted as exf(ex) versus x. A2c F(p.) is then related to the area under this new curve. Whereas t ranged in principle from 0 to oo, x now ranges from -o to o. Hence, we can set up the following integral: F()=j f [f*(x) + f*(-x)] e iX dx (16) or more conveniently, F(g) = O{ [ f*(x)+f* (-x)] cos Ax + i [f*(X)-f*(-x)]sin rx}dx (17)

-10Here we define f*(x) = eXf(eX). This yields real and imaginary parts of F(pL) which we term Fc and Fs respectively. K(.) is similarly composed of real and imaginary parts Kc and Ks. Hence 2 G(u)^ -F(p) Fc + i F (Fc + iFs )(Kc- iK ) K(/i) Kc + iKs K K (18) Next Ir F( ) e-iy d I Jo(Fc+ i Fs)(Kc-iKs) -!K(M) 0 r KC2+ K2 -9o - /0 C S -i sin yg)d/ (19) All odd terms vanish in Equation (19) yielding v I _oFcKc+_Fss K s Fs Kc-FcKs 7 K +K2K KK2+K 2 g (eYa) = fiO( F2K2FK cos yp 2 2sin y/1)d/U C cS S (20) As was pointed out before, if F(Q) is determined for an equidistant set of [i's ranging from 0 to pO in steps of say bA = 0.1, it is then convenient to use tabulated values of the gamma function. (10) In most cases it suffices to use a simple numerical integration scheme such as Simpsonts Rule for Equations (17) and (20), although if the accuracy of the initial data warrants it a more refined procedure such as that of Filon(ll) may be employed. Normally it is not convenient to take experimental data in precisely the right intervals that would best National Bureau of Standards, Applied Mathematics Series 34 (1954). "Tables of the Gamma Function for Complex Arguments." (11) L. N. G. Filon, Proc. Roy. Soc. Edin., 49, n8 (1928-29).

-11suit the method of analysis to be used to evaluate the data. In the present case, two alternatives arise. Either an unequal interval integration scheme could be employed to obtain the Fourier transform function F(l), or else an interpolation procedure might be used to obtain functional values of the initial data at the desired intervals. If a large number of data are available one might choose the former approach, whereas if the initial data are available at only a few points or if the data are badly scattered the latter method might prove the more desirable. V. RESULTS OF THE NUMERICAL EVALUATION To test out the method we have constructed one, two, three, and four component decay curves. In the most accurate of these curves [f(t) = lOOe*.O02t] the data ranged in accuracy from about 1 part in 105 at the beginning of the curve to 3 parts in 10 at the end. The accuracy of the remaining curves ranged from about 0.5 to 1 part in 104 at the beginning to perhaps 5 parts in 10 at the end, except for the curves wherein the data were deliberately distorted in order to study a particular desired effect. To simplify the calculations the data were constructed at equal intervals in the logarithm of t. This yielded many points near the beginning of the decay curve where f(t) is decreasing rapidly, and relatively few points at widely spaced intervals in t near the end of the curve where the longest-lived component is decaying slowly. The data thus produced were more realistically distributed with respect to the variable t than would have been the case had equal intervals in t been chosen, since one would normally tend to take more experimental points where the data were changing rapidly. An example of the decay curve data is given in Table I. Here the initial data for the single component curve f(t) = 100le'02t are

-12TABLE I INPUT DATA FOR THE ACCURATE DECAY CURVE f(t) = 100e-0'02t TOGETHER WITH % DEVIATION FACTORS USED TO SIMULATE A DECAY CURVE WITH SCATTER* t x f(t) Deviation 0 -oo 100.000 1.0000 0 98.020 +4.1 1.2840 0.25 97.467 -8.0 1.6480 0.50 96.758 -6.9 2.1170 0.75 95.854 -6.4 2.7183 1.00 94.706 -1.9 3.4903 1.25 93.259 0 4.4817 1.50 91.426 +0.5 5.7546 1.75 89.128 -1.8 7.3891 2.00 86.261 +3.6 9.4877 2.25 82.717 -2.5 12.182 2.50 78.380 +0.7 15.643 2.75 73.134 -2.1 20.086 3.00 66.915 +2.8 25.790 3.25 59.703 0 33.115 3.50 51.567 +1.9 42.521 3.75 42.724 -1.2 54.598 4.00 33.556 -1.9 70.105 4.25 24.608 -3.6 90.017 4.50 16.524 -.5 115.584 4.75 9.909 0 148.413 5.00 5.139 -1.4 (Continued) * Data points between t=O and t=l not listed.

-13TABLE I (CONT'D) INPUT DATA FOR THE ACCURATE DECAY CURVE f(t) = lOOe-0.02t TOGETHER WITH % DEVIATION FACTORS USED TO SIMULATE A DECAY CURVE WITH SCATTER* t x f(t) Deviation 190.566 5.25 2.212 +4.0 244.692 5.50 0.7493 +5.5 314.191 5.75 0.1866 0 403.429 6.oo 0.0313 -15.9 518.013 6.25 0.0032 +122.0 665.42 6.50 o 0 * Data points between t=O and t=l not listed.

-14shown, together with the factors used to distort the data when the effect of scatter in the original data was studied (see Section V-D). The points between t=O and t=l were omitted from the table since they would normally be obtained by interpolation. The accuracy of these points was, however, about 1 part in 105 for this curve. A. Effect of Cut-Off with Respect to p Using Equation (20) we determine a plot of g(e'Y) versus y by integrating the expression from 0 to.o for each desired value of y. Such a plot is equivalent to a plot of g(X)/x versus X. In Figure 1 we show the effect of increasing the final integration range from \o=2 to o -4, using the data for the single-component curve where X=0.02. It can be seen that the principal peak falls in the same place on each curve, just at the proper X value. The breadth of the principal peak and the smaller ripples appear to be caused by errors, primarily cut-off errors, in the calculation, and errors in the initial data. As the range in. is extended, the resolution of the peak becomes better. In Figure 2 the range in pt has been extended to 0=65 and o=8, and again the resolution increases progressively. It should be noted that the positions of the peaks in the error ripple change as a function of p0, whereas the position of the true peak does not. This fact provides one method for distinguishing small true peaks from error ripples —simply change po and note which peaks do not shift position. Once the positions of the true peaks are found they may be examined more closely by taking smaller intervals in y. The left-hand curve in Figure 3 shows an expanded view of the p =8 peak of Figure 2. Also shown in Figure 5 is the curve obtained for 0o=-9. While the resolution is excellent and the center of the peak falls at the proper place, it can

-15-0.02 t 24 f (t) 100e ^'o = 2.0 lo = 4.0 20 16 16- - - \ 12 4 0 - II I I I 111101 I I I 1111 I l I ii 1i iIII I I 1n I i 1.0 0. 0.0.01 1.0 0.1 0.01 X Fig. 1 Effect of increasing Po from 2.0 to 4.0 in the analysis of a single-component decay curve. X = 0.02, xo = 7.0, A x = 0.25 and A p = 0.1. Ordinate units are arbitrary.

48 f(t) 100e0 o 0=6.0 1 ^o8.0 40... 32 0 24 >116 8 8- ^" ^ -- 10 0.1 0.01 1.0 0.1 0.01 Fig. 2 Effect of increasing p.o from 6.0 to 8.0 in the analysis of a single-component decay curve. X = 0.02, xo = 7.0, A x: 0.25, and A'i = 0.1. Ordinate units are arbitrary.

-1748r f(t ) OOe-0.02t EXPANSION OF PEAK FOR o =8.0 =9.0 40 32 < 24 16 -8 -~8 0.04 0.03 0.02 0.015 1.0 0.1 0.01 x Fig. 3 Detail of the peak shown in Fig. 2 for o0 = 8.0, and the results for the analysis of a single-component decay curve with po = 9.0. X = 0.02, xo = 7.0, A x = 0.25, and A i = 0.1. Ordinate units are arbitrary.

-18be seen that the error ripples are no longer symmetrical with respect to the true peak. In order to obtain finer resolution without increasing the height of the error ripples, it would be necessary to use better initial data and/or perhaps a more accurate integration scheme such as the one due to Filon.(ll) Although it was not checked experimentally, it may only be necessary to extract more points from the initial data by some interpolation procedure so as to be able to extend o to a higher value without introducing more error. Figure 4 illustrates the case of a two-component curve, with Po increasing from 3 to 6. In both curves the principal peaks are of the same height. This is because the coefficient of the second component, which is a factor of 10 smaller than that of the first component, has been divided by a X which is also a factor of 10 smaller than the X for the first component. Furthermore, the breadth of each of the two peaks is the same. This means the resolution is constant over the entire range, due to the fact that the method treats the initial data as a whole and minimizes the error uniformly over the entire curve. This fact is useful in analyzing unknown curves. For example, if two components have X values so close together that the resultant peaks cannot be completely resolved, the peak representing their sum may be wider than would be expected for a single component peak. The symmetry of the peak would give a rough indication of the relative amounts of each component. Consider Figure 5. In this two-component curve we would hope to see the peak for the shorterlived component (X=0.02) rise to one-half the height of the second peak. (11) L. N. G. Filon, Proc. Roy. Soc. Edin., 49, 38 (1928-29).

*'B~~l;Tqa aae sctmn @qxUiTpJ'1T0 = 7 - V puS'g'G = x O'L = Ox'TO'O = 3 T'0 Ty -GAfnD3 -60a9p;uUuoduioo-oA; e JO STSITBUB aqq uT 0'9 o9 0'o oz2 r0 BuOsla 9 ouT jo %so i JJa; *' 100'0 10'0 1'0 01 1000 0 10' 0 01 091 91 ll11 OF li lull I*lo~o'001 +,I~o-@~~l-(XJO l -61-o.Q^Ow^ re" " 0^^'s ________ ^(ro^O~~~l ^ ^o^OO~~l^^_______~ -6-~

-2096- f(t)= 100e-2t + 100e- lt 80- O =6.0 =8.0 64 o<48 532 16-1 16 -16- - - 1.0 0.1 0.01 1.0 0.1 0.01 Fig. 5 Effect of increasing ~o from 6.0 to 8.0 in the analysis of a two-component decay curve. X1 = 0.02, X2 = 0.01, xo = 7.0, A x = 0.25, and A p = 0.1. Ordinate units are arbitrary.

-21With.o=6 the peaks are not separated, but the sum peak is noticeably distorted to the left. With o=8, however, the peaks become separated and appear with about the expected relative heights. Because of the close proximity of the two peaks, the X=0.02 peak has been shifted very slightly to the left giving rise to an apparent X value of slightly more than 0.021. This r 5%o error could be reduced by extending o0 to a slightly higher value or perhaps by an analysis of the profile of the double peak. In Figure 6 the results of the analysis of a three-component decay curve are displayed. Again the increase of io from 6 to 8 greatly improves the accuracy of the final results. Finally, in Figure 7 the results obtained from the analysis of a four-component curve are shown. In this curve, the X values range from 0.5 to 0.01 while the coefficients range from 3750 to 100. In Table II, the relative height of each component peak [g(X)/Xi] and the associated i value as determined graphically from the data in Figure 7 are compared with the actual values. The deviations from the true values give an indication of the accuracy of the method. Here four components were determined from only 29 points on the gross decay curve in the range from t=l to t=1100, together with the associated points in the range from t=0 to t=l which would normally be found by interpolation. By using more points in the initial data and by extending the limit pi slightly, a more accurate determination of the parameters would have been possible. The relatively large error ripple peak, which occurs at a X value of about 0.0035, can be distinguished from a true peak by several means. First, the original data can be checked to see if a component with a halflife roughly three times as long as that of the X=0.01 component is

-22-O.It -0.02t -0.01t 96 f(t) =100Oe + OO + lOOe 96 - Io = 6.0 1 LO 8.0 16 0 48' -16 11111 I til I I I 1.0 0.1 0.01 1.0 0.1 0.01 x Fig. 6 Effect of increasing po from 6.0 to 8.0 in the analysis of a three-component decay curve. k = O.1, X? = 0.02, X3 = 0.01, o = 7.0, A x = 0.25, and A u = 0.1. Ordinate units are arbitrary.

-23-0.5t -Ot -O.-02t -O.01t f(t)=3750e +1000e + OOe +lO00e 96,uo =8.0 80 64 < 48 -16 -,,,,,,,, 1,m,,,,,,, II I I 1.0 0.1 0.01 x Fig. 7 Results of the analysis of a four-component decay curve with ~o = 8.0. xi = 0.5, X2 = 0.1, X3 = 0.02, X4 = 0.01, x0 = 7.0, A x = 0.25, and A [ = 0.1. Ordinate units are arbitrary.

TPBLLE II COMPARISON OF THE RESULTS GRAPHICALLY DETERMINED FROM FIGURE 7 FOR TEHE FOUR-COMPONENT DECAY CURVE WITH THE ACTUAL VALUES HEIGHT OF THE X = 0.1 PEAK NORMALIZED TO 1.000 g(Xi)/Xi IS TERMED Hi IN TABLE H1 H2 H3 Hi Xi1'- 3 Actual Value 0.750 1.000 0.500 1.000 0.500 0.100 0.02 0.010 Graphical Estimates 0.769 1.000 0.466 0.988 0.516 0.0985 0.021 0.0095 (Figure 7) Percent Deviation from Actual Values +2.5 0 -6.8 -1.2 +3.2 -1.5 +5.0 -5.0

-25reasonable. Next the final integration limit uo can be changed and the shift in position of the error ripples noted. Finally, the width of the suspected peak at its base [where g(X)/x = 0] can be compared with the width of true peaks, since error ripples have always appeared narrower than true peaks. The 5 5% shift in the apparent X value of the X=0.01 peak was not encountered when the X=0.02 and X=0.01 two-component curve was analyzed, but did occur in the results of the three-component analysis. Since experience has shown that a 5% shift is relatively large for a major peak, it is possible that the data for the X=0.1 decay curve is slightly in error. The range in X can easily be extended to cover a much greater spread in half-lives provided the coefficients Ni tend to decrease as the Xi values decrease. What is important here is the relative values of the quotient Ni/Xi. Although in Figure 7 the Xi and the Ni values vary by factors of 50 to 37.5 respectively, the quotients Ni/Xi, i.e., the heights of the peaks, differ at most by a factor of 2. The maximum acceptable difference in Ni/Xi will depend on the integration limit po and on the separation in Xi values. For example, with similar i values such as 0.02 and 0.01 in a two-component curve a difference in peak heights of 3 or 4 might be the maximum variation which would still permit a reasonably accurate determination of the parameters Ni and X for a pL of the order of 8. On the other hand with widely separated i values such as 0.5 and 0.02 in a two-component curve, the peak heights might differ by a factor of 10 and still allow a useful solution to be found for a p. of the order of 8. Obviously, the limiting factor is the relative heights of the error ripples compared to the peak heights.

-26B. Effect of Cut-Off with Respect to x In all cases, the experimentally determined function f(t) must be cut off at some finite value of t. Actually, under the change of variables, we are interested in the function exf(ex) versus x. Consider the case of the single-component curve with X=0.02. In Figure 8 we show a plot of the function eXf(ex) and two cases of cut-off. In the first case IXol = 5.25, while IjXl = 6.25 in the second case. It should be kept in mind that equal intervals in x correspond to exponentially increasing intervals in t. While little information is lost on the left, a considerable amount can be lost on the right. Even with a cut-off at jXol = 7 the amount lost is not negligible for high po values (po > 9) since the curve actually extends to x + oo. Figure 9 shows the results of a cut-off at IXol = 5.25. When the final integration is carried out to o=6 the error ripples completely mask the true curve. The dark triangle shows the expected height of the true peak. However, even in this poor case where the initial data have been cut off after about 5.5 half-lives, excellent results can be obtained if we are willing to accept poorer resolution. By restricting po to 4 the error ripples are greatly reduced and the true peak appears at the proper X value. Figure 10 illustrates the case of a cut-off at IXol = 6.25. For po-6 there is no apparent deviation from the results shown in Figure 2 where the data was cut off at xo = 7. Only when Io is increased to 8 do the error ripples appear larger than they were in Figure 2. Even in this case the true peak falls at the proper X value.

-27e f(e )- t [00e-0.0 ] 1000 CUT-OFF AT / CUT-OFF A -' Xo=-5.25 / Xo05.25 00 I I CUT-OFF I AT CUT-OFF 1I I /0^ i X -I6.25 -6 -4-2 0 2 4 6 8 10 X= Int Fig. 8 Plot of the function' exf(ex) versus x for a single component decay curve with 1 = 0.02. Cut-offs at xo = 5.25 and. x = 6.25 shown. curve with X = 0.02. Cut-offs at xo= 5.25 and xo= 6.25 shown.

aJe sqmur a;'uTTpo 0 I G-= r V puu8'ZO - X V z'O - X ~0o o. 0 09 UO~jo SuTseaJsZp o t q;TX^ GAJm Xeop quauoduloo -aiuTs 3 Jo sTsXTSou at; uT' = U~x %;B JJo-;no Jo 0aojj 6'T. 100'0 10'0 10 0O'1 1000 10'0 1' 0'1 ni- "gl - nil n i i i n l" i Sl91 —? -09or —IC V II; VI n iQo 9 5 0?-oz 0 01 — -' 91 5 09 fl- -og | SZ =~X l~t J H1IM) o a0 = lo I 0o- 09 -08 9c? i -i001 -V 0~- -091 g'=~XiV U-flO-ino HUIM) eOOI (eooi( (- O 4 0'0-*

-29- 0.02t f(t)= OOe (WITH CUT-OFF AT X(6.25) 48 F- =6.0 o L 8.0 40 32 24 16 0i -8 1.0 0.1 0.01 1.0 0.1 0.01 Fig. 10 Effect of cut-off at xo = 6.25 in the analysis of a singlecomponent decay curve with o0 increasing from 6.o to 8.0. X 0.02, A x = 0.25, and A A = 0.1. Ordinate units are arbitrary.

-30C. Effect of Poor Extrapolation of the Initial Data Since a large cut-off in the initial data can exhibit such a profound influence on the final results it is naturally of interest to see if it would be advantageous to extrapolate the initial data, even if the extrapolation is somewhat in error. In Figure 11 we show the true decay curve for a single component with X = 0.02 together with two rather poor extrapolations, where in each case the accurate data ended at IXol = 5(t=148.4). The final results obtained for each case appear in Figure 12 and may be compared with Figure 9 where the initial data was cut off at IXol = 5.25. In the case of the high extrapolation the initial data appears to contain an additional component with a X value less than 0.02. Thus, for uo=6 we see a double peak near X = 0.02. The presence of the additional ficticious component has shifted the position of the true peak to a larger value, and has caused large error ripples at the beginning of the curve due to the fact that the initial data in the range from t=O to t=148.4 did not contain the additional component. Although badly distorted, the results are still of much more value than if the extrapolation had not been made, as in the case of the cut-off at IXol = 5.25. At a sacrifice in resolution the high extrapolation data yield excellent results for a -0=4. In the case of the low extrapolation and po=6 we again find that even a poor extrapolation is much better than a large cut-off. With o=6 the error ripples are large but the true peak appears at just about the proper place. Of course, by decreasing So much better results could be obtained. It may be said, then, that even a poor extrapolation will usually yield better results than if the data had not been extrapolated at all, but care must be taken not to try to push the resolution beyond the accuracy inherent in the data.

-31100 ~fit)= looei ^ -0.02t f(t)= 1000e t DATA ENDS HERE 10 HIGH - y EXTRAPOLATION \ / ACTUAL CURVE LOW l EXTRAPOLATION'" F \\ \\' \\\ 0.1 0.001 \ \ curve for X\ = 0.02. f~~~~~' Fig.~~~~~~~~~''4Hg n o xrpltos ftesnl-opnn ea c~~ve for ~. = 0.02.'4'

-32-O. 02t f(t)=1OOe (WITH POOR EXTRAPOLATION) 56 HIGH HIGH LOW EXTRAPOLATION EXTRAPOLATION EXTRAPOLATION 48-.. o "=6.0 o=4.0 o 6.0 1.0 0.1 0.01 1.0 0.1 0.01 1.0 0.1 0.01 24 -- -- Ij I ^ -24.... 1.0 0.1 0.01 1.0 0.1 0.01 1.0 0.1 0.01 X Fig. 12 Effect of poor extrapolation of the initial data in the analysis of a single-component de'cay curve. X = 0.02, xo = 7.0, Ax = 0.25, and A p = 0.1. Ordinate units are arbitrary.

-33D. Effect of Scatter in the Initial Data Since we are dealing with a Fourier-type analysis we expect the method to exert a smoothing effect, in the least square sense, on the data. The greater the smoothing effect the more statistical scatter can be tolerated in the initial data. To investigate this effect we have constructed a curve, using the X=0.02 single component data, which contains a certain amount of more or less random scatter. Table I lists the percent deviation and the sign with which each point of the accurate data was changed to produce the scatter effect. While the scatter so introduced may seem large with respect to certain types of actual experimental data, it is just those cases where accurate data are not available that are in most need of a method of analysis. In Iigure 13 the results obtained with the poor data are shown for No values of 6 and 4. For jo=6 large error ripples are produced but the true peak still appears at the proper X value. By reducing po to 4 relatively excellent results are obtained due to the reduction in the error ripples. It is anticipated that smoothing of the data prior to analysis would have appreciably improved the final results. It should be possible to successfully analyze decay curves containing considerably more scatter provided that the data does not contain a bias and that a sufficient amount of initial data is available to provide a basis for the smoothing effect. E. Effect of Decreased Accuracy in the Numerical Integration In the results presented so far, we have used steps of A x = 0.25 in the numerical integration to obtain F(n ), and steps of A h = 0.1 in the final integration. Several curves have been run, however, with A x=0.5

-34-0.02t f(t)= 100 e (WITH SCATTER) 32 16 <8 X^0=6.0< o=. 24- -, -16 1.0 0. 0.01.0 0.1 0.01 Fig. 13 Effect of scatter in the initial data in the analysis of a singlecomponent decay curve. X = 0.02, xo = 7.0, A x = 0.25, and A = 0.1. Ordinate units are arbitrary.

-35and/or A p = 0.2 or 0.3. As might be expected, for relatively low po values (up to 5 or 6) little if any effect could be detected in the results for single- or well-separated two-component decay curves. For 1io values near 8 or 9, however, the increase in A x did produce a noticeable increase in the error ripple while an increase in A p. still remained relatively ineffective. VI. DISCUSSION As was mentioned in the introduction, it is hoped that the method that has been described here may be of use in the analysis of data from a wide variety of experimental problems in the physical, chemical and biological sciences. The ultimate accuracy will, of course, depend upon the data available. For example, when analyzing a mixture of radioactive isotopes it is usually possible to accumulate a large amount of accurate decay curve data. This is one case where a high level of resolution in the final results should be possible. It is interesting to note that when the usual radioactive decay rate data are analyzed, the heights of the resultant peaks are proportional to the number of atoms of each species. The height of a peak will be proportional to the decay rate of that species at time zero divided by the decay constant, therefore, we have N = (dN/dt)/X since dN/dt = NX for a first-order reaction. Another case where good resolution should be possible occurs when the measurements may be repeated as often as desired and the results from each set of measurements averaged. For example, electrons may be distinguished from heavy particles such as protons or alpha particles by differences in the decay times of the fluorescence produced by each in certain scintillators. Oscilloscope presentations of the pulses may be

-36photographed as often as desired, thus increasing the accuracy of the experiment. Dielectric relaxation properties also lend themselves to repeated measurements. Certain types of chemical kinetics data, on the other hand, or measurements involving the rates at which.injected materials disseminate in a living organism cannot always be obtained with great accuracy. The meager or poor data often, however, represent reactions in which there are only one or two components. Hence, a lower level of resolution may give a satisfactory solution that would not be available using other types of mathematical approaches since the detrimental effects due to the non-orthogonal properties of the exponential series can, to some extent, be avoided using the present method. A few remarks of an empirical nature on the subject of errors are applicable here. With reasonable resolution, the function g(e-Y) versus y will produce symmetrical peak profiles for true components. Furthermore, the breadth of the true peaks will be independent of y. Therefore, an error estimate may be obtained from the peaks themselves, and if a peak appears unsymmetrical or wider than another peak the presence of unresolved components would be indicated. The positions of the peaks in the error ripple depend upon u., while the positions of the true peaks do not. This fact can be used by merely carrying out the final integration to two different maximum Lo values and noting which peaks shift position. Also, the width of the base of a true peak is wider than the base width of an error ripple. Another test depends on the regular damping of the amplitudes of the error ripples. A divergence from the regular trend signals the presence of a true peak. Finally, it may be possible to obtain an error estimate for the integration scheme and also for the cut-off error. When there is doubt about a particular component,

-37calculated decay curves can be constructed, both with and without the suspected component, and the scatter of the experimental points about the calculated curves may be examined statistically with respect to "goodness of fit". The fact that the results of an analysis are available as a functional display particularly at several levels of fineness of resolution, rather than merely a set of values for the parameters, is quite advantageous. Although the determination of the number of components in principle is reduced to counting the number of true peaks, in difficult cases it is desirable to view the solution as a whole, in reference to the mutual interaction of all of the parameters. A study of the curves representing the final results at several levels of resolution can provide a sounder basis for the application of human judgment in cases where there is only a small amount of one of the components or where two components have similar X values. Another advantage is that the initial data are not required to be as accurate as in other methods, and full use is made of the accuracy that is available since the data are treated whole, as opposed to "subtraction-type" methods wherein all but the shortestlived components are determined using fewer points than are actually available. Furthermore, the occurrence of half-lives very close together in magnitude does not endanger the entire solution as in other methods. The major limitation of the method appears to be that information concerning the decay of the longest-lived component must be available for many half-lives. This limitation may not be too serious and several methods for treating the data in this regard have been suggested.