Machine Solutions of Linear Differential Equations A thesis presented by Harvard University Preface The present use of the Mark IV Calculator in the solution of a new dynamic economic model continues the previous work of Dr. H.F. Mitchell in this field. The numerical results obtained should be of interest to economists, and it is believed that a number of the methods employed will prove useful to computers. The writer wishes to express his appreciation to Professor Howard H. Aiken for his helpful guidance and encouragement, and to Professor Wassily W. Leontief who proposed the problem and offered many valuable suggests for its solution. The advice and assistance of numerous members of the Staff of the Computation Laboratory is gratefully acknowledged, in particular that of Dr. Robert C. Minnick and Mr. Anthony G. Oettinger. To Mr. James R. Terrell the writer is indebted for profitable discussions arising out of his use of the matrix subroutines described herein. The writer also wishes to express his thanks
to Miss Jacquelin Sanborn who typed the plates for printing
and to Miss Carmela Ciampa, Miss Joan Boyle, Mr. Robert Burns
and Mr. Paul Donaldson,
all of whom assisted in the preparation of the manuscript.
Synopsis The present study was undertaken to obtain solutions to a dynamic economic model proposed by Leontief, and to develop methods of computation suitable to the Mark IV Computer. The economic model as presented in Chapter 1 is described by a system of linear differential equations with constant coefficients. The data for six systems of order 5, 6, 10, 11, 20 and 21 were supplied by the Bureau of Labor Statistics and the Harvard Economic Research Project. The general solution for each system was sought. In Chapter 3 the complementary function is shown to be expressible in terms of the latent roots and principal vectors of a certain matrix related to the coefficient matrices of the given system. An expression is also given for the particular integrals appropriate to a set of demand functions of interest in economic theory. The socalled Frame method, described in Chapter 4, was chosen for obtaining the latent roots and principal vectors. The Frame process generates the characteristic equation of the matrix and the latent roots of the matrix are then obtained as the roots of the characteristic equation. Methods of solution of polynomial equations of higher degree are therefore discussed in Chapter 5. In Chapter 5 a modification of the quadratic factor method is described which has the advantage of greater uniformity for machine computation. Extensions to a higher order process and to the extraction of quartic factors are considered and a new method for the solution of the quartic is described. A machine method based on a combination of the Graeffe process and the modified quadratic factor method is also proposed. The method is very general in its application and offers several advantages, the most important among them being:
The above method has been used for the solution of several equations of degrees up to twenty. Machine programs for the matrix operations required were organized around a set of matrix subroutines described in Chapter 6. These subroutines are programmed for systems of a general order and are believed to be sufficiently flexible and easy to apply to be of value in general use. They have already been profitably used by another member of the Staff of the Computation Laboratory in an unrelated program. In addition to the standard matrix operations programmed, a new machine method of inversion is described. For machine computation, it has the advantages of uniformity, low storage requirements, and low roundoff error. Several possible extensions to complex matrices are considered. A number of special machine techniques for vector algebra are suggested; among them a process called “floating vector operation” which offers most of the advantages of floating decimal point operations at a fraction of the computing time required by the latter. This process was applied in the calculation of the principal vectors. The requisite details of Mark IV programming are considered in Chapter 2. Particular solutions were obtained for the six systems described.
The results appear in the appendices and are described in Chapter 7.
The complete complementary function was obtained for the systems
of order 5, 6, 10 and 11.
Because of roundoff error accumulation only a partial solution
was obtained for the system of order 20
and a complete solution would necessitate rerunning with greater accuracy.
In view of the solutions of the similar lower order systems already obtained,
the computation required is not warranted by the economic significance
of the missing terms.
Chapter 1 A Dynamic Economic Model Mathematical models have long been used with striking success in the physical sciences and hence attempts have been made to apply them to such diverse fields as economics, biology, and psychology. Early attempts at the application of mathematical models to the field of economics were hampered by (1) the lack of empirical data and (2) the lack of mechanical devices to carry out the extensive calculations required. As a result, economic models were either so oversimplified that they no longer adequately represented real phenomena, or else so complex as to render impossible the calculation of numerical results based upon them. This situation has been radically changed by two factors: (1) the vast amount of statistical data now collected and made available by government and other agencies; and (2) the advent of the large scale automatic calculator which makes feasible the construction of models of much greater complexity than heretofore. The present study was undertaken to obtain general solutions of a dynamic model proposed by Leontief. This is an extension of his earlier static inputoutput model. Both models are fully discussed in references 1 and 2 and only a brief outline will be attempted here. The variables of the static inputoutput model are the outputs of the various “sectors” or “industries” which make up the economy. It is assumed that each sector of the economy will require inputs from some or all of the other sectors, and that a state of equilibrium exists. That is, if x_{1},x_{2},...,x_{n} are the outputs of the several sectors, X_{ij} is the input from the ith sector to the jth sector, and z_{i} is the “final demand” or amount of goods to be made available for use outside the system considered, then
The inputs to each sector are assumed^{[a]} to be linear functions of the output of that sector, i.e. X_{ij} = a_{ij}x_{j}. The “flow” coefficients a_{ij} are assumed to be constants so that the equilibrium condition may now be stated as a set of n linear equations in n unknowns,
For a given final demand z_{i}, the required outputs x_{i} of each sector may now be obtained as a solution of (1.1). The division of the economy into sectors is somewhat arbitrary. In general, the greater the number of sectors considered, the more detailed is the information available in the solution. Since, however, the amount of calculation required varies approximately as the cube of the number of sectors considered, a compromise must be made. Solutions^{3 } have been obtained for systems of various orders, the largest being a system of 192 sectors. Equations (1.1) describe a static model since no dependence on time is implied. Dependence on time may be introduced by making either the coefficient a_{ij} or the demand function z_{i} (or both) functions of time. The simple linear model may also be extended in other ways; e.g. by including terms of higher degree in the variables or by adding terms in the derivatives dx_{i}/dt. For reasons discussed in reference 2 the following dynamic model is chosen^{[b]:}
The flow coefficients a_{ij} and the capital coefficients b_{ij} are assumed to be constants; the outputs x_{i} and the final demands z_{i} are explicit functions of time. The system of firstorder linear differential equations represented by (1.2) may be solved in either of two ways:
The general solution is more difficult to obtain; however, it is of much greater value in studying the behavior of the system than is the particular solution. In the first place, the general solution is expressible as a sum of exponential functions. As in the analogous case of a coupled mechanical system [also described by (1.2)] a knowledge of the socalled “normal modes” or terms of the general solution provides greater insight into the behavior of the system than any tabulated solution, however extensive. Secondly, the particular solution corresponding to any given initial conditions may be obtained fairly readily from the general solution by evaluating the arbitrary parameters. A number of difficulties may arise in connection with the general solution. They may be classified either as mathematical difficulties or as difficulties of interpretation of the model. The mathematical difficulties are problems arising in the actual solution of the system of equations represented by (1.2). They are given consideration in the following chapters. The difficulties of interpretation arise if in certain cases the model does not adequately represent the economic phenomena it purports to describe. These matters will not be discussed here, and again the reader is referred to reference 2. Presumably a final demand can be considered as the input to an additional sector of the economy and by including this sector in the system the equation (1.2) become homogeneous. Such a system is referred to as “closed”. In order to separate the less stable^{[d]} elements of the economy from the more stable, it is convenient to remove the sector classified as “households”from the system and treat it as a final demand on the remaining system. Such a system will be termed an “open” system. The demand function may in general be any function of time, but the functions of most interest are of the form of exponentials, i.e. z_{i} = g_{i}e^{μt} where the g_{i} are constants. Solutions for demand functions of this kind are discussed in section 3F. Flow coefficients for a system of 192 sectors were provided by the Bureau of Labor Statistics and formed the basis for the present calculations. Capital coefficients were supplied by the Harvard Economic Research Project. By a processor of “aggregation” described in Chapter IX of reference 2, the coefficients for several systems of different orders were computed from the original data. The orders of the systems chosen were 6, 11 and 21. Comparison of the solutions for the various systems may be expected to give some indication of the distortion introduced by aggregation. Three “open” systems of order 5, 10 and 20 were derived from the systems of order 6, 11 and 21 by suppressing the “households” section. General solutions for the six systems were obtained using the Harvard Mark IV Calculator. Particular integrals (see Section 3F) were also computed for exponential functions e^{μt} for
At the time this study was undertaken, the Mark IV Calculator had just been completed. The need for solving the extensive economic systems mentioned above therefore provided an excellent opportunity for an investigation of general matrix methods especially suitable for application to Mark IV. Notes Chapter 2 Problem Preparation for an Automatic Computer Section A of this chapter will be devoted to a consideration of the general aspects of problem preparation for an automatic computer as a preliminary to the discussion in Chapters 4 and 5 where possible methods for the solution of our economic model are analyzed and compared on the basis of their suitability for machine computation. Section B will be devoted to programming details of the Harvard Mark IV Calculator as a preliminary to the discussion of Chapters 6 where the programming of the selected methods for the Mark IV Calculator is discussed. Reading of Section B may well be deferred until after Chapter 5. A. General Aspects of Problem Preparation^{[a]} (1) Introduction The solution of virtually any problem of applied mathematics may be reduced to a sequence of elementary operations involving only addition, subtraction, multiplication, division, and reference to previously computed quantities. For this reason it is clear that a truly general purpose calculator may be based on the following components:
The calculator may be made automatic by the addition of a “sequence unit” which determines the sequence of operations to be carried out. The sequence unit operates by reading sequentially a set of instructions presented to it, and carrying out the prescribed operations by determining which numbers are to be transferred, read, or recorded at each step. The sequence instructions or “program” must be coded and recorded in a form suitable to the particular machine. Since each step must be programmed it is clear that little or no labor is saved by the automatic computer unless the operations are made repetitive. Such repetition is made possible by a special instruction referred to as a “call” which interrupts the sequence of instructions being read and cause the sequence unit to begin reading instructions from some point in the program specified by the call. This provides the possibility of a program consisting of a single closed loop which could, for example, be used in the tabulation of a given function as follows: The initial value of the argument is stored in a certain “argument” register and the prescribed function of the number in this register is computed and recorded. The number in the argument register is then incremented and a call is made to the point in the program immediately following the introduction of the initial argument. The computer is made much more versatile
by providing another instruction termed
“conditional call” which effects a call
if and only if the number in a special
“conditional” register is negative.
This allows a choice of programs to be made
depending on results obtained
in the course of the calculation.
For example, a series may be automatically truncated
at the first nonsignificant term as follows:
Suppose that a certain section of the program
is terminated with a call so as to form a closed loop
which repeats indefinitely,
and on each iteration adds a term of the series.
Let the call be replaced by a conditional call
and immediately preceding it let the number
The conditional call is used frequently in programming and, in effect, provides the possibility of making any logical decision necessary in practice. The pattern of the calls and conditional calls occurring in the program divides it into groups of operations called “routines”. (2) Choice of Numerical Method In general, several alternative numerical methods will be available for any given problem; and the first step in the preparation of the problem will be the choice of the most suitable method. This choice will require consideration of the routines necessary to carry out each possible method. The factors which enter into this consideration will now be discussed. The possible error present in computed results must be known if the results are to be of value. Loss of accuracy will be due to two types of errors; those due to certain approximations, and those due to “roundoff”. The first type of error occurs when a certain function is approximated by some more easily computed function as, for example, in the truncation of an infinite series. When such an approximation is used, an estimate of the “truncation error” so committed must be available. Because a computer is limited in the number of digits carried, each number is correct only to a certain number of decimal places. Errors due to this effect are called roundoff errors. The maximum roundoff error introduced at each step is known by the way in which errors from any source propagate depends on the subsequent operations. The total error accumulated is therefore difficult to estimate and depends on the particular method employed. Upper bounds to the possible error can usually be obtained but these bounds are often grossly pessimistic and probable error bounds may have to be considered. Errors of a completely different nature are those introduced by a failure of some part of the computer or by a flaw in the program. Such errors are usually called blunders. They must either be prevented or corrected before the computation is allowed to proceed. For this purpose checks must be included in the program and provision made to rerun the program from some previous point when a blunder occurs. Checks are usually arranged to stop the computer if two quantities which are supposed to agree (because of some mathematical identity) to within a certain tolerance do not agree. Rather than make the machine stop due to an intermittent failure, a rerun may be programmed to occur automatically for a limited number of times before the machine stops. Although a clearout distinction is not always possible it is useful to classify checks as current or final. The purpose of a final check is to insure that the final results recorded are correct. Current checks are provided at frequent intervals, usually at the end of each routine, with the double purpose of reducing the amount of time lost in further computation after an error occurs and of aiding in the detection of a source of error by localizing the trouble to a certain section of the program. Checks should be designed to detect not only machine errors but also any errors in programming. Program errors are by far the most difficult to detect and considerable ingenuity may be required to devise suitable checks. Of two possible methods of solution, the one more general in application is also to be the more complex and the more timeconsuming. Other things being equal, the use of a more general method is to be preferred because of the possibility of using the same program in later problems where the greater generality may be utilized. In some cases, the generality required cannot be completely determined in advance. For example, in Chapter 4 the socalled power method is seen to fail when certain types of multiplicities occur in the roots of the matrix. The amount of computation required is clearly an important factor although it becomes relatively less important as speed and availability of automatic computers [increase]. The time required for a given process may be evaluated roughly by counting the number of multiplications required. This gives a useful measure for comparison. The amount of number storage required by a given method must also be given consideration. The complexity of the program determines the amount of preparation time required of the mathematician, and to a certain extent, the difficulty of running the problem on the machine. When a program is first placed in operation it must be carefully “checked out” to remove any possible programming errors and the more complex the program, the more time this process is likely to take. The probable checkout time should therefore be taken into account in estimating the overall machine time required. The complexity of the program also determines the amount of storage required for sequencing instructions; an important factor in the use of machines with limited storage capacity. The complexity of the program required for a given process depends to a surprising degree on the lack of uniformity of the process. That is to say, an otherwise simple program may become very complex if exceptions must be made for certain special cases. The desirability of uniformity in machine methods is stressed therefore in the discussion of programs in Chapter 6. In summary the major considerations entering into the assessment of a given method are the following:
(3) High Accuracy Operation The digits of a number are stored in successive “columns” of a storage register. In most automatic calculators the decimal point is fixed between a given pair of columns, and the numbers introduced must be suitably scaled so that none of the numbers encountered in the course of the calculation will exceed the capacity of the storage registers. In some calculators such as the Mark IV, the decimal point is fixed during the running of a problem but its location may be chosen at will by the operator. In this case the location of the decimal point must be chosen to prevent the occurrence of overcapacity numbers. In general, the decimal point should be located as far to the left as is consistent with this requirement since the roundoff error is thereby reduced. If very high accuracy is required or if the accumulation roundoff error is very severe, the storage capacity of a single register may not be adequate. In this case each number may be stored in two registers; the high order part of the number in one register, and the low order part in the second. This procedure is know as “double accuracy” or “high accuracy” operation. Provisions are made for the addition and multiplication of such “high accuracy” numbers by suitable programming. It may happen that the number of significant digits required in each number does not exceed the capacity of the registers but that the magnitude of different numbers varies so widely that they cannot be stored at the same decimal point. In this case it is usual to resort to the socalled “floating point” operation in which each number z is represented by two numbers y and p in the semilogarithmic form z = y × 10^{p}. The number y is “normalized”, i.e. it is stored with the first nonzero digit in the highest order column of the register. Some calculators^{[b] } operate directly with numbers of this type and reserve certain columns of each register for representation of the exponent p. In the use of a fixed point machine the corresponding operations must be programmed. (4) Simplification of Programming The speed of modern computers has reached a point where, for many applications, the time required to program a problem may be a more important factor than the actual computing time. Attempts are made, therefore, to simplify the work of programming in various ways. The first step may be made by the programmer in seeking out repetitive operations within the program. Such a repetitive operation need be programmed but once together with instructions to repeat the operation the desired number of times. Similarly an operation which recurs at several different points in a program may be programmed once as a “subroutine” to be called whenever required. Such a subroutine may be arranged to call back automatically to the point in the main program from which it has itself been called. It frequently happens that an operation is repetitive except for a simple permutation of the registers involved. For example, the summation of the numbers in the successive registers 150 is performed by adding successively the registers 1,2,…50 to the previous sum. The only difference in each of the fifty operations is the address of the particular register to be added. Modern machines allow the programmer to make such operations repetitive by making provision for modifying instructions. In the previous example it is only necessary to add one to that part of the instruction specifying the address of the new register to be added each time. In most machines the sequence instructions and numbers are stored in the same set of registers and instructions may therefore be modified in the arithmetic unit. In the Harvard machines the number storage and sequence storage are completely separate and the direct modification of instructions is not possible. Equivalent operation is obtained by the use of certain control registers as follows: Normally each instruction specifies directly the address of the storage register concerned. In Mark IV, a special order is provided which specifies storage address “i_{j}”, where i_{j} is the number stored in a certain control register called the IJ register.^{[c]} Thus, without changing the actual order stored in the sequence storage, the operation effected by it may be modified by changing the number resident in the IJ register. The separation of number storage and sequence storage offers important advantages in making reruns. To make a rerun from any given point in a program all the relevant storage must be restored to the conditions which obtained at that point. For the number storage this may usually be done by so planning the program that the numbers entered into the calculation are still available in the same registers at the time the check is made. Although the number of orders modified in a certain routine may be rather large and the pattern of the modifications rather complex (due perhaps to repetitive loops within the routine), yet the modifications may usually be expressed as simple functions of a very limited number of parameters or “control numbers”. In Mark IV the desired functions of these control numbers are computed and introduced into the IJ register to produce the required modifications. Since the sequence instructions are themselves unchanged, the program may be restored to any point by restoring the new control numbers to their correct values. The rerun procedure is therefore relatively simple and may easily be made to occur automatically. On the other hand, the restoration of the sequence instructions in a machine which modifies them directly is in general so difficult that the simplest procedure is to read the original instructions again from the input device. Automatic reruns become virtually impossible. It should also be noted that easy reruns greatly simplify the process of checking out a program. A disadvantage of the separation of number and sequence storage shows up in the case of problems requiring a large amount of number storage and little sequence storage (or viceversa). In this case the total storage capacity of the machine cannot be utilized. If, however, the storage device used is relatively cheap, e.g. a magnetic drum, the slight increase in total storage capacity required by the separation of numbers and instructions is amply justified by the advantages described above. The algebraic coding machine of Mark IV^{5 } is one approach to the further simplification of programming. Sequencing instructions in essentially algebraic form may be introduced on the keyboard. These instructions are automatically translated into the more complicated program instructions required to sequence the computer and are recorded on tape for later use by the computer. Only the most elementary of the algebraic operations are available by high accuracy and complex operations may also be coded directly. In addition, the elementary functions √x, tan^{1}x, cos x, log_{10}x, and 10^{x}, may be coded directly; e.g. the operation cos x_{0} ⇒ c_{0} computes the cosine of the number in register x_{0} and stores it in register c_{0}. In the simpler operations described above (namely, those requiring less than twelve instructions for their completion) the coding machine records the required sequence instructions directly. In the more complicated operations, e.g. cos x, the coding machine provides only the instructions which are necessary to call a “permanent” subroutine^{6 } which is normally recorded on the sequence drum of the computer at all times. The use of the algebraic coding machine can materially reduce problem preparation time. It also enables a mathematician to use the computer without learning its arbitrary sequence codes. Two minor disadvantages are (1) the coding machine, generally speaking, does not program as efficienctly as is otherwise possible; and (2) in checking out a program and locating errors effectively, the detailed coding must be understood. Even so a mathematician may, with the help of a competent operator, check out a program satisfactorily without knowing the detailed coding. The use of permanent subroutines for elementary functions might well be extended to the use of an extensive library of subroutines. Conceivably such a library could be made to include all the operations likely to be encountered in the majority of problems and programming would then be reduced to the preparation of a relatively simple main program to piece together the required subroutines.^{7} The value of a particular subroutine depends on both the ease and the frequency of its use. On this basis the subroutines for the elementary functions are in general the most valuable and subroutines for linear operations probably take second place. A set of subroutines for matrix operations on Mark IV are discussed in Chapter 6. B. Mark IV Programming^{[d]} This section will be devoted to those details of Mark IV programming^{[e]} which are necessary to an understanding of the routines discussed in Chapter 6. Programs for Mark IV may be prepared either on the algebraic coding machine mentioned in section A or on a “decimal” coding machine, with which the elementary sequence codes are recorded. Since the use of the decimal coding machine leads to a slightly more efficient program, it was chosen for the preparation of the subroutines of Chapter 6. For this reason only the decimal coding will be considered in this section. Orders The units of Mark IV are divided into two groups called β and γ. Under control of a single order a number from any β unit may be transferred to any γ unit or viceversa. Direct transfers between two units of the same group are not possible. An order is made up of three parts called the α, β, and γ codes, which determine the operation to be performed as follows:
A complete list of the codes for Mark IV appear at the end of this section. The calculator has a storage capacity of 10,000 orders, each of which is identified by a “linenumber” in the range 00009999. The function of certain β codes differs depending on whether a readin or readout of the register is involved (i.e., depending on the αcode with which it is combined). Hence it is sometimes necessary to specify the function of a βcode both for readin and for readout. In the coding list this is done by writing “in” or “out” after the appropriate function. Similar remarks apply to the γcodes. Example. The coding list gives for γ code 21: Multiplier in and multiply; low order product out. From this one may deduce the meanings of the following orders:
If neither “in” nor “out” is specified, the function of the code is the same for either, except for direction of transfer. Registers Most of the registers of Mark IV have a capacity of sixteen decimal digits and sign. Such registers are referred to as “regular” in contrast with other “short” registers which do not have a full sixteen columns. All short registers are so arranged that their first (low order) column corresponds to the first column of the regular registers. Readout from any short register to a regular register places zeros in the unoccupied columns. The γ units consist of the γ transfer register and the arithmetic unit (multiplier, divider, and two accumulators). The β units include two hundred fast storage registers and other special registers. The content of any register is not destroyed on readout. Fast Storage Registers The two hundred fast storage registers are referred to by the numbers of the βcodes controlling them and are numbered from 0000 to 0199. It is often convenient to designate them by a letter and subscript, the letter determining the second and third digits as follows:
For example, b_{3} ≡ 0013 and d'_{8} ≡ 0138. Function Registers The function registers are ten in number (β02200229) and are designated as f_{09}. In operation they are identical with the fast storage registers but are normally reserved for use by the permanent subroutines.^{6 } Slow Storage Provision is made for the storage of four thousand numbers in “slow storage” positions numbered 00003999. These numbers are not immediately available but may be transferred to or from fast storage in blocks of ten numbers at a time. Such transfers are made to or from two groups of ten special β registers called S and S' and numbered from S_{0} ≡ 0200 to S'_{9} ≡ 0219. The two groups S and S' are identical in function, but the provision of two groups makes it possible to use one for computation while the other is engaged in a transfer. Transfers may be made in either direction between any block of ten consecutive slow storage positions and either the S or S' group of registers. The slow storage control register (SSCR) is a fourcolumn β register capable of storing a number in the range 00003999. The number storage in the SSCR determines the slow storage position of the first number of the block of ten to be transferred. For example, if the number 0328 stands in the γ transfer register (γ00) the effects of the following orders are:
β codes 1221 and 1231 produce similar results but with the S' registers replacing the S registers. The slow storage control register operates modulo 4,000. Thus the four thousand slow storage positions form a closed ring and may be run through in sequence as many times as desired by simply augmenting the number in the SSCR. The SSCR has ordinary readin and readout (without initiating a slowfast transfer) under β code 3200. The S and S' registers have regular readin and readout and may be used as ordinary fast storage registers. IJ Register The IJ register is a threecolumn β register which can store numbers from 0 to 399. The first (low order) column is referred to separately as the J register, since readin or readout may be made either to the whole IJ register, (β3003) or to the J register alone (β3002). The number standing in the IJ register is referred to as i_{j}. The IJ register controls or modifies the effects of certain β codes as follows: (for this example assume i_{j} = 124)
Shift and Normalize The shift register may be used to shift a number either to the right or to the left by any number of places determined by the number in the shift control register (SCR). The SCR will store any number in the range ±19, a positive number producing a shift to the left, and a negative number producing a shift to the right. The SCR has ordinary readout and readin under β code 0302. The shift register has two types of readout and readin as follows:
Under the first two orders the shift register behaves as an ordinary fast storage register and is often used in this way as a “βtransfer register”; i.e. as intermediary in a transfer between two γ registers. A number is normalized by reading it into the shift register under β code 0301 and out under β code 0301. The amount of shift required for normalization appears automatically in the shift control register. After placing the desired amount of shift in the shift control register a number is shifted by reading it into the shift register under β code 0300 and out under β code 0301. Note that a shift occurs on readout only so that the original number stands in the shift register regardless of which readin or readout is used, and may be read out unshifted under β code 0300 in either case. Decimal Point Switch The decimal point switch is a β register storing a number in the range 016. This number determines the operating decimal point of the machine and represents the number of columns to the right of the decimal point. The decimal point switch is set manually and readout only is possible (β3333). Sign Register The sign register is manually set to store a positive or negative zero and is used to determine whether the decimal point is in the high or the low order part of a number in double accuracy operation. It is also a convenient source of zeros. Readout only is possible (β3332). Readin Switch The readin switch is a manuallyoperated regular register from which it is possible to read out (α07) to either a β or a γ register, or to both simultaneously. Call and Conditional Call The line number register (LNR) is a fourcolumn β register whose content determines the line number of the order to be performed. The line number register is automatically augmented by one as each order is performed so that orders are normally carried out in the sequence determined by their line numbers. The LNR has regular readin (β3010) which is used for calling. If, for example, the number 3255 stands in the γ transfer register then line 3255 may be called by the following order: 10 3010 00. Because the sequence orders are stored on a rotating drum and may be read out only once per drum revolutions (25 cycles) a call will introduce a wait until the called order passes under the reading head. The wait is such that any repetitive program loop consisting of t lines will take r drum revolutions if
Note that the minimum time for a short repetitive loop is one drum revolution and that a loop of twentyfour lines will take twice as long as a loop of twentythree lines. Waits due to calls may be greatly reduced by careful programming. The condition register stores only the sign of a number read into it and is controlled by β code 3012. A conditional call (β3011) performs the same operation as a call if the condition register is negative, and does nothing if it is positive. Note that in the machine x  x = 0. An order to read out of the line number register (β3010) reads out the line number of that same order. γ Transfer Register and External Transfer The γ transfer register is a fast storage register provided in the γ group and may be used as intermediary in any transfer between two β registers. This register also has a special function in the introduction of numbers in the course of a problem. This operation is called external transfer and may be performed in either of two ways:
In this way a complete positive sixteendigit number may be introduced by four orders, the first with α code 04 followed by three with &alpha code 05. Addition Addition is performed in either of two accumulators which operate independently except that they may be used together for double accuracy operations. Accumulator 1 is controlled by γ code 10 and 11 and two different readins are possible. Code 10 is a readin which destroys the number previously stored in the accumulator. Readout in either case is under code 11. The similar operation of accumulator 2 is described in the coding list. Example. The following orders add the contents of β registers 0019 and 0134 and place the result in β register 0022.
Note that any number subtracted from itself yields a negative zero. High Accuracy Addition Using the symbols HO and LO to represent the registers containing respectively the high order and low order parts of a given number, high accuracy addition may be described as follows. The first number to be added is read in with reset, thus:
As many additional numbers as desired may now be added by providing for each of them the two instructions
Finally, the sum may be read out by the following orders:
Multiplication The multiplicand (MC) must be read in first under γ code 20 and remains in the MC register. The multiplier (MP) is read in on the next (or any later) line under γ code 21. The product may be read out on the eighth (or any later) line following the readin of the MP and in the intervening seven lines any orders not involving the multiply unit may be performed. Four different readouts of the product are possible. γ codes 21 and 22 read the product out at the operating decimal point whereas γ codes 23 and 24 read out at decimal point 15. The high order (HO) product is the remainder of the 32digit product formed by the multiply unit and is used in high accuracy operations. The MP is destroyed in the course of multiplication but the MC remains and may be read out at any time except when a multiplication is in progress. Multiplication by a constant multiplicand may be effected by reading the MC in once followed by the successive multiplier. Division The divisor must be read in first, followed by the dividend within fifteen cycles. The quotient may be read out on the twentieth line following readin of the dividend and in the intervening nineteen lines any orders not involving the divide unit may be performed. The divisor and dividend are both destroyed in the process of division and may not be read out. Readout under γ code 34 is not possible following a readout under γ code 32. Vacuous Codes The vacuous order 17 3355 77 performs no operation and is used to provide any desired delay, as for example before readout of a product or quotient when it is impossible to interpose enough useful orders. The vacuous codes 3355 and 77 may also be used separately whenever a β (or γ) register is not involved in the particular order. Coded Stops Two stop codes are provided: β1312 stops the machine on a positive transfer and β1313 on a negative transfer. The two stops are provided to give more reliable operation on a coded check. Example. Suppose that two computed results stored in a_{0} and a_{1} are required to agree to within a certain tolerance stored in a_{2}. The following routine is used to stop the machine if the agreement is not as required:
Lines 14 place a_{0}  a_{1} in accumulator 1. Line 6 adds the tolerance, and so makes the accumulator 1 positive if the agreement is as required. Line 7 then stops the machine if the accumulator is negative. The reason for including line 5 is that the accumulator may be defective and standing in a permanently positive condition, thus passing the check regardless of the value of a_{0}  a_{1}. Lines 5 and 7 together show that the accumulator is negative at one point and positive at another thus indicating it is operative. InputOutput All data are read to and from the machine by means of magnetic tape units which perform the following functions:
β codes 230111 select tape units 111 and with a βin code will record on the appropriate tape the contents of any γ register desired. Read from tape is accomplished by the same β codes combined with a βout code, but the number is read always to the read registers, numbers 1 and 2. The read registers are β registers with normal readout and readin under β codes 2321 and 2322. It is also possible to select “tape unit k”, where k is the number (111) stored in the tape control (K) register. The functions are the same as described above. The K register stores numbers from 0 to 19 and has normal readin and readout under β code 2320. Read to and from the record register without recording is possible under β code 2323.
Printing If final output is to be in printed form, the recorded tapes must be used to control a printer. To accomplish this a “serial” number must be recorded before each number to be printed. This serial number serves to identify the number and to control the printer. The various columns (number right to left) of the serial number are used as follows:
Before printing the printer compares the first column of the following serial number with the preceding one. If the digit has increased by one printing proceeds, but if not the printer passes on to the next number without printing. Thus if any number has be rerecorded (due perhaps to a rerun) only the last recorded value is printed. All numbers are recorded in duplicate on the tape; the printer compares the two recorded values and stops if they do not agree. Waits Certain units which require several cycles to operate are provided with automatic waits to delay any order involving them until the previous operation is completed. This is not true of the multiplier and divider, where the requisite number of orders must be interposed between readin and readout. Even where automatic waits are provided, machine time may be saved by interposing other orders between successive uses of any one of the abovementioned units. The average time required for certain operations is given below in terms of machine cycles. A machine cycle is 1.2 milliseconds, the time required for a single order.
Notes Chapter 3 Linear Differential Equations A. Introduction The general solution of a system of ordinary linear differential equations is well known, as are certain computational methods for evaluating the parameters involved. In this chapter the theory is developed in a manner designed to provide the necessary background for the discussion of computational methods in Chapter 4. In order to make the structure of the solution clear for the case of multiple roots, the whole discussion is based on the Jordan canonical form of the matrix of coefficients. The complementary function is expressed in terms of the system of principal vectors, and expressions for the principal vectors are derived in terms of the adjoint of the characteristic matrix. A knowledge of elementary matrix theory such as contained in the first three chapters of reference 8 will be assumed. In matrix notation the most general form of a system of first order^{[a] } linear differential equations in n variables is
where A and B are square matrices of order n and
are vector functions of the independent variable t. The vector z is a given function and x is the vector to be determined. In the case of constant coefficients, to which this discussion is restricted, the elements of A and B are constants. The general solution of Eq. (3.1) consists of the sum of (1) a particular integral of (3.1), and (2) the complementary function which is any solution involving n arbitrary constants, of the associated homogeneous equation
If B is nonsingular Eq. (3.2) may be reduced to
where D = B^{1}A. The singular case is considered in Section D. The problem of the complementary function will be considered first and to this end certain results in the theory of matrices are first developed. In this chapter the following conventions are adopted:
The proof of the following theorem, omitted here,
may be found on page 64 of
B. The Complementary Function Equation (3.3) is repeated here for reference
Substitution of Eq. (3.9) gives
Using the fact that SIS^{1} = I and then premultiplying by S^{1} we have
Let y = S^{1}x. Then = S^{1} and finally Eq. (3.3) becomes
Substitution of the expression
in the lefthand side of Eq. (3.10) gives
The bracketed quantity is identically zero by virtue of (3.7) and therefore (3.11) is a solution of (3.10) for an arbitrary vector c. This is the general solution sought since it contains n arbitrary constants c_{1},c_{2},…,c_{n}. The expanded form of (3.11) is exhibited below for the particular
matrix^{[b]}
In terms of components:
The general solution of Eq. (3.3) is now given by
Since S is nonsingular x also contains n arbitrary constants. Let be the vector of initial conditions; i.e. at
C. Reduction to Canonical Form The theorem cited at the end of section A assures the existence of a nonsingular matrix S such that the similarity transformation S^{1}DS reduces D to canonical form. A process for obtaining a suitable matrix S for any given matrix D will now be derived. (1) Principal Vectors and Latent Roots A principal vector of order q is defined as nontrivial vector x satisfying the conditions
The corresponding value of λ is called the latent root, or more briefly the root, associated with the principal vector x.
Proof: In order to relate the unit vectors to the structure of J_{n}, the following notation is adopted:
Thus x_{kj} is a column vector with zero components except for a unit element in the jth position of the group of rows occupied by the kth Jordan box J_{nk}(λ). As k and j range over the integers 1 to m and 1 to n_{k} respectively, the vectors x_{kj} range over the set of unit vectors e_{1},e_{2},…,e_{n}. Clearly then,
For example, if k = j = 2, and J is the matrix of section 3A11,
If λ = λ_{k} then
Clearly then
Thus x_{kj} is a principal vector of order j associated with the latent root λ_{k} and the theorem is proved.
Proof: Let J_{n} be the canonical form of D. Then
Therefore
Therefore Sx_{kj} is a principal vector of order j of the matrix D and is associated with the root λ_{k}. The principal vectors of A are therefore the columns of S and since S is nonsingular, they are linearly independent. In other words, the matrix S required for the canonical reduction of D and hence for the solution of (3.3) is made up of the principal vectors of D arranged in descending order for each latent root as indicated by the relations (3.15). These relations carry over to the matrix D in the form
where S_{kj} is the rth column of S and r = j + n_{k}. Eigenvectors. A principal vector of order one is called an eigenvector. If in the canonical form of D, n_{k} = 1, k = 1,2,…,m, the principal vectors are all eigenvectors and they form a linearly independent set. The n_{k} are necessarily equal to unity if
The proof of (2) rests on the fact that if D is symmetric so is J. Therefore J has no offdiagonal elements. The significance of a linearly independent set of eigenvectors is that the solution of (3.3) then involves only exponentials and no terms of the form t^{k}e^{λt}. (2) The Adjoint Matrix It is easily verified that the pth power of a canonical matrix J may be expressed as
Hence for any polynomial function of J
Further, if J_{nk}(λ_{k}) is a Jordan box then
For example,
Formally applying the above results to the function J_{n}^{1} we have for each Jordan box
Premultiplication by J_{nk} gives
Therefore J_{n}^{1} as formally define by (3.20) is indeed the inverse of J_{n} and the preceding results may therefore be extended to negative powers of J. To avoid difficulties due to the possible singularity of the matrix J we consider instead the adjoint of J, to be denoted by Δ(J) and defined by
Since the determinant J is equal to
and since the largest negative power of λ occurring in J^{1} is λ_{k}^{nk} it follows that the adjoint is finite even for singular matrices. Since J and J^{1} commute
(3) Principal Vectors Derived from the Adjoint Since J_{n}(λ_{i},n_{i}) – λI = J_{n}(λ_{i}–λ,n_{i}) is a canonical matrix, substitution in (3.22) gives
If λ is chosen equal to any latent root λ_{k} the determinant on the right vanishes and therefore any nontrivial column of the adjoint Δ[J_{n}(λ_{i}–λ,n_{i})] is an eigenvector of J_{n} corresponding to the root λ_{k} (see Eqs. (3.14) and the definition of an eigenvector in section 3C1). The existence of such nontrivial columns will now be demonstrated. The form of a single box of the adjoint of (J – λI) is shown below for n_{k} = 3.
In the definition of J_{n}(λ_{i},n_{i}) in section 3A11 the symbol λ_{k} is used to denote the diagonal element of the kth Jordan box. This notation is not meant to imply that the elements λ_{1},λ_{2},…,λ_{m} are necessarily distinct. To this point the occurrence of two or more equal values of λ_{k} has not entered the discussion but a distinction must now be made. Let m_{k} be the multiplicity of the root λ_{k} in J. Thus if T is the set of latent roots equal to λ_{k} the multiplicity of λ_{k} is
In particular, if the λ_{k} are distinct m_{k} = n_{k}; and if a root λ_{k} is repeated in more than one Jordan box then m_{k} > n_{k}. The symbol will now be adopted for the product of the λ’s. The prime indicates that the factors are chosen only from a set of distinct λ_{k}. Then
The product in (3.23) may now be written as
Factoring out (λ_{j} – λ)^{mk} and multiplying it into the matrix, Eq. (2.23) appears as follows:
The same relation for a general value of n_{k} may be written concisely in terms of the matrices I_{k} as follows:
Suppose now that λ_{k} is a latent root
occurring in a single Jordan box so that
is some nonzero constant c_{0}. The kth box is also null except for an element equal to
in the upper right corner. This may be verified by referring either to the example of (3.24) or to the general expression (3.25). Therefore the only nonzero column occurring in the adjoint is the vector (1)^{nk1}c_{0}x_{k1}. Eigenvectors (and principal vectors) are defined by a homogeneous equation and therefore any scalar multiple of an eigenvector is also an eigenvector. For convenience, the nonzero scalar multiplier (1)^{nk1}c_{0} in the column of the adjoint may be dropped to give x_{k1} as the eigenvector corresponding to λ_{k}. This is the same result previously obtained from Eqs. (3.16). The remaining principal vectors associated with λ_{k} may be obtained by differentiating the adjoint with respect to λ. For convenience the scalar constants c_{p} are defined as
Differentiating (3.24) or (3.25) gives
or
Since c_{0} ≠ 0 the derived adjoint contains a nonzero column proportional to the eigenvector. If c_{1} = 0 it contains in addition a linearly independent column with a component proportional [to] x_{k2}, which is the principal vector of order 2 associated with λ_{k}. If c_{1} ≠ 0 this column contains also a component in the direction of x_{k1}. This brings out a point already implicit in (3.16); namely that a principal vector of order greater than 1 is nonunique to the extent that it may have added to it any linear combination of the associated principal vectors of lower order. The pattern should now be clear. The final column in the kth Jordan box of
is a principal vector of order p+1 associated with the eigenvalue λ_{k}. At most p other columns are nonzero and they are principal vectors of lower order. Note also that for p ≥ n_{k} the derived adjoint is null. The case where λ_{k} is repeated
in two or more boxes may now be disposed of briefly.
The element of lowest degree in (λ_{k} – λ)
which occurs in the adjoint is
(λ_{k} – λ)^{mknk}
as shown by (3.24) or (3.25).
In this case m_{k} > n_{k}
so that the adjoint evaluated at
(4) Principal Vectors of a General Matrix The results obtained above for a canonical matrix may now be extended to a general matrix A by the use of the following train of relations: If
and
Then
Since S is nonsingular Eq. (3.26) shows that Δ(A – λI) possesses the same number of independent columns as does Δ(J – λI). Also from (3.27) the functions
are seen to have the same zeros. From the definition of a principal vector and the results obtained for the canonical matrix, it follows that the derived adjoints of A contain columns proportional to the principal vectors of A. D. The Singular Case In passing from Eq. (3.2) to Eq. (3.3) it is necessary to assume that the matrix B is nonsingular. The singular case will now be considered. If the matrix B is of rank (nr) it is possible by rearranging and taking linear sums of the component equations of (3.2) to reduce them to an equivalent system in which the coefficients of the derivatives in the last r equations are zero. The last r equations may therefore be solved directly for r of the variables in terms of the remaining variables. These r variables may then be eliminated from the first (nr) equations. We are then left with a reduced system of equations in (nr) variables of the form
where B_{1} is nonsingular. Solutions of the reduced system may now be obtained by the methods previously described and the r eliminated variables are expressible in terms of these solutions. The general solution now contains only (nr) arbitrary constants and the system may be said to have lost r degrees of freedom. The initial values of the variables are therefore subject to r linear restraints as determined by those r equations which are free of terms in derivatives. If the matrix A is nonsingular an alternative method is available. Premultiplication of (3.2) by A^{1} gives:
where
Let J = S^{1}CS be the canonical form of C. Then Eq. (3.29) may be transformed to
where
Then
is the general solution of (3.30). [Compare Eq. (3.11).] If B is nonsingular so are C and J and hence (3.31) is a valid result. If B is singular J^{1} does not exist but the solution may still be made meaningful in the following way. Because C is singular it has at least one latent root equal to zero. Suppose it has exactly one zero root and let S be so chosen that this root occurs as the last element of J. Formally inverting J by the use of Eq. (3.20) the last latent root of J^{1} becomes infinite. The solution y = e^{J1t}c therefore has a discontinuity at t = 0 unless the vector c (i.e., the vector of initial conditions) is so chosen that it has no component in the direction of the eigenvalue of J corresponding to the zero latent root. A single linear restraint is thus imposed on the variables and at the same time one of the exponential terms in the general solution is dropped. This eigenvector of the canonical matrix J is of course the vector e_{n} (see Eq. (3.16)). The restraint is therefore simply c_{n} = 0 and hence y_{n} = 0. This is the same linear restriction arrived at by the previous method as may be verified as follows: The transformation
is one of the possible ways of reducing Eq. (3.2) to the form (3.28), for we obtain (3.30) namely,
Since the last row of the matrix J is null the system may now be reduced to the form (3.28) by dropping the last equation and using it as an auxiliary restraint. This equation is simply y_{n} = 0 as was to be demonstrated. The same reasoning may be used to extend this result to the case of multiple zero roots of B. We may therefore lay down the following procedure for the case when B is singular.
This result is very useful in the present problem since the matrix A is known to be nonsingular and the matrix B is in some cases singular. Also the inverse of the matrix A happens to be of interest in itself since A^{1} is the solution of the static system (1.1). E. The Particular Integral The particular integral of (3.1) is now easily obtained.
A particular integral is then given by
as may be verified by substitution in (3.32). Finally,
F. The Economic Model The equations of motion (1.2) of the dynamic economic model proposed in Chapter 1 may be expressed in matrix notation as
where A = [a_{ij}] is the matrix of flow coefficients and B = [b_{ij}] is the matrix of capital coefficients. As indicated in section 3D the associated homogeneous equation may be put in the form
where D = (I – A)^{1}B. The general solution is then expressible in terms of the latent roots and principal vectors of D. As indicated in Chapter 1 the final demand functions of interest in the present problem are of the form
where g is an arbitrary vector and μ is a constant. A particular integral of (3.34) then takes the simple form
This may be verified by substitution in (3.34) as follows:
Notes Chapter 4. Latent Roots and Principal Vectors In Chapter 3 the solution of linear differential equations with constant coefficients was reduced to the problem of evaluating the latent roots and principal vectors of a certain matrix. A number of possible computational methods will now be considered with a view to choosing those methods best suited to the solution of the problem outlined in Chapter 1. Methods such as that of Jacobi which apply only to symmetric matrices will not be considered. The complete solution of this problem requires the calculation of all the latent roots and all the principal vectors of the matrix (I – A)^{1}B. Since the matrix is not symmetric, complex roots may be expected to occur and indeed these complex roots are of particular interest in the theory of the business cycle.^{[a] } The method chosen must therefore be capable of extracting all roots, real or complex, and the associated principal vectors. A. Bounds on the Roots The rate of accumulation of roundoff error in some of the methods to be described depends critically on the magnitude of the latent roots of the matrix concerned. The latent roots also determine the magnitude of the largest number to be encountered in the course of the calculation.^{[b] } For two reasons then it is useful to have a preliminary estimate of the latent roots. The best such estimate available with a small amount of computation is furnished by Gergšgorin’s theorem^{11} which follows.
Proof: Let x be an eigenvector of A corresponding to the latent root λ. Then Ax = λx and
Let x_{k} = max x_{j} ≠ 0 since x is nontrivial. Then
Dividing by x_{k} we have
Thus λ lies in one of the circles C_{i} and the theorem is proved. Since A^{T} has the same latent roots as A a similar theorem holds for circles D_{i} with centers a_{ii} and radii
Clearly then the two bounds to the roots are given by
B. The Power Method (1) Single Dominant Root The power method in its simplest form can be used only to obtain the dominant root of a matrix and owes its popularity to the large number of problems in which a single dominant root is of interest. The method is based on the fact that for any vector x with a nonzero component in the direction of the largest eigenvector^{[c] } the sequence of vectors Ax,A^{2}x,…,A^{p}x approaches the largest eigenvector of the matrix A. The ratios of corresponding components of two successive vectors of the sequence approach the largest latent root of the matrix A. The proof of these statements is easily established with the use of the Jordan canonical form of A. For if
If λ_{1} is the dominant root of A then J^{p} has the form^{[d]}
where each submatrix is of the form
Since λ_{1} is dominant
Hence
where α is the first component of (S^{1}x) and e_{1} is a unit vector. Since the eigenvector of A corresponding to λ_{1} is Se_{1} by Eq. (3.17) and since x has a component in the direction of this eigenvector it follows that α is not zero. Thus
The expression on the right is an eigenvector of A. Furthermore,
The convergence is linear and depends on the ratios
For this reason troubles arise when two or more roots of equal modulus^{[e]} occur. This problem will be deferred and for the present the roots are assumed to occur singly. (2) Smaller Roots The simple power method may be extended to obtain roots other than the dominant root in either of two ways, which will be referred to as the modified power method and the method of deflation. The Modified Power Method. Let the roots of matrix A be denoted by λ_{1},λ_{2},…,λ_{n} where λ_{1} > λ_{2} > … > λ_{n}. Equations (3.18) and (3.19) show that if F is any polynomial function then the roots μ_{i} of the matrix F(A) are given by
In particular if
Clearly any one of the μ_{i} may be made dominant by choosing α close enough to the corresponding root α_{i}. Thus each value of μ_{i} (and from it λ_{i}) may be obtained from the matrix B by a suitable choice of α. In the absence of any knowledge of the location of λ_{i}, the required values of α must be obtained by repeated trials. The Method of Deflation. This method, like the preceding, modifies the matrix A so as to make a secondary root dominant, but differs in that the dominant root only is changed and is replaced by zero. If the roots of A are distinct and if c_{i} and r_{i} are respectively the eigenvectors of A and A^{T} corresponding to the root λ_{i} and with length so chosen that r_{i}^{T}c_{i} = 1, then
For
If i ≠ j then λ_{i} ≠ λ_{j} and hence r_{i}^{T}c_{j} = 0. Equation (4.2) is therefore established. If will now be shown that
Note that each product c_{i}r_{i}^{T} is a square matrix of order n. To establish the relation (4.3) consider the matrix product
Since this is true for each of the n linearly independent vectors c_{j}, it follows that the matrix in brackets above is null and (4.3) is established. Suppose the dominant root and corresponding eigenvector of A have been obtained. Then the corresponding eigenvector of A^{T} may be obtained by reapplying the power method, or more simply (since λ_{1} is now known), by solving the singular homogeneous system
to obtain a nontrivial vector x. The “deflated” matrix
is now formed. Clearly B has the same roots as A except that the dominant root λ_{1} has been replaced by zero. The root λ_{2} of B may now be obtained by a further application of the power method. The process may be repeated for successive roots but the error accumulation is rather rapid. Other methods of carrying out the deflation process are discussed by Feller and Forsythe.^{14} If A is symmetric r_{i} = c_{i} and a separate calculation of r_{1} is not required. (3) Multiple Roots If multiple roots occur two cases must be distinguished:
If the matrix A is nondefective there are two or more independent eigenvectors corresponding to the same root. Furthermore any linear combination of these independent eigenvectors is itself an eigenvector. Thus the power method will converge to the dominant latent root but the eigenvector obtained will not be unique. By starting from different initial vectors the process will in general converge to different possible eigenvectors associated with the dominant root so that finally all of the independent eigenvectors are obtained. If the matrix A is defective the form of a power of the Jordan box given in (3.19) shows that a single limit of the sequence Ax,A^{2}x,…,A^{p}x exists. For example, if λ_{1} is dominant and n_{1} = 3 then
Therefore
Thus the sequence Ax,A^{2}x,…,A^{p}x converges always to the eigenvector associated with the largest root. The power method therefore fails to obtain the principal vectors of order greater than one. C. Methods Employing the Characteristic Polynomial The characteristic polynomial of a matrix A is defined as
The equation P(λ) = 0 is called the characteristic equation of the matrix A. P(λ) is a polynomial of degree n and for a canonical matrix the expansion of the determinant clearly yields
Thus the zeros of P(λ) are the latent roots of J and each root occurs with the same multiplicity in P(λ) as in J. This also holds true for a general matrix A since similar matrices have the same characteristic polynomial, as will now be demonstrated. Suppose A and B are similar matrices, then
and since
Clearly then the latent roots of a matrix may be obtained in the following two steps:
Since the degree of P(λ) is the same as the order of the matrix the numerical solution of the characteristic equation may be very difficult. The problem of solving polynomial equations is taken up in Chapter 5. We now proceed to consider two closely related models for generating the characteristic polynomial of a matrix. (1) The Method of Bingham ^{15} Let
be a polynomial with zeros λ_{1},λ_{2},…,λ_{n}. Defining
it is possible with the help of Newton’s identities (reference 16, p. 147) to express the coefficients of P(λ) in terms of the quantities s_{k}, k = 1,2,…,n as follows:
The quantities s_{k} corresponding to the characteristic polynomial of a canonical matrix J are easily obtained, for clearly the sum of the diagonal elements of J is equal to s_{1}. Furthermore, Eqs. (3.18) and (3.19) show that the sum of the diagonal elements of J^{k} is equal to s_{k}. The sum of the diagonal elements of a matrix A is usually called the trace of A or more briefly Tr A. Using this term we may write
Moreover, the trace is invariant under a similar transformation for if
then in terms of components
Consequently if (4.6) is true for the canonical form of a matrix A it is also true for A. Hence
To summarize, the coefficients of the characteristic polynomial of the matrix A may be obtained in the following steps:
(2) The Method of Frame ^{17} The Latent Roots. Let c_{0},c_{1},…,c_{n} be the coefficients of the characteristic polynomial of the matrix A and let F_{j}(A) be a sequence of polynomial functions of A defined by
Now
and in general
by virtue of (4.5).
Using this expression for c_{k}, Eqs. (4.8) may be written
Equations (4.9) define the sequence of operations in the Frame method. The Bingham process computes the sequence of matrices A,A^{2},A^{3},…,A^{n} and from the traces s_{1},s_{2},…,s_{n} of these matrices computes the coefficients by the use of Newton’s identities. The Frame method also uses repeated multiplication by A but at each step modifies the matrix so as to derive the coefficients directly without explicit use of Newton’s identities. The Frame method has the further advantage that
This relation serves as a useful overall check on the accuracy of the computation. The Principal Vectors. The adjoint of (A – λI) may be expressed as a polynomial in λ with coefficients formed from the sequence of matrices (4.8) as follows:
by (4.8) and since F_{n} = 0 then
by the definition of the characteristic equation.
Thus Δ(A – λI) as given by (4.11)
is indeed the adjoint of D. Error Analysis The power method is an iterative process and error accumulation is therefore not a problem. In the use of the Frame method however the error accumulation may be expected to be rather rapid and a bound on the error will now be established. Assume that A is nondefective and consider the series of matrices F_{k}(A) defined by (4.8). Because of roundoff error the computed value _{k}(A) will differ from F_{k}(A) and we may write
The propagation of the error introduced at the first step will be investigated first. Let the length of the longest column vector of E_{k} be denoted by ε_{k}. Then
Since for any vector x,
where λ is the dominant root of A, then
The error computed above is due to the roundoff committed at the first step. Since a similar roundoff occurs at each step the total error bound ε will be a sum of such terms as follows:
For large values of n the error accumulation may therefore be very severe but in work with matrices the maximum error bounds are usually found to be too pessimistic. E. Improving an Approximate Solution Let A be a nondefective matrix and let S be the matrix which diagonalizes A. If S_{0} is an approximation to S then the method of Jahn^{18} may be used to improve this approximation. The object is to compute S where AS = SΛ and Λ is a diagonal matrix. Since S_{0} ≐ S we have
where Λ_{1} is diagonal and B is a matrix
whose elements are small;
in particular the diagonal elements are zero.
Let the diagonal elements of Λ_{1} be
λ_{1},λ_{2},…,λ_{n}
and let
Then S_{1} = S_{0} (I + C) is an improved
approximation to S correct to second order.
For clearly
Hence
F. Comparison of Methods The amount of computation required to generate the characteristic equation by either the Frame or Bingham method is of the order n^{4} multiplications. In addition the solution of the equation so obtained will require considerable labor. Consequently if only a few of the roots are wanted, use of the power method is indicated. If a complete solution is desired, the power method has serious shortcomings. Firstly, the process is not uniform since special modifications are required to care for special cases such as complex or multiple roots. Furthermore in the case of a defective matrix it fails to obtain the principal vectors of order > 1. The attractive features of the Frame method are its complete generality and uniformity. Unlike the power method the process proceeds in a purely routine manner independent of the type of roots encountered. The completion of the method of course assumes the solution of a polynomial equation of degree equal to the order of the matrix and some of the essential difficulties may reappear in this phase of the problem. In any case the solution of polynomial equations is an important problem in its own right and the development of suitable methods is worth considerable effort. It therefore seems reasonable to base the solution of the problem of latent roots on the assumption of the solvability of the characteristic equation. If the roundoff error in the Frame method should be so severe that the desired accuracy is not attained, the method of Jahn may be used to improve the accuracy. Each application of the Jahn method requires the equivalent of about four matrix multiplications. Notes 5. Zeros of a Polynomial A. Methods Depending on an Initial Approximation The problem of locating the zeros of a polynomial of high degree is greatly simplified when good initial approximations to the zeros are known. This situation is often met in practice because of the frequent occurrence of families of polynomial equations. If, for example, the solution of a polynomial equation is required within some larger iterative process, the coefficients may be expected to vary gradually from one step to the next and the set of roots obtained at one step may be expected to provide good approximations to the roots of the next equation. Under these conditions, the BirgeVieta^{19} method and the quadratic factor method are probably the best, since they converge quadratically near a root. Before proceeding to a discussion of these methods certain useful processes will be developed. 1. Synthetic Division^{16 } Let P_{n}(x) = a_{0}x^{n} + a_{1}x^{n1} + … + a_{n1}x + a_{n} be a polynomial of degree n with real coefficients. Let
where the coefficients b_{i} are defined by
Then
Therefore
Hence b_{n} is the value of the polynomial evaluated at p and when b_{n} is zero, Q(x) is the socalled reduced equation. The effect of the algorithm (5.1) is to divide the polynomial P(x) by the linear factor (x–p) to obtain the quotient polynomial Q(x) and the remainder b_{n}. The algorithm is therefore called synthetic division. Adopting the convention that all coefficients with negative subscripts are zero, Eq. (5.1) may be written more concisely as
Synthetic division by factors of any degree is easily defined but the case of a quadratic factor is of particular interest. In this case let
Then
as may be verified by expanding the lefthand side. 2. Change of Origin A translation or change of origin to the point
x = p is defined by the substitution
The coefficients b_{k} are most conveniently computed by the process known as Horner’s method.^{16} This method is nothing more than a repeated application of synthetic division to the successive quotient polynomials. For n = 4 the scheme of computation is as follows:
The polynomial in y is then
To prove that the polynomial in y may be so obtained let the successive quotient polynomials Q_{k}(x) be defined by
Then
Repeating the differentiation of the relation
yields the values of the successive derivatives of P_{n}(x) evaluated at the point p. Thus
3. The BirgeVieta Process The socalled BirgeVieta process for locating a real root of a polynomial was first proposed by Vieta in 1600 and was discarded by later mathematicians as requiring too much computation. It was not until the advent of the modern desk calculator that the method was reestablished by Birge and came into wide use. The method is perhaps the most satisfactory known for obtaining real roots with a desk calculator. The BirgeVieta method is based on the wellknown NewtonRaphson process (reference 20, page 84) which is defined as follows: Let be a real roots of the equation
Then for any value x_{0} sufficiently close to the quantities x_{i} defined by
will converge to the root .
The convergence is quadratic; i.e., the difference
The process is then repeated using the quantity
and in general
Equation (5.8) defines the BirgeVieta process. This process is usually employed only for real roots, but may be extended to complex roots by simply allowing the variables x_{0},x_{1},… to be complex. Since a polynomial is an analytic function^{21} the derivative P_{n}'(x) exists for x complex, and the process carries over without change to complex numbers. (Note, however, that if the initial approximation x_{0} is real, the successive x_{i} remain real.) Every operation is now complex and hence the amount of computation is increased by a factor of about four. The complex BirgeVieta method can be applied equally well to the solution of polynomial equations with complex coefficients, but for real polynomials the method to be described below is generally considered to be superior. 4. The Quadratic Factor Method Since the complex zeros of a real polynomial occur in conjugate pairs it is possible to factor any real polynomial into a set of real linear and quadratic factors. Using a method similar to the BirgeVieta process it is possible to improve an approximate quadratic factor by repeated synthetic divisions of the type described by Eq. (5.4). The goal is to so adjust the coefficients p and q as to make both remainders b_{n1} and b_{n} zero. Treating b_{n1} and b_{n} as functions of p and q and expanding in Taylor’s series we have
A generalization of the NewtonRaphson process to two variables then gives the desired corrections δp and δq as solutions of the pair of linear equations
The quadratic factor method as first proposed by Bairstow^{22} and as presented in standard works^{23, 24} on numerical analysis requires the use of a second synthetic division to evaluate the derivatives in (5.9) and the application of the corrections δp and δq so derived. A modification of the method, which provides a more uniform computational scheme^{[a]} and simplifies the extension to a higher order process, will now be described. The proof is also simplified. We defined a modified algorithm for synthetic division, similar to that given by (5.4), namely,
This process is completely uniform, without exception for the computation of c_{n}. Comparing with (5.4)
letting
Eq (5.5) becomes
The modification of the quadratic factor method now consists in choosing δp and δq so as to reduce to zero the quantities c_{n1} and c_{n} rather than b_{n1} and b_{n}. Note that the vanishing of c_{n1} and c_{n} is the necessary and sufficient condition for the vanishing of b_{n1} and b_{n}. Treating x, p and q as independent variables and differentiating (5.11) partially with respect to p we have
Thus ∂c_{n1}/∂p and ∂c_{n}/∂p
are respectively the first and second
remainders on applying the modified algorithm
to the polynomial
Differentiating (5.11) with respect to q we have
and the first and second remainders on applying the modified algorithm to Q(x) are ∂c_{n1}/∂q and ∂c_{n}/∂q respectively. Hence
Applying Eq. (5.9) to the quantities c_{n} and c_{n1} and substituting the expressions obtained above for the derivatives we have
Thus:
The complete uniformity of the method is clearly
show in the following example for
δp and δq are then given by (5.14). In comparison with the complex BirgeVieta method the quadratic factor method is seen to require only half as many multiplications per iteration since it requires only two real multiplications per step while the complex BirgeVieta method requires one complex^{[b]} multiplication per step. Because the coefficients computed in the complex BirgeVieta process are in general complex, two registers are required to store the real and complex parts. The storage required is therefore double that required for the quadratic factor method. On the other hand, some additional calculation is required to resolve the quadratic factors obtained to yield the actual roots. The time required for this additional calculation is, however, rather insignificant. 5. Higher Order Processes The NewtonRaphson process for obtaining
the solution of the
equation
after the second term. In a similar way higher order processes may be derived by considering further terms of the expansion. Thus a third order process is obtained by solving
for . Then
where the smaller of the two possible correction terms is chosen. In the BirgeVieta method the convergence is already quadratic near a single root and little is gained by a higher order process. In the case of a double root the convergence of the BirgeVieta method is linear and may be made quadratic by the use of (5.15). Use of the higher order process may therefore be advisable if multiple roots are encountered. Similar remarks apply to the application of a higher order process based on the quadratic factor method. The value of f "(x_{0}) is given by (5.7) and may be obtained by a third synthetic division. To extend the quadratic factor method to a higher order process we adopt a new notation for the remainders as indicated in the scheme below.
Successive lines of the scheme correspond to successive applications of the modified algorithm. Using the notation
Eqs. (5.12) and (5.13) become
Because of the complete uniformity of the algorithm
the relation between any two quantities in the scheme above
is determined by their relative position.
For example, since G and C stand in the same relation
as do D and A and since
Extending the Eqs. (5.9) to include second order terms of the Taylor’s series and using the present notation we have
Hence
This pair of nonlinear equations may be solved by iteration providing that δp and δq are not too large. Each approximation to the values of δp and δq are substituted in the righthand side and the resulting linear equations are solved for a new approximation. The initial values of δp and δq are taken as zero. 6. Factors of Higher Degree Aitken^{25} has studied the extraction of real factors of higher degree than second by means of methods called “penultimate remaindering”. The modified quadratic factor method may also be extended to extract factors of higher degree. Unless all roots are known to be real the real factor must be of even degree. The case of a quartic factor will be developed briefly. As for the quadratic factor the synthetic division algorithm
yields the quotient polynomial
such that
Again the modified algorithm
yields the same quotient polynomial Q(x). Adopting the notation
we have
Note that the vanishing of the quantities A, B, C, and D is the necessary and sufficient condition for the vanishing of the remainders b_{n}, b_{n1}, b_{n2} and b_{n3}. Differentiating Eq. (5.17) we have
Rearranging the last equation to the form
shows that the application of the modified algorithm to Q(x) leaves D_{s} as the first remainder. Thus if the scheme of two successive applications of the algorithm is shown as
we have D_{s} = K. The coefficient of x_{2} in (5.21) is
But from the form of the remainders in (5.17) it is clear
that the coefficient of x^{2} is
Similarly B_{s} = I and A_{s} = H. Rearranging (5.18) to the form
shows the lefthand side to be the polynomial formed of the coefficients c_{0},c_{1},…,c_{n4}, D, C, B, and therefore the derivatives may be obtained as before but with K, J, I, and H replaced respectively by H, G, F, and E. Thus D_{p} = H, C_{p} = G, B_{p} = F, A_{p} = E. In a similar way derivatives with respect to q and r may be obtained from Eqs. (5.19) and (5.20) and in matrix notation we have finally
Expanding A, B, C, and D in Taylor’s series and retaining only linear terms we have (compare Eq. (5.9))
Solutions of the fourth order linear system (5.22) yield the necessary corrections. This method removed four roots at a time but because of the greater complexity the quadratic factor method is probably to be preferred. When finally obtained, the quartic factor must itself be resolved. A method of solution is considered below. 7. Solution of the Quartic^{[c]} The following method is easily programmed for an automatic computer and could be used in conjunction with the method discussed above. However, it is probably of greatest value as a method for desk calculators and is believed to be superior for this purpose to existing methods. The method is capable of extracting either real or complex roots but works essentially well if all four roots are complex. Given a real quartic
two real quadratic factors are assumed to be
Then
Setting
satisfies all but the second of the equations above. Imposing the further condition
satisfies (5.24) and the values of a, b, α, β so obtained furnish the desired solution of the quartic in terms of its quadratic factors. Equation (5.28) may be solved for β by the method of false position, using the fact that β is real and less than √D in absolute value. If all the roots are complex than β is also positive and (5.28) has a single real root in (0,√D). The computation of ø(β) is carried out by evaluating b, α and a according to (5.27) and substituting in (5.28). The method fails only if b  β = 0.
This case can arise only if
Example.
A useful modification of the method consists in first making the substitution x= y + h with h = A/4 so that the resulting equation in y contains on term in y^{3}. Considering the equation in y we now have A = 0 and Eqs. (5.27) and (5.28) simplify to
The special case mentioned above now exhibits itself
immediately by the fact that An attempt to use a NewtonRaphson process to solve Eq. (5.28a) leads to a difficult expression for dø(β)/dβ. If, however, ø is expressed as a function of the new variable
Although the method works for either real or complex roots the behavior of ø(γ) is greatly simplified by the assumption that all roots are complex, for then ø(γ) has exactly one real zero and since b and β are both positive γ > 2D^{½}. Also ø'(γ) < 0 and ø"(γ) > 0 for γ > 2D^{½}. Thus ø(γ) has one real zero and no minimum or inflection point in the interval (2√D,∞). Because of these properties the zero of ø(γ) is very easy to locate. Using a machine which will store d as a constant multiplier
the only quantities necessary to record are γ and ø.
However since the final value of Example. The substitution x = y  0.6376 in the quartic of the previous example gives
Computing C^{2}, C^{2}/2 and 4D we have
A few rough mental calculations show that ø(50) > 0 and ø(52) < 0.
Hence b  β = √416.2704 = 20.4027,
Thus
The quartic formed from these roots agrees with the given quartic to four decimal places in each coefficient. Since
is easily computed it is possible to use the higher order correction
where the sign is chosen to give the smaller correction.
Application of this in the preceding example gives
8. Solution of the Sextic An extension of the method of real factors to the sextic yields a practical computational method as follows: Assume that the sextic
has real factors x^{2} + αx + β and x^{4} + ax^{3} + bx^{2} + cx + d. Again assuming that A is made zero by translation we obtain the following relations
Then if d = F/β, γ = d  β^{2},
and β is chosen as any real root of
the values of α and β so obtained satisfy Eqs. (5.29). The remaining coefficients are then furnished by the relations
At least one root of (5.31) lies in the interval (F^{1/3},F^{1/3}). If all the zeros of P_{6}(x) are complex the zeros of (5.31) are all positive and exactly three in number. Hence either the interval (0,F^{1/3}) or the interval (F^{1/3},∞) contains one real zero only. The two values of α given by Eq. (5.30) introduce a complication but it is usually easy to recognize which value should be chosen. The method fails only if d  β^{2} = 0,
in which case β = F^{1/3} and
(since α is finite)
the expression for α being obtained as the limiting form of (5.30). Example:
Since 0 ≤ β ≤ F^{1/3} ≅ 22 choose β_{0} = 10. Linear prediction is used even if two successive values of ø do not differ in sign. R^{2} represents the square of the radical.
B. Initial Approximations 1. Method of Search Any point in the complex plane from which one of the methods of section A will converge to a particular zero of given polynomial will be considered a satisfactory initial approximation to that zero. If the complex plane is covered with a set of points on a sufficiently fine mesh then at least one point of the set will be sufficiently close to each zero of a given polynomial so as to serve as a suitable initial approximation. If bounds on the roots are known, the number of mesh points to be considered is finite and a program which will choose in succession the points of this finite set may therefore be used to provide initial approximations. Since the fineness of the mesh required is not generally known in advance it is advisable to cover the region of interest first on a coarse mesh and repeat on successively finer meshes until all zeros are located. This approach was employed in locating the zeros of the truncated series for the exponential^{26} namely
for values of n ≤ 23. Approximations were chosen on a rectangular mesh and each was introduced as an initial approximation to (a) the BirgeVieta process if the point lay on the real axis, or (b) the quadratic factor method if the point was complex. This method was found to be extremely slow for values of n > 10. In order to study the behavior of the process, an optional routine was included in the program to print out each successive approximation if desired. This record of successive approximation showed two reasons for the slowness of the process. (1) In general a very close first approximation was required for convergence. This meant tat the search had to be conducted on a very fine mesh. (2) The process would often follow blind leads with the remainders continuing to decrease for many iterations before finally diverging. Because of the computing time required this method was discarded as impractical. 2. Graeffe’s Method^{27 } Let r_{1},r_{2},…,r_{n} be the zeros of the polynomial
The coefficients are expressible in terms of the zeros as follows:^{16}
If the zeros are real and widely separated it is clear that
and in general
Thus a polynomial with widely separated roots breaks into a set of linear equations whose zeros are approximately the zeros of the original polynomial. More generally if a set of s zeros of equal magnitude occur and the rest are widely separated the polynomial breaks into (n–s) linear equations and a cluster of (s+1) successive coefficients defining a polynomial of degree s whose zeros are (approximately) the s zeros of equal modulus. The Graeffe process consists in generating a sequence of polynomials such that the zeros of any one of the sequence are the negative squares of the zeros of the preceding polynomial. Applications of the process k times raises the original zeros to the 2^{k} power. Hence if all zeros are distinct a few iterations serve to separate them widely. If b_{0},b_{1},…,b_{n} are coefficients of the polynomials so generated, we have
Hence r_{1},r_{2},…,r_{n} may be obtained. More generally if r_{p+1},r_{p+2},…,r_{p+s} are a set of s roots with equal moduli then the polynomial
has the s zeros, r_{p}^{2k},…,r_{p+s}^{2k}. Therefore
Thus the moduli of all the roots are obtainable. Certain refinements of Graeffe’s method allow the zeros r_{1},r_{1},…,r_{n} themselves to be evaluated. If, however, some of the zeros are complex this involves essentially the solution of the polynomial equation
and so leads back to the original problem of solving a polynomial equation of high degree if s is large. The rootsquaring process is obtained by rearranging the polynomial equation
so that odd powers of x appear on one side of the equation and even powers on the other. Squaring this relation gives
The substitution y = x^{2} then gives
where
Equation (5.35) defines the Graeffe rootsquaring process and the roots of (5.34) are the negative squares of the roots of (5.33). C. A General Method of Solution The problem of finding all the zeros of a polynomial of high degree is of very frequent occurrence and a satisfactory general method of solution is greatly to be desired. The standard methods are unsatisfactory however because of their failure in certain special cases. For example, the BirgeVieta and quadratic factor methods may fail to converge unless a very good first approximation is available; the method of false position and Sturm’s functions fail to locate complex roots; and methods such as that of Lin^{28} which are included under the process of penultimate remaindering converge only under certain conditions as established by Aitken^{25}. The method of Graeffe also fails in certain cases as pointed out in the previous section. A mathematician using a desk computer is not seriously troubled by these special cases since they may be fairly easily recognized and an intelligent choice of a different method usually resolves the difficulty. Conceivably a computer may also be programmed to recognize the special cases and choose a new method accordingly. However, decisions of this type which may appear simple to a human are exceedingly difficult to program and a single general method is to be preferred. Moore has considered this situation in some detail^{29} and has suggested a possible general method of solution. Moore’s method is considered impractical because of the amount of computation required. The number of multiplications required is of the order n^{3} (where n is the degree of the polynomial) and each involves complex numbers represented with a floating decimal point. A general method of solution is now proposed consisting of the following five steps:
Since the zeros are known to lie on each of the two families of circles described in step 4 it is clear that the first four steps alone could provide solutions. Such a method would however suffer from the following disadvantages:
If, however, the method is used only to supply initial approximations for step 5 the above disadvantages cause no difficulty. For clearly, the intersections not corresponding to zeros will be rejected and the genuine approximations will be improved to the accuracy desired, by a method requiring only 2n or 4n multiplications per iteration and possessing quadratic convergence. In general the Graeffe process should be repeated only enough times that the given approximations will converge to the desired zero. No general rule can be given for attaining this^{[d]} but suppose the accuracy required is such that two zeros differing by less than one part in hundred need not be separated, i.e. they will appear as a pair of roots with moduli equal to the geometric mean of their true moduli. Since (0.99)^{29} ≐ 0.005 it is clear that two roots with moduli in the ratio 0.99 will after 9 iterations have the ratio 0.005. In practice fewer than nine iterations may well be used, the number being increased if difficulties are encountered. The choice of the parameter d in shifting the origin
is arbitrary but should be made a function
of the moduli of the roots.
The choice of Since one iteration of the Graeffe process
and a change of origin each require n^{2}/2
multiplications,
the whole process described
(assuming k applications of the Graeffe process
and neglecting the computation of the intersections
of the circles) requires in all
Because of the wide range in magnitude of the numbers derived from Graeffe’s process it is clear that floating point operation will be required. Floating point operation is probably to be desired in any case since careful scaling of the coefficients to avoid overcapacity numbers is rendered unnecessary. The advantages offered by the proposed method may now be summarized:
Further discussion of the method will be found in Chapter 6 where suitable checks and other programming details are considered. Notes Chapter 6 Machine Programs A. Summary of Calculations The general solution of the model proposed in Chapter 1
has been shown in Chapters 3 and 4 to be expressible
in terms of the latent roots and principal vectors
of the matrix The solution was obtained by the following steps:
All programs were made to handle a system of general order and in general the calculations for the six systems of order 5, 6, 10, 11, 20, 21 were run concurrently. This mode of operation had the advantage that the programs could be checked out on the low order system. The programs employed will be examined in some detail in the following sections. The desirability of uniformity in machine programs discussed in section 2A will be stressed throughout. B. Matrix Subroutines 1. Introduction Since it was not clear at the outset what matrix methods might be required, it was decided to organize the programming around a number of basic subroutines which were designed to carry out all of the more important matrix operations likely to be encountered. In this way any method to be employed would consist of a relatively simple main program piecing together the subroutines required. So far as possible the use of the subroutines should reduce the programming of matrix operations to the simplicity of the programming of the corresponding algebraic operations. In designing the subroutines a number of requirements were laid down. These requirements are listed below:
Since the fast storage registers can accommodate a matrix of order fifteen at most, the first requirement makes the use of slow storage registers necessary. A vector x = [x_{1},x_{2},…,x_{n}] will be said to be stored in location b if the components x_{1},x_{2},…,x_{n} are stored respectively in the registers b+1,b+2,…,b+n. The negative sum of the components will be designated by x_{0}, and will be stored in register b. In this way the location of the vector may be completely specified by the order n and the number b. Also each vector carries with it its own check since the sum of the elements x_{0},x_{1},x_{2},…,x_{n} must be zero. The use of the negative sum makes the process of summing the check more uniform. Each row of a matrix A = [a_{ij}] may be considered as a vector and will be stored as indicated above. The collection of negative row sums then forms a “zeroth” column of the matrix. In the same way negative column sums are included as a zero row and finally the negative sum of this row completes the (n+1)×(n+1) square array:
where
Note that
by either definition above. Sum checks prove extremely useful in linear operations and for this reason all matrices (including and output on magnetic tape) are stored as indicated above. The simplest use of the sumcheck is for checking the components of a vector each time it enters into the computation. As an example of other ways in which the sumcheck may be used, consider the matrix multiplication
which is normally defined as
Since the matrix A is stored with a zeroth row
this definition may be extended to
The quantity in parentheses is zero for each j so that finally
This last is an identity which checks the whole matrix multiplication. Note the complete uniformity of the operations. In machine calculation the roundoff error may cause a discrpancy of as much as (n+1) in the last column in the identity (6.1). Hence a tolerance of (n+1) should be allowed in checking. It should, however, be standard practice to then “correct” the sum y_{0} so that
In this way the only error to be allowed at each step is the maximum roundoff error of that step. Otherwise the error accumulates and it becomes difficult to estimate what tolerance should be allowed. A matrix may be stored most compactly by assigning the first element to the register following the last element of the proceeding row. It is necessary, however, to space the rows in slow storage so that the transfer of a given row to slow storage will not destroy numbers of the following row.^{[a]} The rth row of a matrix A beginning at location b (which stores a_{00}) is therefore stored at location b+rs where s is the minimum number of blocks of then registers required to store (n+2)^{[b]} numbers. The number s is computed from n and is stored in register 198. To illustrate the foregoing consider the matrix A stored at slow storage location 1830 where
Then n = 3, s = 1 and the complete matrix is stored as shown below where the numerical quantities on the left of the arrow are stored in the corresponding registers indicated on the right.
The fast storgage registers 190199 are utilized as follows:
In order to meet the requirement of flexibility without increasing the number of orders required in calling the subroutine, the following device is used. If a certain control number (e.g., the row number r) is essentially positive then the positive absolute value may be used where necessary in the routine. Thus in calling the subroutine this control number may be entered either positivelly or negatively, and at some point in the subroutine the sign of the number may be used to make a choice between two alternative operations. Examples will be found in the discussion of the individual subroutines. Although checks are prepared in the subroutine (e.g., the components of a vector transferred from slow storage will be summed), the final check is left to the main program if possible. This makes it possible to program reruns instead of checkstops, if desired, and also makes it easier to determine at what point in the main program a failure has occurred. 2. Floating Vector Operation Floating point operation (see section 2A) offers decided advantages in reducing roundoff error and in avoiding overflow.^{[c]} … 3. Matrix Programs C. Matrix Inversions 1. A New Machine Method 2. Inversion of a Complex Matrix D. The Method of Frame E. Solution of the Characteristic Equation 1. The Graeffe Process 2. The Quadratic Factor Method F. The Principal Vectors If the Frame method (section 4C) is applied to a matrix A to obtain the matrices F_{0},F_{1},…,F_{n}, then Eq. (4.11) may be used to evaluate the adjoint Δ(A  λI) of the matrix A  λI. Denoting the ith column of the matrix F_{j} by f_{ij}, the ith column ν_{i}(λ) of the adjoint is given by
If λ_{k} is a latent root of the matrix A, any nontrivial columns of the adjoint Δ(A  λI) are proportional to the associated eigenvector. Hence it usually suffices to compute only a single column ν_{i}(λ_{k}) for each distinct root λ_{k}. If multiple roots occur the required number of independent columns of the derived adjoint may be computed by differentiating Eq. (6.10). The polynomial (6.10) may be evaluated most efficiently by grouping the operations as follows:
This is simply an extension of the synthetic division process to a polynomial with vector coefficients. Because of the extreme range of magnitudes to be encountered in the quantities f_{ij}^{[s]} the calculation requires floating point operation. The roots λ_{k} in general are complex so that floating point complex operation is required. Such operation is extremely slow and may be avoided by the use of floating vector operation as follows: The root λ_{k} is first put in the form λ_{k} = r_{k}μ_{k} where r_{k} is real and μ_{k} = 1. Equation (6.11) may now be written as
The powers of r_{k} are formed using
floating point operation and the vectors
r_{k}^{n1j}f_{ij},
Notes Chapter 7 Results and Conclusions The tables of initial data and final results appear in Appendices IVII and are described in section A of this chapter. The chapter is concluded with remarks of a general nature in section B. A. Table of Results Equation (1.2) of Chapter 1 which describes the economic model under consideration is repeated here for reference:
As described in Chapter 1 the variables x_{1},x_{2},…x_{n} represent the outputs of each of n industries making up a given economy. Table 27 of Appendix VII shows the industry classification^{[a] } adopted for three systems of order 21, 11 and 6. For the 21 industry classification the numbers 121 indicate which of the variables x_{1}x_{21} are assigned to represent each industry; e.g. x_{5} represents the output of the chemical industry. For the 11 and 6 industry classifications the converging lines indicate which industries of the 21 country classification are grouped together to be represented as a single industry in the lower order classification. Three open^{[b] } systems of order 20, 10 and 5 are derived from the foregoing systems by suppressing in each the variable representing households and considering it instead as a final demand on the remaining system. Tables 226 contain square matrices of various orders. Except in Tables 813 of Appendix III where the row and column number of each element is shown directly, the columns of the table represent the columns of the corresponding matrix. Because sum checks prove as useful in hand computation as in machine computation, each matrix is recorded with a zeroth row containing the negative sums of the columns and a zeroth column containing the negative sums of the rows as described in section 6B1. The matrix A = [a_{ij}] of flow coefficients derived^{[a] } from data supplied by the Bureau of Labor Statistics is shown in Appendix I. The corresponding matrix of capital coefficients B = [b_{ij}], supplied by the Harvard Economic Research Project, appears in Appendix II. The broken lines in each of these tables (27) indicate the submatrix of coefficients which is retained in the open systems considered. The column and row sums must of course be adjusted for these systems. Four decimal places were given in the initial data (Tables 27) and all results are rounded to four decimal places. Twelve decimal places were carried in the calculation of Tables 819. A final check on the results shown in Tables 2024 and the agreement obtained are described at the end of this section. To avoid repetition, the use of the tables in the solution of the economic model will be described only in terms of the system of order five and all references to tables will be to the fifth order system. Expressed in matrix notation and specialized to the particular final demand functions to be considered, Eq. (1.2) appears as
The general solution of Eq. (7.1) is a sum of
(See Eqs. (3.31) and (3.12).) The general solution is therefore
where g is a given vector determining the demand function and c is a vector to be determined by the initial conditions = x(0). Thus
The matrix^{[c] } S
is given in Table 20, the inverse S^{1}
in Table 25, and the matrix The only quantity now left undefined in Eq. (7.2)
is the matrix J.
Let us denote the kth eigenvalue of the matrix
Finally e^{J1t}c is the column vector
If one of the latent roots λ_{k} is zero then the results must be interpreted according to section 3D as follows: The corresponding component c_{k} of the vector c is made zero. This in turn impose a linear restraint on the initial condition according to Eq. (7.3). For if S^{1}_{k} is the kth row of the matrix S^{1} then
The complete system of eigenvectors and latent roots is given for the systems of order 5, 6, 10 and 11. To each complex conjugate pair of latent roots there corresponds a conjugate pair of eigenvectors. These pairs appear in Appendix II as double columns separated by the symbol ±i. Partial solutions were also obtained for the system of order 20 but due to roundoff error accumulation the complete set was not obtained to satisfactory accuracy. The final check used on the eigenvectors ν_{k} was the identity [A–λ_{k}]ν_{k} = 0. For all eigenvectors show in Appendix V agreement to six or more places was obtained. B. Conclusions The results of the present study show the feasibility of obtaining general solutions of large order systems of linear differential equations. The rapid accumulation of roundoff error indicates the need for maintaining rather extreme accuracy in the calculations and suggests that the rate of error accumulation in a given process is of major importance. The method of solution of polynomial equations (proposed in section 5C) has proven satisfactory for polynomials of degree twenty and less. In spite of the rather large amount of calculation required by the rootsquaring process it is believed that in the absence of initial approximations to the roots this is perhaps the most satisfactory method available for equations of high degree. The matrix subroutines described in section 6B have proven very valuable in reducing problem preparation time in the present work. However, this does not necessarily imply that the subroutines possess practical value for general use. There are two major reasons for this situation:
As evidence of the practicability of the matrix subroutines for general use the present routines have been used with profit and a minimum of difficulty by another member of the Staff of the Computation Laboratory in the solution of a vibration problem.^{34 } Indeed, many of the difficulties in the use of subroutines arise in possible conflicts in the use of sequence or number storage and these in turn stem in part from the rather limited storage capacity of most computers. The advantage of large storage capacity and of separate sequence and number storage combine to make the efficient use of a library of subroutines a practical possibility for Mark IV. As a result of the experience gained in their use a number of changes in the subroutines are now proposed.
The machine computations were begun during the final testing of Mark IV when the machine was first coming into operation. Hence conditions were far from ideal and the lack of experience in programming and operation combined with machine failures to cause a considerable loss of time. However it is under just such adverse conditions that the relative ease of programming reruns on Mark IV proves especially valuable. The reliability of the computer has now been greatly improved and at present writing it has been operating with more than 85% good running time during the month of November 1953. The possibilities of the Leontief dynamic economic model have not been exhausted by the present study. For example, a stepbystep solution might well be attempted which could take into account the irreversibility of certain economic processes such as the investment in capital goods. It is believed that the economic interpretation of the present results will prove valuable in pointing the direction to be taken by further research. Notes Bibliography
