Machine Computation of Mass Spectrometer Analyses Triangular Inverse Method HARRY F. HOPP and RICHARD WERTZLER Sinclair Research laboratories, Inc., Harvey, 111.
b The calculation of percentage composition of gas samples run on the mass spectrometer i s a suitable application for electronic data-computing machines. Of the methods in common use, calculation with an inverse matrix is the fastest. The conventional inverse method, however, leads to certain errors inherent in the use of physical data, which are reduced b y using a modification called the triangular inverse method. This method solves for each successive component b y means of successively fewer elements. At each stage of the calculation spurious negative quantities are set equal to zero and d o not adversely influence subsequent results. The large relative errors arising from subtraction of one large term from another are substantially reduced. Total machine time, including reading, calculation, and punching, for a 21-component mixture i s 50 seconds on the IBM 650.
T
case for applying rapid dataprocessing machines to the calculation of mass spectrometer analyses was eloquently expressed by Dudenbostel and Priestley (1) in their discussion of the card programed calculator. McTeer and Morgan ( 2 ) have described in detail a method for applying the I B M 650 stored program computer to the simultaneous calculation of 33 components in gas samples run on the mass spectrometer. The speed and versatility of the I B N 650 permit a choice among several good methods of calculation: Each method in common use has advantages for particular cases, but any method chosen must be capable of solving simultaneous linear equations. These equations arise in the following manner. The output of an analytical mass spectrometer is commonly manifested as a sequence of discrete peaks observed as transient deflections from a fixed base line. Each peak can be uniquely identified with an integer, i. The height or magnitude, P,, of any such peak, i, produced by a substance identified by the integer, j , is proportional to the quantity, xi, of t h a t substance introduced into the mass spectrometer. The proportionality constant is known as the sensitivity coefficient and may be represented by a$,. HE
Experiment has shown that each substance in a mixture makes its contribution to a peak independently of the others present; these contributions are additive. Thus, if a mixture contains n components, n-e have the equations:
Quantitative analysis of the mixture requires evaluation of the x’s and can be accomplished if n independent equations of the above type can be written. In these equations, the a’s are known from calibration runs on pure samples of the suspected components, and the P’s are obtained by running the unkno\m mixture on the mass Spectrometer. The process of solving these equations is lengthy, when each sample is treated as an individual case without the aid of a high speed computer. Therefore, it has been common practice t o use matrix inversion to simplify and reduce the work of sample calculation done by hand. The aiiJs can be arranged in a matrix, [ i l l . Then, the inverse matrix, [ B ] , consists of the elements b,;, which are those needed in the following equations explicit for z:
Each x can be calculated as the sum of n terms involving the P’s and b’s, all of which are known. The inversion process, though lengthy, is required only as a consequence of new calibration runs. Calculation of individual samples is reduced to a process requiring but little skill and minimum time. LIMITATIONS OF SQUARE INVERSE SOLUTION
This process is easily adapted to machine calculation. However, it has a disadvantage in any application to actual physical data where the a’s and P’s are subject to ordinary errors of measurement. Furthermore, the a’s in the mass spectrometer problem are not truly constant. They show random fluctuations and a long-term drift between calibrations. These variations can cause appreciably larger errors in
the 2 s than other methods which do not use the inverse matrix. In the inverses normally used, many of the elements are necessarily negative, and the negative elements frequently approach or exceed the positive elements in magnitude. When one large term is subtracted from another, the relative error becomes magnified. Therefore, a small instability or error in one of these terms can produce a large relative error in the 1: calculated from these terms. One manifestation of this effect is the appearance of negative z’s n hich have no physical meaning in our problem, A negative z can be reported as zero, which is nearer the truth. However, its influence is felt in other components as well. This adverse influence can be appreciably diminished by use of a triangular inverse. I n calculating the z’s by using the square inverse, each .z is calculated independently of the others, as the sum of n products, some of vihich may be zero. The same peak heights are used for every component. At the end of the calculation any negative zJs are set equal to zero before normalization, but there remain compensating positive errors that are undetected. I n using the triangular inverse, homever, the 2’s are calculated in a definite sequence. As each z is calculated, its effect, as expressed by the original array of sensitivity coefficients, is subtracted from the peak heights before the next 2: is calculated. A negative value of z is set equal to zero as soon as it appears, and its effect is not transmitted to the 2)s to be calculated later. Furthermore, as the peak heights are corrected, if one becomes negative, its Talue is taken as zero for all subsequent calculations. These two procedures give the triangular inverse its superiority. It is somen hat inferior in speed. That the negative coefficients in the triangular inverse are less in both total numbers and magnitude is apparent from Table I. The mass spectrometer gas analysis problem is particularly amenable to this treatment; thus fewer positive coefficients are required in the triangular solution. OBTAINING TRIANGULAR INVERSE
The coefficients of the triangular inverse may be obtained during the VOL. 30, NO. 5, MAY 1958
877
EXAMPLE
Table 1.
Number and Magnitudes o f Elements in Square Inverse and in Triangular Inverse
K‘umber of Elements Positive Segative TriTriSquare angular Square angular 8 0 9 1 10 0 7 1 10 0 6 1
11? 3-* CH,
co .
8 8 9 11
C2-
C2 CO‘
Air
6
Cq-
6 6 7 7 8 r
C; C4i-Cn n-Ca i-Ci
8
n-Cs Cs=C6-
3 1 1
n-C6
118
:
1 1
8 8 r
1 1 1
5 r
2
1
0 2 0
0
10 10
4
0 4 0 4 6
r
1 6
; 8
5 6
-
9
r
6 G
i
3
3
5 1
1 1 -
0
1 33
1 ~
44
130
Sum of Elements Positive Negative TriTriSquare angular Square angular 3788 3582 355 0 7333 1598 45090 0 2160 1902 1351 0 29976 28289 668h 8% __. 5i10 2823 7296 0 5624 5392 728 212 5418 1535 2071 0 983.1 9821 658 0 21510 2800 15821 0 7382 1526 14067 1631 3830 3258 7338 0 1922 486 26241 291 16414 12598 860 439 3870 3293 28700 2i0 25321 14610 1762 1718 3772 3T72 346 48 4522 4522 6 0 l04!)9 10499 25 25
>
c
All
L.
Figure matrix
1.
Augmented
’
third-order
process of matrix inversion. The way in which the coefficients appear will depend on the particular scheme of matrix inversion used. The method used here is a form of Gaussian elimination (S), in which the pivotal elements are chosen in sequence and no back solution is necessary. One row of the triangular inverse appears after each transformatioii. If the original array of calibration coefficients is augmented, as s1ion.n for a third-order matrix in Figure 1, a single formula can be used to compute each element in each transformation.
(3) I n Figure 2, which represents the results of the first transformation, the circled element is taken as the first row of the triangular matrix. It is used to calculate the last conipoiient in the analysis, after corrections have been made on the corresponding peak height for all other components in the sample. The elements shoiyn in Figure 3 are calculated from those of Figure 2 by a second transformation using Formula 3, but replacing A by A’ and A’ by A ” . The circled elenients of Figure 3
878 *
I
*;3
ANALYTICAL CHEMISTRY
Figure 2. Augmented matrix after firs+transformation
To evaluate the extent to which the use of the triangular inverse might iniprove the accuracy of the results, a set of ten analyses \\-as made up, for n hich peak heights nere calculated. Each peak height n-as then perturbed by adding to it a random normal deviate, drawn from a population with zero mean and unit variance. Each analysis was so treated with eight different sets of nornial deviates, resulting in 80 “samples.” These perturbed samples n ere then calculated using the square inverse, and the triangular inverse. The accuracy of the results was evaluated by taking the differences of each calculated percentage from the ‘[true” percentage, adding the squares of theqe deviations, and dividing by the degrees of freedom. This result is the variance. The square root of the result is the standard deriation. The total variance using the square inverse was 0.346; using the triangular inverse it n a s 0.219. This reduction in variance. tested by the F ratio, is significant a t the 1% level. There are large differences in the effect for different components. Table I1 lists the coniponents in order of variance using the square matriu, and gives the standard deriation using the square matrix and the triangular matrix. -1s a matter of interest. a standard deviation of 0.04 represents the level of rouiidoff error in this calculation. PROGRAM
-I
0
0
01
Figure 3. Augmented matrix after second transformation
, Ai’,
A32
Figure 4. Triangular inverse matrix
become the second roil (next to last component) of the triangular inverse. Similarly, succeeding transformations are carried out-equal in number to the number of equations in the original set. The final triangular inverse is assembled in Figure 4. I n the original array, columns correspond to coniponents; in the inverse, rows correspond to components.
Data input t’o the 650 for each sample is represented by a deck of I B l I cards. Peak heights are punched one to a card (Spectro-Sadie or keypunch) nith attenuation, peak number, and sample number. Peak cards not required in the calculation are automatically ignored by the 650. X “Header” card is merged n-ith each sample deck before re:idin to provide additional indicative information, such as date and saniple pressure, and t o signal the type of calculation to be performed. The latter function is necessary because several arrays are available in the program deck. To calculate a group of samples, the sample cards are stacked in any order behind the program deck and the necessary snitches set to initiate the calculations. S o further attention is required until the output cards are removed from the hopper to be listed for reporting. \Then a Header card is reachedidentified by 2-punch in column 6the observed pressure is calculated and punched on an identification card along iritli the other identifying data such as complete sample number and date. Also, certain sn-itches are set n ithin the program to cause the par-
ticular type of sample to be calculated properly. Immediately thereafter, the peakheight input cards are read. As each is read, it is niultiplied by its attenuation factor, corrected for backgrountl, ant1 storcd in tn-o appropriate locat’ions, one for calculation and t h e other for readout. Data cards may be in any order within a given saiiiplc deck. \Then the Header of the nest sample is sensed as the nest card to he read, c:tlculation of the composition of the sample begins. The calculation begins with the summing of the appropriate products of invcree clemcrits by peak heights t’o oht’ain the partial pressurc of the first c~~mpoiicntto lie calculated, zls. If x18 is positive, all peaks are diniinished liy tht. terms x38a,.18 to yield P’%. .cli is calculntecl in the same way’ cxccpt that the P’,’s are used in place of the Pi’s and onlj- I I - 1 ternis are required. -411 pPaks are then corrected for T17 if positive. If any corrected pcak beconies negative, it is set equal to zero for further calculation. This process is caont,inucd until all the 5‘s have heen calculated. Sormalixation of the partial pressure yields mole or gas volume per cents. If w i g h t a d liquid volume per cents are required, the appropriate inolccular weight and molecular volume arc applied, follo\\-ed again by normalization. .It the same time residual peak heights are availnhlc for all peaks used in the calculation as 11-ell as for threci cherk 1)eaks. These residual peak h(,ight?,
~~
Table II. Standard Deviations of Test Analyses b y Square Inverse and b y Triangular Inverse
Component Isopentane Isobutane n-Pentane Sitrogen Butenes Pentenes Propene n-Hexane n-Butane Air Carbon monoxide Propane Ethane
Ethene Hydrogen Hesenes Met hane Carbon dioxide -1verage
Standard Ileviation Square Triangular inverse inverse 1.25 1.16 0 84 0.76 0 . 2 0 65 0 58 0 46 0 46 0.31 0.24 (1.22 0 is 0 lti 0 15 0.12 0.1“ 0.04 0.50
1.04 0 95 0.88 0.23 0.42 0 33 0 41 0 -17 0 42 0.29 0.26 0.21 0 18 0 10 0 12 0.12 0.08 0 05 0,4i
nliicli may be positive or negative, are available in determining the validity of the separate percentages and in identifying the reasons for poor analyses. The r e d u a l peak heights may he calculated and punched out by using either the square inverse or tlie triangular inver-e method. These values arc very useful in spotting erroneous aiialy.es, and often in detecting a d corrwtirig a n crroneous peak. But RS red u a l peak heights must be calculated us a n additional operation after determination of partial preisures 11 hen :I q u a r e inverse is u w l , niuch of tlie
potential speed advantage is lost. No additional time i b required to obtain residual peaks I\ hen the triangular inverse is used. At this point, all the desiicd figures being available, a data card is punched for each component. Each data card contains identification, three kinds of percrntagw, aiitl a comparison of residual peak height n ith the original peak height of the coniponciit’s principal peak. -4 totals card is punched for each saiiiple. The G50 then proceeds to the nest sample. Each sample ir calculated in the proper iiianner until the “Read-Feed” is enipty. Cntire time to calculate 21 conipoiients and three additional check peaks is 50 sccondq. ACKNOWLEDGMENT
The subbt:iiiti:il coriti ibution. iiiadc hy Grace 11. Grant and George W. C‘leary in applying this method to niabb spectrometer calculation^ are gratefully acknon ledgcd by the authors. LITERATURE CITED
(1) Dudenbostcl, B.
Jr.,
.hAL.
F.,Jr., l’riestleJ-j W., CHEJI.
1275-8
26,
(1954).
(2) McTeer, 11. I,., l l o r g a n , G . S., Jr.,
Union Carbide Chemic:& Co., unpuhlished report. ( 3 ) Satl. Bur. Standards, “Contributions to the Solution of Syytems of Linear Equations and t h e Lktermination of Eigenvalues,” .4ppl. ?\lath. Series 39, pp. 12-14.
RECEIVED for review December 14, 1957. Accepted March ti, 1958.
Automatic Computer Program for the Reduction of Routine Emission Spectrographic Data FRANK W. ANDERSON and JAMES H. MOSER Martinez Research laboratory, Shell Oil Co., Martinez, Calif. T ,o relieve the burden of manual computation, a general-purpose, medium speed, digital computer program has been prepared for the conversion of emission spectrographic film line transmittances to element concentrations. Estimates of the precision of each analysis are also computed, the kind of estimate being dependent on whether the input data are from a semiquantitative or a quantitative analytical procedure. The program also detects and signals the occurrence of certain types of errors in input and rejects such data without stopping the computer. Both speed and reliability of spectrographic calculations have been improved b y this program.
for an I B l I Type G50 computer with 2000 ten-digit words of memory is in use for conversion of emission spectrographic film line transmittances to element concentrations in samples. Such conversions have been performed previously by use of graphical methods or a “calculating board’’ (a special-purpose slide rulc). Harvey ( 1 ) has discussed thcse methods and has given instructions for the preparation of the graphs and the acales for use with the calculating board. The present digital computer program follows the principles of the mechanical method of calculation t o perform the desired data transformations. PROGRAM
PRINCIPLES OF CALCULATION
Given the “ganinia” curve of the film relating per cent transmittance of a n image to relative intensity (Figure I), the transiiiittance of the line of the element in question. T E , and that of the internal standard, TR,are converted to the relative intensities, RE and RR, respectively. The relative intensity ratio I = KE/RR
is computed. Referring to Figure 2, the relative intensity ratio is then converted to the estimated concentration of the element in the sparked sample. Finally, the average coacenVOL. 30, NO. 5,
M A Y 1958
a
879
