Showing papers in "Mathematics of Computation in 1980"

Updating Quasi-Newton Matrices With Limited Storage

Abstract: We study how to use the BFGS quasi-Newton matrices to precondition minimization methods for problems where the storage is critical. We give an update formula which generates matrices using information from the last m iterations, where m is any number supplied by the user. The quasi-Newton matrix is updated at every iteration by dropping the oldest information and replacing it by the newest informa- tion. It is shown that the matrices generated have some desirable properties. The resulting algorithms are tested numerically and compared with several well- known methods. 1. Introduction. For the problem of minimizing an unconstrained function / of n variables, quasi-Newton methods are widely employed (4). They construct a se- quence of matrices which in some way approximate the hessian of /(or its inverse). These matrices are symmetric; therefore, it is necessary to have n(n + l)/2 storage locations for each one. For large dimensional problems it will not be possible to re- tain the matrices in the high speed storage of a computer, and one has to resort to other kinds of algorithms. For example, one could use the methods (Toint (15), Shanno (12)) which preserve the sparsity structure of the hessian, or conjugate gradient methods (CG) which only have to store 3 or 4 vectors. Recently, some CG algorithms have been developed which use a variable amount of storage and which do not require knowledge about the sparsity structure of the problem (2), (7), (8). A disadvantage of these methods is that after a certain number of iterations the quasi-Newton matrix is discarded, and the algorithm is restarted using an initial matrix (usually a diagonal matrix). We describe an algorithm which uses a limited amount of storage and where the quasi-Newton matrix is updated continuously. At every step the oldest information contained in the matrix is discarded and replaced by new one. In this way we hope to have a more up to date model of our function. We will concentrate on the BFGS method since it is considered to be the most efficient. We believe that similar algo- rithms cannot be developed for the other members of the Broyden 0-class (1). Let / be the function to be nnnimized, g its gradient and h its hessian. We define

LINPACK User's Guide.

Abstract: We provide further discussion of the use of least-squares techniques for this and G. W. Stewart, LINPACK Users' Guide (Society of Industrial and Applied. User's Guide for Intel® Math Kernel Library 11.3 for Linux* OS. Revision: Benchmark your cluster with Intel® Optimized MP LINPACK Benchmark for Clusters. Running Linpack on Linux High performance computer. No problem. I need a working guide for hpl if someone can do help me at this. The errors are.

Polynomial approximation of functions in Sobolev spaces

Abstract: Constructive proofs and several generalizations of approximation results of J. H. Bramble and S. R. Hilbert are presented. Using an averaged Taylor series, we represent a function as a polynomial plus a remainder. The remainder can be manipulated in many ways to give different types of bounds. Approximation of functions in fractional order Sobolev spaces is treated as well as the usual integer order spaces and several nonstandard Sobolev-like spaces.

An incomplete factorization technique for positive definite linear systems

Abstract: This paper describes a technique for solving the large sparse symmetric linear systems that arise from the application of finite element methods. The technique combines an incomplete factorization method called the shifted incomplete Cholesky factorization with the method of generalized conjugate gradients. The shifted incomplete Cholesky factorization produces a splitting of the matrix A that is dependent upon a parameter ..cap alpha... It is shown that if A is positive definite, then there is some ..cap alpha.. for which this splitting is possible and that this splitting is at least as good as the Jacobi splitting. The method is shown to be more efficient on a set of test problems than either direct methods or explicit iteration schemes.

Monotone difference approximations for scalar conservation laws

Abstract: : A complete self-contained treatment of the stability and convergence properties of conservation form, monotone difference approximations to scalar conservation laws in several space variables is developed. In particular, the authors prove that general monotone difference schemes always converge and that they converge to the physical weak solution satisfying the entropy condition. Rigorous convergence results follow for dimensional splitting algorithms when each step is approximated by a monotone difference scheme. The results are general enough to include, for instance, Godunov's scheme, the upwind scheme (Differenced through stagnation points), and the Lax-Friedrichs scheme together with appropriate multi-dimensional generalizations. (Author)

On the coupling of boundary integral and finite element methods

Abstract: We prove some error estimates for a procedure obtained by combining the boundary integral method and the usual finite element method.

The spectral transformation Lánczos method for the numerical solution of large sparse generalized symmetric eigenvalue problems

