Machine Solutions of Linear Differential Equations
Applications to a Dynamic Economic Model

 

A thesis presented by

Kenneth E. Iverson

to

The Division of Applied Science
in partial fulfullment of the requirements
for the degree of
Doctor of Philosophy
in the subject of Applied Mathematics
 

Harvard University
Cambridge, Massachusetts, January, 1954
 

Preface

The present use of the Mark IV Calculator in the solution of a new dynamic economic model continues the previous work of Dr. H.F. Mitchell in this field. The numerical results obtained should be of interest to economists, and it is believed that a number of the methods employed will prove useful to computers.

The writer wishes to express his appreciation to Professor Howard H. Aiken for his helpful guidance and encouragement, and to Professor Wassily W. Leontief who proposed the problem and offered many valuable suggests for its solution. The advice and assistance of numerous members of the Staff of the Computation Laboratory is gratefully acknowledged, in particular that of Dr. Robert C. Minnick and Mr. Anthony G. Oettinger. To Mr. James R. Terrell the writer is indebted for profitable discussions arising out of his use of the matrix subroutines described herein.

The writer also wishes to express his thanks to Miss Jacquelin Sanborn who typed the plates for printing and to Miss Carmela Ciampa, Miss Joan Boyle, Mr. Robert Burns and Mr. Paul Donaldson, all of whom assisted in the preparation of the manuscript.
 
 

Synopsis

The present study was undertaken to obtain solutions to a dynamic economic model proposed by Leontief, and to develop methods of computation suitable to the Mark IV Computer.

The economic model as presented in Chapter 1 is described by a system of linear differential equations with constant coefficients. The data for six systems of order 5, 6, 10, 11, 20 and 21 were supplied by the Bureau of Labor Statistics and the Harvard Economic Research Project.

The general solution for each system was sought. In Chapter 3 the complementary function is shown to be expressible in terms of the latent roots and principal vectors of a certain matrix related to the coefficient matrices of the given system. An expression is also given for the particular integrals appropriate to a set of demand functions of interest in economic theory.

The 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 x₁,x₂,...,x_n are the outputs of the several sectors, X_ij is the input from the ith sector to the jth sector, and z_i is the “final demand” or amount of goods to be made available for use outside the system considered, then

xi – Xij = zi .

The inputs to each sector are assumed^[a] to be linear functions of the output of that sector, i.e. X_ij = a_ijx_j. The “flow” coefficients a_ij are assumed to be constants so that the equilibrium condition may now be stated as a set of n linear equations in n unknowns,

xi – aijxj = zi ,     i = 1,2,…,n.                             (1.1)

For a given final demand z_i, the required outputs x_i of each sector may now be obtained as a solution of (1.1).

The division of the economy into sectors is somewhat arbitrary. In general, the greater the number of sectors considered, the more detailed is the information available in the solution. Since, however, the amount of calculation required varies approximately as the cube of the number of sectors considered, a compromise must be made. Solutions³  have been obtained for systems of various orders, the largest being a system of 192 sectors.

Equations (1.1) describe a static model since no dependence on time is implied. Dependence on time may be introduced by making either the coefficient a_ij or the demand function z_i (or both) functions of time. The simple linear model may also be extended in other ways; e.g. by including terms of higher degree in the variables or by adding terms in the derivatives dx_i/dt. For reasons discussed in reference 2 the following dynamic model is chosen^[b]:

xi – aijxj – bij

dxj dt

= zi ;     i = 1,2,…,n.                   (1.2)

The flow coefficients a_ij and the capital coefficients b_ij are assumed to be constants; the outputs x_i and the final demands z_i are explicit functions of time.

The system of 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. z_i = g_ie^μt where the g_i are constants. Solutions for demand functions of this kind are discussed in section 3F.

Flow coefficients for a system of 192 sectors were provided by the Bureau of Labor Statistics and formed the basis for the present calculations. Capital coefficients were supplied by the Harvard Economic Research Project. By a processor of “aggregation” described in Chapter IX of reference 2, the coefficients for several systems of different orders were computed from the original data. The orders of the systems chosen were 6, 11 and 21. Comparison of the solutions for the various systems may be expected to give some indication of the distortion introduced by aggregation. Three “open” systems of order 5, 10 and 20 were derived from the systems of order 6, 11 and 21 by suppressing the “households” section.

General solutions for the six systems were obtained using the Harvard Mark IV Calculator. Particular integrals (see Section 3F) were also computed for exponential functions e^μt for

At the time this study was undertaken, the Mark IV Calculator had just been completed. The need for solving the extensive economic systems mentioned above therefore provided an excellent opportunity for an investigation of general matrix methods especially suitable for application to Mark IV.

Notes

Chapter 2   Problem Preparation for an Automatic Computer

Section A of this chapter will be devoted to a consideration of the general aspects of problem preparation for an automatic computer as a preliminary to the discussion in Chapters 4 and 5 where possible methods for the solution of our economic model are analyzed and compared on the basis of their suitability for machine computation. Section B will be devoted to programming details of the Harvard Mark IV Calculator as a preliminary to the discussion of Chapters 6 where the programming of the selected methods for the Mark IV Calculator is discussed. Reading of Section B may well be deferred until after Chapter 5.

A. General Aspects of Problem Preparation^[a]

(1) Introduction

The solution of virtually any problem of applied mathematics may be reduced to a sequence of elementary operations involving only addition, subtraction, multiplication, division, and reference to previously computed quantities. For this reason it is clear that a truly general purpose calculator may be based on the following components:

The calculator may be made automatic by the addition of a “sequence unit” which determines the sequence of operations to be carried out. The sequence unit operates by reading sequentially a set of instructions presented to it, and carrying out the prescribed operations by determining which numbers are to be transferred, read, or recorded at each step. The sequence instructions or “program” must be coded and recorded in a form suitable to the particular machine.

Since each step must be programmed it is clear that little or no labor is saved by the automatic computer unless the operations are made repetitive. Such repetition is made possible by a special instruction referred to as a “call” which interrupts the sequence of instructions being read and cause the sequence unit to begin reading instructions from some point in the program specified by the call. This provides the possibility of a program consisting of a single closed loop which could, for example, be used in the tabulation of a given function as follows: The initial value of the argument is stored in a certain “argument” register and the prescribed function of the number in this register is computed and recorded. The number in the argument register is then incremented and a call is made to the point in the program immediately following the introduction of the initial argument.

The computer is made much more versatile by providing another instruction termed “conditional call” which effects a call if and only if the number in a special “conditional” register is negative. This allows a choice of programs to be made depending on results obtained in the course of the calculation. For example, a series may be automatically truncated at the first 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 |T| - R be read into the condition register. If T is the last term added and R is a certain tolerance, the conditional call will fail when |T| becomes less than the tolerance and the computer will proceed to the next section of the program.

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 × 10^p. 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 “i_j”, where i_j is the number stored in a certain control register called the IJ register.^[c] Thus, without changing the actual order stored in the sequence storage, the operation effected by it may be modified by changing the number resident in the IJ register.

