ENGINEERING, DESIGN, AND EQUIPMENT

Refractive Index Curve. Small weighing bottles of about 30- calibration points from the line was within the limits of accuracy
ml. capacity, with female-joint tops, were found convenient for (0.05 O C. ) of t h e temperature-measuring equipment. Subse-
weighing out the liquids in the determination of t h e refractive quent voltage readings were corrected b y adding t h e AE for
index curve. T h e fact t h a t the grinding on the bottle proper each When the boiling point of the carefully puSified
was on t h e outside prevented accidental wetting of t h e ground toluene was measured, it was 110.59",as compared t o a National
portion as t h e liquids were poured in from small graduated Bureau of Standards value of 110.62'. This serves as a further
cylinders. The smallest quantity of a n y one liquid weighed in indication of the accuracy of the temperature measurements.
these determinations was 1 gram, and because the balance was
accurate to 0.0002 gram, t h e combined error for any sample could Literature cited
not be over 0.04%.
Pressure-Regulating System. The barostat used was that Bromiley, E. C., and Quiggle, D., IND.ENG.CHEX.,25, 1136
designed by Willingham and Rossini ( 2 0 ) a t t h e National Bureau (1933).
of Standards. It was thoroughly cleaned, then heated under Brown, J., and Ewald, A . H., Australian J . Sci. Research, 3A,
high vacuum and filled with clean, redistilled mercury, a vacuum 306-23 (1950); 4A, 198-212 (1951).
being obtained b y carefully boiling out dissolved air under vac- Cornell, L. W., and Montonna, R. E., IKD.ENQ.CHEM.,25,
uum IThile the filling was in progress. The barostat was installed 1331-35 (1933).
in an air thermostat, the temperature of which was controlled t o Fowler, R. T., J. SOC.Chem. Ind. (London), 68, 131-2 (1949);
u-ithin 3" C. b y a Fenwall thermoswitch, which was surrounded Ind. Chemist, 24, 717, 824 (1948).
by the circulating air stream. T h e temperature of the mercury Garner, F. H., Ibid., 25, 238 (1948).
lagged behind t h a t of the air, and consequently did not change Gillespie, D. T. C., IND.ENG.CHEM.,ANAL.ED.,18, 575 (1940).
over as a i d e a range. 4 dibutyl phthalate manometer regis- Jones, C. A., Schoenborn, E. >I., and Colburn, A. P., IND.
tered the difference between the system pressure and t h a t of t h e ENG.CHEM.,35, 666-72 (1943).
atmosphere. B y reading this manometer and the barometer, Kortiim, G., Moegling, D., and Woerner, F., Chem. Eng. Tech.,
it was possible t o calculate the pressure within the apparatus. 22, 453-7 (1950).
Kitrogen was supplied t o the barostat through a diaphragm Natl. Bur. Standards, Circ. C461, 39, 43, 123, 126 (1947).
valve, from a standard nitrogen cylinder. When the mercury Othmer, D. F., IND.ENQ.CHEM.,20,743-6 (1928); Anat. Chem.,
column reached the upper contact of the barostat, a circuit was 20, 763 (1948).
completed which activated a relay t h a t opened a solenoid valve Regnault, Mem. mad. inst. France, 26, 727 (1862).
and allowed escape of nitrogen through a needle valve. All runs Rose, Arthur, Williams, E. T., Sanders, W. W., Heiny, R. L.,
were made a t 760 mm. of mercury T h e prevailing pressure in and Ryan, J. F., IND.ENG.CHEW,45, 1568-72 (1963).
State College is about 730 mm. The pressure variation due t o Rose-Innes and Young, Phil. Mag., (V) 47, 353 (1899).
the variation of the temperature of the mercury was about 0.1 Saddington, A. W., and Krase, N. W., J . Am. Chem. Soc., 56,
mm. This far exceeds pressure variation due t o the opening and 353-61 (1934).
closing of the solenoid valve. It was easy to adjust t h e dia- Sage, B. H., and Lacey, W. N., Trans. Am. Inst. Mining Met.
phragm valve on the nitrogen cylinder and the final outlet needle Engrs., 136, 13G-57 (1940).
valve, so that the pressure variation could not be detected on the Sieg, L., Chem. Eng. Tech., 22, 322-6 (1950).
dibutyl phthalate manometer. Steinhauser, H. H., and White, R. R., IND.ENG.CHEY.,41,
Temperature Measurement. The temperature-measuring sys- 2912-20 (1949).
tem consisted of a calibrated single-junction copper-constantan Swietoslawski. W., J . Chem. Educ., 5, 469 (1928).
thermocouple, a Rubicon galvanometer with a sensitivity of Williams, E. T., Ph.D. thesis, Pennsylvania State College,
1.1 pv. per mm., and a Leeds & Northrup Type K-2 potenti- University Park, Pa., 1952.
ometer, which could be read to 0.1 mv. Willingham, C. B., and Rossini, F. D., API Research Project
The calibration was accomplished b y reading the electromotive 6, December 1945.
force with carefully purified n-heptane, and with pure water Zawideki, J., 2. physik. Chem., 35, 129 (1900).
boiling in the apparatus. The difference between the observed
voltage and the voltage of a standard sample was taken a t each R E C E I V Efor
D review April 21, 1953. ACCEPTED February 2 , 1955.
point and plotted against the observed voltage. A straight Presented before the Division of Industrial end Engineering Chemistry a t
line was drawn from the zero point ( A E = 0, EOb4 = 0) the 123rd Meeting of the AMERICAN CHEMICAL SOCIETY,Los .4ngeles,
which passed between the plotted points. The distance of t h e Calif,

Applications of Matrix Mathematics to

