A P P L I E D M A T H E M A T I C S
H. L. DRYDEN J. M. LESSELLS
E D I T E D B Y
T. C. FRY W. PRAGER J. L. SYNGE
TH. v. KARMAN i. S. SOKOLNIKOFF
H. BATEMAN J. P. DEN HARTOG K. O. FRIEDRICHS F. D. MURNAGHAN R. V. SOUTHWELL
W IT H T H E C O L L A B O R A T IO N O F
M. A. BIOT H. W. EMMONS J. N. GOODIER S. A, SCHELKUNOFF G, I. TAYLOR
L. N. BRILLQUIN W. FELLER
P. LE CORBEILLER W. R. SEARS S. P. TIMOSHENKO
V o l u m e I I
JANUARY • 194S
N u m b e r 4OF
A P P L I E D M A T H E M A T I C S
T his periodical is published quarterly under the sponsorship of Brown Uni
versity, Providence, R ,I. For its support, an operational fund is being set up to which industrial organizations m ay contribute. T o date, contributions of the following industrial companies are gratefully acknowledged:
Be ll Te l e p h o n e La b o r a t o r ie s, In c.; Ne w Yo r k, N . Y ., Th e Br ist o l Co m p a n y; Wa t e r b u r y, Co n n.,
Cu r t iss Wright Co r p o r a t io n; Air p l a n e Di v is io n; Bu f f a l o, N . Y ., Eastm a n Kod a k Co m p a n y; Ro c h este r, N , Y .
Ge n e r a l Electric Co m p a n y; Sc h e n e c t a d y, N . Y.,
Gu l f Rese a r c h a n d De v e l o p m e n t Co m p a n y; Pit t s b u r g h, Pa., L e e d s & N o r t h r u p C om pany; P h ila d e lp h ia . P a .,
Pratt & Wh it n e y, Div is io n Nil e s-Be m e n t-Pond Co m p a n y; West Ha r t
f o r d, Co n n.,
Re p u b l ic Av ia t io n Co r p o r a t io n; Fa r m in g d a l e, Long Is l a n d, N . Y ., Un it e d Aircr a ft Co r p o r a t io n; Ea st Ha r tfo r d, Co n n.,
We s t in g h o u s e Electric a n d Ma n u f a c t u r in g Co m p a n y; Pit t s b u r g h, Pa. T he Qu a r ter ly prints original papers in applied m athem atics which have an intim ate connection with application in industry or practical science. I t is ex
pected th at each paper will be of a high scientific standard; th at the presentation will be of such character th at the paper can be easily read by those to whom it would be of interest; and th at the mathem atical argument, judged by the standard of the field of application, will be of an advanced character.
Manuscripts submitted for publication in the Q u a r t e r l y o f A p p l i e d M a t h e m a t i c s should be sent to the Managing Editor, Professor W, Prager, Quarterly of Applied Mathematics, Brown University, Providence 12, R. I., either directly or through any one of the Editors or Collaborators, In accordance with their general policy, the Editors welcome particularly contributions which wili be of interest both to mathematicians and to engineers. Authors will receive galley proofs only.
Seveuty-five reprints without covers will be furnished free: additional reprints and covers will be supplied at cost.
The subscription price for the Q u a r t e r l y is $6.00 per volume (Aprii-January), single copies
$2.00. Subscriptions and orders for single copies may be addressed to: Quarterly of Applied Mathe
matics, Brown University, Providence 12, R.I., or to 450 Ahnaip St., Menasha, Wisconsin.
Entered as second class matter March 14,1944, at the post office at Providence, Rhode Island, under the act of March 3, 1379. Additional entry at Menasha, Wisconsin.
GEOSCE SANTA PUBLISHING COMPANY, MENASHA, WISCONSIN
Q U A R T E R L Y
OF
A P P L I E D M A T H E M A T I C S
EDITED BY
H. L. DRYDEN T. C. FRY TH. v. KARMAN
J. M. LESSELLS W. PRAGER I. S. SOKOLNIKOFF
J. L. SYNGE
H. BATEMAN J. P. DEN HARTOG K. O. FRIEDRICHS F. D. MURNAGHAN R. V. SOUTHWELL
WITH THE COLLABORATION OF M. A. BIOT
H. W. EMMONS J. N. GOODIER S. A. SCHELKUNOFF G. I. TAYLOR
L. N. BRILLOUIN W. FELLER
P. LE CORBEILLER W. R. SEARS
S. P. TIMOSHENKO
Vo l u m e I I 1944
G e o r g e B a n t a P u b l i s h i n g C o m p a n y
Menasha, Wisconsin
T. Alfrey: Non-hom ogeneous stresses in visco-elastic m e d ia ... 113 A. Barjansky: T he distortion of the Boussinesq field due to a circular hole . . 16 H . J. Barten: On the deflection of a cantilever b e a m ...168 G. K. Batchelor: Power series expansions of the velocity potential in com
pressible f l o w ... 318 A. M . Binnie: Stresses in the diaphragms of d iap h ragm -p u m p s... 37 K. E. B isshopp: Lateral bending of sym m etrically loaded conical discs . . . 205 B. S. Cain: On the treatm ent of discontinuities in beam deflection problems . 353 G. F. Carrier: T he thermal-stress and body-force problems of the infinite
orthotropic s o l i d ... ' 31 W. Z. Chien: T he intrinsic theory of thin shells and plates:
Part II. A pplication to thin p l a t e s ...4-3 Part III. Application to thin s h e l l s ...120 A. P. Cowgill: T he m athem atics of weir f o r m s ... 142 H . B. Curry: T he method of steepest descent for non-linear minimization
p r o b l e m s ... 258 H . W. Em m ons: T he numerical solution of partial differential equations . . 173 L. Fox: Solution by relaxation methods of plane potential problems with
mixed boundary c o n d it io n s ...251 J. N . Goodier: On combined flexure and torsion, and the flexural buckling of a
tw isted b a r ...93 M . G reenspan: Effect of a small hole on the stresses in a uniformly loaded
p l a t e ... 60 M. H erzberger: Studies in optics:
I. General coordinates for optical system s with central or axial sym m etry 196 II. A nalysis of a given system with the help of the characteristic function,
using the direct method of a n a ly s is ... 336 N . J. Hoff: Corrections to my paper “A strain energy derivation of the tor
sional-flexural buckling loads of straight columns of thin-walled open sections” ... 172 B. H offm ann: K ron’s method of s u b s p a c e s ...218 W . H . Jurney: N o te on flow in c a n a l s ... 342 E. G. K eller: Some present nonlinear problems of the electrical and aero
nautical i n d u s t r i e s ... 72 E. K osko: On the treatm ent of discontinuities in beam deflection problems . 271 O. L aporte: Rigorous solutions for the spanwise lift distribution of a certain
class of a i r f o i l s ... 232 K. Levenberg: A method for the solution of certain non-linear problems in
least s q u a r e s ...164 A. N. Lowan: N ote on the problem of heat conduction in a sem i-infinite hollow
c y l i n d e r ... 348 A. N. Lowan and H . E. Salzer: Formulas for complex interpolation . . 272 W. R. M acLean: T he resonator action t h e o r e m ...329
---Effect of a small hole on the stresses in a uniformly loaded plate . . 350 D . M oskovitz: T he numerical solution of Laplace’s and Poisson’s equations . 148 H. Putm an: T he dynam ics of a diffusing g a s ... 267 D. H . Rock (See A . W e in s te in)
H. E. Salzer (See A . N . L o w a n)
S. A. Schelkunoff: Impedance concept in w ave g u i d e s ... 1 ---Proposed sym bols for the modified cosine and exponential integrals . 90 ---On w aves in bent p i p e s ...171 F. Steinhardt: N ote on the elliptic w i n g ...346 H . J. Stew art: T he aerodynam ics of a ring a ir f o il... 136 J. L. S y n g e : A geom etrical interpretation of the relaxation method . . . . 87 T he problem of Saint V enant for a cylinder with free sides . . . . 307 H . W ayland: Expansion of determ inantal equations into polynomial form . 277 A. W einstein and D. H . Rock: On the bending of a clamped plate . . . . 262 Book R e v i e w s ... 9 1 ,2 7 5 ,3 5 4
Q U A R T E R L Y O F A P P L I E D M A T H E M A T I C S
V o l. I I J A N U A R Y , 19 4 5 N o . 4
EXPANSION OF DETERMINANTAL EQUATIONS INTO POLYNOMIAL FORM*
H A R O L D W A Y L A N D C alifornia In stitu te o f Technology**
I. I N T R O D U C T I O N
Determ inantal equations of the form
| A 0X " + A iAn_1 -f- ■ ■ • + A n \ = 0 ,
where the coefficients A { are square matrices of order m , arise in quantum mechanics, electrical circuit theory, the theory of small vibrations, and in m any other branches of physics and engineering. In the process of solving such equations, it is frequently desirable to expand them into polynomial form. Since most of the techniques for such expansion are not discussed in the standard works on numerical com putation, it has seemed advisable to make a critical comparison of the methods available. T he most im portant methods will be described in sufficient detail to aid the non-professional computer.
1. B asic techniques of solution of determinantal equations. Three basic techniques for solution of determinantal equations m ust be considered: (1) Direct solution of the equation by numerical methods [ l] . (2) Direct solution by the method of matrix m ultiplication [2 -6 ] (directly applicable only to the case \ A — /X | = 0 ) . (3) Expansion into polynomial form [6-14] and solution of the polynomial equation by standard m ethods.
In spite of the fact th at much effort has been spent in developing techniques to avoid expansion into polynomial form, th at very technique frequently proves to be m ost economical of effort. Let us consider the numerical solution of the sim plest case:
D{ A) A - JX =
Ö11 — X ÖJ2 ' ' ' O-ln
&21 Û22 — X ' • • 32«
3|i2
( 1 )
If w e assign a value to X and evaluate D(K) using Chib’s expansion [15 ] of th e numeri
cal determinant, we shall make the following number of operations:
* R eceived M a rc h 24, 1944.
** O n leav e fro m th e U n iv e rsity of R ed lan d s, R ed la n d s, C a lifo rn ia.
(1 /3 )(n 3 + 2n — 3) M-D (multiplications and divisions), («/6)(m + l)(2 n + 1) A-S (additions and subtractions).
T he number of additions and subtractions includes the n subtractions required to evaluate the diagonal terms of (1). If w e have D(A) in polynom ial form, we can evalu
ate it by syn th etic division with only n m ultiplications and n additions and subtrac
tions. W e can obtain the polynomial expansion of (1) by D anielew sky’s method (§2c) with n 3 — 2w + l m ultiplications and divisions and n ( n — l ) 2 additions and sub
tractions; hence to obtain k different values of D(K) by first obtaining the polynomial expansion, and then evaluating the polynomial, w e need
(w3 — 2 n + 1) + k n M-D,
n ( n — l ) 2 -f- kn A-S.
To obtain k values of D(K) from the determ inant will require (k / Z ) { n 3 + I n - 3) M-D, ( k n / 6 ) ( n + 1)(2m + 1) A-S.
A comparison of these results shows that it will be quicker to obtain the polynomial expansion first if w e need more than three values of D(K).
If w e use iterative methods to solve (1), w e form the sequence of matrices
A X , A - X , A 3X , ■ • • , (2)
where X is an arbitrary column matrix,*
X = { X lt ■ ■ • , X „ |.
To form each member of the sequence (2) requires n* multiplications and divisions and n ( n — 1) additions and subtractions; hence for k iterations we need
kn* M-D,
k n ( n — 1) A-S.
If we have the polynomial form, each iteration requires only n m ultiplications and divisions and n — 1 additions and subtractions. Adding these to the operations required to expand (1) into polynomial form by th e modified Danielew sky method (§2c) w e need
(« 3 - 2n + 1) + k n M-D,
»(« — l ) 2 + k ( n — 1) A-S.
Hence if £ ^ « + 1 the expansion to polynomial form represents a net saving, still using the powerful iterative method [3, 4, 16, 17, 18].
I I . E X P A N S IO N O F A D E T E R M I N A N T A L E Q U A T IO N IN T O P O L Y N O M IA L F O R M
2. M ethods applicable to the case | A — 7A[ = 0 . Direct expansion is tedious except for the very lowest orders, although som etim es it is desirable because all the elem ents need not be given numerically. Purely numerical m ethods, such as the method of un
determined coefficients or the use of an interpolation formula, will be described further
* I n th is p a p e r, a row of q u a n titie s enclosed in b ra ce s w ill d e n o te a co lu m n m a trix .
on, as th ey are applicable to the most general case. Of the various methods particu
larly applicable to the present case (Eq. 1), five will be discussed. These will be taken up in chronological order of discovery.
a) T h e method o f Leverrier [7, 12, 19]. Until recently, Leverrier’s method was probably the best general method for obtaining the polynomial expansion of the char
acteristic equation of the matrix A . Let the characteristic values of this matrix, i.e., th e roots of the equation | A — IX| = 0 , be Xi, X2, • • ■ , An. Then from well known rela
tions between the coefficients of a polynomial equation and its roots [20], w e have
| A - IX | = ( - 1)’-A" + ( - l) - i ( X i + X2 + • ■ • + XJX"-1
+ ( - 1 ) ’! -2(XxX2 + XxX3 + • • • + X„_xXn)Xn _ 2 + • • •
+ ( X i X 2 ■ • • X „ ) . ( 3 )
B ut direct expansion of the determ inant gives
I A — 7X | = ( — l ) nA” + ( — l ) n_I(an + 022 + • • • + 0nn)Xn_1 + ■ • ■ , whence we have the relation
« 1 1 + 0 2 2 + ■ • • + 0 n n = X l + X 2 + • ' • + X „ . ( 4 ) If we consider the ¿th power of the matrix A , the characteristic values of the matrix
A k will be Ai, Xf, • • - , X* [21]. From this we obtain
| A — In | = ( — 1 ) fi + ( — 1) (X i + X2 + ■ • • + X„)m + • • •
= ( - 1)V + ( — 1)" \ a £ + 0 2 2’ + • ■ • + o in V ” + • • • . which yields the relation
(k) , (k) (*> k k k
Sk = a, n + 022 + • • • + a„ „ = X i + X2 + • • • + X„, (5)
where the terms a[f are taken from the principal diagonal of A k.
If we write our original equation as
| A - 7X | = ( - 1 )» [X » + jf>iXn_I + + • • . + Pn_iX + p„) = 0.
we have, from (3), (4) and (5),
Pl — — ( X l + X 2 + • • ■ + X „ ) = — ( 0 1 1 + 022 + • • • + an n ) = — 5 l , p i — — S\p\ = ( X i + X 2 + • • • + X „ )
= Xi + X2 + • ■ • + Xn + 2(X iX2 + X1X3 + • • • + X„_iX„) = s2 + 2pi,
or
p i — — s 1, P2 = — i ( s i p ! + i 2) , • • • ,
which are of a set of n sim ultaneous linear equations for th e coefficients p k (¿ = 1,
2
Leverrier’s method can be summarized as follows: W e form the powers of thé matrix A , i.e., A k (k — 1, 2, • • • , n ), and add the diagonal terms of each matrix to obtain
(*> , <*0 Sk — «11 + Û22 +
(*) + Ann •
W e then set up the set of n simultaneous linear equations
pk — — ( 1 / k ) (s ip k -i + Sipk-i + • • • + Sk-ipi + Sk) (k — 1, 2, • • • , »), (6) which can be solved for pk (k — 1, 2, ■ • • , n ). W e can then write the equation
D ( \ ) = | A - I \ | - ( - l)"(Xn + />iXn_1 + p i \ n~2 + • • • + p n).
Leverrier’s method is rather tedious because of the labor of forming the powers of the matrix A . Each elem ent of a product matrix will require n m ultiplications and n — 1 additions and subtractions, and each matrix will have w2 terms, except th at one need form only th e n diagonal terms for the last matrix. T he fastest w ay to solve the sim ultaneous equations is to solve them successively, i.e., to use the solution of the first to solve the second, the solutions of the first tw o to solve the third, and so on.
This will require $ ( n 2+ n — 2) m ultiplications and divisions and $ n (n — 1) additions and subtractions. Consequently, Leverrier’s method will require, in the general case,
\ ( n — l ) ( 2 w 3 — 2 m2 + n + 2 )
$11(11 — 1)(2 h2 — 411 + 3) E x a m p le . Let us consider the matrix
A =
M-D,
A-S. (7 )
“ 6 - 3 4 1“
4 2 4 0
4 - 2 3 1
1-----
2 3 1 _
T he sums of the elem ents of the principal diagonals of A , A 2, A3, A 4 are found to be
S i = 12, 52 = 56, 5 3 = 2 8 8, 54 = 1504. Thus Eqs. (6) take the form
P1 = - 1 2 ,
p 2 = ~ ( i / 2 ) ( 1 2 p 1 + 56) = 44,
p 3 = - (l/3)(12/>2 + 56pi + 288) = - 48, p \ — (l/4)(12/>3 + 5 6 ^ + 288/>i + 1504) = 16, and w e have
7X = X4 - 12X3 + 44X2 - 48X + 16.
b) T h e method o f K r y lo v [7 ]. Even as originally formulated by K rylov, this method represents a considerable saving of effort over the method of Leverrier when the order of th e determ inant is greater than four. As modified by Fraser, Duncan and Collar
[22], th e saving is even greater. T he modified form of this method is as follows.
T h e Cayley-H am ilton theorem [23] states th at a square matrix satisfies its own characteristic equation when interpreted as a matrix equation, i.e., if
X" + M " - 1 + M n~ 2 + ' ‘ ' + Pn = 0
is the characteristic equation of the matrix A , then
A n + p i A n~' + P2A" - 2 H + p j = 0. (8 )
If we postm ultiply (8) by an arbitrary column matrix
C ( 0 ) = {cio , C20, • • • , Cno},
and define the sequence
C ( l ) = A C ( 0 ) = ( c u , c 2i, ■ • • , c „ i } ,
C ( 2 ) = A C { 1) = A 2C ( 0 ) = {C12, c M , ■ ■ • , cn2 },
w e obtain the matrix equation
C ( 0 ) p n + C ( l ) p „ - i + • • • + C ( n — 2 ) p i + C ( n — l ) p i = — C ( n ) .
This is equivalent to the set of n sim ultaneous linear equations in the 11 unknowns Pu Pi, • • ' 1Pn,
n - 1
y / Cikpn—k — . ( i = 1, 2, • • • , w). (9)
* -0
Solving this set of equations for the p k ’s, we can readily write down the polynomial equation.
T he formation of each of the column matrices C ( k ) requires n2 multiplications and n ( n — 1) additions and subtractions. T he solution of Eqs. (9) by A itken’s method [28]
requires ( l / 2 ) « 2(w + 3 ) m ultiplications and divisions and ( « /2 ) ( « 2 —1) additions and subtractions. Consequently, w e require
(3 /2 )h2(w + 1 ) M-D,
( 1 0 )
(» /2 )(m - 1)(3«■+ 1) A-S.
K rylov’s original method can be shown to require
(1/3) (« 4 + 4m3 + 2»2 — n — 3) M-D,
( 1 1 )
(1/6 )» (» - 1)(2«2 + 7n - 1 ) A-S.
E x a m p l e . Let us again consider the matrix of the previous example. We have C(0) = [1, 1, 1, 1}, C (l) = [8, 10, 6, 10], C(2) = [52, 76, 40, 8 0 ], C(3) = [324, 520, 256, 560}, C(4) = {1968, 3360, 1584, 3664}.
Thus (9) becomes
Pa + 8 / > 3 + 5 2 ^ + 3 2 4 Pi + 1 9 6 8 = 0 , p i + 1 0 / > 3 + 7 6 p t + 5 2 0 p l + 3 3 6 0 = 0 , . p i + 6 p 3 + 4 0 ^ 2 - f - 25 6 p i + 1 5 8 4 = 0 , p i + 1 0 / > 3 + 8 0 / > 2 + 5 6 0 p i + 3 6 6 4 = 0 , whence
Pi = 16, p 3 = - 48, p 3 = 44, Pk - - 12.
The characteristic equation of A is then
16 - 48X + 44X2 - 12X3 + X4 = 0.
p i ~ x p2 P 3- • • Pn
1 - X 0 • • • 0
II 0 1 - X • • • 0
0 0 0 • • • - X
c) The method o f D a n ie le w s k y [8]. The essence of D anielew sky’s method is the transformation of the expression D ( \ ) = \ A — IX| to the Frobenius standard form
T his gives the polynomial expansion directly.
Danielew sky starts with the (w — l)th elem ent of the nth row (a constant term), reduces it to unity, and then uses this to elim inate the constant terms from the other elem ents of the nth row. T his process introduces extraneous terms in X in the (n — l)th row, which can then be removed by m ultiplying the other rows by appropriate con
stants and adding to the (n —l)th row. A similar procedure is then followed with the (« — 2)th elem ent of th e (n —l)th row, and the reduction is continued until the stand
ard form is reached.
T his process of elimination can readily be carried out by matrix multiplications.
Let us consider a matrix of order 6 which has already had two rows reduced. It is then of the form
C =
~ Cli C12 ¿13 C14 Cu Cl6 ~ C 21 Cn C23 C 24 C26 C26 C31 C 32 c 33 C 34 c 36 C36 Cil C42 C43 C44 C45 c 46
0 0 0 1 0 0
_ 0 0 0 0 1 0 _
and the determ inant has been reduced to the form Z?(X) = | C - A | . W e now postm ultiply the matrix C by th e matrix E ,
E
- 1 0 0 0 0 0
0 1 0 0 0 0
— C4 1/C43 — C4 2/C43 1/C43 — C4 4/C43 — c u / c a — C46/C43
0 0 0 1 0 0
0 0 0 0 1 0
_ 0 0 0 0 0 1
to obtain
__ r / / ! / r —, c 11 C12 Cj3 C\\ Cl5 Cl6
/ / / / / /
C 21 C22 ^23 Cu C 25 C 26
f / / > 1 r
Czi C32 C33 C 34 C 35 £36
0 0 1 0 0 0
0 0 0 1 0 0
0 0 0 0 1 0 _
After this transformation w e have
D ( \ ) E = | C' - E \ | .
If this expression is now premultiplied by E ~ l, we return to a form which is e q u a l\o our original expression, but one step closer to the standard form,
Z>(X) = | E ~ 'C ' - E ~ ' E \ | = | C " - 7X | . Fortunately, £ _1 can be written down directly,
“ i 0 0 0 0 0 "
0 1 0 0 0 0
Ci\ C\ 2 C43 Cu Gi 5 C46
0 0 0 1 0 0
0 0 0 0 1 0
_ 0 0 0 0 0 1 _
T his particular prem ultiplication changes only the third row of C', hence the determi
nant has been transformed to the form
/ . / / 1 / 1
Cu — X C12 G\3 Cu Cl6 Cl6
/ / . 1 / f /
£21 G 22 — ^ Giz C 24 C 26 C 26
// It " \ 11 // a
£31 C32 C33 — A C 34 c 36 c 36
0 0 1 - X 0 0
0 0 0 1 - X 0
0 0 0 0 1 - X
A continuation of this process will eventually yield th e normal form.
If the matrix method is followed strictly, it involves an undue am ount of writing.
W ith only a slight increase in th e number of operations, it can be abridged to give greater ease of calculation with a calculating machine, and to permit checking at every stage of the com putation. A numerical example will best illustrate the method and the check.
E x a m p le . For the matrix of the two previous examples, the schem e of calculations will run as follows:
(1) (2) (3) (4)
(4 ')
(5) (6) (7) (8)
(9)
(9 ')
(10) (11) (12) (13) (14) (1 4 ')
(15) (16) (17) (18) (19)
Hence
"4"
2 3
_1J
’ 0 ‘ - 3 6
12 . - 4 .
48 11.3 3 3 3 -4 4
4 5 .3 3 3 3 J
£ £ '
6 - 3 4 1 8
4 2 4 0 10
4 - 2 3 1 6
4 2 3 1 10
(1 .3 3 3 3 0 .6 6 6 7 1 0 .3 3 3 3 3 .3 3 3 3 )
0 .6 6 6 7 - 5 . 6 6 6 7 1 .3333 - 0 .3 3 3 3 - 4 - 5 .3 3 3 3
- 1 . 3 3 3 3 - 0 . 6 6 6 7 1.3333 - 1 .3 3 3 3 - 2 - 3 .3 3 3 3
0 - 4 1 0 - 3 - 4
0 0 1 0 1
0 - 3 6 12 - 4 - 2 8
(0 1 - 0 . 3 3 3 3 0.1 1 1 1 0 .7 7 7 8 )
0 .6 6 6 7 0 .1 5 7 4 - 0 . 5 5 5 6 0 .2 9 6 3 0 .5 6 4 8 0 .4 0 7 4 - 1 . 3 3 3 3 0 .0 1 8 5 2 1.1111 - 1 . 2 5 9 3 - 1 . 4 6 3 0 - 1 .4 8 1 5
0 1 0 0 1
0 0 1 0 1
48 11.3333 - 4 4 4 5 .3 3 3 3 6 0 .6 6 6 7
a 0 .2 3 6 1 - 0 .9 1 6 7 0 .9 4 4 4 1 .2 6 3 8 )
0 .0 1 3 8 9 0 0 .0 5 5 5 6 - 0 . 3 3 3 3 - 0 . 2 6 3 9 - 0 . 2 7 7 8
1 0 0 0 1
0 1 0 0 1
0 0 1 0 1
12 - 4 4 48 - 1 6 0
D (\) = - 12X3 + 44X2 - 48X + 16.
T he explanations of this schem e are as follows. W e first postm ultiply the matrix whose elem ents are given in lines 1, 2, 3 and 4, by the matrix
’ 1 0
0 1
- 4 / 3 - 2 / 3
0 0
0 0 •
0 0
1/3 - 1 / 3
0 1
This is accomplished by dividing the elem ents of row 4 by the elem ent in the third column, 3, yielding row 4'. T he second-order minors of the unit elem ent in th e third column of row 4' are then formed with rows 1, 2, and 3, the unit elem ent ahvays taking the leading position. T his is easily done by writing row 4' on a card, and form
ing the cross products with the various rows. These minors are entered in rows 5, 6 and 7, under the column corresponding to the elem ents with which the cross product is formed. The elem ents in column 3 are formed by dividing the corresponding ele
m ents of column 3, rows 1, 2 and 3 by the italicized elem ent in row 4. Row 8 is im m ediately written as shown. Thus, the elem ent —5.6667 in row 5, column 2, comes from 1( —3) — (0.6667)(4), and the elem ent 1.333 in row 5 column 3 comes from divid
ing the elem ent 4 in row 1 column 3 by 3.
W e next prem ultiply the matrix with th e elem ents given in rows 5, 6, 7 and 8, by the matrix
" 1 0 0 0 "
0 1 0 0
4 2 3 1
_ 0 0 0 1 _
T his is accomplished by writing the elem ents of row 4 in a column on a card (here shown written to the left of rows 5, 6, 7 and 8), and forming the sum of the products of these numbers with the elem ents of the columns of the rows 5, 6, 7, and 8. This yields row 9, which is the transformed form of row 7 after the matrix multiplication.
Since the rest of the matrix is unchanged, it is unnecessary to rewrite it.
T he whole process is now repeated starting with row 9, dividing each element of th at row by the italicized elem ent ( — 36) to obtain row 9'. T his is written on a card and used to form the cross products with rows 5 and 6, giving the elements in the first, third and fourth columns of rows 10 and 11. T he elem ents in the second column are obtained by dividing the corresponding elem ents of rows 5 and 6 by —36. Rows 12 and 13 can be w ritten down im m ediately as shown. T he process is continued until row 19 is reached. A t this stage it is unnecessary to rewrite the entire matrix, since the desired coefficients appear in only the first row, i.e., row 19. T he polynomial can now be written down as shown.
T he columns labelled X ) and X ^ are used for checking the work. T he elements in the column labelled X] are obtained by summing the elem ents of the first four col
umns of th at row, while those in ^ 2 ' come from only three columns, om itting that column which contains the elem ent used as the pivot for the previous set of cross m ultiplications. T he cross products formed with th eX ) columns should equal the ele
ments of the X ^ column at the next stage of the transformation, e.g., the elem ent
— 3.333 in row 6, column 222' comes either from adding the elem ents ( —1.3333) + ( — 0 .6 6 6 7 )+ ( — 1.3333) of row 6, or from the cross product (1) (10) — (4) (3.3333).
Since these give equal results, the com putation of row' 6 is probably correct. This check is not applicable to the row ju st reduced, so there is no point in forming X ^ for that row, or the rows already in standard form. T he accuracy of row 9, and similar rows, is checked by forming the sum of the products of th e column (4, 2, 3, 1) and the elem ents of column X}- Since this product-sum is equal to the sum of the elem ents of row 9, the accuracy of that row is checked. Compensating errors can occur, so the check is not absolute, but it is a great help in avoiding an accum ulation of errors.
W e m ust next consider the exceptional case in which a zero appears for th e ele
m ent with which w e expect to divide in making th e next reduction, i.e., th e elem ent one place to the left of the diagonal. The following two cases arise: (1) There is at least one elem ent in the row under consideration which does not have a vanishing constant term. (2) All of the constant terms in the row under consideration vanish.
Case (1) can be decomposed into two subcases, according as the non-vanishing elem ent is (a) to the left of the diagonal, (b) to th e right of the diagonal. In subcase (a), we add the elements of the column containing the non-vanishing elem ent to the column in which we wish to introduce a constant term. (This technique can also be used if the elem ent im m ediately to the left of the diagonal has a fairly large tabular error, and another term farther to the left is more certain.) T his will not only intro
duce the desired constant term, but in some row it will introduce an unwanted term in X off the diagonal. T his can be removed, however, by subtracting the appropriate row from the row containing the extraneous X. T he reduction can then go ahead as usual. In subcase (b), the determinant is im m ediately factorable into the product of two determinants, one of which is already in standard form.
T hese subcases can best be illustrated by examples. In subcase (a), Ipt us suppose that after two reductions w e reach the form
4 - X 3 - 2 5 3
1 2 -- X - 1 4 1
2 0 4 - X - 1 6
0 0 1 - X 0
0 0 0 1 - X
If we add column 1 to column 2, w e obtain
4 - X 7 - X - 2 5 3
1 3 - X - 1 4 1
2 2 4 - X - 1 6
0 0 1 - X 0
0 0 0 1 - X
This has an extraneous X in row 1, which we can elim inate by subtracting row 2 from row 1, to obtain
3 - X 4 - 1 1 2
1 3 - X - 1 4 1
2 2 4 - X - 1 6 .
0 0 1 - X 0
0 0 0 1 - X
This is now in a form capable of treatment b y the general method. In subcase (b), we might reach the form
4 - X 3 - 2 5 3
1 2 - X - 1 4 1
0 0 4 - X - 1 6
0 0 1 - X 0
0 0 0 1 - X
which can be factored into the product
The determinant on the right is already in the Frobenius standard form, while that on the left can be expanded im m ediately, although in the general case it would have to be reduced further by the genera! method.
In case (2), the vanishing of all constant terms in a given row indicates that A is a factor of D ( \ ) . If the rth row from the bottom has vanishing constant terms, it means that Ar is a factor of D(K). T he determ inant for the lower degree polynomial which we have y et to determine can readily be constructed from the elem ents above the vanishing row. As an example, let us consider
This is equal to
4 - X 3 - 2 5 3
1 2 - X - 1 4 1
0 0 - X . 0 0
0 0 1 - X 0
0 0 0 1 - X
4 - X 3
1
( - X ) '
In its original form, D anielew sky’s method requires
( „ 2 _ 2)(M - 1) M-D,
n ( n — l ) 2 A-S,
and in the modified form given in detail above, it requires (* - l)(n 2 + n - 1) M-D,
n (n - l ) 2 A-S.
( 1 2 )
(13)
In spite of the extra operations required, the modified form is to be preferred, because it is better adapted to routine com putation -with a calculating machine, and because it can be checked a t each stage of the com putation.
d) Reierstfl’s method. Reiers^l [14] bases his method of obtaining the coefficients of the determinantal equation
D (X ) A I = ( - l ) n( An - M - M ”- Pn)
on the fact that the coefficients p t can be calculated as ( — I)**-1 times the sum of all
¿-rowed principal minors of the matrix A . T he method is powerful for low values of n, but for large n th e labor is considerable.
Reiers^l uses a method for com puting the principal minors of A based on the same pivotal method used by Aitken in various numerical processes dealing with de
terminants [28]. In the process of evaluating a determ inant by Chio’s method [15], simple quotients of various minors are obtained, and the method is easily extended to give all of the principal minors of the matrix.
This form of calculation is used here since it is uniform with most of th e other methods described in this paper, *and requires essentially the same number of opera
tions as the application of Reiers^l’s recursion formulae. In fact, it is merely a schema- tization of his method.
In th e process of evaluating the determ inant \ A | by Chio’s method [15], w e use the identity
A =
011 012 013 - • 01 n / 1 r
022 023 - • 02n
021 0 22 023 ‘ • 02n / / r
032 033 ’ ’ 03n 031 032 033 ' • 03n — 011
/ 1 /
0n2 0n3 ’ 0nn Onl 0n2 0n3 • 0nn
where 0y = 0 , 7 —a,iai,'/au, (i, j — 2, 3, • • • , n ) . If we define our principal minors (m i
nors obtained by striking out the same rows as columns) as B r„ .. .w (r < s < • ■ • < w ) , where r, s, • • • , w are the indices of the rows (and columns) used to form the minor, then we have
B i 0110«*»/
Carrying the reduction one stage farther, w e obtain
A — 011022
033n 03411
n //
043 044
/ / //
0n3 0n4
03n//
04n//
where, as before a'tj — alj Then
73 ' "
2>121 — ^11^22 d ll,
and so on through the reduction. This gives all the second order principal minors of the form Bi„ the third order minors of form B u t , th e fourth order minors of th e form -Bi23u, and so on, including the value of th e determinant itself.
For the minors in which th e first two subscripts are 1 and 3, we start with the determinant
a n
033 • • • a3 n at3t
a n3 •
and carry through th e pivotal reduction as before. For minors of the form B 2. . . we start with
a22 a23 ■ * • a 2, 032 033 • • • 03.
0n2 0n3
In this w ay all of the principal minors can be built up. Since we always start with a term on the principal diagonal, the algebraic sign presents no problem.
T he number of operations required by this method is:
5-2" — (m2 + A n + 5) M-D, 4-2" - (»2 + 3 w + 4) A-S.
(14)
E x a m p le . Let us again consider the expression
6 - X - 3 4 1
4 2 - X 4 0
D ( X) = .
4 - 2 3 - X 1
4 2 3 1 - X
T he schem e of calculations will run as follows:
flu = 6 1 —1/2 2 /3 1 /6 B x = 6
4 2 4 0 b2 = 2
4 - 2 3 1 B z = 3
4 2 3 1 B< = 1
022 ” ^ 1 1/3 •1 /6 B n = 6 -4 = 24
0 1/3 1/3 B n = 6 1 / 3 = 2
4 1/3 1/3 B 14 = 6 -1 /3 = 2
^ = 1/3 1 1 B 123 = 6 - 4 1 / 3 = 8
- 1 1 B124 = 6 -4 -1 = 24
a £ = 2 1 B 1234 = 6 -4 - (1 /3 ) - 2 = 16
an = 6
Ö33 = 1/3 1 1
1/3 1/3
0 B m = 6- (1 /3 ) 0 = 0.
T his com pletes all terms with 1 as the first index. Starting with the third order de
terminant w e obtain by striking out th e first row and column, a 22
bL = 7
b " = 8 /7
1 - 2
2
2
3 3
0 1 1
1/7
1
Bza —
3 1 3 1
1
= 0 .
B n = 2 -7 = 14 B u = 2 • 1 = 2 23234 = 2 - 7 -8 /7 = 16
Hence we have
p 1 = ( - 1)*(5! + B 2 + B 3 + B 4) = 12,
p i — ( — l ) 3(Bi2 + B u + B u + B 23 -f- B u + .8 3 4) = — 44, Pi = ( ~
l)4(Tl23
+^124
+ B\2\ + B 234) = 48,p 4 = ( - l ) 5- 16 = - 16,
DPS) = X4 - 12X3 + 44X2 - 48X + 16.
e) S a m u e ls o n 's method. Samuelson [13 ] has devised one of the fastest m ethods yet developed. H is method requires a few more operations than D anielew sky’s, but the routine involved is extrem ely simple.
T o get the polynomial expansion of
D ( \ ) = \ A - JX| = ( - 1)"(X" - M ”- 1 “ M n_2 />"), we consider the differential equation
D ( d / d t ) Xl(t) = ( - l ) n[*in)(/) - p i x i ',~ 1\ l ) - p 2x lr * \ t )---p nXl(t) = 0, (15) where the superscripts in brackets denote derivatives with respect to t. Equation (15) can be written as a set of n sim ultaneous first order differential equations
»1
x i (/) = Y u bijXj(l) (i = 1, 2 , ■ • • , «),
; - i
where
- p i p2 P3 • • pn- 1P n ~
1 0 0 • • 0 0
M = 0 1 0 • • 0 0
_ 0 0 0 • • 1 0 _
is th e “companion m atrix” to the polynomial in question.
A ctually we need a schem e to go from a system in many variables to a high order system in one variable. Samuelson accomplishes this in the following manner.
L et us consider the system
A x ( t ) = x'{t), (16)
where x{t) is the column matrix
x(t) = { * i(0 . *2(0 . • • • . *n(/)}.
Equation (16) gives us n equations in the 2n variables X\, x 2, ■ ■ ■ , x n, x i , x 2 , ■ ■ • , x i . There are insufficient equations to elim inate all of the variables except those carrying the subscript 1. However we can differentiate (16) (« —1) times with respect to t, obtaining the n1 equations
A £ (n-1)(/) = x M (l), A x < ~ r* > (0 = * < - » ( 0 ,
(17) A x(l) = x ’{t).
W e now have n 2 linear equations in the « 2+ w variables (rci, ar2, • ■ • , x{ , x i, • • • ,
A n) An)
1 y 2 } ■ , *“ ). We can use all but one equation to elim inate the m2 —1 variables not involving the subscript 1, and substitute in the remaining equa
tion to get the desired high order equation in Xi and its derivatives.
L et us consider
A =
a n a n • • a
®21 Û22 ' • 02n
a n2 • Onn
an R
5 M
If w e transfer the variables with subscript 1 to the right of (17), w e can rearrange and rewrite it in the form
W =
- I M 0 • • ■ 0 0 0 - 5 0 • • • 0
0 - I M ■ ■ 0 0 0 0 —S • • • 0
0 0 0 • • • - I M 0 0 0 • • • - 5
0 0 0 • • • 0 R 0 0 0 • • • —an
0 0 0 • • • R 0 0 0 0 ■ • • 0
0 0 R ■ ■■ 0 0 0 1 — an • • • 0
0 R 0 • • • 0 0 1 — an 0 • • • 0
(18)
T he elimination of the m 2 — 1 unwanted variables from the first m 2 — 1 equations and the subsequent substitution in the remaining equation can be performed by pivotal reduction [28], always using elem ents of the matrix on th e left of (18) until a single row remains on the right.
Reduction down to the first row containing R can be made in the general form, yielding
R R M R M2
L R M n~ l
0 0 0 0 • ■• 0 0 1 — a n
0 0 0 0 • • • 0 1 - a n —R S
0 0 0 0 • • • 1 — a ii — R S - R M S
1 - a n - R S - R M S . . . - R M ^ S J In practice this is th e point at which to start the reduction.
To set up the matrix (19) requires
m ( m - l ) 2 M-D,
« (» — 1)(» — 2) A-S.
Pivotal reduction of (19) will require
4n2 - 13m + 12 M-D,
4m2 - 1 3m + 12 A-S,
(19)
making totals for the method of
n 3 + 2«2 - 12« + 12 M-D,
„3 + n t _ n M + 1 2 A-S. ( 2 0 )
Samuelson uses a method due to Crout for the reduction of his equations. Crout's method involves forming exactly the sam e products and sums as are formed in A itken’s method used above, although Crout’s formulation involves som ew hat less writing than the above method, but also requires keeping in mind som ew hat more com plicated formulae. For the average engineer or physicist, ease is fully as im portant as speed.
E x a m p le . W e again consider the matrix
] C\ - 3 4 1 "
4 2 4 0
4 - 2 3 1
- 4 2 3 1 _
«u R
5 M
We find th at R M = [ -12, 3, 5], 2?M2= [ - 20, - 2 4 , 8], R M * = [ 2 4 , - 1 2 8 , - 1 6 ] , R S = [8], R M S — [ — 16], R A P S = [ — 144]. T he matrix to be reduced then becomes
~ - 3 4 1 0 0 0 1 - 6 ~
- 1 2 3 5 0 0 1 - 6 - 8
- 2 0 - 2 4 8 0 1 - 6 - 8 16
_ 24 - 1 2 8 - 1 6 1 - 6 - 8 16 144 _
T he reduction then proceeds as follows:
13
1300 507
E
1 4
3
1
3 0 0 0 1
3 2 1
- 1 2 3 5 0 0
✓
1 - 6 - 8 - 1 7
- 2 0 - 2 4 8 0 1 - 6 - 8 16 - 3 3
24 - 1 2 8 - 1 6 1 - 6 - 8 16 144 27
1 1
0 0 1 10 16 5
13 13 13 13 13
4
7
0 1 - 644
~ 3 56 - 1 3
- 9 6 - 8 1 - 6 - 8 24 96 3
1 0 507 5018 12324 3224 3289
1300 1300 1300 1300 1300
200 1 - 6 200 1272 * 288 519
13 13 13 13 13
1 - 1 2 44 - 4 8 16 1
T he pivotal elem ent for each succeeding reduction is made equal to unity by dividing th at row by the value of th at elem ent. T h e column marked X) usec^ as a check. For any stage of the reduction, the cross products are formed using th eX ] column as if it were any other column, and the values entered as usual. These values should equal the sum of the elem ents in the row in which they appear. T he check is not absolute, but it is very useful.
3. M ethods applicable to the case \ A — 2?X| = 0 . -a) T h e method o f reciprocation.
T he equation
\ A - B \ I = 0 has the sam e roots as the equation [24]
I B ~ 'A - I \ I = 0, ( 2 1 )
where B _1 is the reciprocal of B , provided only th at B is not singular. T he matrix product B ~ lA can readily be formed by A itk en ’s method [28], and the determ inant (21), which contains X’s only along the principal diagonal, can then be expanded by one of the methods of §2.
T h e formation of the product B ~ yA requires ( « 7 2 ) ( 3 » _ i)
(« /2 )(« — 1)(3k 1)
M-D, A-S.
Using the modified D anielew sky method to obtain the polynomial form, we shall need all told
(1/2) (» + 1)(5«2 - 6» - f 2) (m/2 ) (m - 1)(5« - 3) E x a m p le . L et us consider the determ inant
D(X) =
M-D, A-S.
( 2 2 )
We have
A =
B ~ 'A =
- 9 + 2 X - 8 + 3 X — 7+X - 7 + 2 X 15-3 X 1 6-5X 1 3 -2 X 1 5 - 4X
— 8+X - 8 + 2 X - 7 + 2 X - 8 + 3 X 23-3 X 24-5 X 1 9 -3 X 2 2 - 6X
- 9 - 8 - 7 - 7 2 3 1 2
15 16 13 15 3 - 5 - 2 - 4
B = )
- 8 - 8 - 7 - 8 1 2 2 3
_ 23 24 19 22 _ 3 - 5 - 3 — 6 _
8 0 - 2 3 " 8—X 0 - 2 3
5 0 1 - 2 5 - X 1 - 2
4 8 4 7 /(X) =
i 8 4 - X 7
6 - 8 - 5 —7 _ —6 - 8 — 5 — 7 —X
B ( B ~ 1A ) =
~ 2 3 1 2 “ - - 8 0 - 2 3“ ~ - 9 - 8 - 7 - 7 “
- 3 - 5 - 2 - 4 5 0 1 — 2 15 16 13 15
1 2 2 3 4 8 4 7 - 8 - 8 - 7 - 8
3 - 5 - 3 — 6 _ _ —6 - 8 - 5 — 7_ _ 23 24 19 2 2 _
= A .
Thus
/(X) = X4 + 11X3 + 33X2 + 8X + 8.
T he calculation of B ~ l can be carried through conveniently :by means of A itken’s method [28] of obtaining the reciprocal of a matrix. T he product B { B ~ lA ) can be formed as a check.
b) T h e D a n ie le w s k y - M a s n y a m a method. An extension of D anielew sky’s method of transforming a X determ inant to the Frobenius standard form has been made by M asuyam a [6] for the case \ A — B \ \ = 0 . This method requires
(l/2 4 )« (« - 1)(7«2 + 13« + 66) M-D, (1 /2 4 )(« - 1)(7«3 + 5«2 + 58« - 48) A-S. (23)
For low orders (through the fifth) this represents a small saving in the number of operations over th e method given in §3a, but it is not as well adapted to machine com putation. It is applicable, however, to the case in which the determ inant of the matrix B vanishes, but in this case the method of §4b, using N ew ton ’s interpolation formula, is to be preferred. For these reasons, we shall not consider the method in more detail here.
4. M ethods applicable to the case | J4 0Xn+^4iXn_1+ • • • + j4 „ | = 0 . -a) T r a n s fo r m a tio n to the f o r m \ A — 7X| = 0 . I t is possible to transform an ?«th order determinant, the terms of which are polynom ials at most of degree n in X, into a determ inant of order m n with terms linear in X, provided that the matrix of the coefficients of X" is not singular. This can be done in more than one w ay, but the following seem s most convenient [25].
If the determ inant we wish to transform is
Z)(X) = | AoXn + ^iX”-1 + • • • + A„_iX + A„ | = 0 , (24) we consider the related set of sim ultaneous linear differential equations
(A oD» + A i D " - 1 + • • • + A„_iD + A n) x = 0 (D = d/d t), (25) where * is the column matrix
{ #1» #2, * , } ,
U o j ^ O . (26)-
and
Because of th e condition (26) th e reciprocal matrix Ao“ 1 exists. Consequently we can prem ultiply Eq. (25) by A r l, obtaining
{ I mD - + Ao-MiZ?"-1 + ••'• + Ao-Vl^-D + A r ' A J x = 0, (27) where I m is a square matrix of order m with units on the principal diagonal and zeroes everywhere else. W e can write (27) in the form
I mD x (n~ v + ¿o-b li*'" -1) + A f hl2* ("-2) + • • • + Ao-1An- 1* fl) + A f ' A nx(0) 0, (28)
