Showing papers in "Mathematics of Computation in 1983"

Numerical approximations to nonlinear conservation laws with locally varying time and space grids

Abstract: Numerical approximations to the initial value problem for nonlinear systems of conservation laws are considered The considered system is said to be hyperbolic when all eigenvalues of every real linear combination of the Jacobian matrices are real Solutions may develop discontinuities in finite time, even when the initial data are smooth In the investigation, explicit finite difference methods which use locally varying time grids are considered The global CFL restriction is replaced by a local restriction The numerical flux function is studied from a finite volume viewpoint, and a differencing technique is developed at interface points between regions of distinct time increments

On the asymptotic convergence of collocation methods

Abstract: We prove quasioptimal and optimal order estimates in various Sobolev norms for the approximation of linear strongly elliptic pseudodifferential equations in one independent variable by the method of nodal collocation by odd degree polynomial splines. The analysis pertains in particular to many of the boundary element methods used for numerical computation in engineering applications. Equations to which the analysis is applied include Fredholm integral equations of the second kind, certain first kind Fredholm equations, singular integral equations involving Cauchy kernels, a variety of integro-differential equations, and two-point boundary value problems for ordinary differential equations. The error analysis is based on an equivalence which we establish between the collocation methods and certain nonstandard Galerkin methods. We compare the collocation method with a standard Galerkin method using splines of the same degree, showing that the Galerkin method is quasioptimal in a Sobolev space of lower index and furnishes optimal order approximation for a range of Sobolev indices containing and extending below that for the collocation method, and so the standard Galerkin method achieves higher rates of convergence.

Preconditioning and two-level multigrid methods of arbitrary degree of approximation

Abstract: Let h be a mesh parameter corresponding to a finite element mesh for an elliptic problem. We describe preconditioning methods for two-level meshes which, for most prob- lems solved in practice, behave as methods of optimal order in both storage and computa- tional complexity. Namely, per mesh point, these numbers are bounded above by relatively small constants for all h - ho, where ho is small enough to cover all but excessively fine meshes. We note that, in practice, multigrid methods are actually solved on a finite, often even a fixed number of grid levels, in which case also these methods are not asymptotically optimal as h - 0. Numerical tests indicate that the new methods are about as fast as the best implementations of multigrid methods applied on general problems (variable coefficients, general domains and boundary conditions) for all but excessively fine meshes. Furthermore, most of the latter methods have been implemented only for difference schemes of second order of accuracy, whereas our methods are applicable to higher order approximations. We claim that our scheme could be added fairly easily to many existing finite element codes. 1. Introduction. Consider the numerical solution of elliptic boundary value prob- lems discretized by finite element methods. We assume that the boundary is polygonal or consists of planes. We note that in practical problems one often has a fine enough grid already after the definition of the boundary and the minimal number of vertices needed for a first (coarse) triangulation. Anyhow, if not so, in most cases one makes only a few steps of mesh refinement. Hence the power of multigrid methods-their optimal order of computational complexity-is most often not achieved fully, because optimality requires a large number of recursively defined meshes (for details see, e.g., (4) and for further references see (7)). Hence one might as well consider other methods, perhaps simpler and more effective on a fixed mesh, but which are not asymptotically optimal. Here we shall describe a method which uses only a fixed mesh, but for which one nevertheless achieves a low order of computational complexity and of seemingly optimal order except for, from a practical viewpoint, excessively small meshes. To be more precise, the computational cost per mesh point is bounded by clog N for N < No, where N is the number of mesh points, No is large enough to cover most applications and c is small enough that the method is competitive with multigrid methods. As is well known, the latter need recursion and the usual smoothing followed by corrections of the solutions on the different mesh levels. We claim that the new method is more suitable for implementation in existing finite element packages. In fact most packages for the multigrid methods are only for second order difference methods.

Some extensions of W. Gautschi’s inequalities for the gamma function