Chemical Engineering Problems
ANDREAS ACRlVOSl AND NEAL R. AMUNDSON2
Departmenf of Chemical Engineering, University of Minnesofa, Minneapolis 14, Minn.

R ECENTLY, chemic81 engineers have become increasingly

aware of the usefulness of the calculus of finite differences
as a mathematical tool for solving problems in stagewise opera-
by Tiller and Tour ( 2 5 ) , Tiller (ZC), Mason and Piret ( l 7 ) ,
Lapidus and Amundson (16), and Amundson ( I ).
This article presents a method t h a t will enable the chemical
tions. This particular section of mathematics has remained engineer t o deal with more involved problems t h a t arise in con-
dormant for some years, because of the restricted use it has nection with the unsteady-state behavior of stagewise opera-
found in most branches of applied mathematics, where mostly tions; there has been, of late, considerable interest in these
continuous, rather than discrete, phenomena are studied. There chemical engineering problems. This new approach will make
are, however, many textbooks on the calculus of finite differences, extensive use of matrix algebra, the salient features of which are
including those of Boole ( 3 ) , Jordan ( I S ) , Milne-Thompson described.
( I $ ) , Wallenberg and Guldberg (271, and Norlund (21). Ap- It is surprising t o learn t h a t most results and theorems of prac-
plications of the calculus of finite differences t o a variety of rela- tical interest concerning matrices were discovered almost a
tively simple chemical engineering problems have been made century ago. However, their usefulness in many branches
1 Present address, Department of Chemistry and Chemical Engineering,
of applied mathematics was not appreciated until very recently.
Univeisity of California, Berkeley, Calif. The present interest in matrices was first stimulated by the work
* Present address (on leave), Cambridge University, Cambridge, England. of Born, Heisenberg and P. Jordan in quantum mechanics, a n d

August 1955 INDUSTRIAL AND ENGINEERING CHEMISTRY 1533

ENGINEERING, DESIGN, AND EQUIPMENT

it seems t h a t the aeronautical and electrical engineers were the where t h e ai's are scalars, is, of course, another matrix, P(A).
first among t h e engineers t o have employed matrices extensively Then
i n their researches. At present, of course, this interest is wide-
spread among applied mathematicians, as revealed b y the ap-
propriate literature.
The present article will, i t is hoped, serve two purposes-to
present t h e essential features of matrix algebra, and thus make is said t o b e a polynomial of A of t h e n t h degree. This, of course,
them known to a larger circle of chemical engineers; and t o is a n extension of the definition of polynomials for scalars. Note
illustrate the usefulness of this new approach by means of some that A", A"-',. . ., A2, A, Zare all of the same order.
rather interesting examples. N o general theorems will be proved, Power Series in Matrices. A natural extension of the above
except in some special cases, and constant reference t o the text- definition of a polynomial of a matrix t o include infinite series
books on the subject will be made. Some of the recommended in matrices can now be put forth. Let there be a sequence, Ao,
textbooks in English are those of Perlis (22), Ferrar (6),and AI, A2,. . ., of matrices t h a t are all of t h e same order, and let
the first chapter of Courant and Hilbert ( 4 ) , which are primarily P m

devoted t o the basic principles; those of Michal (18),Ferrar S, = Ai. Then the infinite series, Ai, is said t o be con-
0 0
(?'), and especially Dwyer (6) and Frazer, Duncan, and Collar
vergent if every element in S, converges t o a bounded limit as p
(9),which contain material of considerable practical interest ;
approaches infinity.
and t h a t of MacDuffee (16), in which most of the finer points
Of particular interest in applications are the infinite power
are examined. m

series in a matrix, A-that is, matrices of the form, Xa,Al,

Properties of matrices 0
where the a 1 ' s are scalars. The discussion of the convergence
Some of the simple rules obeyed by matrices, namely addition of such power series and of their summation is presented with t h e
and multiplicatiori, have already been stated by Amundson formula of Sylvester.
( I ) . I n the present article, use is made only of square matrices, One of the most important matrix power series is the matrix
which are denoted by capital letters, and column matrices, called exponential function, eA, which is defined, in complete analogy
vectors, which are represented by lower ca8e letters. Boldface with the scalar exponential function, b y the power series
type is used for matrices t o dist,inguish them from scalars. Ma- m
trices with real elements will be studied.
The following definitions apply t o a matrix:
eA=Z+Affl-+...
A2
=c- A"
m!

(211 (212 (213 . . @;) It can easily be seen from the above definition t h a t
(221 (222 (223 . .
A = [ u c ~=] . e*e-A =
(Un1 an2 an3 . 1 ann and
A square matrix is of order n if it has n rows and n columns. eAeB -. eA+B

The transposed A' of a matrix, A, is defined t o be the matrix

A' only if A and B are commutative matrices-that is, only if
which has its rows identical with t h e columns of A-Le.,
AB = BA. Other examples of matrix power series are the func-
= [aii].
The adjoint of A, adj A, is equal t o [ A i < ]where
, Aii is the tions, sin A and cos A, defined by
cofactor of aii in A . A3 A5
The matrix, I = [&,I, where 869 is t h e Kronecker symbol
sinA=A-
3!
~ +% ...
+
( S i j = 0 if i j and 6 ~ =
i l),is called the unit matrix, since
A2
A I = ZA = A
COS A = I - -
2!
+ 4! A4
- (4)

for an arbitrary A of the same order as I . If A is a square matrix of order n,the determinantal equation
T h e matrix, 0,all the elements of which are equal to zero, is
/ X I - AI = 0 (5)
called the null matrix, since
is called the characteristic (or secular) equation of A. X is a
A 0 = OA = 0 scalar, and b y expanding the above determinant it can be re-
again for a n arbitrary A. duced t o a polynominal in X o€ degree n. Thus,

++
The inverse matrix, A-1, of a matrix A has the property t h a t 1XI- A l ~ P ( X ) ~ X " + p i X " - ' + p z X n - 2 + ,..
AA-1 = A-'A = I pn-1X pa 0 (6)
where the p ' s are scalar coefficients. The algebraic equation,
A necessary and sufficient condition for the existence of A-' is ( k = 1,2,. . which are called t h e
Equation 6, has n roots,
t h a t the determinant of A, denoted b y 1 A I or det A, be different
characteristic, or latent, roots of the square matrix, A. The
from zero. I n t h a t case
Cagley-Hamilton theorem states t h a t P( A ) is the null matrix,
or, in other words, t h a t
P(A) E A" + PIA"..' + pzAn-2 + . . . + pn-jA + p n I = 0 (7)
A matrix, A, for which I A 1 = 0, is said to be singular. If A' =
where t h e p l ' s are identical with the corresponding coefficients
A, A is said t o be symmetric. If A' = A-I, A is said to be or-
of the characteristic equation. ilnother formulation of the
thogonal.
The matrix, B = C-'AC, where C is any nonsingular matrix, above theorem is t h a t any square matrix, A, satisfies its own
characteristic equation.
is said t o be obtained from A by a similarity transformation.
Sylvester Formula. The Sylvester formula, a n exceedingly
Matrix Polynomials. The finite sum
useful expression, can readily be derived from the Cayley-Hamil-
wA" + oriA"-l + . . . 4- or,-~A + ornI ton theorem and the Lagrange interpolation formula. It states

1534 INDUSTRIAL AND ENGINEERING CHEMISTRY Vol. 47, No. 8

 ENGINEERING, DESIGN, AND EQUIPMENT

that, if the characteristic roots of a square matrix, A, of order in the vapor leaving the nth plate, and z n ( i j , the mole fraction of
n are all distinct, and F( A ) is any matrix polynomial, or any con- the ith component in the liquid leaving the n t h plate, is given by
vergent power series, in A , Ydi) = Kn(i)zn(i) (12)
where the function, K n ( i ) has
, to be estimated a priori. The state
of the rectification column can then be described b y the following
equation, which can be obtained b y combining a material balance
around a n y plate in t h e tower with Equation 12:
Thus
Ln-izn-~(i)- [Ln + V n K n ( i ) l ~ n (+i )
(9) Vn+1K,+1(i)sn+.(i)= 0 (13)
for 0 5 n iN , b u t n # f, where f denotes the feed tray. I n
Equation 13, L, and V, denote the total moles of liquid and
from which it follows t h a t e A , is defined b y the infinite series,
vapor, respectively, leaving plate n. If z~(i) denotes the mole
exists for all square matrices, A, since ex exists for all scalars,
fraction of the i t h component in the feed, when n = f,
A. A similar argument can be used also for any power series in
A, and it follows t h a t a power series in A , F ( A j , converges only Lj-l~j-,(i) - [ L , + V / K / ( i ) l ~ t (+i )
as long as all the characteristic roots of A lie within the radius Vj+iKj+i(ihj+i(i) + Z F ( ~ )= 0 (14)
of convergence of F( A). Equations 13 and 14 can next be combined into a single matrix
The usefulness of the Sylvester theorem cannot be over-
equation b y defining t h e vectors
estimated. Fractional powers of a matrix can now be defined
without difficulty, and it can be stated that, under suitable con-
ditions, a n y analytic function of a square matrix, A , of order n
can be expressed as a polynominal in A of degree n - 1. -4s
a matter of fact, many functions of a matrix A can be formu-
lated satisfactorily only by means of the Sylvester formula.
When multiplicities in the latent roots occur,
the Sylvester formula has to be modified. Thus
if, for example, a set of roots, X1, b, . . . , A,, are ~ ( iE)
0 0 LZ
equal, Sylvester's formula can be stated in its con-
fluent form ( 9 )
where D is the nusnber of moles of overhead product. The matrix
equation, which takes the place of Equations 13 and 14, is

where &(A) = ( A -A, + 1) (X - A, + 5). . .(X - A,) and the sum A(i)x(ij + z(ij = 0 (15)
is taken over all the distinct roots of A . It is worth remembering and its solution is
t h a t the following identity exists
~ ( i= )-A-'(i)Z(ij (16)
adj ( A i l - A ) = (-l)n-l ?r (Ail - A) (11)
i #j I n practice it is best t o obtain zo(i)from Equation 16 and then
to calculate z,(i) for n > 0 b y iterating Equation 13. Example
Differentiation and Integration of Matrices. If u ( t ) is a square I11 will demonstrate the method of evaluating the determinant
matrix or a vector, the elements of which are fractions of a param- and the appropriate cofactor of a matrix of same form as A .
eter, t, the derivative of u ( t ) with respect to t is defined as the Suppose that, as a result of the preceding calculations, the
matrix obtained by differentiating each element of u(t). estimate in the function, K n ( i ) ,has t o be revised slightly so t h a t
the new matrix thus obtained is
thus
A i = A + B
where A is identical with t h a t defined previously and B is an-
other square matrix, the elements of which are much smaller in
similarly magnitude than t h e corresponding elements of A .
(A +
B ) - l = A-I(Z +
BA-l)-' =
A-'[I - BA-' +
(BA-')Z - . . . j
or, approximately,
I n addition, there exists a definition for the derivative of a
matrix with respect t o a matrix, but this concept will not be ~ i ( i=) - A i - ' ~ ( i ) ( I - A-'Bjx(i) (17)
used in this discussion. where ~ ( iis )given by Equation 16. This method appears at-
tractive, in t h a t once a good approximation to the answer has
Example 1. Inversion of a matrix by iteration been arrived at-that is, once A and A-I have been determined-
an iteration can be set up b y varying Bin Equation 17.
It often happens t h a t the investigator is able to discover a n
approximate solution t o a given problem, and then is faced with
Example II. Analytical solution for
t h e task of devising a simple iterative procedure which will equations of multicomponent rectifscation
converge and which, moreover, will improve the answer rapidly.
Such a n example is presented in the design of fractionating Some years ago Underwood ($6) and Murdoch (60),among
columns, if use is made of the Geddes and Thiele method (25). others, arrived a t an analytic solution for the equations of multi-
It is assumed t h a t the column has a total of N plates, including component rectification. Their solution can be derived b y the
the reboiler, and a condenser, the plate 0, and that, moreover, use of matrix algebra in a n elegant way.
t h e relation between yn(i), the mole fraction of the i t h component If constant molal overflow is assumed, constant relative

August 1955 INDUSTRIAL AND ENGINEERING CHEMISTRY 1535

ENGINEERING, DESIGN, AND EQUIPMENT

volatilities and a relation between x , ( i ) and y,(i) of the form 1 +1 a1

1 +11 an 1
1
A = 1 1 f a a .

i=l

th*en a material balance around an arbitrary section of the en-

= ?r
i
ai j l + 7d;i
riching column results in (23) Therefore,

i= 1
and the roots of A are identical with those of
where R is the reflux ratio, x o ( i ) denotes the composition of the
overhead product, m the total number of components in the mix-
ture and p ( i ) 2 [ a ( i ) ]-I, where a ( i ) is the relative volatility of
t h e i t h component. Equation 18 can now be linearized b y a
change in the dependent variable-that is, by defining a new It is, of course, very easy to establish t h a t the roots of Equation
function, X n ( i ) ,as 26 are all positive and distinct (20). Moreover, it can be shown
that

j= 1

in which case Equation 18 can be transformed into

where the representative term, bi, is equal t o

together with the condit,ion t h a t

and that, therefore,
The substitution given b y Equation 20 follows as a logical
consequence of the plate-by-plate method and will be described
in detail in a subsequent publication.
However, Equation 21 can be expressed as

where

and
Equation 28 can be shown t o be identical to the formulas derived
by Underwood (26) and Murdoch ( 2 0 ) . It is also apparent that,
in general, if n is small, Equation 24 is t h e more convenient form
of the solution, while if n is large, Equation 27 is preferable.

Example 111. Transient behavior

of some stagewise systems
X , is a vector and A is a square matrix. The solution of Equa-
tion 23 is Mathematical Model and Its Solution. The following section
considers, in some detail, a simple mathematical model t h a t ap-
pears wherever the unsteady-state behavior of absorption,
which expresses the result of t h e iterative method of solution, extraction, and distillation columns, or other kinds of stagewise
the so-called plate-by-plate method, in the convenient matrix systems with reflux, is investigated. For the sake of illustration,
notation. However, Equation 24 can be transformed by means a distillation column is considered, but, of course, this restric-
of the Sylvester formula so t h a t tion should not be interpreted literally.
If a material balance is made around a representative plate,
n, in the column the following is obtained:

i f it is assumed, for the moment, t h a t the characteristic roots of

A are all distinct. This equation can be simplified considerably.
By subtracting the first column from each of the others and then Again, as in Example I, i is the component-designating symbol;
expanding the resulting determinant, the determinant z n ( i , t ) and yn(i,t) are, respectively, the mole fraction of the ith

1536 INDUSTRIAL AND ENGINEERING CHEMISTRY voi. 47, NO. a

 ENGINEERING, DESIGN. AND EQUIPMENT

component in the liquid and the vapor stream leaving the n t h which is but a product of c square matrices. Similarly,
plate; L,(t) and V,(t) are, respectively, the moles of liquid and
vapor leaving the nth plate; while h, is the liquid holdup on the
nth plate. In constructing the above model, any time lags that
[M;(A)]-I = &1 [I- A(t,-,) :] (38)

are always present have been neglected because, of course, the

so that, in this way, Equation 33 can be evaluated numericaily
various streams are moving across the column a t a finite
with the aid of a sufficiently powerful digital computer, by per-
speed. These time lags would have to be taken into account
forming 2c multiplications of a vector b y a square matrix.
in any more exact calculations, because of the pronounced
However, the problem under consideration has certain simpli-
effects they have on the response of the column to an applied
fying features, because the vector, z ( t ) , has all of its elements b u t
change in the operating conditions. Another simplification
one equal to zero, unless side streams are introduced into or re-
occurred when the vapor holdup was neglected. The vapor
moved from the column, and A ( t ) is o€ the form
holdup may be included and more general models may be
constructed but the solution, in any event, follows t h a t presented
in this discussion.
It is assumed that
a12

as2 -aar
0
a23
0
0a4
a :) (39)

where K,(i) is a function which can be determined a priori. I n where

particular, if the solution t o Equation 29 is t o be used to predict
the transient characteristics of the column, following a small
fluctuation of the various parameters of the system, such as the
composition of the feed, the heat content of the feed, the reflux
ratios, etc., then y,(i) and x,(i) can be made to denote the devia-
tions from the corresponding steady-state values in the vapor
and liquid, respectively, and Equation 30 may be written
It is important to note that the form of A ( t ) remains unchanged
if side streams are introduced into or removed from the column.
I n the following discussion the restriction will be to problems
By combining Equat,ions 29 and 30, and by restriction to a where the square matrix, A , is of the form given by Equation 39.
representative component, i, a single matrix equation can be Special Case. I n the special problem A ( t ) E A, if i t is
obtained. assumed that the various flows in the column are independent of
the time, so t h a t A ( t ) is a constant square matrix and if, in addi-
tion, it is supposed that
z ( t ) = z(0) for t <0
A problem similar to the present one, but simpler, has already
been examined b y Lapidus and Amundson ( 1 5 ) ,and an equation and
similar t o Equation 29 was solved by the use of the Laplace z ( t ) = z(0) + S(t) for t > 0
transform. More general results may be obtained by making
use of the Peano-Baker method ( 9 ) ,which is more straightforward Equation 33 can readily b e simplified to
than t h a t of the Laplace transform, and which can be employed
even in cases where the use of the Laplace transform would fail x(t) = x(0) + J' e A ( l - a ) i ( s )ds
to give a solution.
Equation 32 can be solved by iteration to give
Equation 41 is equivalent, b y the Sylvester formula, to

where M ( A ) is the matrizant of A .

Properties of the Matrizant. By definition, or to

where N is the total number of plates, including the still pot.

It was shown b y Frazer and coworkers (9) t h a t M ( A ) is always Finally, if Z = z, a constant, then the solution is
convergent.
-
x ( t ) = ~ ( 0 ) A-lz + A-'eAtz (44)
M h ( A ) = M ; , ( A )34:1(A) (35)
or, by the Sylvester formula,
Actually, the form of the matrizant given by Equation 34 is
useless for computational purposes. However, i t is equivalent t o

Stability of the Solution. It is of some interest t o examine the

solution as given by Equation 45. The following become
and, therefore, if c is fairly large,
obvious:

(37) 1. T h e elements of x ( t ) will increase indefinitely if even one of

the roots of A is positive.

August 1955 INDUSTRIAL AND ENGINEERING CHEMISTRY 1537

ENGINEERING, DESIGN, AND EQUIPMENT

2. T h e elements of x ( t ) will exhibit a damped harmonic

motion if all the roots of A have negative real parts and if a t
least two of them are complex. Table I. Numerical Values for Coefficients bSi for
3. T h e elements of x ( t ) will be monotonic increasing functions
of t, but bounded, for m(i)> 0, if all the roots of A are real and
1 < 5,j 6 12
Values of 1.j
negative. s = 1 s = 2 8 = 3 5=4
1 1
It will be shown t h a t for any matrix, A, of t h e form of Equation 2 - 0.500000 0.500000
39, the elements of which are given by Equation 40, all the roots 3 0.333333 - 0.500000 0.166667
4 - 0.250000 0.468333 -0.250000 0.041667
of the characteristic equation 5 0.200000 -0.416667 0.291667 - 0.083333
6 -0.166667 0.380566 - 0.312500 0.118056
IXZ- AI = 0 7
8
0.142857
-0.125000
- 0.380000 0.322222 -0.145833
0.324100 -0.325694 0.167882
9 0.111111 -0.301984 0.325518 -0.185417
0.199427
are real and negative. T o prove this assertion, however, the 10 - 0,100000 0.282897 -0,323165
following important theorems must be stated: 11 0.090909 -0.271360 0.319504 -0.210676
12 - 0,083333 0.251656 -0,315068 0.219745
1. T h e matrix, S-lAS-ihat is, the matrix obtained from s = 5 s = 6 5=7 5=8
A b y a similarity transformation, where S is any lionsingular ( X 10') ( X 102) ( X 108) ( X 109
matrix-has t h e same characteristic roots as A. 5 0,083333 ...
2. All the roots of a symmetric matrix are real. 6 - 0.20833 0 .' i i s s s g ...
7 0.347222 -0.416667 o.'iQs413
3. A sufficient condition for all the roots of a symmetric 8 - 0.486111 0.79861 1 - 0.694444 0.248016
matrix, H, to be negative, is t h a t quadratic form, H ( z , z ) , be 9 0,618634 - 1.250000 -2.604167
1.504630 -0.992063
negative definite. 10 - 0,742188 1.743634 2.397487
11
12
0.856013 - 2.269838 3.952546 - 4.546958
- 0,960242 2.784862 - 5.606366 7.461833
The first step in the proof would be to reduce A by means of
a similarity transformation to a symmetric form. A is assumed s = 11 s = 12
(x 107) ( X 10s)
t o be of the form given by Equation 39, where the elements, 9 0.275573 ...
, positive real members. The matrices, S I , SB,. . .,S, + I ,
a S j are 10 - 1.240079 0 .'275573 ..,
11 - 3.306877
6.820461 -1.377866 0 .'260521
are defined as follows: 12 4.018784 - 1.377869 0 .'208768

desk calculator, use will be made, for the determination of the

roots of A, of a method proposed b y Hicks ( 1 2 ) . First it is noted,
and so on. however, t h a t if
These matrices, S,, are all diagonal, and since, if Ki is the matrix
derived from Z by replacing the 1 in the i t h diagonal position by A-1 0 , Ao = I, Ai --ail,
IC, then K,A is the result of multiplying the ith row of A by k ,
and, similarly, AK, is the result of multiplying the i t h column then A, = - a n n A n - l - an-l,,&,n-lAn-~ (47)
of A b y k , then it can be shown without difficulty t h a t
an iterative formula which enables t h e calculation of
SGil . . . S;'S;'ASiSz . . S Y + ~= H
A ~ + i IAl
z
where H , a symmetric matrix, is
without difficulty. The method proposed b y Hicks consists of
the reduction of the characteristic equation t o a polynomial in
X by means of Kewton's interpolation formula.
A'+ 1
If A(A) = a,(X - A,)" (48)
(46) 0

Therefore, b y means of a similarity transformation, A has been and if A(A) for X = A0 +

r dA is computed, where 0 < < N + 1,
1'
reduced t o a symmetric form. A sufficient condition for all and dX' is arbitrary, it is found that
the roots of H to be negative is t h a t the quadratic form, H ( z , z ) ,
b e negative definite, which means t h a t ( 8 ) uo = A(Ao)

A, = -dll < 0, Az > 0,A3 < 0,etc.

and (49)
But J = 8

A2 = -CLBZAI- a l u ~ z i=
1
- [LiLz + G1Lz + GIG21 >0 where AlA(A0) = A(Ao + dA) - A@,), etc.
hihz
The coefficients, b,j, can be computed once and for all, as follows:
where G, = V,K,
and, in general,

where

c j ( j ) = 1 and ci(j + 1) = ci-i(j) - j c % ( j ) ; co(j) = 0

Since H ( z , z ) is clearly negative definite, the roots, Ah, of A are

all real and negative, and, as a consequence, t h e elements of the
vector, x ( t ) , given by Equation 45, are monotonic increasing Numerical values for a few coefficients, b,,, are shown in
functions of t, for Z F ( ~ )> 0, and bounded. Table I.
I 1
Evaluation of Roots. Siricr the evaluation of A where A is a Numerical Example. This problem has already been solved
square matrix of the form of Equation 39, can be carried out with by Lapidus and iimundson (16). T t will be treated in this dis-
ease on a sufficiently powerful digital computer, or even on a cussion b y the matrix method. From Equations 29, 30, and 32,

1538 INDUSTRIAL AND ENGINEERING CHEMISTRY Vol. 47, No. 8

 ENGINEERING, DESIGN, AND EQUIPMENT

$ = -Ax + z(l) Example IV. Transient behavior of a system

of N continuously agitated tank reactors in
series in which a chemical reaction i s occurring
with

A=(-;
a + B -a
a-6+ p
0
o rfo
-01 -01
0
0
i ) Let the chemical reaction

take place in a system of continuously stirred tank reactors con-

nected in series. T h e following symbols are used:
~ ( n=) concentration of t h e i t h substance in t h e n t h reactor,
l < i < 3 and 1 < n < N
where 01 and p, assumed for simplicity t o he constant from plate ~ ( 0=
) concentration of t h e i t h substance in the feed entering
t o plate, are t h e first tank
8 = holding time, assumed equal in each reactor

The continuity equation for the first tank reactor results in

If the system is a t steady state for 1 < 0, the matrices, x and z,

can be made to measure changes in t h e composition of t h e various
streams in t h e column and in t h e two feeds, for t > 0. Moreover,
if i t is assumed t h a t z is independent of the time, t h e solution and similarly for t h e other two substances. However, t h e three
t o Equation 51 is given by continuity equations thus obtained can be expressed in the single
relation
x ( t ) = (I - e--46)A-k (52)
(55)
Or, hecause of t h e Sylvester theorem

(
and
where Akr 1 < k 6 N, are the characteristic roots of A . (kl + ;) -k:

As a n example, consider the following: B = -kl (k: + + 8'>

k2

a = 0.634, p = 0.539, S = 6 plates, zo = 0, x1 =0.0881

0 -kz (k: + f>
Because of t h e simplicity of A , a n expression for t h e characteristic
roots, Ah, in a closed form can easily b e found. For illustrative If it is supposed, for simplicity, t h a t a t t = 0 the concentration
purposes, however, the characteristic equation will be obtained of all the substances in all the tanks is equal to zero, and that,
by the method recommended above. moreover, c(0) is independent of t , the solution t o Equation 55 is
given b y
A(A) E 1AI - Ai

and by using Equation 47,

x A(X) A'A(X) A2A(X) A'A(h) A4A(X) A5A(X) ASA(h) The more general problem, for which the con-
0 0.29430 -0.27446 0.23418 -0.15676 -0.01452 1.73248 11.25001 centrations in all t h e tanks are not all equal t o
0.5 0,01984 -0.04028 0.07742 -0.17128 1.71796 12.98249 zero a t t = 0, can easily be reduced t o the one under
1.0 -0.02044 0.03714 -0.09386 1.54668 14.70045
1.0 0,01670 -0,05672 1.45282 16.24713 consideration b y changing t h e dependent varia-
2.0 -0.04002 1.39610 17.69995
2.5 1.35608 19,09605 bles.
3.0 20.45213 For t h e second tank reactor

Because of Equations 48 and 49, the characteristic equation is 1

-d dt = ; c(1) - Bc(2)
A(A) X6 - 7.0380A' + 18.9304A' - 24.2625Xa + so t h a t
14.9928Xe - 3.9374A + 0.2943 = 0 1
42) = - [ I - e-Bt - e - B f ( B t ) ] B - 2 ~ ( 0 )
the roots of which are 02

Similarly, for any tank n,

A I = 0.1196, Xz = 0.4441, 0.9129,
Aa = 1.4333, h j = 1.9021, Xe = 2.2266

Finally, substitution of the above into Equation 53 yields

z g = 0.0397(1 - e - a . 1 1 g 6 f )+ 0.0347(1 - e-0.44416) + which, because of t h e Sylvester theorem, can be expressed as
+
0.0262(1 - e - 0 . 9 1 2 8 t ) 0.0167(1 - e--1.4333i) +
0.0081(1 - e - 1 . 9 0 2 1 t )-I- 0.0021(1 - e--2.22662)
k= 1 0 afk
which is the composition of t h e fat oil leaving the sixth plate. (60)

August 1955 INDUSTRIAL AND ENGINEERING CHEMISTRY 1539

 ENGINEERING, DESIGN, AND EQUIPMENT

where Xk, 1 k 6 63, are the characteristic roots of B . Of course, wise process can be described by the matrix Equation 36, with
Equation 59 holds no matter what t h e order of the matrix, B; the matrix, A ( t ) , given by
. Equation
- 39. There is one difference,
hence any number of consecutive reactions may be considered however, in t h a t the elements of A ( t ) are not numbers any longer
As an example, consider the following: but are themselves square matrices of order m. Thus, the ele-
ments, A l j , of A are
e = 1, kl = 0.20, k: = 0.05, k z = 0.10, k: = 0.05
C i ( 0 ) = 1, C2(0) = ($0) = 0
Then
1.20 -0.05
B = (-0.20 1.15 -0.05
0 -0.10 1.05
O )
T h e characteristic equation for B is
/XI - BI Xa - 3.4000X' + 3.8325X - 1.4325 = 0 \.
which has a s roots
XI = 1.2866, XZ = 1.1133, and X 3 = 1.0000
and, therefore, from Equation 60,
and
1 - (1.2866t)i
n - 1 ___ ](-~ : ~ ~ ~ $ )
e-1.2868t

c
0 i! 0.4027

+ (0.8982)" { 1-
n-1 (1.1133t);} ( ::::::)
e-l.1133<
C 0 T - 1.0186 A ( t ) is a partitioned matrix.

+
{ I + e-t 2 '.
0
0.0769
:}(0.3080)
0.6159
From this point the treatment parallels exactly the one given
in the general solution to Equation 32, and, if A ( t ) is independent
of t , the answer is given by Equations 41 and 42. However, the
This one formula gives t h e concentration in the n t h reactor matrix, A , is now of order m ( N +
I ) , something which adds,
of all three components. Additional reactions would merely naturally, t o the complexity of the problem. I n addition, since
increase the order of the vectors. A now does not possess the simplifying features of Equation 39,
and since the evaluation of det A would, as a consequence, pre-
sent some rather formidable computational difficulties, the
Example V. More complicated
method of Hicks, described in t,he previous example, for reducing
systems-partitioned matrices
the characteristic equation of A t o a polynomial, cannot, in gen-
I n examining the transient behavior of stagewise processes, eral, be recommended.
the simplifying assumption has been made t h a t the characteristics B s a matter of fact, the problem of how t o calculate the
of such phenomena could be described b y m independent equa- latent roots of a matrix is a very important one, and has therefore
tions, Equation 32, where A ( t ) is a square matrix of the form of received some attention. There are two main directions along
Equation 39, and where m denotes the total number of compo- which one can proceed. One can attempt to reduce the charac-
nents in the system. However, if the changes t h a t occur in the teristic equation t o a polynomial or t r y to obtain the roots of a
process, following a change in the operating variables, are not matrix by some iterative procedure.
small, Equation 30 does not usually remain valid, and a better Probably the more promising method for reducing the char-
approximation is obtained if i t is assumed t h a t acteristic equation t o a polynomial is t h a t developed b y Frame
m. ( 5 ) . It has the disadvantage, however, like all the methods of a
similar nature, t h a t the accuracy t o which all the calculations
y ( i , n, 1 ) = K ( i , k , n ) z ( k , n, t )
k=l have t o be performed is usually very extensive. One could also
refer t o the article by Fettis ( 8 ) and t o the very interesting and
It is not difficult t o see, then, t h a t the state of the system under complete paper b y Wayland (68). On the other hand, the itera-
consideration, the distillation column in particular, will not be tive methods ( 2 , 6, 9, 19, 1 4 ) for calculating characteristic roots
described any longer by m, ( N +
1) dimensional vectors, xi, invariably result in rather good estimates for the larger roots
but by a single m(N +
1)dimensional vector, x-that is, a column with progressively poorer estimates for the smaller roots.
matrix with m(N + 1) elements. T h e investigation of such a However, from examination of Equations 41 and 42 i t is ap-
system, however, becomes simpler if use is made of partitioned parent t h a t the Sylvester formula becomes useful for large values
matrices. of the time, t , and that, consequently, one is primarily interested
Again, let i be the component-designating symbol, 1 6 i m, 6 in the smaller roots of A . This particular point is open t o debate,
and n be the plate designating symbol, 0 6 n 6 N and, as before, but it appears that, if a sufficiently powerful computer is available,
consider, for the sake of illustration, the transient features of it is best t o use, in general, Equation 41 in preference t o Equation
rectification columns. Define the vectors 42, and then t o estimate the smallest root, Amin., of A from the
behavior of x ( 1 ) for large values of t , since as t + m ,
x(t) - x( m ) - eXmlnJc
where c is a column vector.

T h e vector, x(n), has m elements which are scalars, while the

vector x has N +I elements which are, themselves, vectors;
Summary
x is said t o be a partitioned vector. This article presents a rather brief survey of certain features
With these conventions, the unsteady state of such a stage- of matrix algebra, which are of importance in applied work, and

1540 INDUSTRIAL AND ENGINEERING CHEMISTRY Vol. 47, No. 8

 ENGINEERINGr DESIGN, AND EQUIPMENT

a number of illustrative examples of interest t o chemical engineers. literature cited

Its purpose is to familiarize a wider circle of chemical engineers (1) Amundson, N. R., Trans. Am. Inst. Chem. Engrs., 42,939 (1946).
with t h e field of matrix algebra, which is becoming more and (2) Arnoldi, N. E., Quart. A p p l . Math., 9, 17 (1951).
more useful in many branches of applied mathematics, and t o (3) Boole, G., “Treatise on Calculus of Finite Differences,” Mac-
demonstrate the power of this new approach by examining some millan, London, 1880.
typical chemical engineering problems. (4) Courant, R., and Hilbert, D., “Methods of Mathematical
Physics,” Interscience, New York, 1953.
It is felt t h a t the use of matrix methods will develop rapidly (5) Dwyer, P. S., “Linear Computations,” Wiley, Mew York, 1951.
in chemical engineering problems since i t seems t o be a n ideal (6) Ferrar, W. L., “Algebra,” Oxford University Press, Oxford, 1941.
tool for application t o multicomponent and multistage systems. (7) Ferrar, W. L., “Finite Matrices,” Oxford Univ. Press,
Oxford, 1951.
( 8 ) Fettis, H. E., Quart. A p p l . Math., 8 , 206 (1950).
Nomenclature (9) Frazer, R. A., Duncan, W. J., and Collar, A. R., “Elementary
Matrices,” Cambridge Univ. Press, Cambridge, England.
A = square matrix 1938.
A-’ = inverse matrix (10) Fry, T. C., Quart. A p p l . Math., 3, 89 (1945).
A’ = transposed matrix (11) Hicks, B. L., J. Chem. Phys., 8 , 569 (1940).
= element in i t h row and j t h column of A (12) Jahn, H. A., Quart. J . Mech. and A p p l . Math., 1 , 131 (1948).
= cofactor of a,j (13) Jordan, Ch., “Calculus of Finite Differences,” Chelsea, New
= adjoint matrix of A York, 1949.
= det A = determinant of A (14) Kincaid, W. M., Quart. A p p l . Math., 5, 320 (1947).
= column matrix or vector (15) Lapidus, L., and Amundson, N. R., IND.ENG.CHEM.,42, 1071
= unit matrix (1950).
= null matrix (16) LlacDuffee, C. C., “Theory of Matrices,” Chelsea, New York,
= concentration of i t h substance in n t h reactor 1950.
= reaction velocity constants (17) Mason, D. R., and Piret, E. L., IND.ENG.CHEM.,43,1210 (1951).
= vapor liquid equilibrium constant of i t h component in (18) Michal, A. D., “Matrix and Tensor Calculus,” Wiley, New
liquid mixture leaving n t h stage York, 1950.
= liquid holdup on n t h plate, moles (19) Milne-Thompson, L. M., “Calculus of Finite Differences,”
= liquid leaving n t h plate, moles/unit time Macmillan, London, 1951.
= reciprocal relative volatility of i t h component (20) Murdoch. P. G.. Chern. Eno. Proar.. 44. 855 (1948).
(21) Norlund, N. E., “Vorlesungen uder Differenzenrechnungen,” J .
= reflux ratio
= time
Springer, Berlin, 1924.
(22) Perlis, S.,“Theory of Matrices,” Sddison-Wesley Press, Boston,
= vapor leaving n t h plate, moles/unit time 1952.
= mole fraction of i t h component in overhead (23) Robinson, C. S., and Gilliland, E. R., “Elements of Fractional
= mole fraction of i t h component in feed Distillation,” McGraw-Hill, New York, 1951.
= mole fraction of i t h component in liquid stream leaving (24) Tiller, F. M., Chem. Eng. Progr., 44, 299 (1948).
feed tray (25) Tiller, F. M., and Tour, P. S., Trans. Am. Inst. Chem. Engrs.,
= mole fraction of i t h component in liquid leaving n t h 40, 319 (1944).
tray (26) Underwood, A. J. V., Chem. Eng. Progr., 44, 603 (1948).
= mole fraction of i t h component in vapor leaving n t h (27) Wallenberg, G., and Guldberg, A., “Theorie der Linearen
tray Differenzengleichungen,” B. G. Teubner, Berlin, 1911.
= a determinant
(28) Wayland, H., Quart. A p p l . Math., 2, 277 (1945).
= a characteristic root of A RECEIVED
for review September 7, 1954. ACCEPTED
M a r c h 1 , 1955

A Test for the Validity

WALTER P. REID
U. S. Naval Ordnance Ted S f a f i o n , China Lake, Calif.

T HE person performing a sedimentation experiment wishes

t o know how large a concentration he may safely use. I n the
literature one can find many statements t o the effect t h a t inter-
matically. T h e results obtained in this way serve to help one
guess what is happening during sedimentation.
Throughout this paper it is assumed t h a t the particles fall
actions between particles may be neglected provided t h a t the independently, and hence that they pass freely through each
concentration of sedimentary material is below about 1 t o 5%. other. T h e expected number of times, h, t h a t a particle OC ran-
T h e analysis t h a t follows indicates t h a t one must also take into dom initial position passes through other particles, and the ex-
account the shape of the particle size distribution curve. pected number of times, H , t h a t other particles pass through it
T h e sedimentation method of determining particle size dis- are calculated. These results are g k e n in Equations 9 and 6.
tributions is based on the assumption t h a t the particles fall inde- Both h and H depend on the particle size distribution, g ( r ) , and
pendently of each other. For this assumption t o be completely each is proportional, not to the concentration, but t o the
valid, i t would be necessary, among other things, for the particle,. amount of sedimentary material per unit of cross-sectional area
t o pass freely through each other. This, of course, they cannot of the (cylindrical) settling vessel. T h e expected fraction, f,
do. Xevertheless, it is possible to assume t h a t the particles can of its settling time t h a t a given size particle spends passing
pass through each other and t h a t is what is done in this paper. through other particles is given in Equation 14. T h e expected
-4ny attempt t o take into account the intricacies of particle inter- fraction, F , of its settling time during which other particles are
actions during sedimentation would lead t o a very difficult prob- passing through it is given in Equation 12. Bothf and F depend
lem. However, by assuming t h a t the particles fall independently on the particle size distribution, and each is directly proportional
one obtains a simple problem which is easily handled mathe- t o the concentration.

August 1955 INDUSTRIAL AND ENGINEERING CHEMISTRY 1541