The separation of number storage and sequence storage offers important advantages in making reruns. To make a rerun from any given point in a program all the relevant storage must be restored to the conditions which obtained at that point. For the number storage this may usually be done by so planning the program that the numbers entered into the calculation are still available in the same registers at the time the check is made. Although the number of orders modified in a certain routine may be rather large and the pattern of the modifications rather complex (due perhaps to repetitive loops within the routine), yet the modifications may usually be expressed as simple functions of a very limited number of parameters or “control numbers”. In Mark IV the desired functions of these control numbers are computed and introduced into the IJ register to produce the required modifications. Since the sequence instructions are themselves unchanged, the program may be restored to any point by restoring the new control numbers to their correct values. The rerun procedure is therefore relatively simple and may easily be made to occur automatically.

On the other hand, the restoration of the sequence instructions in a machine which modifies them directly is in general so difficult that the simplest procedure is to read the original instructions again from the input device. Automatic reruns become virtually impossible. It should also be noted that easy reruns greatly simplify the process of checking out a program.

A disadvantage of the separation of number and sequence storage shows up in the case of problems requiring a large amount of number storage and little sequence storage (or 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 IV⁵  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, log₁₀x, and 10^x, may be coded directly; e.g. the operation cos x₀ ⇒ c₀ computes the cosine of the number in register x₀ and stores it in register c₀.

In the simpler operations described above (namely, those requiring less than twelve instructions for their completion) the coding machine records the required sequence instructions directly. In the more complicated operations, e.g. cos x, the coding machine provides only the instructions which are necessary to call a “permanent” subroutine⁶  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.⁷ The value of a particular subroutine depends on both the ease and the frequency of its use. On this basis the subroutines for the elementary functions are in general the most valuable and subroutines for linear operations probably take second place. A set of subroutines for matrix operations on Mark IV are discussed in Chapter 6.

B. Mark IV Programming^[d]

This section will be devoted to those details of Mark IV programming^[e] which are necessary to an understanding of the routines discussed in Chapter 6. Programs for Mark IV may be prepared either on the algebraic coding machine mentioned in section A or on a “decimal” coding machine, with which the elementary sequence codes are recorded. Since the use of the decimal coding machine leads to a slightly more efficient program, it was chosen for the preparation of the subroutines of Chapter 6. For this reason only the decimal coding will be considered in this section.

Orders

The units of Mark IV are divided into two groups called β and γ. Under control of a single order a number from any β unit may be transferred to any γ unit or 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.

Registers

Most of the registers of Mark IV have a capacity of sixteen decimal digits and sign. Such registers are referred to as “regular” in contrast with other “short” registers which do not have a full sixteen columns. All short registers are so arranged that their first (low order) column corresponds to the first column of the regular registers. 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, b₃ ≡ 0013 and d'₈ ≡ 0138.

Function Registers

The function registers are ten in number (β0220-0229) and are designated as f_0-9. In operation they are identical with the fast storage registers but are normally reserved for use by the permanent subroutines.⁶ 

Slow Storage

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 S₀ ≡ 0200 to S'₉ ≡ 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.

IJ Register

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 i_j. The IJ register controls or modifies the effects of certain β codes as follows: (for this example assume i_j = 124)

Shift and Normalize

The shift register may be used to shift a number either to the right or to the left by any number of places determined by the number in the shift control register (SCR). The SCR will store any number in the range ±19, a positive number producing a shift to the left, and a negative number producing a shift to the right. The SCR has ordinary 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).

Sign Register

The sign register is manually set to store a positive or negative zero and is used to determine whether the decimal point is in the high or the low order part of a number in double accuracy operation. It is also a convenient source of zeros. Read-out only is possible (β3332).

Read-in Switch

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

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:

Multiplication

The multiplicand (MC) must be read in first under γ code 20 and remains in the MC register. The multiplier (MP) is read in on the next (or any later) line under γ code 21. The product may be read out on the eighth (or any later) line following the 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.

Division

The divisor must be read in first, followed by the dividend within fifteen cycles. The quotient may be read out on the twentieth line following 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.

Vacuous Codes

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.

Coded Stops

Two stop codes are provided: β1312 stops the machine on a positive transfer and β1313 on a negative transfer. The two stops are provided to give more reliable operation on a coded check.

Example. Suppose that two computed results stored in a₀ and a₁ are required to agree to within a certain tolerance stored in a₂. The following routine is used to stop the machine if the agreement is not as required:

Lines 1-4 place -|a₀ - a₁| in accumulator 1. Line 6 adds the tolerance, and so makes the accumulator 1 positive if the agreement is as required. Line 7 then stops the machine if the accumulator is negative. The reason for including line 5 is that the accumulator may be defective and standing in a permanently positive condition, thus passing the check regardless of the value of |a₀ - a₁|. Lines 5 and 7 together show that the accumulator is negative at one point and positive at another thus indicating it is operative.

Input-Output

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.

Printing

If final output is to be in printed form, the recorded tapes must be used to control a printer. To accomplish this a “serial” number must be recorded before each number to be printed. This serial number serves to identify the number and to control the printer. The various columns (number right to left) of the serial number are used as follows:

Before printing the printer compares the first column of the following serial number with the preceding one. If the digit has increased by one printing proceeds, but if not the printer passes on to the next number without printing. Thus if any number has be rerecorded (due perhaps to a rerun) only the last recorded value is printed. All numbers are recorded in duplicate on the tape; the printer compares the two recorded values and stops if they do not agree.

Waits

Certain units which require several cycles to operate are provided with automatic waits to delay any order involving them until the previous operation is completed. This is not true of the multiplier and divider, where the requisite number of orders must be interposed between 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.

Notes

Chapter 3   Linear Differential Equations

A. Introduction

The general solution of a system of ordinary linear differential equations is well known, as are certain computational methods for evaluating the parameters involved. In this chapter the theory is developed in a manner designed to provide the necessary background for the discussion of computational methods in Chapter 4. In order to make the structure of the solution clear for the case of multiple roots, the whole discussion is based on the Jordan canonical form of the matrix of coefficients. The complementary function is expressed in terms of the system of principal vectors, and expressions for the principal vectors are derived in terms of the adjoint of the characteristic matrix. A knowledge of elementary matrix theory such as contained in the first three chapters of reference 8 will be assumed.

In matrix notation the most general form of a system of first order^[a]  linear differential equations in n variables is

Ax – B = z                                             (3.1)

where A and B are square matrices of order n and

	x = (x1,x2,…,xn),
	= (	dx1 dt	,	dx2 dt	,…,	dxn dt	),
and
	z = (z1,z2,…,zn)

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

Ax – B = 0.                                             (3.2)

If B is non-singular Eq. (3.2) may be reduced to

Dx – = 0                                               (3.3)

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:

