Machine Solutions of Linear Differential Equations
A thesis presented by
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.
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 so-called 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 round-off 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 round-off 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 over-simplified 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 input-output 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 input-output 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 x1,x2,...,xn are the outputs of the several sectors, Xij is the input from the ith sector to the jth sector, and zi 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. Xij = aijxj. The “flow” coefficients aij 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 zi, the required outputs xi 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. Solutions3 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 aij or the demand function zi (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 dxi/dt. For reasons discussed in reference 2 the following dynamic model is chosen[b]:
The flow coefficients aij and the capital coefficients bij are assumed to be constants; the outputs xi and the final demands zi are explicit functions of time.
The system of first-order 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 so-called “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. zi = gieμt where the gi 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.
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]
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 non-significant 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 “round-off”.
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 round-off errors. The maximum round-off 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 clear-out 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 time-consuming. 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 so-called 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 check-out time should therefore be taken into account in estimating the over-all 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 over-capacity numbers. In general, the decimal point should be located as far to the left as is consistent with this requirement since the round-off error is thereby reduced.
If very high accuracy is required or if the accumulation round-off 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 so-called “floating point” operation in which each number z is represented by two numbers y and p in the semi-logarithmic form z = y × 10p. The number y is “normalized”, i.e. it is stored with the first non-zero 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 1-50 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 “ij”, where ij 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 vice-versa). 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 IV5 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-1x, cos x, log10x, and 10x, may be coded directly; e.g. the operation cos x0 ⇒ c0 computes the cosine of the number in register x0 and stores it in register c0.
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” subroutine6 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.
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 vice-versa. 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 “line-number” in the range 0000-9999.
The function of certain β codes differs depending on whether a read-in or read-out 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 read-in and for read-out. 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.
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. Read-out 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 read-out.
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, b3 ≡ 0013 and d'8 ≡ 0138.
The function registers are ten in number (β0220-0229) and are designated as f0-9. In operation they are identical with the fast storage registers but are normally reserved for use by the permanent subroutines.6
Provision is made for the storage of four thousand numbers in “slow storage” positions numbered 0000-3999. 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 S0 ≡ 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 four-column β register capable of storing a number in the range 0000-3999. 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 read-in and read-out (without initiating a slow-fast transfer) under β code 3200. The S and S' registers have regular read-in and read-out and may be used as ordinary fast storage registers.
The IJ register is a three-column β register which can store numbers from 0 to 399. The first (low order) column is referred to separately as the J register, since read-in or read-out 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 ij. The IJ register controls or modifies the effects of certain β codes as follows: (for this example assume ij = 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 read-out and read-in under β code 0302. The shift register has two types of read-out and read-in 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 read-out only so that the original number stands in the shift register regardless of which read-in or read-out 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 0-16. 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 read-out only is possible (β3333).
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. Read-out only is possible (β3332).
The read-in switch is a manually-operated 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 four-column β 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 read-in (β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 twenty-four lines will take twice as long as a loop of twenty-three 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 sixteen-digit number may be introduced by four orders, the first with α code 04 followed by three with &alpha code 05.
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 read-ins are possible. Code 10 is a read-in which destroys the number previously stored in the accumulator. Read-out 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:
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 read-in of the MP and in the intervening seven lines any orders not involving the multiply unit may be performed.
Four different read-outs 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 32-digit 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.
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 read-in 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. Read-out under γ code 34 is not possible following a read-out under γ code 32.
The vacuous order 17 3355 77 performs no operation and is used to provide any desired delay, as for example before read-out 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.
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 a0 and a1 are required to agree to within a certain tolerance stored in a2. The following routine is used to stop the machine if the agreement is not as required:
Lines 1-4 place -|a0 - a1| 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 |a0 - a1|. Lines 5 and 7 together show that the accumulator is negative at one point and positive at another thus indicating it is operative.
All data are read to and from the machine by means of magnetic tape units which perform the following functions:
β codes 2301-11 select tape units 1-11 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 read-out and read-in under β codes 2321 and 2322.
It is also possible to select “tape unit k”, where k is the number (1-11) 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 read-in and read-out under β code 2320. Read to and from the record register without recording is possible under β code 2323.
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.
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 read-in and read-out.
Even where automatic waits are provided, machine time may be saved by interposing other orders between successive uses of any one of the above-mentioned 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.
Chapter 3 Linear Differential Equations
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 non-singular Eq. (3.2) may be reduced to
where D = B-1A. 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-1x. Then = S-1 and finally Eq. (3.3) becomes
Substitution of the expression
in the left-hand 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 c1,c2, ,cn.
The expanded form of (3.11) is exhibited below for the particular
In terms of components:
The general solution of Eq. (3.3) is now given by
Since S is non-singular 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 non-singular matrix S such that the similarity transformation S-1DS 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 non-trivial 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.
Thus xkj 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 Jnk(λ). As k and j range over the integers 1 to m and 1 to nk respectively, the vectors xkj range over the set of unit vectors e1,e2, ,en. Clearly then,
For example, if k = j = 2, and J is the matrix of section 3A11,
If λ = λk then
Thus xkj is a principal vector of order j associated with the latent root λk and the theorem is proved.
Proof: Let Jn be the canonical form of D. Then
Therefore Sxkj 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 non-singular, 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 Skj is the r-th column of S and r = j + nk.
Eigenvectors. A principal vector of order one is called an eigenvector. If in the canonical form of D, nk = 1, k = 1,2, ,m, the principal vectors are all eigenvectors and they form a linearly independent set. The nk 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 off-diagonal 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 tkeλ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 Jnk(λk) is a Jordan box then
Formally applying the above results to the function Jn-1 we have for each Jordan box
Premultiplication by Jnk gives
Therefore Jn-1 as formally define by (3.20) is indeed the inverse of Jn 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 Jn(λi,ni) – λI = Jn(λi–λ,ni) 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 non-trivial column of the adjoint Δ[Jn(λi–λ,ni)] is an eigenvector of Jn corresponding to the root λk (see Eqs. (3.14) and the definition of an eigenvector in section 3C1). The existence of such non-trivial columns will now be demonstrated.
The form of a single box of the adjoint of (J – λI) is shown below for nk = 3.
In the definition of Jn(λi,ni) 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 mk 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 mk = nk; and if a root λk is repeated in more than one Jordan box then mk > nk.
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 nk may be written concisely in terms of the matrices Ik as follows:
Suppose now that λk is a latent root
occurring in a single Jordan box so that
is some non-zero constant c0. 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 non-zero column occurring in the adjoint is the vector (-1)nk-1c0xk1.
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 non-zero scalar multiplier (-1)nk-1c0 in the column of the adjoint may be dropped to give xk1 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 cp are defined as
Differentiating (3.24) or (3.25) gives
Since c0 ≠ 0 the derived adjoint contains a non-zero column proportional to the eigenvector. If c1 = 0 it contains in addition a linearly independent column with a component proportional [to] xk2, which is the principal vector of order 2 associated with λk. If c1 ≠ 0 this column contains also a component in the direction of xk1. This brings out a point already implicit in (3.16); namely that a principal vector of order greater than 1 is non-unique 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 non-zero and they are principal vectors of lower order. Note also that for p ≥ nk 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 – λ)mk-nk
as shown by (3.24) or (3.25).
In this case mk > nk
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
Since S is non-singular 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 non-singular. The singular case will now be considered.
If the matrix B is of rank (n-r) 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 (n-r) equations. We are then left with a reduced system of equations in (n-r) variables of the form
where B1 is non-singular. 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 (n-r) 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 non-singular an alternative method is available. Premultiplication of (3.2) by A-1 gives:
Let J = S-1CS be the canonical form of C. Then Eq. (3.29) may be transformed to
is the general solution of (3.30). [Compare Eq. (3.11).]
If B is non-singular 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 = eJ-1tc 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 en (see Eq. (3.16)). The restraint is therefore simply cn = 0 and hence yn = 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 yn = 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 non-singular 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 = [aij] is the matrix of flow coefficients and B = [bij] 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)-1B. 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:
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)-1B. 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 round-off 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 theorem11 which follows.
Proof: Let x be an eigenvector of A corresponding to the latent root λ. Then Ax = λx and
Let xk = max |xj| ≠ 0 since x is non-trivial. Then
Dividing by |xk| we have
Thus λ lies in one of the circles Ci and the theorem is proved.
Since AT has the same latent roots as A a similar theorem holds for circles Di with centers aii 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 non-zero component in the direction of the largest eigenvector[c] the sequence of vectors Ax,A2x, ,Apx 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 Jp has the form[d]
where each submatrix is of the form
Since λ1 is dominant
where α is the first component of (S-1x) and e1 is a unit vector. Since the eigenvector of A corresponding to λ1 is Se1 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 ci and ri are respectively the eigenvectors of A and AT corresponding to the root λi and with length so chosen that riTci = 1, then
If i ≠ j then λi ≠ λj and hence riTcj = 0. Equation (4.2) is therefore established. If will now be shown that
Note that each product ciriT 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 cj, 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 AT may be obtained by re-applying the power method, or more simply (since λ1 is now known), by solving the singular homogeneous system
to obtain a non-trivial 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 ri = ci and a separate calculation of r1 is not required.
(3) Multiple Roots
If multiple roots occur two cases must be distinguished:
If the matrix A is non-defective 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,A2x, ,Apx exists. For example, if λ1 is dominant and n1 = 3 then
Thus the sequence Ax,A2x, ,Apx 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
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
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 sk, k = 1,2, ,n as follows:
The quantities sk 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 s1. Furthermore, Eqs. (3.18) and (3.19) show that the sum of the diagonal elements of Jk is equal to sk.
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 c0,c1, ,cn be the coefficients of the characteristic polynomial of the matrix A and let Fj(A) be a sequence of polynomial functions of A defined by
and in general
by virtue of (4.5).
Using this expression for ck, 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,A2,A3, ,An and from the traces s1,s2, ,sn 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 Fn = 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 non-defective and consider the series of matrices Fk(A) defined by (4.8). Because of round-off error the computed value k(A) will differ from Fk(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 Ek 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 round-off committed at the first step. Since a similar round-off 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 non-defective matrix and let S be the matrix which diagonalizes A. If S0 is an approximation to S then the method of Jahn18 may be used to improve this approximation. The object is to compute S where AS = SΛ and Λ is a diagonal matrix. Since S0 ≐ 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
Then S1 = S0 (I + C) is an improved
approximation to S correct to second order.
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 n4 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 re-appear 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 round-off 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.
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 Birge-Vieta19 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 Division16
Let Pn(x) = a0xn + a1xn-1 + + an-1x + an be a polynomial of degree n with real coefficients. Let
where the coefficients bi are defined by
Hence bn is the value of the polynomial evaluated at p and when bn is zero, Q(x) is the so-called 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 bn. 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
as may be verified by expanding the left-hand side.
2. Change of Origin
A translation or change of origin to the point
x = p is defined by the substitution
The coefficients bk 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 Qk(x) be defined by
Repeating the differentiation of the relation
yields the values of the successive derivatives of Pn(x) evaluated at the point p. Thus
3. The Birge-Vieta Process
The so-called Birge-Vieta 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 re-established by Birge and came into wide use. The method is perhaps the most satisfactory known for obtaining real roots with a desk calculator.
The Birge-Vieta method is based on the well-known Newton-Raphson process (reference 20, page 84) which is defined as follows: Let be a real roots of the equation
Then for any value x0 sufficiently close to the quantities xi 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 Birge-Vieta process. This process is usually employed only for real roots, but may be extended to complex roots by simply allowing the variables x0,x1, to be complex. Since a polynomial is an analytic function21 the derivative Pn'(x) exists for x complex, and the process carries over without change to complex numbers. (Note, however, that if the initial approximation x0 is real, the successive xi remain real.) Every operation is now complex and hence the amount of computation is increased by a factor of about four.
The complex Birge-Vieta 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 Birge-Vieta 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 bn-1 and bn zero. Treating bn-1 and bn as functions of p and q and expanding in Taylor’s series we have
A generalization of the Newton-Raphson 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 Bairstow22 and as presented in standard works23, 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 cn. Comparing with (5.4)
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 cn-1 and cn rather than bn-1 and bn. Note that the vanishing of cn-1 and cn is the necessary and sufficient condition for the vanishing of bn-1 and bn.
Treating x, p and q as independent variables and differentiating (5.11) partially with respect to p we have
Thus ∂cn-1/∂p and ∂cn/∂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 ∂cn-1/∂q and ∂cn/∂q respectively. Hence
Applying Eq. (5.9) to the quantities cn and cn-1 and substituting the expressions obtained above for the derivatives we have
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 Birge-Vieta 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 Birge-Vieta method requires one complex[b] multiplication per step. Because the coefficients computed in the complex Birge-Vieta 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 Newton-Raphson process for obtaining
the solution of the
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 Birge-Vieta 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 Birge-Vieta 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 "(x0) 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
This pair of non-linear 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 right-hand 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
Aitken25 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
Again the modified algorithm
yields the same quotient polynomial Q(x). Adopting the notation
Note that the vanishing of the quantities A, B, C, and D is the necessary and sufficient condition for the vanishing of the remainders bn, bn-1, bn-2 and bn-3.
Differentiating Eq. (5.17) we have
Re-arranging the last equation to the form
shows that the application of the modified algorithm to Q(x) leaves Ds as the first remainder. Thus if the scheme of two successive applications of the algorithm is shown as
we have Ds = K. The coefficient of x2 in (5.21) is
But from the form of the remainders in (5.17) it is clear
that the coefficient of x2 is
Similarly Bs = I and As = H. Rearranging (5.18) to the form
shows the left-hand side to be the polynomial formed of the coefficients c0,c1, ,cn-4, 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 Dp = H, Cp = G, Bp = F, Ap = 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
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
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 y3. 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 Newton-Raphson 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 C2, C2/2 and 4D we have
A few rough mental calculations show that ø(50) > 0 and ø(52) < 0.
Hence b - β = √416.2704 = 20.4027,
The quartic formed from these roots agrees with the given quartic to four decimal places in each coefficient.
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 x2 + αx + β and x4 + ax3 + bx2 + 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 (-F1/3,F1/3). If all the zeros of P6(x) are complex the zeros of (5.31) are all positive and exactly three in number. Hence either the interval (0,F1/3) or the interval (F1/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 β = F1/3 and
(since α is finite)
the expression for α being obtained as the limiting form of (5.30).
Since 0 ≤ β ≤ F1/3 ≅ 22 choose β0 = 10. Linear prediction is used even if two successive values of ø do not differ in sign. R2 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 exponential26 namely
for values of n ≤ 23. Approximations were chosen on a rectangular mesh and each was introduced as an initial approximation to (a) the Birge-Vieta 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 Method27
Let r1,r2, ,rn 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 2k power. Hence if all zeros are distinct a few iterations serve to separate them widely. If b0,b1, ,bn are coefficients of the polynomials so generated, we have
Hence |r1|,|r2|, ,|rn| may be obtained. More generally if rp+1,rp+2, ,rp+s are a set of s roots with equal moduli then the polynomial
has the s zeros, rp2k, ,rp+s2k. Therefore
Thus the moduli of all the roots are obtainable.
Certain refinements of Graeffe’s method allow the zeros r1,r1, ,rn 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 root-squaring 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 = -x2 then gives
Equation (5.35) defines the Graeffe root-squaring 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 Birge-Vieta 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 Lin28 which are included under the process of penultimate remaindering converge only under certain conditions as established by Aitken25. 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 detail29 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 n3 (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 n2/2
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 over-capacity 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.
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
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 = [x1,x2, ,xn] will be said to be stored in location b if the components x1,x2, ,xn are stored respectively in the registers b+1,b+2, ,b+n. The negative sum of the components will be designated by x0, 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 x0,x1,x2, ,xn must be zero. The use of the negative sum makes the process of summing the check more uniform.
Each row of a matrix A = [aij] 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:
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 sum-check is for checking the components of a vector each time it enters into the computation. As an example of other ways in which the sum-check 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 round-off 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 y0 so that
In this way the only error to be allowed at each step is the maximum round-off 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 a00) 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 190-199 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 check-stops, 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 round-off 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 F0,F1, ,Fn, 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 Fj by fij, the ith column νi(λ) of the adjoint is given by
If λk is a latent root of the matrix A, any non-trivial 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 fij[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 = rkμk where rk is real and |μk| = 1. Equation (6.11) may now be written as
The powers of rk are formed using
floating point operation and the vectors
Chapter 7 Results and Conclusions
The tables of initial data and final results appear in Appendices I-VII 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 x1,x2, xn 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 1-21 indicate which of the variables x1-x21 are assigned to represent each industry; e.g. x5 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 2-26 contain square matrices of various orders. Except in Tables 8-13 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 = [aij] 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 = [bij], supplied by the Harvard Economic Research Project, appears in Appendix II. The broken lines in each of these tables (2-7) 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 2-7) and all results are rounded to four decimal places. Twelve decimal places were carried in the calculation of Tables 8-19. A final check on the results shown in Tables 20-24 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 eJ-1tc 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 ck 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-1k 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 round-off 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.
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 round-off 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 root-squaring 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 step-by-step 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.