Abstract: A new algorithm is developed which computes a specified number of eigenvalues in any part of the spectrum of a generalized symmetric matrix eigenvalue problem. It uses a linear system routine (factorization and solution) as a tool for applying the Lanczos algorithm to a shifted and inverted problem. The algorithm determines a sequence of shifts and checks that all eigenvalues get computed in the intervals between them. It is shown that for each shift several eigenvectors will converge after very few steps of the Lanczos algorithm, and the most effective combination of shifts and Lanczos runs is determined for different sizes and sparsity properties of the matrices. For large problems the operation counts are about five times smaller than for traditional subspace iteration methods. Tests on a numerical example, arising from a finite element computation of a nuclear power piping system, are reported, and it is shown how the performance predicted bears out in a practical situation.

Stable and entropy satisfying approximations for transonic flow calculations

Abstract: Finite difference approximations for the small disturbance equation of transonic flow are developed and analyzed. New schemes of the Cole-Murman type are presented fpr which nonlinear stability is proved. The Cole-Murman scheme may have entropy violating expansion shocks as solutions. In the new schemes the switch between the subsonic and supersonic domains is designed such that these nonphysical shocks are guaranteed not to occur. Results from numercial calculations are given which illustrate these conclusions

Analysis of mixed methods using mesh dependent norms

Abstract: : This paper presents a new approach to the analysis of mixed methods for the approximate solution of 4th order elliptic boundary value problems. In this approach one introduces a pair of mesh dependent norms and proves the approximation method is stable with respect to these norms. The error estimates then follow in a direct manner. In a mixed method, one introduces an auxiliary variable, usually representing another physically important quantity, and writes the differential equation as a lower order system. One then considers Ritz-Galerkin approximation schemes based on a variational formulation of this lower order system, thereby obtaining direct approximations to both the original and auxiliary variables. Three particular mixed methods for the approximate solution of the biharmonic problem are examined in detail.

Iterative refinement implies numerical stability for Gaussian elimination

Abstract: Because of scaling problems, Gaussian elimination with pivoting is not always as accurate as one might reasonably expect. It is shown that even a single iteration of iterative refinement in single precision is enough to make Gaussian elimination stable in a very strong sense. Also, it is shown that without iterative refinement row pivoting is inferior to column pivoting in situations where the norm of the residual is important.

Convergence of multigrid iterations applied to difference equations

Abstract: Convergence proofs for the multi-grid iteration are known for the case of finite element equations and for the case of some difference schemes discretizing boundary value problems in a rectangular region. In the present paper we give criteria of convergence that apply to general difference schemes for boundary value problems in Lipschitzian regions. Furthermore, convergence is proved for the multi-grid algorithm with Gauss-Seidel's iteration as smoothing procedure.

Generalized OCI schemes for boundary layer problems

Abstract: A family of tridiagonal formally fourth-order difference schemes is developed for a class of singular perturbation problems. These schemes have no cell Reynolds number limitation and satisfy a discrete maximum principle. Error estimates and numerical results for this family of methods are given, and are compared with those for several other schemes.

A weak discrete maximum principle and stability of the finite element method in _{∞} on plane polygonal domains. I

Abstract: Let Q2 be a polygonal domain in the plane and Shr(92) denote the finite element space of continuous piecewise polynomials of degree 2) defined on a quasi-uniform triangulation of Q2 (with triangles roughly of size h) It is shown that if uh E Sh(n) is a "discrete harmonic function" then an a priori estimate (a weak maximum principle) of the form 11uhIllL 42) This says that (modulo a logarithm for r = 2) the finite element method is bounded in Lo0, on plane polygonal domains 0 Introduction and Statement of Results The purpose of this paper is to discuss some estimates for the finite element method on polygonal domains In particular, we shall consider the validity of (for want of a better terminology) a "discrete weak maximum principle" for discrete harmonic functions and then use this result to discuss the boundedness in Loo of the finite element projection In this part we shall discuss the case of a quasi-uniform mesh In Part II we shall concern ourselves with meshes which are refined near points Let us first formulate the problems we wish to consider and state our results References to other work in the literature which are relevant to our considerations will be given as we go along For simplicity let Q2 be a simply connected (this is not essential) polygonal domain in R2 with boundary 3Q2 and maximal interior angle a, 0 < a < 2ir, where we emphasize that in general Q2 is not convex On Q2 we define a -family of finite element spaces For simplicity of presentation we shall restrict ourselves to a special but important class of piecewise polynomials For each 0 < h < 1, let Th denote a triangulation of Q2 with triangles having straight edges We shall assume that each triangle r is contained in a sphere of radius h and contains a sphere of radius yh for some positive constant y We shall also assume that the family {Tn } of triangulations Received September 19, 1978 AMS (MOS) subject classifications (1970) Primary 65N30, 65N15 *This work was supported in part by the National Science Foundation ? 1980 American Mathematical Society 0025-571 8/80/0000-0004/$0475 77 This content downloaded from 1575539253 on Wed, 08 Jun 2016 05:26:01 UTC All use subject to http://aboutjstororg/terms