1.	The letter n refers exclusively to the order of the matrix under consideration. The order of a submatrix will be denoted by the letter n with a subscript.
2.	Matrices are denoted by capital letters A, B, or A(λ), B(λ), etc., where the elements are functions of a parameter λ. A subscript n may be added to specify the order of the matrix.
3.	Vectors are denoted by lower case Roman letters x, y, etc.
4.	All equations are to be understood as matrix or vector relations unless otherwise specified.
5.	The Kronecker delta is defined as:   δij  = 0 if i ≠ j,   = 1 if i = j.
6.	The symbols e1,e2,…en will be used to denote the so-called unit vectors defined by   ei = (ei1,ei2,…ein) where eij = δij.
7.	The unit matrix is denoted by I and the matrices Ik are defined by   Ik = [δi,j-k],     k ≥ 0. Thus Ik is a matrix whose elements are zero except for a diagonal of unit elements located k places above the main diagonal. For example if n = 3   I0 = 1 0 0 ;       I1 = 0 1 0 ;       I2 = 0 0 1 .   0 1 0 0 0 1 0 0 0   0 0 1 0 0 0 0 0 0 If   X = [Xpq] = IkIm, then   Xpq = δp,j-kδj,q-m = δp,q-(m+k) . Therefore   X = IkIm = Ik+m. Similarly,   ImIk = Ik+m. Thus,   I0 = I .                                 (3.4)   Ikm = Imk = Ik+m   Im = 0 if m ≥ n
8.	Since differentiation is a linear operation it is easily extended to matrices and vectors. Thus if   A(λ) = [aij(λ)] is a matrix whose elements are functions of a scalar λ, the derivative dA/dλ is defined as the matrix   dA/dλ = [daij/dλ].
9.	Two matrices A and B are said to be “similar” if there exists a non-singular matrix S such that   B = S-1AS .
10.	A matrix Jn(λ) of the form   Jn(λ) = λI + I1                                         (3.5) will be called a Jordan box of order n, e.g.,   J3(λ) = λ 1 0 .   0 λ 1   0 0 λ
11.	A matrix J(λi,ni) made up of Jordan boxes arranged along the main diagonal with zeros elsewhere is said to be in Jordan canonical form and will be referred to briefly as a canonical matrix. Thus:   Jn(λi,ni) = Jn1(λ1)                           (3.6)     Jn2(λ2)       …         Jnm(λm) where   ni = n. For example, if n1 = 3, n2 = 2, n3 = 1   J6(λi,3;λ2,2;λ3,1) = λ1 1  0  0  0  0  .   0  λ1 1 0  0  0    0  0 λ1 0  0  0    0  0  0  λ2 1 0    0  0  0  0 λ2 0    0  0  0  0  0  λ3
12.	The exponential function of a matrix is defined formally by the following series, convergent for any infinite matrix A:   eA = I + A +  A2 2!  +  A3 3!  + … . This function may be shown to have the usual properties of the exponential function (reference 9, p. 41). In particular   deAt dt = AeAt .                                       (3.7) If A and B commute the following relation also holds:   eAeB = eBeA = eA+B . Consequently   eAe-A = I .                                           (3.8) For example:   eJ4(λ)t = e(λI+I1)t = eλteI1t     = eλt (I0 + tI1 + t2 2! I2 +  t3 3! I3).  Therefore   eJ4(λ)t = eλt 1  t  t2 2! t3 3! .     1  t  t2 2!       1  t          1

The proof of the following theorem, omitted here, may be found on page 64 of reference 10.

B. The Complementary Function

Equation (3.3) is repeated here for reference

Dx – = 0                                               (3.3)

Substitution of Eq. (3.9) gives

SJS-1 – = 0.

Using the fact that SIS^-1 = I and then premultiplying by S^-1 we have

JS-1x – S-1 = 0.

Let y = S^-1x. Then = S^-1 and finally Eq. (3.3) becomes

J(λi,ni)y – = 0.                                          (3.10)

Substitution of the expression

in the left-hand side of Eq. (3.10) gives

[ J(λi,ni)eJ(λi,ni)tc – eJ(λi,ni)t ] c.

The bracketed quantity is identically zero by virtue of (3.7) and therefore (3.11) is a solution of (3.10) for an arbitrary vector c. This is the general solution sought since it contains n arbitrary constants c₁,c₂,…,c_n.

The expanded form of (3.11) is exhibited below for the particular matrix^[b] J = J₆(λ₁,3;λ₂,2;λ₃,1) appearing in section 3A11:

y =

eλ1t

1

t

t2/2!

c1

1

t

c2

1

c3

eλ2t

1

t

c4

1

c5

eλ3t

c6

In terms of components:

	y1 = (c1 + c2t + c3	t2 2!	)eλ1t
	y2 = (c2 + c3t) eλ1t
	y3 = c3eλ1t
	y4 = (c4 + c5t) eλ2t
	y5 = c5eλ2t
	y6 = c6eλ3t .

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

t = 0,   x = .

= Sc.

c = S-1.                                               (3.13)

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

	(A – λI)q-1x	≠ 0,		(3.14)
	(A – λI)qx	= 0.

The corresponding value of λ is called the latent root, or more briefly the root, associated with the principal vector x.

Proof: In order to relate the unit vectors to the structure of J_n, the following notation is adopted:

	xkj = er,
where
	r = j + np ;
	k = 1,2,…,m;
	j = 1,2,…,nk.

Thus x_kj is a column vector with zero components except for a unit element in the jth position of the group of rows occupied by the kth Jordan box J_nk(λ). As k and j range over the integers 1 to m and 1 to n_k respectively, the vectors x_kj range over the set of unit vectors e₁,e₂,…,e_n. Clearly then,

For example, if k = j = 2, and J is the matrix of section 3A11,

	[Jn(λi,ni) – λI] x22 =		λ1–λ	1							0		=		0		.
				λ1–λ	1						0				0
					λ1–λ						0				0
						λ2–λ	1				0				1
							λ2–λ				1				λ2–λ
								λ3–λ			0				0

If λ = λ_k then

	(Jn – λkI)xkj	= xk,j-1,	j ≠ 1;		(3.15)
		= 0,	j = 1.

Clearly then

	(Jn – λkI)j-1xkj	= xk,1 ≠ 0		(3.16)
and
	(Jn – λkI)jxkj	= 0.

Thus x_kj is a principal vector of order j associated with the latent root λ_k and the theorem is proved.

Proof: Let J_n be the canonical form of D. Then

Therefore

	S-1(D – λkI)j-1Sxkj = (J – λkI)j-1xkj ≠ 0		(3.17)
and
	S-1(D – λkI)jSxkj = (J – λkI)j-1xkj = 0.

Therefore Sx_kj is a principal vector of order j of the matrix D and is associated with the root λ_k.

The principal vectors of A are therefore the columns of S and since S is 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 S_kj is the r-th column of S and r = j + n_k.

Eigenvectors. A principal vector of order one is called an eigenvector. If in the canonical form of D, n_k = 1, k = 1,2,…,m, the principal vectors are all eigenvectors and they form a linearly independent set. The n_k are necessarily equal to unity if