Abstract: It has been shown by W. Gautschi that if 0 I Xi-s < F(x ) < exp[(I s)x + 1)]. The following closer bounds are proved: exp[(I s)4(x + 12)] < F + ) < exp[(I s) (x + s I)] F(x ? s)2 and [x + 2]

On convergence of monotone finite difference schemes with variable spatial differencing

Abstract: Monotone finite difference schemes used to approximate solutions of scalar conservation laws have the advantage that these approximations can be proved to converge to the proper solution as the mesh size tends to zero. The greatest disadvantage in using such approximating schemes is the computational expense encountered since monotone schemes can have at best first order accuracy. Computation savings and effective accuracy could be gained if the spatial mesh were refined in regions of expected rapid solution variation. In this paper we prove that standard monotone difference schemes, (satisfying a fairly unrestrictive CFL condition), converge to the "correct" physical solution even in the case when a nonuniform spatial mesh is employed.

Numerical methods for a model for compressible miscible displacement in porous media

Abstract: A nonlinear parabolic system is derived to describe compressible miscible displace- ment in a porous medium. The system is consistent with the usual model for incompressible miscible displacement. Two finite element procedures are introduced to approximate the concentration of one of the fluids and the pressure of the mixture. The concentration is treated by a Galerkin method in both procedures. while the pressure is treated by either a (lalerkin method or by a parabolic mixed finite element method. Optimal order estimates in L2 and essentially optimal order estimates in Lx are derived for the errors in the approximate solutions for both methods. Introduction. We shall consider the single-phase, miscible displacement of one compressible fluid by another in a porous medium under the assumptions that no volume change results from the mixing of the components and that a pressure-den- sity relation exists for each component in a form that is independent of the mixing. These equations of state will imply that the fluids are in the liquid state. Our model will represent a direct generalization of the model (3), (4), (7) that has been treated extensively for incompressible miscible displacement. The reservoir S will be taken to be of unit thickness and will be identified with a bounded domain in R2. We shall omit gravitational terms for simplicity of exposi- tion; no significant mathematical questions arise when the lower order terms are included.

A method for interpolating scattered data based upon a minimum norm network

Abstract: A method for interpolating scattered data is described. CTiven (xi, yi, zi), i 1. N, a bivariate function S with continuous first order partial derivatives is defined which has the property that S(xi, Yi) zi, i 1,...,N.-The method is based upon a triangulation of the domain and a curve network which has certain minimum pseudonorm properties. Algorithms and examples are included.

Runge-Kutta theory for Volterra and Abel integral equations of the second kind

Abstract: The present paper develops the local theory of general Runge-Kutta methods for a broad class of weakly singular and regular Volterra integral equations of the second kind. Further, the smoothness properties of the exact solutions of such equations are investigated.

On the finite element method for singularly perturbed reaction-diffusion problems in two and one dimensions

Abstract: Second order elliptic boundary value problems which are allowed to degenerate into zero order equations are considered. The behavior of the ordinary Galerkin finite element method without special arrangements to treat singularities is studied as the problem ranges from true second order to singularly perturbed.

Error estimates for the numerical identification of a variable coefficient

Abstract: Error estimates are derived for the approximate identification of an unknown transmissivity coefficient in a partial differential equation describing a model problem in groundwater now. The approximation scheme considered determines the coefficient by least squares fitting of the observed pressure data. 1. Introduction. In this paper we wish to present an error analysis of a common numerical scheme used in the identification of parameters in distributed systems. Specifically, we shall concern ourselves here with a model problem in groundwater flow. The problem is to identify a spatially varying transmissivity coefficient a(x) from observations of the piezometric head u(x) in a two-dimensional static aquifer s, where a and u are related by the equations (1) -div(avw) =/ inS2,

Numerical Methods for Stiff Equations and Singular Perturbation Problems.

Abstract: 1. Introduction.- Summary.- 1.1. Stiffness and Singular Perturbations.- 1.1.1. Motivation.- 1.1.2. Stiffness.- 1.1.3. Singular Perturbations.- 1.1.4. Applications.- 1.2. Review of the Classical Linear Multistep Theory.- 1.2.1. Motivation.- 1.2.2. The Initial Value Problem.- 1.2.3. Linear Multistep Operators.- 1.2.4. Approximate Solutions.- 1.2.5. Examples of Linear Multistep Methods.- 1.2.6. Stability, Consistency and Convergence.- 2. Methods of Absolute Stability.- Summary.- 2.1. Stiff Systems and A-stability.- 2.1.1. Motivation.- 2.1.2. A-stability.- 2.1.3. Examples of A-stable Methods.- 2.1.4. Properties of A-stable Methods.- 2.1.5. A Sufficient Condition for A-stability.- 2.1.6. Applications.- 2.2. Notions of Diminished Absolute Stability.- 2.2.1. A (?)-stability.- 2.2.2. Properties of A(?)-stable Methods.- 2.2.3. Stiff Stability.- 2.3. Solution of the Associated Equations.- 2.3.1. The Problem.- 2.3.2. Conjugate Gradients and Dichotomy.- 2.3.3. Computational Experiments.- 3. Nonlinear Methods.- Summary.- 3.1. Interpolatory Methods.- 3.1.1. Certaine's Method.- 3.1.2. Jain's Method.- 3.2. Runge-Kutta Methods and Rosenbrock Methods.- 3.2.1. Runge-Kutta Methods with v-levels.- 3.2.2. Determination of the Coefficients.- 3.2.3. An Example.- 3.2.4. Semi-explicit Processes and the Method of Rosenbrock.- 3.2.5. A-stability.- 4 Exponential Fitting.- Summary.- 4.1. Exponential Fitting for Linear Multistep Methods.- 4.1.1. Motivation and Examples.- 4.1.2. Minimax fitting.- 4.1.3. An Error Analysis for an Exponentially Fitted F1.- 4.2. Fitting in the Matricial Case.- 4.2.1. The Matricial Multistep Method.- 4.2.2. The Error Equation.- 4.2.3. Solution of the Error Equation.- 4.2.4. Estimate of the Global Error.- 4.2.5. Specification of P.- 4.2.6. Specification of L and R.- 4.2.7. An Example.- 4.3. Exponential Fitting in the Oscillatory Case.- 4.3.1. Failure of the Previous Methods.- 4.3.2. Aliasing.- 4.3.3. An Example of Aliasing.- 4.3.4. Application to Highly Oscillatory Systems.- 4.4. Fitting in the Case of Partial Differential Equations.- 4.4.1. The Problem Treated.- 4.4.2. The Minimization Problem.- 4.4.3. Highly Oscillatory Data.- 4.4.4. Systems.- 4.4.5. Discontinuous Data.- 4.4.6. Computational Experiments.- 5. Methods of Boundary Layer Type.- Summary.- 5.1. The Boundary Layer Numerical Method.- 5.1.1. The Boundary Layer Formalism.- 5.1.2. The Numerical Method.- 5.1.3. An Example.- 5.2. The ?-independent Method.- 5.2.1. Derivation of the Method.- 5.2.2. Computational Experiments.- 5.3. The Extrapolation Method.- 5.3.1. Derivation of the Relaxed Equations.- 5.3.2. Computational Experiments.- 6. The Highly Oscillatory Problem.- Summary.- 6.1. A Two-time Method for the Oscillatory Problem.- 6.1.1. The Model Problem.- 6.1.2. Numerical Solution Concept.- 6.1.3. The Two-time Expansion.- 6.1.4. Formal Expansion Procedure.- 6.1.5. Existence of the Averages and Estimates of the Remainder.- 6.1.6. The Numerical Algorithm.- 6.1.7. Computational Experiments.- 6.2. Algebraic Methods for the Averaging Process.- 6.2.1. Algebraic Characterization of Averaging.- 6.2.2. An Example.- 6.2.3. Preconditioning.- 6.3. Accelerated Computation of Averages and an Extrapolation Method.- 6.3.1. The Multi-time Expansion in the Nonlinear Case.- 6.3.2. Accelerated Computation of $$\bar f$$.- 6.3.3. The Extrapolation Method.- 6.3.4. Computational Experiments: A Linear System.- 6.3.5. Discussion.- 6.4. A Method of Averaging.- 6.4.1. Motivation: Stable Functionals.- 6.4.2. The Problem Treated.- 6.4.3. Choice of Functionals.- 6.4.4. Representers.- 6.4.5. Local Error and Generalized Moment Conditions.- 6.4.6. Stability and Global Error Analysis.- 6.4.7. Examples.- 6.4.8. Computational Experiments.- 4.6.9. The Nonlinear Case and the Case of Systems.- 7. Other Singularly Perturbed Problems.- Summary.- 7.1. Singularly Perturbed Recurrences.- 7.1.1. Introduction and Motivation.- 7.1.2. The Two-time Formalism for Recurrences.- 7.1.3. The Averaging Procedure.- 7.1.4. The Linear Case.- 7.1.5. Additional Applications.- 7.2. Singularly Perturbed Boundary Value Problems.- 7.2.1. Introduction.- 7.2.2. Numerically Exploitable Form of the Connection Theory.- 7.2.3. Description of the Algorithm.- 7.2.4. Computational Experiments.- References.

The numerical solution of equality constrained quadratic programming problems

Abstract: This paper proves that a large class of iterative schemes can be used to solve a certain constrained minimization problem. The constrained minimization problem considered involves the minimization of a quadratic functional subject to linear equality constraints. Among this class of convergent iterative schemes are generalizations of the relaxed Jacobi, Gauss-Seidel, and symmetric Gauss-Seidel schemes.

Applications of a computer implementation of Poincaré’s theorem on fundamental polyhedra

Abstract: PoincarCs Theorem asserts that a group F of isometries of hyperbolic space H is discrete if its generators act suitably on the boundary of some polyhedron in H, and when this happens a presentation of F can be derived from this action. We explain methods for deducing the precise hypotheses of the theorem from calculation in F when F is "algorithmi- cally defined", and we describe a file of Fortran programs that use these methods for groups F acting on the upper half space model of hyperbolic 3-space H. We exhibit one modest example of the application of these programs, and we summarize computations of repesenta- tions of groups PSL(2, 0) where 6 is an order in a complex quadratic number field. In the early 1880's H. Poincare discovered a general theorem allowing one to deduce the discreteness of, and a presentation for, a group G of isometries of hyperbolic space from its action on a hyperbolic polyhedron under certain condi- tions. Theorems of this sort are part of the foundations of his theories of Fuchsian and Kleinian groups that have become very popular again, and H. Seifert has recently given us a modern proof of a fairly general version of Poincare's Theorem in (12); see also (7). This theorem has been little used over the past century, perhaps partly because its hypotheses have seemed very difficult to verify for a given group G except in very special circumstances. One reason for doubting that Poincare's Theorem is unreasonably difficult to apply to fairly general discrete groups is that no general alternative method for accomplishing its tasks has been proposed. The present paper is devoted to demonstrating that Poincare's Theorem can indeed be applied to given groups in apparently difficult cases and that much of the work can be done by a computer. Our experience suggests that the theorem is really very helpful in guiding the user to an understanding of the details of the action of G starting from a state of near ignorance. An outline of this paper is as follows. In Section 1 we begin by stating Seifert's version of Poincare's Theorem and explaining how we would apply it to an " algorithmically defined" group G of isometries of hyperbolic space H 'o. In Section 2 we specialize to the situation, 00 , for which we wrote our file, Poincare', of Fortran

Computation of Faber series with application to numerical polynomial approximation in the complex plane

Abstract: Kovari and Pommerenke [19], and Elliott [8], have shown that the truncated Faber series gives a polynomial approximation which (for practical values of the degree of the polynomial) is very close to the best approximation. In this paper we discuss efficient Fast Fourier Transform (FFT) and recursive methods for the computation of Faber polynomials, and point out that the FFT method described by Geddes [13], for computing Chebyshev coefficients can be generalized to compute Faber coefficients. We also give a corrected bound for the norm of the Faber projection (that given in Elliott [8], being unfortunately slightly in error) and very briefly discuss a possible extension of the method to the case when the mapping function, which is required to compute the Faber series, is not known explicitly.

The determination of the value of Rado’s noncomputable function Σ() for four-state Turing machines

Abstract: The well-defined but noncomputable functions E(k) and S(k) given by T. Rado as the "score" and "shift number" for the k-state Turing machine "Busy Beaver Game" were previously known only for k 13 and S(4) a 107, reported by this author, supported the conjecture that these lower bounds are the actual particular values of the functions for k 4. The four-state case has previously been reduced to solving the blank input tape halting problem of only 5,820 individual machines. In this final stage of the k = 4 case, one appears to move into a heuristic level of higher order where it is necessary to treat each machine as representing a distinct theorem. The remaining set consists of two primary classes in which a machine and its tape are viewed as the representation of a growing string of cellular automata. The proof techniques, embodied in programs, are entirely heuristic, while the inductive proofs, once established by the computer, are completely rigorous and become the key to the proof of the new and original mathematical results: l:(4) 13 and S(4) = 107. " In any case, even though skilled mathematicians and experienced programmers attempted to evaluate E(3) and S(3), there is no evidence that any known approach will yield the answer, even if we avail ourselves of high-speed computers and elaborate programs. As regards L(4), S(4), the situation seems to be entirely hopeless at present."-Tibor Rado, 1963. 1. Background and Introduction. The "Busy Beaver Game" was devised by Tibor Rado [8] for the purpose of illustrating the notion of noncomputability. Given the set of Turing machines of exactly k states which operate with the minimum alphabet of two symbols (a space and a mark or 0 and 1) one considers the problem of behavior of these machines on a tape which is initially all blank (all O's). This is a finite set of machines, there being exactly (4k + 1)2k distinct machines, where, with two table entries per state, each table entry may consist of printing 0 or 1, moving right or left, branching to state 1, 2,... ,k or simply an entry to declare a halt.* Since the input tape is blank, each machine faces one of two possible fates: either it eventually halts, or else it continues running forever. After a particular machine Received July 22, 1981; revised June 28, 1982 and September 7, 1982. 1980 Mathematics Subject Classification. Primary 03D10, 03B35; Secondary 68C30, 68D20.

Additive Runge-Kutta methods for stiff ordinary differential equations

Abstract: Certain pairs of Runge-Kutta methods may be used additively to solve a system of n differential equations x' = J(t)x + g(t, x). Pairs of methods, of order p < 4, where one method is semiexplicit and A-stable and the other method is explicit, are obtained. These methods require the LU factorization of one n X n matrix, and p evaluations of g, in each step. It is shown that such methods have a stability property which is similar to a stability property of perturbed linear differential equations. 1. Introduction. In a recent article (2) the authors showed that certain pairs of methods may be used in an additive fashion to solve an initial value problem for a system of n differential equations x' = f(t, x), x(a) = xo, a - t - b,

Divisors of Mersenne numbers

Abstract: We add to the heuristic and empirical evidence for a conjecture of Gillies about the distribution of the prime divisors of Mersenne numbers We list some large prime divisors of Mersenne numbers Mp in the range 17000

Convergence of Galerkin Approximations for the Korteweg-de Vries Equation,

Abstract: Standard Galerkin approximations, using smooth splines on a uniform mesh, to 1-periodic solutions of the Korteweg-de Vries equation are analyzed. Optimal rate of convergence estimates are obtained for both semidiscrete and second order in time fully discrete schemes. At each time level, the resulting system of nonlinear equations can be solved by Newton's method. It is shown that if a proper extrapolation is used as a starting value, then only one step of the Newton iteration is required.

A general orthogonalization technique with applications to time series analysis and signal processing

Abstract: A new orthogonalization technique is presented for computing the QR factorization of a general n X p matrix of full rank p (n 2 p). The method is based on the use of projections to solve increasingly larger subproblems recursively and has an O(np2) operation count for general matrices. The technique is readily adaptable to solving linear least-squares problems. If the initial matrix has a circulant structure the algorithm simplifies significantly and gives the so-called lattice algorithm for solving linear prediction problems. From this point of view it is seen that the lattice algorithm is really an efficient way of solving specially structured least-squares problems by orthogonalization as opposed to solving the normal equations by fast Toeplitz algorithms.

Numerical methods for flows through porous media. I

Abstract: The degenerate parabolic equation aut = v ( Oulvvu), v >1 has been used to model the flow of gas through a porous medium. Error estimates for continuous and discrete time finite element procedures to approximate the solution of this equation are proved, and several new regularity results are given. 1. A Porous Medium Equation. Introduction. We shall study the porous medium equation (1.1) au/at = V (I u I"vu) on 2 X (0, T], (1.2) au/an = 0 on aa X [O, T], (1.3) u(x, 0)= uO(x) on 2, where v > 1 is a parameter and Q is a bounded domain in RN, N ? 3, with a smooth boundary. The initial function u0 is assumed to be nonnegative and four times continuously differentiable on U. Notice that the compatibility condition auo/an = 0 holds on asa. Our main result is the derivation of error estimates for numerical approximations to the problem (1. 1)-(1.3), which we shall refer to as "the porous medium equation" or "IPME". The PME does not, in general, admit classical solutions. Existence and uniqueness of weak solutions was proved in one space dimension by Oleinik, Kalashnikov, and Czou [15], [16] and in several space dimensions by Lions [12]. These proofs concern the PME with different boundary conditions, but the arguments carry over to the PME (1.1)-(1.3). The maximum principle implies that, since u0 is nonnegative on 92, u(x, t) is nonnegative for all (x, t) E Q X [0, T]; see [15], [16]. If u0 is nonzero, the Neumann boundary condition implies that u will eventually become strictly positive and (1.1) will become nondegenerate for all time t > To, To sufficiently large. We can rewrite (1.1) in the form (1.4) au/at=AK(u) on QX(0,T], Received May 11, 1979; revised July 24, 1981. 1980 Mathematics Subject Classification. Primary 35K65, 65N30, 65N15. *Current address: Reservoir Engineering Branch, Petroleum Engineering Division, BP Exploration (DOS), Britannic House, Moore Lane, London EC2Y 9BU, United Kingdom. (D1983 American Mathematical Society 0025-5718/82/0000-1 063/$08.25 435 This content downloaded from 157.55.39.55 on Tue, 23 Aug 2016 04:12:54 UTC All use subject to http://about.jstor.org/terms

A spline-trigonometric Galerkin method and an exponentially convergent boundary integral method

Abstract: We consider a Galerkin method for functional equations in one space variable which uses periodic cardinal splines as trial functions and trigonometric polynomials as test functions. We analyze the method applied to the integral equation of the first kind arising from a single layer potential formulation of the Dirichlet problem in the interior or exterior of an analytic plane curve. In constrast to ordinary spline Galerkin methods, we show that the method is stable, and so provides quasioptimal approximation, in a large family of Hilbert spaces including all the Sobolev spaces of negative order. As a consequence we prove that the approximate solution to the Dirichlet problem and all its derivatives converge pointwise with exponential rate.