The construction of Jacobi and periodic Jacobi matrices with prescribed spectra

Abstract: : The spectral properties of Jacobi and periodic Jacobi matrices are analyzed and algorithms for the construction of Jacobi and periodic Jacobi matrices with prescribed spectra are presented. Numerical evidence demonstrates that these algorithms are of practical utility. These algorithms have been used in studies of the periodic Toda lattice, and might also be used in studies of inverse eigenvalue problems for Sturm-Louiville equations and Hill's equation. (Author)

On an algorithm for finding a base and a strong generating set for a group given by generating permutations

Abstract: This paper deals with the problem of finding a base and strong generating set for the group generated by a given set of permutations. The concepts of base and strong generating set were introduced by Sims (51, (61 and provide the most effective tool for computing with permutation groups of high degree. One algorithm, originally proposed by Sims [71, is described in detail; its behavior on a number of groups is studied, and the influence of certain parameters on its performance is investigated. Another algorithm, developed by the author, is given, and it is shown how the two algorithms may be combined to yield an exceptionally fast and

Reciprocal polynomials having small measure. II

Abstract: The measure of a monic polynomial is the product of the absolute value of the roots which lie outside and on the unit circle. We describe an algorithm, based on the root-squaring method of Graeffe, for finding all polynomials with integer coefficients whose measures and degrees are smaller than some previously given bounds. Using the algorithm, we find all such polynomials of degree at most 16 whose measures are at most 1.3. We also find all polynomials of height 1 and degree at most 26 whose measures satisfy this bound. Our results lend some support to Lehmer's conjecture. In particular, we find no noncyclotomic polynomial whose measure is less than the degree 10 example given by Lehmer in 1933.

Negative norm estimates and superconvergence in Galerkin methods for parabolic problems

Abstract: Negative norm error estimates for semidiscrete Galerkin-finite element methods for parabolic problems are derived from known such estimates for elliptic problems and applied to prove superconvergence of certain procedures for evaluating point values of the exact solution and its derivatives. Our first purpose in this paper is to show how known negative norm error estimates for Galerkin-finite element type methods applied to the Dirichlet problem for second order elliptic equations can be carried over to initial-boundary value problems for nonhomogeneous parabolic equations. We then want to describe how such estimates may be used to prove superconvergence of a number of procedures for evaluating point values of the exact solution and its derivatives. These applications include in particular the case of one space dimension with continuous, piecewise polynomial approximating subspaces, where we analyze methods proposed by Douglas, Dupont and Wheeler [3]. Further, in higher dimensions we discuss the application of an averaging procedure by Bramble and Schatz [11 for elements which are uniform in the interior and in the nonuniform case a method employing a local Green's function considered by Louis and Natterer [4]. The error analysis of this paper takes place in the general framework introduced in Bramble, Schatz, Thome'e and Wahlbin [21 allowing approximating subspaces which do not necessarily satisfy the homogeneous boundary conditions of the exact solution. These subspaces are assumed to permit approximation to order O(h') in L2 (r > 2) and to yield O(h2r2) error estimates for the elliptic problem in norms of order -(r 2). The superconvergent order error estimates which we aim for in the parabolic problem are then of this higher order. In [2], estimates of the type considered here were obtained for homogeneous parabolic equations by spectral representation; our basic results in this paper are derived by the energy method. 1. Preliminaries. We shall be concerned with the approximate solution of the initial-boundary value problem (ut = au/at, R+ = {t; t > O}) Lu-ut +Au-f in Q x R+, ;( .I) u(x, t) O on a x R , u(x, 0) =v(x) on Q2. Received September 12, 1978. AMS (MOS) subject classifications (1970). Primary 65N15, 65N30. ? 1980 American Mathematical Society 0025-571 8/80/0000-0005/$06.25 93 This content downloaded from 157.55.39.58 on Tue, 15 Nov 2016 03:53:30 UTC All use subject to http://about.jstor.org/terms 94 VIDAR THOMEE Here Q2 is a bounded domain in RN with sufficiently smooth boundary M, Au N (aj auk + au, j, k=l 1 i k) with alk and ao sufficiently smooth time-independent functions, the matrix (alk) symmetric and uniformly positive definite and ao nonnegative in Q2. In order to introduce some notation, we consider first the corresponding elliptic problem (1.2) Au =f in 2, u = 0 on a2, and denote by T: L2(Q) Hol) n H2(Q) its solution operator, defined by u = Tf. Notice that by the symmetry of A, T is selfadjoint and positive definite in L2(Q). Recall also the elliptic regularity estimate 1I Tflls+2 6 Cllf 11, for s > 0, where 11 Il denotes the norm in Hs(Q). Set now for s a nonnegative integer and v, w E L2(Q), with (,) the inner product in L2(Q), (1.3) (v, w)_S = (TSv, w), IIvIL_ = (TSv, v)l 2. Since T is positive definite, (, *) is an inner product. One can show that 11 IIis equivalent to the norm

A high-order difference method for differential equations

Abstract: This paper analyzes a high-accuracy approximation to the mth-order linear ordinary differential equation Mu = f. At mesh points, U is the estimate of u; and U satisfies MnU = Inf, where MnU is a linear combination of values of U at m + 1 stencil points (adjacent mesh points) and Inf is a linear combination of values of f at J auxiliary points, which are between the first and last stencil points. The coefficients of Mn, In are obtained "locally" by solving a small linear system for each group of stencil points in order to make the approximation exact on a linear space S of dimension L + 1. For separated two-point boundary value problems, U is the solution of an nby-n linear system with full bandwidth m + 1. For S a space of polynomials, existence and uniqueness are established, and the discretization error is O(h L+l ); the first m 1 divided differences of U tend to those of u at this rate. For a general set of auxiliary points one has L = J + m; but special auxiliary points, which depend upon M and the stencil points, allow larger L, up to L = 2J + m. Comparison of operation counts for this method and five other common schemes shows that this method is among the most efficient for given convergence rate. A brief selection from extensive experiments is presented which supports the theoretical results and the practicality of the method.

A conservative finite element method for the Korteweg-de Vries equation

Abstract: . A finite element method for the 1-periodic Korteweg-de Vries equation "t + 2uux + "xxx = ° is analyzed. We consider first a semidiscrete method (i.e., discretization only in the space variable), and then we analyze some unconditionally stable fully discrete methods.In a special case, the fully discrete methods reduce to twelve point finite differenceschemes (three time levels) which have second order accuracy both in the space andtime variable. 1. Introduction. The purpose of this paper is to study a Galerkin-type methodfor the 1-periodic Korteweg-de Vries equation (1.1) "f + 2""* + "xxx = °. w(x, 0) = u0(x),for 0 0 is a fixed real number. This equation arises for exampleas a model equation for unidirectional long waves in nonlinear dispersive media. For adiscussion of this equation we refer the reader to Whitham [9] and references giventhere.We derive the numerical method by writing the equation (1.1) in the conserva-tive form(1.2) ut-wx = 0,where the flux w is given by

The rate of convergence of Hermite function series

Abstract: Let a > 0 be the least upper bound of -y for which f(z) O(ell'y) for some positive constant q as Izi oon the real axis. It is then proved that at least an infinite subsequence of the coefficients {an} in

Numerical integrators for stiff and highly oscillatory differential equations

Abstract: Some L-stable fourth-order explicit one-step numerical integration formulas which requiare no matrix inversion are proposed to cope effictively with systems of ordinary differential equations with large Lipschitz constants (including those having highly oscillatory solutions). The implicit integration procedure proposed in Fatunla (11) is further developed to handle a larger class of stiff systems as well as those with highly oscillatory solutions. The same pair of nonlinear equations as in (11) is solved for the stiffness/oscillatory parameters. However, the nonlinear systems are transformed into linear forms and an efficient computational procedure is developed to obtain these parameters. The new schemes compare favorably with the backward differentiation formula (DIFSUB) of Gear (13), (14), and the blended linear multistep methods of Skeel and Kong (24), and the symmetric multistep methods of Lambert and Watson (17).

Spectral and pseudospectral methods for advection equations

Abstract: Spectral and pseudo spectral methods for advection equations are investigated. A basic framework is given which allows the application of techniques used in finite element analysis to spectral methods with trigonometric polynomials. Error estimates for semidiscrete spectral and pseudo spectral as well as fully discrete explicit pseudo spectral methods are given. The approximation schemes are shown to converge with infinite order.

Analysis of convergence of the $T$-transformation for power series

Abstract: Recently the present author has given some convergence theorems of general nature for Levin's nonlinear sequence transformations. In this work these theorems are extended and sharpened to cover the case of power series, both inside and on their circle of convergence. It is shown that one of the two limiting processes considered in the previous work can be used for analytic continuation and a realistic estimate of its rate of convergence is given. Three illustrative examples are also appended.