Showing papers in "Mathematics of Computation in 1971"

Numerical methods for computing angles between linear subspaces

Abstract: Assume that two subspaces F and G of unitary space are defined as the ranges (or nullspaces) of given rectangular matrices A and B. Accurate numerical methods are developed for computing the principal angles $\theta_k (F,G)$ and orthogonal sets of principal vectors $u_k\ \epsilon\ F$ and $v_k\ \epsilon\ G$, k = 1,2,..., q = dim(G) $\leq$ dim(F). An important application in statistics is computing the canonical correlations $\sigma_k\ = cos \theta_k$ between two sets of variates. A perturbation analysis shows that the condition number for $\theta_k$ essentially is max($\kappa (A),\kappa (B)$), where $\kappa$ denotes the condition number of a matrix. The algorithms are based on a preliminary QR-factorization of A and B (or $A^H$ and $B^H$), for which either the method of Householder transformations (HT) or the modified Gram-Schmidt method (MGS) is used. Then cos $\theta_k$ and sin $\theta_k$ are computed as the singular values of certain related matrices. Experimental results are given, which indicates that MGS gives $\theta_k$ with equal precision and fewer arithmetic operations than HT. However, HT gives principal vectors, which are orthogonal to working accuracy, which is not in general true for MGS. Finally the case when A and/or B are rank deficient is discussed.

The fast Fourier transform in a finite field

Abstract: A transform analogous to the discrete Fourier transform may be defined in a finite field, and may be calculated efficiently by the 'fast Fourier transform' algorithm. The transform may be applied to the problem of calculating convolutions of long integer sequences by means of integer arithmetic.

On descent from local minima

Abstract: When a local minimum of a function of several variables has been found by use of an algorithm for finding such minima numerically, one often runs the same algorithm many times with different starting values in the hopes of finding a lower minimum. Here, under the assumption that a local minimum is known, a process with analytical criteria is described which sometimes finds smaller local minima in an algorithmic manner. Methods of descent are useful for minimizing functions of several variables. Generally, one can always obtain points (if such exist) for which the gradient vanishes, and moreover, points which are local minima. At saddle points. one can continue descent with second derivative information. A point which is a local minimum for a function may or may not be a global minimum. At this juncture one resorts to search techniques to attempt to further decrease the function. The process to be described sometimes finds smaller local minima in an algorithmic manner with analytical criteria. One has no general test, of course, for a global minimum. Consider first the problem of finding the global minimum for a 2nth degree polynomial P1(x) in one variable. The coefficient of x24 will be positive. Let xl be a local minimizer of P1. Then one may write

On Jacobi and Jacobi-like algorithms for a parallel computer

Abstract: Many existing algorithms for obtaining the eigenvalues and eigenvectors of matrices would make poor use of such a powerful parallel computer as the ILLIAC IV. In this paper, Jacobi's algorithm for real symmetric or complex Hermitian matrices, and a Jacobi-like algorithm for real nonsymmetric matrices developed by P. J. Eberlein, are modified so as to achieve maximum efficiency for the parallel computations. 1. Introduction. With the advent of parallel computers, the study of compu- tationally massive problems became economically possible. Such problems include, for example, solution of sets of partial differential equations over sizable grids, and multiplication, inversion, or determination of eigenvalues and eigenvectors of large matrices. An example of a parallel computer is the ILLIAC IV.* This computer is es- sentially an array of coupled arithmetic units driven by instructions from a common control unit. Each of the arithmetic units, called processing elements (PE's), have 2048 words of 64-bit memory with an access time under 420 nanoseconds. Each PE is capable of 64-bit floating-point multiplication in about 550 nanoseconds. Two 32-bit floating-point operations may be performed in each PE in approximately the same times. The PE instruction set is similar to that of conventional machines with two exceptions. First, the PE's are capable of communicating data to four neigh- boring PE's by means of routing instructions. Second, the PE's are able to set their own mode registers to effectively disable or enable themselves. For a more detailed