The proof of (2) rests on the fact that if D is symmetric so is J. Therefore J has no 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 t^ke^λt.

(2) The Adjoint Matrix

It is easily verified that the pth power of a canonical matrix J may be expressed as

	Jnp(λi,ni) =		Jn1p(λ1)					.
				Jn2p(λ2)
					…
						Jnmp(λm)

Hence for any polynomial function of J

	F(Jn) =		F[Jn1(λ1)]					.                 (3.18)
				F[Jn2(λ2)]
					…
						F[Jnm(λm)]

Further, if J_nk(λ_k) is a Jordan box then

F[Jnk(λk)] = F(λkI0 + I1) = F(λk)I0 + F'(λk)I1 + (λk)I2 + … .      (3.19)

For example,

	J34 =		λk4	4λk3	6λk2		.
				λk4	4λk3
					λk4

Formally applying the above results to the function J_n^-1 we have for each Jordan box

	F[Jnk(λk)]	= Jnk-1(λk) = λk-1I0 – λk-2I1 + λk-3I2 – λk-4I3 + …
		= (-1)pλk-(p+1)Ip.                             (3.20)

Premultiplication by J_nk gives

	Jnk(λk) Jnk-1(λk)	= (λkI0 + I1) (-1)pλk-(p+1)Ip
		= (-1)pλk-pIp + (-1)pλk-(p+1)Ip+1
		= (-1)pλk-pIp – (-1)pλk-pIp
		= I0.

Therefore J_n^-1 as formally define by (3.20) is indeed the inverse of J_n and the preceding results may therefore be extended to negative powers of J.

To avoid difficulties due to the possible singularity of the matrix J we consider instead the adjoint of J, to be denoted by Δ(J) and defined by

Since the determinant |J| is equal to

(λk)nk

and since the largest negative power of λ occurring in J^-1 is λ_k^-n_k it follows that the adjoint is finite even for singular matrices. Since J and J^-1 commute

(3) Principal Vectors Derived from the Adjoint

Since J_n(λ_i,n_i) – λI = J_n(λ_i–λ,n_i) is a canonical matrix, substitution in (3.22) gives

If λ is chosen equal to any latent root λ_k the determinant on the right vanishes and therefore any non-trivial column of the adjoint Δ[J_n(λ_i–λ,n_i)] is an eigenvector of J_n corresponding to the root λ_k (see Eqs. (3.14) and the definition of an eigenvector in section 3C1). The existence of such non-trivial columns will now be demonstrated.

The form of a single box of the adjoint of (J – λI) is shown below for n_k = 3.

	Δ[Jnk(λk – λ)] = (λj – λ)nj		(λk – λ)-1	-(λk – λ)-2	(λk – λ)-3		(3.23)
				(λk – λ)-1	-(λk – λ)-2
					(λk – λ)-1

In the definition of J_n(λ_i,n_i) in section 3A11 the symbol λ_k is used to denote the diagonal element of the kth Jordan box. This notation is not meant to imply that the elements λ₁,λ₂,…,λ_m are necessarily distinct. To this point the occurrence of two or more equal values of λ_k has not entered the discussion but a distinction must now be made.

Let m_k be the multiplicity of the root λ_k in J. Thus if T is the set of latent roots equal to λ_k the multiplicity of λ_k is

mk = nk.

In particular, if the λ_k are distinct m_k = n_k; and if a root λ_k is repeated in more than one Jordan box then m_k > n_k.

The symbol will now be adopted for the product of the λ’s. The prime indicates that the factors are chosen only from a set of distinct λ_k. Then

λknk = λkmk.

The product in (3.23) may now be written as

(λj – λ)mj.

Factoring out (λ_j – λ)^m_k and multiplying it into the matrix, Eq. (2.23) appears as follows:

	Δ[Jnk(λk – λ)] = (λj – λ)mj		(λk – λ)mk-1	-(λk – λ)mk-2	(λk – λ)mk-3
				(λk – λ)mk-1	-(λk – λ)mk-2
					(λk – λ)mk-1
					(3.24)

The same relation for a general value of n_k may be written concisely in terms of the matrices I_k as follows:

Δ[(Jnk(λk – λ)] = (λj – λ)mj (-1)p(λk – λ)mk-(p+1) Ip .       (3.25)

Suppose now that λ_k is a latent root occurring in a single Jordan box so that m_k = n_k. If the adjoint is evaluated at λ = λ_k all Jordan boxes other than the kth are null due to the occurrence of the factors (λ_k – λ_k)^m_k in the product . For the kth box the product

(λk – λj)mj

is some non-zero constant c₀. 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)^n_k-1c₀x_k1.

Eigenvectors (and principal vectors) are defined by a homogeneous equation and therefore any scalar multiple of an eigenvector is also an eigenvector. For convenience, the non-zero scalar multiplier (-1)^n_k-1c₀ in the column of the adjoint may be dropped to give x_k1 as the eigenvector corresponding to λ_k. This is the same result previously obtained from Eqs. (3.16).

The remaining principal vectors associated with λ_k may be obtained by differentiating the adjoint with respect to λ. For convenience the scalar constants c_p are defined as

cp = (λk – λ)mj .

Differentiating (3.24) or (3.25) gives

Δ[Jnk(λk – λ)] = c1

(-1)nk-1

– c0

(-1)nk-2

0

(-1)nk-2

or