The convergence of an algorithm for solving sparse nonlinear systems

Abstract: A new algorithm for solving systems of nonlinear equations where the Jacobian is known to be sparse is shown to converge locally if a sufficiently good initial estimate of the solution is available and if the Jacobian satisfies a Lipschitz condition. The results of numerical experiments are quoted in which systems of up to 600 equations have been solved

Uniform computation of the error function and other related functions

Abstract: Uniform methods of computation, to any required degree of accuracy, for the error and other closely related functions are given. from which the above quadrature formula was derived, were included in the contour c. This means that t was restricted to the range t > h2/47r2. In this paper, we extend the method to include the cases when h is such that the poles at z = -iw lie on and outside the contour c. The respective results are

Algorithms for Hermite and Smith normal matrices and linear Diophantine equations

Abstract: New algorithms for constructing the Hermite normal form (triangular) and Smith normal form (diagonal) of an integer matrix are presented. A new algorithm for determining the set of solutions to a system of linear diophantine equations is presented. A modification of the Hermite algorithm gives an integer-preserving algorithm for solving linear equations with real-valued variables. Rough bounds for the number of operations are cubic polynomials involving the order of the matrix and the determinant of the matrix. The algorithms are valid if the elements of the matrix are in a principal ideal domain.

Least Squares Methods for 2mth Order Elliptic Boundary-Value Problems

Abstract: In this paper we consider a general class of boundary-value problems for 2mth order elliptic equations including nonhomogeneous essential boundary conditions and nonselfadjoint problems. Approximation methods involving least squares approximation of the data are presented and corresponding error estimates are proved. These methods can be considered in the category of Rayleigh-Ritz-Galerkin methods and have the special feature that the trial functions need not satisfy the boundary conditions. A special case of the trial functions which is studied are spline functions defined on a uniform mesh of width h (or more generally piecewise polynomial functions). For a given \"well set\" boundary-value problem for a 2mth order operator the theory presented will provide a method with any prescribed order of accuracy r which is optimal in the sense that the best approximation in the underlying subspace is of order of accuracy r.

Whittaker's Cardinal Function in Retrospect*

Abstract: This paper exposes properties of the Whittaker cardinal function and illustrates the use of this function as a mathematical tool. The cardinal function is derived using the Paley-Wiener theorem. The cardinal function and the central-difference expansions arelinked through their similarities. A bound is obtained on the difference between the cardinal func- tion and the function which it interpolates. Several cardinal functions of a number of special functions are examined. It is shown how the cardinal function provides a link between Fourier series and Fourier transforms, and how the cardinal function may be used to solve integral equations.

The method of Christopherson for solving free boundary problems for infinite journal bearings by means of finite differences

Abstract: A method for solving free boundary problems for journal bearings by means of finite differences has been proposed by Christopherson We analyse Christopherson's method in detail for the case of an infinite journal bearing where the free boundary problem is as follows: Given T > 0 and /(t) find r E (0, T) and p(t) such that (i) (h3p')' = h' for t E (0, r), (ii) p(O) = 0, (iii) p(t) = 0 for t (T, Ir7, and (iv) p'(r - 0) = 0 First, it is shown that the discrete approximation is accurate to O((At)2) where At is the step size Next, it is shown that the discrete problem is equivalent to a quadratic pro- gramming problem Then, the iterative method for computing the discrete approximation is analysed Finally, some numerical results are given 1 Introduction A journal bearing consists of a rotating cylinder which is separated from a "bearing surface" by a thin film of lubricating fluid (see Fig 1) The fluid is fed in at A and flows out at B The width of the film is smallest at C, and we set 1 = O/0, where 0 is as shown in Fig 1 Between C and B, the width of the film increases so that the pressure in the lubricating fluid may be expected to decrease We assume that for t - r the pressure becomes so low that the fluid vaporizes The point t T the interface between the two phases of the fluid, is called the free boundary The mathematical problem can now be formulated (see Pinkus and Sternlicht

Least squares methods for 2th order elliptic boundary-value problems

Abstract: In this paper we consider a general class of boundary-value problems for 2mth order elliptic equations including nonhomogeneous essential boundary conditions and nonselfadjoint problems. Approximation methods involving least squares approximation of the data are presented and corresponding error estimates are proved. These methods can be considered in the category of Rayleigh-Ritz-Galerkin methods and have the special feature that the trial functions need not satisfy the boundary conditions. A special case of the trial functions which is studied are spline functions defined on a uniform mesh of width h (or more generally piecewise polynomial functions). For a given \"well set\" boundary-value problem for a 2mth order operator the theory presented will provide a method with any prescribed order of accuracy r which is optimal in the sense that the best approximation in the underlying subspace is of order of accuracy r.

Uniform asymptotic expansions of the Jacobi polynomials and an associated function

Abstract: Asymptotic expansions have been obtained using two theorems due to Olver for the Jacobi polynomials and an associated function. These expansions are uniformly valid for complex arguments over certain regions, for large values of the order.

On the Convergence of Broyden's Method for Nonlinear Systems of Equations

Abstract: This paper uses majorant techniques to study the convergence of Broyden's single-rank update method for nonlinear systems of equations. It also contains a very elementary proof of the local convergence of the method. The heart of the method is a procedure for generating an approximation to the Jacobian of the system using only in- formation on hand and not requiring partial derivatives. 1. Introduction. C. G. Broyden (2) suggested an algorithm for iterating to a solution of a system of nonlinear equations which has shown its mettle in dealing with practical problems. The purpose of this paper is to provide a Kantorovich-type analysis and an elementary local convergence proof for this method. In fact, the analysis is applicable to the entire class of 'single-rank update' methods just to the extent that it seems to justify heuristically the generally superior performance of Broyden's method over the rest of the class. These methods generate sequences {xn}, {//"}, one consisting of approximate roots and the other of the corresponding approximate inverse Jacobian matrices. At the nth step, one obtains xn+i from xn and H" by setting (1) *»+l = Xn — JnHnFiXn), where F(xn) = (ji(x"), • ■ • , fn(xn))T is the residual vector at the point x" and yn is a real number about which more will be said later. Hn+i is obtained by using an idea due to Davidon (5). It is required to satisfy the equation

Homomorphisms of knot groups on finite groups

Abstract: We describe trial and error computer programs for finding certain homomorphisms of a knot group on a special projective group LF(2, p), p prime, and programs to evaluate H1(=; Z) where M is a finitely sheeted branched covering space of S3 associated with such a homomorphism. These programs have been applied to several collections of examples, in particular to the Kinoshita-Terasaka knots, and we state numerous conjectures based on these experiments. About forty years ago a universal method for obtaining algebraic invariants of knot type was proposed and became standard. The method, as applied to a knot k of type K and having group 7rK = 7r1(S3 -k; *), begins with the determination of the homomorphisms of wrK on a given group G. These homomorphisms fall into equivalence classes under the action of the automorphisms of G, and a crude preliminary invariant of K is the number of homomorphism classes. In the next stage of the method, we fix a transitive permutation representation of G, perhaps of infinite degree. Each homomorphism class of 7rK on G is associated with a covering space c% of S3 k such that the number of sheets in the covering is the degree of the permutation representation, and the group H&(Ut; Z) is an algebraic invariant of the knot type K. In this paper, we shall discuss the means and results of implementing the universal method on a computer when the group G is chosen to be one of the special projective groups L, = LF(2, p) = PSL(2, p), where p is a prime integer. The universal method has been most thoroughly examined in the case where the group G is cyclic. The determination of the homomorphism classes becomes completely trivial, and all the homology invariants can be deduced from a single matrix, the Alexander matrix. These "cyclic invariants" have been applied with good effect to just about every problem in knot theory, but, alas, when the Alexander polynomial A(x) of the knot reduces to the constant 1 these invariants degenerate and are worthless. It is no good choosing a solvable group for G in such a case; to get useful results from the universal method, we must use nonsolvable groups. Of course, simple groups receive first consideration in this context, and of the families of classical finite simple groups, the family IL,4 is the most manageable. As further encouragement for this study, the group L4 is isomorphic to the alternating group AS, and R. H. Fox has shown in several papers ([5] is a good example) that the homomorphisms of a knot group on A5 have some interesting applications. There are, however, very few clues to suggest what the best method for finding homomorphisms on the special projective groups could be. In addition, almost nothing is known about the homology invariants of homomorphisms on noncyclic groups G. It therefore seems reasonable Received July 17, 1970, revised December 2, 1970. AMS 1969 subject classifications. Primary 5520.

A fourth-order finite-difference approximation for the fixed membrane eigenproblem

Abstract: The fixed membrane problem All + Xu = 0 in Q, u = 0 on aQ, for a bounded region Q of the plane, is approximated by a finite-difference scheme whose matrix is monotone. By an extension of previous methods for schemes with matrices of positive type, 0(h 4) convergence is shown for the approximating eigenvalues and eigenfunctions, where h is the mesh width. An application to an approximation of the forced vibration problem Au + qu = J in Q, u = 0 in aQ, is also given.

Comparing Error Estimators for Runge-Kutta Methods

Abstract: A way is proposed to compare local error estimators. This is applied to the major estimators for fourth-order Runge-Kutta procedures. An estimator which leads to a production code 18% more efficient than a code using the standard one is recommended and supported by numerical tests. 1. Introduction. The design of a production code for the numerical solution of ordinary differential equations requires three major choices (at least): the method of advancing the solution one step, the method of estimating the local error incurred in one step, and the strategy for choosing the next step length. The first choice has been intensively studied in the literature of numerical analysis but the latter two have been rather neglected in spite of their great practical significance. The aim of this paper is to compare alternatives for estimating the local error. We restrict ourselves to estima- tors appropriate to fourth-order Runge-Kutta methods for two reasons. Production codes based on fourth-order Runge-Kutta processes are quite common and we are able to recommend improvements based on this study. Secondly, texts frequently mention how easy it is to change the step size when using Runge-Kutta procedures but they often fail to point out the associated practical disadvantages. Namely, it is rather difficult and expensive to estimate the local error and so decide when to change the step size. Furthermore, the standard estimator forces one to use a relatively crude step size strategy. It is, therefore, important to consider the reliability and efficiency of error estimating schemes. In Section 2, we shall propose a way of comparing numerically local error estima- tors which seems to be of broad applicability. Section 3 describes the most useful error estimators known to the authors, which are then compared in the fourth section.

The lattice structure of multiplicative congruential pseudo-random vectors

Abstract: The lattice structure of points in an n-dimensional space produced by an ap- propriate grouping of pseudo-random numbers obtained from multiplicative congruential generators is discussed. Examples are given for 2 < n < 6. The work is based on the theory of the reduction of positive quadratic forms in n variables.

Calculation of the Gamma Function by Stirling's Formula

Abstract: In this paper, we derive a simple error estimate for the Stirling formula and also give numerical coefficients. Stirling's formula is: log r(s) - (s - 2) log s - s + 2 log 27r

An implementation of Christoffel’s theorem in the theory of orthogonal polynomials

Abstract: An algorithm for the construction of the polynomials associated with the weight function w(t)P(t) from those associated with w(t) is given for the case when P(t) is a polynomial which is nonnegative in the interval of orthogonality. The relation of the algorithm to the LR algorithm is also discussed. Introduction. In several problems of numerical analysis, particularly in the construction of Gaussian quadrature rules with preassigned nodes, the following problem arises. Given the orthogonal polynomials {pj(t)j associated with a weight function w(t) on the interval (a, b) and a polynomial P(t) of degree m which is nonnegative on the interval (a, b), construct the orthogonal polynomials { q j(t)} associated with the weight function P(t)w(t) on the same interval. A theorem of Christoffel [1] gives an explicit expression for the polynomial qQ(t) in the form Pn(t) Pn+ 1(t) ..P., .(t) Pn(zl) Pn+l(Zl) ..P.+. (z1) (l1) q.(t)P(t) = P.(Z2) P.+l(Z2) P.+. (Z2) Pn(Zm) Pn+i(Zm) ... Pn+m(Zm) where the Zk, k = I(I)m, are the roots of P(t). If some root, zi, is of multiplicity j, then the corresponding rows of (1) are replaced by the derivatives of order 0, 1, ... * j 1 of the polynomials pr(t), r n(l)n + m, at t = zi. For numerical calculations, Eq. (1) is very clumsy to use, even for simple evaluation of the polynomial qn(t) at a point, unless m is small. Often, the three-term recurrence relation (2) pj(t) = (tbj)pj_(t) -gjpj2(t), j = 1, 2, * with p0(t) = 1 and p1(t) = 0, is known because it is more convenient to obtain [4], [5] and to use [5], [6]. The main result of this paper is to prove a theorem, equivalent to Christoffel's, which states an explicit construction of the three-term recurrence relation (3) qj(t) = (t Bj)q1_1(t) Gjqq2(t), j = 1, 2, . . Received March 9, 1970, revised May 13, 1970. AMS 1969 subject classifications. Primary 6510; Secondary 6525, 6550, 3345.

The triangular decomposition of Hankel matrices

Abstract: An algorithm for determining the triangular decomposition H R*DR of a Hankel matrix H using O(n') operations is derived. The derivation is based on the Lanczos algorithm and the relation between orthogonalization of vectors and the triangular decomposition of moment matrices. The algorithm can be used to compute the three-term recurrence relation for orthogonal polynomials from a moment matrix.

Pseudo-random numbers. The exact distribution of pairs

Abstract: 4bstract. Pseudo-random numbers are usually generated by linear congruential methods. Starting with an integer yo, a sequence (y1 j is constructed by yi+1 _ ay; + r (mod m), m, a, r being integers. The derived fractions xi _ y,/mn are taken as samples from the uniform distribution on [0, 1). In this paper it is shown that the joint probability distribution of pairs xi, xi+, cani be calculated exactly. Explicit calculations slhow that this distribution is surprisingly near to the uniform distribution for most 'reasonable' generators. The best approximation to the uniform distribution on the unit-square is achieved if the continued frawtion for at and m (or at and mr/I) is long.

Adjusted forms of the Fourier coefficient asymptotic expansion and applications in numerical quadrature

Abstract: The conventional Fourier coefficient asymptotic expansion is derived by means of a specific contour integration. An adjusted expansion is obtained by deforming this contour. A corresponding adjustment to the Euler-Maclaurin expansion exists. The effect of this adjustment in the error functional for a general quadrature rule is investigated. It is the same as the effect of subtracting out a pair of complex poles from the integrand, using an unconventional subtraction function. In certain applications, the use of this subtraction function is of practical value. An incidental result is a direct proof of Erdelyi's formula for the Fourier coefficient asymptotic expansion, valid when f(x) has algebraic or logarithmic singularities, but is otherwise analytic. 1. The Fourier Coefficient Asymptotic Expansion. The Fourier coefficient asymptotic expansion (F.C.A.E.) (1.3) below is a classical formula which is elementary to derive using a standard application of the formula for integration by parts, namely, rl~~~~~~~~(v 1 (-)/2)

An Example of Ill-Conditioning in the Numerical Solution of Singular Perturbation Problems*

Abstract: The use of finite-difference methods is considered for solving a singular perturba- tion problem for a linear ordinary differential equation with an interior turning point. Computational results demonstrate that such problems can lead to very ill-conditioned matrix equations.