Δ[(Jnk(λk – λ)] = c1(-I1)nk-1 – c0(-I1)nk-2 .

Since c₀ ≠ 0 the derived adjoint contains a non-zero column proportional to the eigenvector. If c₁ = 0 it contains in addition a linearly independent column with a component proportional [to] x_k2, which is the principal vector of order 2 associated with λ_k. If c₁ ≠ 0 this column contains also a component in the direction of x_k1. This brings out a point already implicit in (3.16); namely that a principal vector of order greater than 1 is 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 ≥ n_k the derived adjoint is null.

The case where λ_k is repeated in two or more boxes may now be disposed of briefly. The element of lowest degree in (λ_k – λ) which occurs in the adjoint is (λ_k – λ)^m_k-n_k as shown by (3.24) or (3.25). In this case m_k > n_k so that the adjoint evaluated at λ = λ_k is null, as are all derived adjoints up to but not including the (m_k–n_k)th. The (m_k–n_k)th derived adjoint has a single non-zero column proportional to the eigenvector associated with the kth box. The successive derived adjoints up to the (m_k–1)st repeat the pattern established above and yield in all the n_k principal vectors associated with the kth box. Note, however, that this same process is going on simultaneously for all boxes with diagonal elements equal to λ_k. Thus more than one new independent column may be introduced by each differentiation. In particular, the (m_k–1)st derived adjoint introduces all the principal vectors of highest order associated with each occurrence of λ_k.

(4) Principal Vectors of a General Matrix

The results obtained above for a canonical matrix may now be extended to a general matrix A by the use of the following train of relations: If

and

Δ(A – λI) = S [ Δ(J – λI)] S-1,     p = 0,1,2,…

Then

(A – λI)q Δ(A – λI) = S [(J – λI)q Δ(J – λI)] S-1.            (3.27)

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

(A – λI)q Δ(A – λI)   and   (J – λI)q Δ(J – λI)

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

A1x – B1 = 0                                           (3.28)

where B₁ 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:

Ix – C = 0                                             (3.29)

where

Let J = S^-1CS be the canonical form of C. Then Eq. (3.29) may be transformed to

y – J = 0                                              (3.30)

where

Then

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 = e^J-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 e_n (see Eq. (3.16)). The restraint is therefore simply c_n = 0 and hence y_n = 0. This is the same linear restriction arrived at by the previous method as may be verified as follows: The transformation

	S-1A-1[Ax – B]
	x = Sy

is one of the possible ways of reducing Eq. (3.2) to the form (3.28), for we obtain (3.30) namely,

y – J = 0.

Since the last row of the matrix J is null the system may now be reduced to the form (3.28) by dropping the last equation and using it as an auxiliary restraint. This equation is simply y_n = 0 as was to be demonstrated. The same reasoning may be used to extend this result to the case of multiple zero roots of B. We may therefore lay down the following procedure for the case when B is singular.

This result is very useful in the present problem since the matrix A is known to be 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.

Ax – B = z.                                             (3.1)

Again assuming B non-singular,
	B-1Ax – = B-1z.
Let
	B-1A = SJS-1.
Then
	S-1JSx – S-1IS = B-1z.
Letting y = Sx yields
	Jy – = SB-1z.                                           (3.32)

A particular integral is then given by

y = -eJt e-Jt SB-1z dt                                     (3.33)

as may be verified by substitution in (3.32). Finally,

x = S-1y = -S-1eJt e-Jt SB-1z dt .

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

(I – A)x – B = z                                        (3.34)

where A = [a_ij] is the matrix of flow coefficients and B = [b_ij] is the matrix of capital coefficients. As indicated in section 3D the associated homogeneous equation may be put in the form

x – D = 0                                             (3.35)

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:

Notes

Chapter 4.   Latent Roots and Principal Vectors

In Chapter 3 the solution of linear differential equations with constant coefficients was reduced to the problem of evaluating the latent roots and principal vectors of a certain matrix. A number of possible computational methods will now be considered with a view to choosing those methods best suited to the solution of the problem outlined in Chapter 1. Methods such as that of Jacobi which apply only to symmetric matrices will not be considered.

The complete solution of this problem requires the calculation of all the latent roots and all the principal vectors of the matrix (I – A)^-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 theorem¹¹ which follows.

Theorem: Let Ci be the closed region of the complex plane enclosed by a circle with center at aii and radius   ri = |aij| . Then all the latent roots of the matrix A = [aij] lie in the region formed by the union of the circles Ci, i = 1,2,…,n.

Proof: Let x be an eigenvector of A corresponding to the latent root λ. Then Ax = λx and

	aijxj = λxi ,		i = 1,2,…,n;
∴	aijxj = (λ – aii)xi
∴	|aij||xj| ≥ |λ – aii||xi| .

Let x_k = max |x_j| ≠ 0 since x is non-trivial. Then

|aij||xk| ≥ |λ – aii||xk| .

Dividing by |x_k| we have

|λ – aii| ≤ |aij|

Thus λ lies in one of the circles C_i and the theorem is proved.

Since A^T has the same latent roots as A a similar theorem holds for circles D_i with centers a_ii and radii

|aij| .

Clearly then the two bounds to the roots are given by

		|λmax| ≤ |aij|
and
		|λmax| ≤ |aij| .

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,A²x,…,A^px 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 λ₁ is the dominant root of A then J^p has the form^[d]

	Jp =		λ1p
				Jn2p(λ2)
					…
						Jnmp(λm)

where each submatrix is of the form

	Jkp(λk) =		λkp	nλkp-1	λkp-2		…		.
				λkp	nλkp-1		…
					λkp		…
						…
							λkp

Since λ₁ is dominant

	Jp = λ1p		1					.
				0
					…
						0

Hence

Jp(S-1x) = λ1pαe1 ,

where α is the first component of (S^-1x) and e₁ is a unit vector. Since the eigenvector of A corresponding to λ₁ is Se₁ by Eq. (3.17) and since x has a component in the direction of this eigenvector it follows that α is not zero. Thus

Apx = SJp(S-1x) = λ1pαSe1 .

The expression on the right is an eigenvector of A. Furthermore,

		(Apx)j	=	λ1pα(Se1)j	= λ1 .
		(Ap-1x)j		λ1p-1α(Se1)j

The convergence is linear and depends on the ratios

|

λ1 λk

|.

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 λ₁,λ₂,…,λ_n where λ₁ > λ₂ > … > λ_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 B = A – αI then the roots of B are λ_i – α. Thus if α ≥ λ₁ the dominant root of B becomes μ_n = λ_n – α. It is then possible to determine μ_n by applying the power method to B and so obtain λ_n. Furthermore, Eq. (3.20) shows that the roots of A^-1 are the reciprocals of the roots of A. Thus the roots of the matrix B = (A – αI)^-1 are the quantities

Clearly any one of the μ_i may be made dominant by choosing α close enough to the corresponding root α_i. Thus each value of μ_i (and from it λ_i) may be obtained from the matrix B by a suitable choice of α. In the absence of any knowledge of the location of λ_i, the required values of α must be obtained by repeated trials.

The Method of Deflation. This method, like the preceding, modifies the matrix A so as to make a secondary root dominant, but differs in that the dominant root only is changed and is replaced by zero. If the roots of A are distinct and if c_i and r_i are respectively the eigenvectors of A and A^T corresponding to the root λ_i and with length so chosen that r_i^Tc_i = 1, then

For

If i ≠ j then λ_i ≠ λ_j and hence r_i^Tc_j = 0. Equation (4.2) is therefore established. If will now be shown that

A = λiciriT .                                         (4.3)

Note that each product c_ir_i^T is a square matrix of order n. To establish the relation (4.3) consider the matrix product

	(A – λiciriT)cj	= Acj – λiciriTcj
		= λjcj – λiciδij
		= (λj – λj)cj = 0.

Since this is true for each of the n linearly independent vectors c_j, it follows that the matrix in brackets above is null and (4.3) is established.

Suppose the dominant root and corresponding eigenvector of A have been obtained. Then the corresponding eigenvector of A^T may be obtained by re-applying the power method, or more simply (since λ₁ is now known), by solving the singular homogeneous system

to obtain a non-trivial vector x. The “deflated” matrix

B = A – λ1c1r1T = 0·c1r1T + λiciriT

is now formed. Clearly B has the same roots as A except that the dominant root λ₁ has been replaced by zero. The root λ₂ 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.¹⁴ If A is symmetric r_i = c_i and a separate calculation of r₁ 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,A²x,…,A^px exists. For example, if λ₁ is dominant and n₁ = 3 then

	J3p(λ1)	=		λ1p	pλ1p-1	½ p(p–1)λ1p-2
					λ1p	pλ1p-1
						λ1p

	= ½ p(p–1)λ1p-2		2λ12		2λ1		1		.
			p(p–1)		(p–1)
					2λ12		2λ1
					p(p–1)		(p–1)
							2λ12
							p(p–1)

Therefore

	J3p(λ1) = ½ p(p–1)λ1p-2		0		0		1		.
					0		0
							0

Thus the sequence Ax,A²x,…,A^px 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

P(λ) = | J – λI | = (λk – λ)nk .

Thus the zeros of P(λ) are the latent roots of J and each root occurs with the same multiplicity in P(λ) as in J. This also holds true for a general matrix A since similar matrices have the same characteristic polynomial, as will now be demonstrated.

Suppose A and B are similar matrices, then

and since

Clearly then the latent roots of a matrix may be obtained in the following two steps:

Since the degree of P(λ) is the same as the order of the matrix the numerical solution of the characteristic equation may be very difficult. The problem of solving polynomial equations is taken up in Chapter 5. We now proceed to consider two closely related models for generating the characteristic polynomial of a matrix.

(1) The Method of Bingham ¹⁵

Let

be a polynomial with zeros λ₁,λ₂,…,λ_n. Defining

sk = λik ,                                             (4.4)

it is possible with the help of Newton’s identities (reference 16, p. 147) to express the coefficients of P(λ) in terms of the quantities s_k, k = 1,2,…,n as follows:

The quantities s_k corresponding to the characteristic polynomial of a canonical matrix J are easily obtained, for clearly the sum of the diagonal elements of J is equal to s₁. Furthermore, Eqs. (3.18) and (3.19) show that the sum of the diagonal elements of J^k is equal to s_k.

The sum of the diagonal elements of a matrix A is usually called the trace of A or more briefly Tr A. Using this term we may write

Moreover, the trace is invariant under a similar transformation for if

then in terms of components

	bij    = s-1ikakmsmj
and
	Tr B = bii = s-1ikakmsmi
	= akm (smis-1ik)
	= akmδmk = akk = Tr A.

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 ¹⁷

The Latent Roots. Let c₀,c₁,…,c_n be the coefficients of the characteristic polynomial of the matrix A and let F_j(A) be a sequence of polynomial functions of A defined by

	F0(A) = I,		(4.8)
	Fk(A) = AFk-1(A) + ckI,     k = 1,2,…,n.

Now

and in general

by virtue of (4.5).

Using this expression for c_k, Eqs. (4.8) may be written

	F0(A) = I,		(4.9)
	Fk(A) = AFk-1(A) – 1/k Tr[AFk-1(A)] I.

Equations (4.9) define the sequence of operations in the Frame method.

The Bingham process computes the sequence of matrices A,A²,A³,…,Aⁿ and from the traces s₁,s₂,…,s_n of these matrices computes the coefficients by the use of Newton’s identities. The Frame method also uses repeated multiplication by A but at each step modifies the matrix so as to derive the coefficients directly without explicit use of Newton’s identities. The Frame method has the further advantage that

This relation serves as a useful overall check on the accuracy of the computation.

The Principal Vectors. The adjoint of (A – λI) may be expressed as a polynomial in λ with coefficients formed from the sequence of matrices (4.8) as follows:

by (4.8) and since F_n = 0 then

by the definition of the characteristic equation. Thus Δ(A – λI) as given by (4.11) is indeed the adjoint of (A – λI). When evaluated at a value of λ equal to a latent root the columns are proportional to the eigenvectors associated with that root. Furthermore since λ is general in the expression (4.11) the derived adjoints may be obtained by differentiation. The derived adjoints furnish all the principal vectors as demonstrated in section 3C.

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 F_k(A) defined by (4.8). Because of round-off error the computed value _k(A) will differ from F_k(A) and we may write

k(A) – Fk(A) = Ek.

The propagation of the error introduced at the first step will be investigated first. Let the length of the longest column vector of E_k be denoted by ε_k. Then

	k+1(A) – Fk+1(A)	= A[k(A) – Fk(A)] – 1/k Tr{A[k(A) – Fk(A)]}I
		= AEk – 1/k Tr(AEk)I.

Since for any vector x,

where λ is the dominant root of A, then

	εk+1	≤ | λ |εk + n/k | λ |εk
		= | λ | n+k k εk
∴	εn-1	≤ | λ |n-1 (2n-1)!  [(n-1)!]2 ε0 .

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:

ε < ε0

[2(n–k)+1]!  [(n–k)!]2

| λ |n-k .

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 S₀ is an approximation to S then the method of Jahn¹⁸ may be used to improve this approximation. The object is to compute S where AS = SΛ and Λ is a diagonal matrix. Since S₀ ≐ S we have

where Λ₁ is diagonal and B is a matrix whose elements are small; in particular the diagonal elements are zero. Let the diagonal elements of Λ₁ be λ₁,λ₂,…,λ_n and let C = [c_rs] where

crs =

brs λr – λs

.

Then S₁ = S₀ (I + C) is an improved approximation to S correct to second order. For clearly B = CΛ₁ – Λ₁C

∴	S0-1AS0 = Λ1 + CΛ1 – Λ1C     and     AS1 = AS0 (I + C).

Hence

F. Comparison of Methods

The amount of computation required to generate the characteristic equation by either the Frame or Bingham method is of the order n⁴ 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.

Notes

5.   Zeros of a Polynomial

A. Methods Depending on an Initial Approximation

The problem of locating the zeros of a polynomial of high degree is greatly simplified when good initial approximations to the zeros are known. This situation is often met in practice because of the frequent occurrence of families of polynomial equations. If, for example, the solution of a polynomial equation is required within some larger iterative process, the coefficients may be expected to vary gradually from one step to the next and the set of roots obtained at one step may be expected to provide good approximations to the roots of the next equation. Under these conditions, the Birge-Vieta¹⁹ method and the quadratic factor method are probably the best, since they converge quadratically near a root. Before proceeding to a discussion of these methods certain useful processes will be developed.

1. Synthetic Division¹⁶

Let P_n(x) = a₀xⁿ + a₁x^n-1 + … + a_n-1x + a_n be a polynomial of degree n with real coefficients. Let

where the coefficients b_i are defined by

	b0 = a0 ,		(5.1)
	bi = ai + pbi-1 ,       i = 1,2,…,n.

Then

Therefore

Hence b_n is the value of the polynomial evaluated at p and when b_n 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 b_n. The algorithm is therefore called synthetic division. Adopting the convention that all coefficients with negative subscripts are zero, Eq. (5.1) may be written more concisely as

Synthetic division by factors of any degree is easily defined but the case of a quadratic factor is of particular interest. In this case let

	bi = ai + pbi-1 + qbi-2,     i = 0,1,2,…,n–1,		(5.4)
	bn = an + qbn-2.

Then

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 y = x – p. The change in the coefficients of P_n(x) induced by such a translation may be derived from the relation

The coefficients b_k are most conveniently computed by the process known as Horner’s method.¹⁶ This method is nothing more than a repeated application of synthetic division to the successive quotient polynomials. For n = 4 the scheme of computation is as follows:

The polynomial in y is then

To prove that the polynomial in y may be so obtained let the successive quotient polynomials Q_k(x) be defined by

Then

	Pn(x)	= (x–p)[(x–p)Qn-2 + bn-1] + bn		(5.6)
		= (x–p){(x–p)[((x–p)Qn-3 + bn-2] + bn-1} + bn
		= … = (x–p)nb0 + (x–p)n-1b1 + … + (x–p)bn-1 + bn
		= b0yn + b1yn-1 + … + bn.

Repeating the differentiation of the relation

yields the values of the successive derivatives of P_n(x) evaluated at the point p. Thus

dkPn(x)    p dxk

= k!bn-k,     k = 0,1,2,…,n.                         (5.7)

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 x₀ sufficiently close to the quantities x_i defined by

xi+1 = xi –

f(xi) f '(xi)

,     i = 0,1,2,… ,

will converge to the root . The convergence is quadratic; i.e., the difference x_i – is approximately squared at each step. When the function f(x) is a polynomial the discussion of Horner’s method and Eq. (5.7) make it clear that two synthetic divisions with p = x₀ will furnish the quantities

The process is then repeated using the quantity

x1 = x0 –

bn bn-1

,

and in general

xi+1 = xi –

bn(xi) bn-1(xi)

.                                       (5.8)

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 x₀,x₁,… to be complex. Since a polynomial is an analytic function²¹ the derivative P_n'(x) exists for x complex, and the process carries over without change to complex numbers. (Note, however, that if the initial approximation x₀ is real, the successive x_i remain real.) Every operation is now complex and hence the amount of computation is increased by a factor of about four.

The complex 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 b_n-1 and b_n zero. Treating b_n-1 and b_n as functions of p and q and expanding in Taylor’s series we have

	bn(p+δp, q+δq)	= bn	+	∂bn ∂p	δp +	∂bn ∂q	δq + …
	bn-1(p+δp, q+δq)	= bn-1	+	∂bn-1 ∂p	δp +	∂bn-1 ∂q	δq + … .

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

	bn	+	∂bn ∂p	δp +	∂bn ∂q	δq = 0		(5.9)
	bn-1	+	∂bn-1 ∂p	δp +	∂bn-1 ∂q	δq = 0.

The quadratic factor method as first proposed by Bairstow²² and as presented in standard works^{23, 24} on numerical analysis requires the use of a second synthetic division to evaluate the derivatives in (5.9) and the application of the corrections δp and δq so derived. A modification of the method, which provides a more uniform computational scheme^[a] and simplifies the extension to a higher order process, will now be described. The proof is also simplified.

We defined a modified algorithm for synthetic division, similar to that given by (5.4), namely,

This process is completely uniform, without exception for the computation of c_n. Comparing with (5.4)

letting

Eq (5.5) becomes

The modification of the quadratic factor method now consists in choosing δp and δq so as to reduce to zero the quantities c_n-1 and c_n rather than b_n-1 and b_n. Note that the vanishing of c_n-1 and c_n is the necessary and sufficient condition for the vanishing of b_n-1 and b_n.

Treating x, p and q as independent variables and differentiating (5.11) partially with respect to p we have

0 = (x2 – px – q)

∂Q(x) ∂p

– xQ(x) + x

∂cn-1 ∂p

+

∂cn ∂p

– p

∂cn-1 ∂p

– cn-1

∴

xQ(x) + cn-1 = (x2 – px – q)

∂Q(x) ∂p

+ x

∂cn-1 ∂p

+ (

∂cn ∂p

– p

∂cn-1 ∂p

).

Thus ∂c_n-1/∂p and ∂c_n/∂p are respectively the first and second remainders on applying the modified algorithm to the polynomial xQ(x) + c_n-1, i.e. to the coefficients c₀…c_n-1. Let the coefficients so obtained be d₀…d_n-1, then

dn-2 =

∂cn-1 ∂p

and     dn-1 =

∂cn ∂p

.                          (5.12)

Differentiating (5.11) with respect to q we have

0 = (x2 – px – q)

∂Q(x) ∂q

– Q(x) + x

∂cn-1 ∂q

+ (

∂cn ∂q

– p

∂cn-1 ∂q

)

and the first and second remainders on applying the modified algorithm to Q(x) are ∂c_n-1/∂q and ∂c_n/∂q respectively. Hence

dn-3 =

∂cn-1 ∂q

and     dn-2 =

∂cn ∂q

.                          (5.13)

Applying Eq. (5.9) to the quantities c_n and c_n-1 and substituting the expressions obtained above for the derivatives we have

Thus:

	δp	=	-(cn-1dn-2 – cndn-3) D			(5.14)
	δq	=	(cn-1dn-1 – cndn-2) D
	D	=	dn-22 – dn-1dn-3.

The complete uniformity of the method is clearly show in the following example for n = 5. The initial approximations are p₀ and q₀. Coefficients with negative subscripts are by convention assumed to be zero.

δ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 equation f(x) = 0 may be derived by truncating the Taylor’s series expansion

f(x) = f(x0) + (   x0)f '(x0) + (x – x0)2

f "(x0) 2!

+ …

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

0 = f() ≅ f(x0) + ( – x0)f '(x0) + ( – x0)2

f "(x0) 2!

for . Then

≅ x0 –

f '(x0) ± √ {[f '(x0)]2 – 2f(x0)f "(x0)} f "(x0)

(5.15)

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 "(x₀) 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

Ap =

∂A ∂p

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 D = A_q, then G = C_q = A_pq. Clearly then

Extending the Eqs. (5.9) to include second order terms of the Taylor’s series and using the present notation we have

	0 = A + Apδp + Aqδq +	App 2	δp2 + Apqδpδq +	Aqq 2	δq2
	0 = B + Bpδp + Bqδq +	Bpp 2	δp2 + Bpqδpδq +	Bqq 2	δq2.

Hence

	Cδp + Dδq	= – (A +	F 2	δp2 + Gδpδq +	H 2	δpδq2)
	Dδp + Eδq	= – (B +	G 2	δp2 + Hδpδq +	I 2	δpδq2).

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

Aitken²⁵ has studied the extraction of real factors of higher degree than second by means of methods called “penultimate remaindering”. The modified quadratic factor method may also be extended to extract factors of higher degree. Unless all roots are known to be real the real factor must be of even degree. The case of a quartic factor will be developed briefly.

As for the quadratic factor the synthetic division algorithm

yields the quotient polynomial

such that

Again the modified algorithm

yields the same quotient polynomial Q(x). Adopting the notation

we have

Note that the vanishing of the quantities A, B, C, and D is the necessary and sufficient condition for the vanishing of the remainders b_n, b_n-1, b_n-2 and b_n-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 D_s as the first remainder. Thus if the scheme of two successive applications of the algorithm is shown as

we have D_s = K. The coefficient of x₂ in (5.21) is

But from the form of the remainders in (5.17) it is clear that the coefficient of x² is J - pK.

Similarly B_s = I and A_s = H. Rearranging (5.18) to the form

shows the left-hand side to be the polynomial formed of the coefficients c₀,c₁,…,c_n-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 D_p = H, C_p = G, B_p = F, A_p = E. In a similar way derivatives with respect to q and r may be obtained from Eqs. (5.19) and (5.20) and in matrix notation we have finally

Ap

Aq

Ar

As

=

E

F

G

H

.

Bp

Bq

Br

Bs

F

G

H

I

Cp

Cq

Cr

Cs

G

H

I

J

Dp

Dq

Dr

Ds

H

I

J

K

Expanding A, B, C, and D in Taylor’s series and retaining only linear terms we have (compare Eq. (5.9))

E

F

G

H

δp

=

-A

.

(5.22)

F

G

H

I

δq

-B

G

H

I

J

δr

-C

H

I

J

K

δs

-D

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 x² + ax + b and x² + αx + β.

Then

Setting

	b =	D	,	α =	C-Aβ	;		a = A - α	(5.27)
		β			b-β

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 C - AD^½ = 0 and a solution may then be obtained directly from the relations

Example.

	x4 + 2.5504x3 + 37.1185x2 - 38.4650x + 520.3597 = 0
	b =	520.3597 β	;		α =	-38.4650 - 2.5504β b - β	;
	a = 2.5504 - α;				ø = 37.1185 - b - β - aα.

Choose β0 ≅

D½ 2

.

A useful modification of the method consists in first making the substitution x= y + h with h = -A/4 so that the resulting equation in y contains on term in y³. Considering the equation in y we now have A = 0 and Eqs. (5.27) and (5.28) simplify to

	b = D/β,		α =	C b - β	,		a = - α,	(5.27a)
	ø(β) = B - b + α2 - β = 0.							(5.28a)

The special case mentioned above now exhibits itself immediately by the fact that C = 0.

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

	γ	=	β + b
then
	ø(γ)	=	B - γ + d
	ø'(γ)	=	- 1 + γd2 C2/2
where
	d	=	C2 . γ2 - 4D

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 γ² - 4D is required to compute b - β = √γ² - 4D and since ø' need not be recomputed in the final stages it is convenient to record both of these quantities as well.

Example. The substitution x = y - 0.6376 in the quartic of the previous example gives

Computing C², C²/2 and 4D we have

	d =	7009.8421 γ2 - 2237.9164	,	ø = 34.6783 - γ + d
	ø' = -	1 + γd2 . 3504.9210

A few rough mental calculations show that ø(50) > 0 and ø(52) < 0.

Hence b - β = √416.2704 = 20.4027,

Thus

	y =		-2.0518	± 5.6348i
			2.0518	± 3.3687i
	x =		-2.6894	± 5.6348i
			1.4142	± 3.3687i.

The quartic formed from these roots agrees with the given quartic to four decimal places in each coefficient.

Since

	ø"(γ) =	2C2(3γ2 + 4D)
		(γ2 - 4D)3

is easily computed it is possible to use the higher order correction

	- ø' ± √ ø' 2 - 2øø"
	ø"

where the sign is chosen to give the smaller correction. Application of this in the preceding example gives ø"(51) = 2.9410 and a correction of 0.5288 which is much better than that given by the Newton-Raphson method. However the simpler method is probably to be preferred.

8. Solution of the Sextic

An extension of the method of real factors to the sextic yields a practical computational method as follows: Assume that the sextic

has real factors x² + αx + β and x⁴ + ax³ + bx² + cx + d. Again assuming that A is made zero by translation we obtain the following relations

	a		=		-α			(5.29)
	b + β + α2		=		B
	c + bα - βα		=		C
	d + cα - bβ		=		D
	dα - cβ		=		E
	dβ		=		F.

Then if d = F/β, γ = d - β²,

	α =		E/2	±	√  (	E/2	)2 +	F + β [-D + β(B-β)]	,     α real          (5.30)
			γ			γ		γ

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

	b = B - β + α2,       c =	E - dα	.
		β

At least one root of (5.31) lies in the interval (-F^1/3,F^1/3). If all the zeros of P₆(x) are complex the zeros of (5.31) are all positive and exactly three in number. Hence either the interval (0,F^1/3) or the interval (F^1/3,∞) contains one real zero only. The two values of α given by Eq. (5.30) introduce a complication but it is usually easy to recognize which value should be chosen.

The method fails only if d - β² = 0, in which case β = F^1/3 and (since α is finite) FB³ - D³ = 0. The solution is then given immediately by

	β = F3,                   α =	1	(3β - 2B +	D	),
		3E		β

the expression for α being obtained as the limiting form of (5.30).

Example:

	A	=	0		B	=	46.5813		C	=	-89.2555
	D	=	1355.5763		E	=	-2198.2332		F	=	10076.5517
	d	=	10076.5517 β		γ	=	d - β2		δ	=	-1099.1166 γ
	α	=	δ ± √ δ2 + 100076.5517 + β [β(46.5813 - β) - 1355.5763] d - β2
	ø	=	-2198.2332 + 89.2555 + α [46.5813 - 2β + α2 - d ]. β β

Since 0 ≤ β ≤ F^1/3 ≅ 22 choose β₀ = 10. Linear prediction is used even if two successive values of ø do not differ in sign. R² represents the square of the radical.

B. Initial Approximations

1. Method of Search

Any point in the complex plane from which one of the methods of section A will converge to a particular zero of given polynomial will be considered a satisfactory initial approximation to that zero. If the complex plane is covered with a set of points on a sufficiently fine mesh then at least one point of the set will be sufficiently close to each zero of a given polynomial so as to serve as a suitable initial approximation. If bounds on the roots are known, the number of mesh points to be considered is finite and a program which will choose in succession the points of this finite set may therefore be used to provide initial approximations. Since the fineness of the mesh required is not generally known in advance it is advisable to cover the region of interest first on a coarse mesh and repeat on successively finer meshes until all zeros are located.

This approach was employed in locating the zeros of the truncated series for the exponential²⁶ namely

Sn(z) = zk/k!

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 Method²⁷

Let r₁,r₂,…,r_n be the zeros of the polynomial

The coefficients are expressible in terms of the zeros as follows:¹⁶

If the zeros are real and widely separated it is clear that

	- a1/a0 ≐ r1 ;
	- a2/a1 = (a2/a0)(- a0/a1) =	r1r2 (1 + r1r3/r1r2 + … + rn-1rn/r1r2) r1 (1 + r2/r1 + r3/r1 + … + rn/r1)	≐ r2 ;

and in general