Showing papers in "Mathematics of Computation in 1997"

On some inequalities for the gamma and psi functions

Abstract: We present new inequalities for the gamma and psi functions, and we provide new classes of completely monotonic, star-shaped, and superadditive functions which are related to Γ and?.

On some inequalities for the incomplete gamma function

Abstract: Let p ¬= 1 be a positive real number. We determine all real numbers α = α(p) and β = β(p) such that the inequalities formula math. formula math. are valid for all x > 0. And, we determine all real numbers a and b such that - log(1 - e -ax ) ≤ √ x ∞ e-t/t ≤ - log(1 - e -bx ) hold for all > 0.

A globally convergent Lagrangian barrier algorithm for optimization with general inequality constraints and simple bounds

Abstract: We consider the global and local convergence properties of a class of Lagrangian barrier methods for solving nonlinear programming problems. In such methods, simple bound constraints may be treated separately from more general constraints. The objective and general constraint functions are combined in a Lagrangian barrier function. A sequence of such functions are approximately minimized within the domain defined by the simple bounds. Global convergence of the sequence of generated iterates to a first-order stationary point for the original problem is established. Furthermore, possible numerical difficulties associated with barrier function methods are avoided as it is shown that a potentially troublesome penalty parameter is bounded away from zero. This paper is a companion to previous work of ours on augmented Lagrangian methods.

On the rapid computation of various polylogarithmic constants

Abstract: We give algorithms for the computation of the d-th digit of certain transcendental numbers in various bases. These algorithms can be easily implemented (multiple precision arithmetic is not needed), require virtually no memory, and feature run times that scale nearly linearly with the order of the digit desired. They make it feasible to compute, for example, the billionth binary digit of log(2) or π on a modest work station in a few hours run time.

A posteriori error estimate for the mixed finite element method

Abstract: A computable error bound for mixed finite element methods is established in the model case of the Poisson-problem to control the error in the H(div,Ω) x L 2 (Ω)-norm. The reliable and efficient a posteriori error estimate applies, e.g., to Raviart-Thomas, Brezzi-Douglas-Marini, and Brezzi-Douglas-Fortin-Marini elements.

A stochastic particle method for the McKean-Vlasov and the Burgers equation

Abstract: In this paper we introduce and analyze a stochastic particle method for the McKean-Vlasov and the Burgers equation; the construction and error analysis are based upon the theory of the propagation of chaos for interacting particle systems. Our objective is three-fold. First, we consider a McKean-Vlasov equation in [0,T] x R with sufficiently smooth kernels, and the PDEs giving the distribution function and the density of the measure μ t , the solution to the McKean-Vlasov equation. The simulation of the stochastic system with N particles provides a discrete measure which approximates μkΔt for each time kΔt (where Δt is a discretization step of the time interval [0, T]). An integration (resp. smoothing) of this discrete measure provides approximations of the distribution function (resp. density) of μkΔt. We show that the convergence rate is O (1/√N+ for the approximation in L i (Ω x R) of the cumulative distribution function at time T, and of order O (e 2 +1/e) (1/√N + √Δt for the approximation in L l (Ω x R) of the density at time T (Ω is the underlying probability space, e is a smoothing parameter). Our second objective is to show that our particle method can be modified to solve the Burgers equation with a nonmonotonic initial condition, without modifying the convergence rate O (1/√N+√Δt). This part extends earlier work of ours, where we have limited ourselves to monotonic initial conditions. Finally, we present numerical experiments which confirm our theoretical estimates and illustrate the numerical efficiency of the method when the viscosity coefficient is very small.

A least-squares approach based on a discrete minus one inner product for first order systems

Abstract: The purpose of this paper is to develop and analyze a least-squares approximation to a first order system. The first order system represents a reformulation of a second order elliptic boundary value problem which may be indefinite and/or nonsymmetric. The approach taken here is novel in that the least-squares functional employed involves a discrete inner product which is related to the inner product in H -1 (Ω) (the Sobolev space of order minus one on Ω). The use of this inner product results in a method of approximation which is optimal with respect to the required regularity as well as the order of approximation even when applied to problems with low regularity solutions. In addition, the discrete system of equations which needs to be solved in order to compute the resulting approximation is easily preconditioned, thus providing an efficient method for solving the algebraic equations. The preconditioner for this discrete system only requires the construction of preconditioners for standard second order problems, a task which is well understood.

Preconditioning in H (div) and applications

Abstract: We consider the solution of the system of linear algebraic equations which arises from the finite element discretization of boundary value problems associated to the differential operator I-grad div. The natural setting for such problems is in the Hilbert space H(div) and the variational formulation is based on the inner product in H(div). We show how to construct preconditioners for these equations using both domain decomposition and multigrid techniques. These preconditioners are shown to be spectrally equivalent to the inverse of the operator. As a consequence, they may be used to precondition iterative methods so that any given error reduction may be achieved in a finite number of iterations, with the number independent of the mesh discretization. We describe applications of these results to the efficient solution of mixed and least squares finite element approximations of elliptic boundary value problems.

Calculation of Gauss-Kronrod quadrature rules

Abstract: The Jacobi matrix of the (2n+1)-point Gauss-Kronrod quadrature rule for a given measure is calculated efficiently by a five-term recurrence relation. The algorithm uses only rational operations and is therefore also useful for obtaining the Jacobi-Kronrod matrix analytically. The nodes and weights can then be computed directly by standard software for Gaussian quadrature formulas.

Generators and irreducible polynomials over finite fields

Abstract: Weil's character sum estimate is used to study the problem of constructing generators for the multiplicative group of a finite field. An application to the distribution of irreducible polynomials is given, which confirms an asymptotic version of a conjecture of Hansen-Mullen.

A search for Wieferich and Wilson primes

Abstract: An odd prime p is called a Wieferich prime if 2 P-1 = 1 (mod p 2 ) alternatively, a Wilson prime if (p - 1)|= -1 (mod p 2 ). To date, the only known Wieferich primes are p = 1093 and 3511, while the only known Wilson primes are p = 5,13, and 563. We report that there exist no new Wieferich primes p < 4 x 10 12 , and no new Wilson primes p < 5x 10 8 . It is elementary that both defining congruences above hold merely (mod p), and it is sometimes estimated on heuristic grounds that the probability that p is Wieferich (independently: that p is Wilson) is about 1/p. We provide some statistical data relevant to occurrences of small values of the pertinent Fermat and Wilson quotients (mod p).

Numerical solution of isospectral flows

Abstract: In this paper we are concerned with the problem of solving numerically isospectral flows. These flows are characterized by the differential equation L' = [B(L), L]; L(0) = L 0 , where L 0 is a d x d symmetric matrix, B(L) is a skew-symmetric matrix function of L and [B, L] is the Lie bracket operator. We show that standard Runge-Kutta schemes fail in recovering the main qualitative feature of these flows, that is isospectrality, since they cannot recover arbitrary cubic conservation laws. This failure motivates us to introduce an alternative approach and establish a framework for generation of isospectral methods of arbitrarily high order.

Estimating the number of asymptotic degrees of freedom for nonlinear dissipative systems

Abstract: We show that the long-time behavior of the projection of the exact solutions to the Navier-Stokes equations and other dissipative evolution equations on the finite-dimensional space of interpolant polynomials determines the long-time behavior of the solution itself provided that the spatial mesh is fine enough. We also provide an explicit estimate on the size of the mesh. Moreover, we show that if the evolution equation has an inertial manifold, then the dynamics of the evolution equation is equivalent to the dynamics of the projection of the solutions on the finite-dimensional space spanned by the approximating polynomials. Our results suggest that certain numerical schemes may capture the essential dynamics of the underlying evolution equation.

Analysis and convergence of a covolume method for the generalized Stokes problem

Abstract: We introduce a covolume or MAC-like method for approximating the generalized Stokes problem. Two grids are needed in the discretization; a triangular one for the contimlity equation and a quadrilateral one for the momentum equation. The velocity is approximated using nonconforming piecewise linears and the pressure piecewise constants. Error in the L 2 norm for the pressure and error in a mesh dependent H 1 norm as well as in the L 2 norm for the velocity are shown to be of first order, provided that the exact velocity is in H 2 and the true pressure in H I . We also introduce the concept of a network model into the discretized linear system so that an efficient pressure-recovering technique can be used to simplify a great deal the computational work involved in the augmented Lagrangian method. Given is a very general decomposition condition under which this technique is applicable to other fluid problems that can be formulated as a saddle-point problem.

Composition constants for raising the orders of unconventional schemes for ordinary differential equations

Abstract: Many models of physical and chemical processes give rise to ordinary differential equations with special structural properties that go unexploited by general-purpose software designed to solve numerically a wide range of differential equations. If those properties are to be exploited fully for the sake of better numerical stability, accuracy and/or speed, the differential equations may have to be solved by unconventional methods. This short paper is to publish composition constants obtained by the authors to increase efficiency of a family of mostly unconventional methods, called reflexive.

Numerical solution of the scalar double-well problem allowing microstructure

Abstract: The direct numerical solution of a non-convex variational problem (P) typically faces the difficulty of the finite element approximation of rapid oscillations. Although the oscillatory discrete minimisers are properly related to corresponding Young measures and describe real physical phenomena, they are costly and difficult to compute. In this work, we treat the scalar double-well problem by numerical solution of the relaxed problem (RP) leading to a (degenerate) convex minimisation problem. The problem (RP) has a minimiser u and a related stress field σ = DW**(⊇u) which is known to coincide with the stress field obtained by solving (P) in a generalised sense involving Young measures. If u h is a finite element solution, σ h := DW**(⊇u h ) is the related discrete stress field. We prove a priori and a posteriori estimates for σ-σ h in L 4/3 (Ω) and weaker weighted estimates for ⊇u - ⊇u h . The a posteriori estimate indicates an adaptive scheme for automatic mesh refinements as illustrated in numerical experiments.

An a posteriori error estimate for a first-kind integral equation

Abstract: In this paper we present a new a posteriori error estimate for the boundary element method applied to an integral equation of the first kind. The estimate is local and sharp for quasi-uniform meshes and so improves earlier work of ours. The mesh-dependence of the constants is analyzed and shown to be weaker than expected from our previous work. Besides the Galerkin boundary element method, the collocation method and the qualocation method are considered. A numerical example is given involving an adaptive feedback algorithm.

A fast algorithm to compute cubic fields

Abstract: We present a very fast algorithm to build up tables of cubic fields. Real cubic fields with discriminant up to 10 11 and complex cubic fields down to -10 11 have been computed.

Locking-free finite element methods for shells

Abstract: We propose a new family of finite element methods for the Naghdi shell model, one method associated with each nonnegative integer k. The methods are based on a nonstandard mixed formulation, and the kth method employs triangular Lagrange finite elements of degree k+2 augmented by bubble functions of degree k + 3 for both the displacement and rotation variables, and discontinuous piecewise polynomials of degree k for the shear and membrane stresses. This method can be implemented in terms of the displacement and rotation variables alone; as the minimization of an altered energy functional over the space mentioned. The alteration consists of the introduction of a weighted local projection into part, but not all, of the shear and membrane energy terms of the usual Naghdi energy. The relative error in the method, measured in a norm which combines the H I norm of the displacement and ro tation fields and an appropriate norm of the shear and membrane stress fields, converges to zero with order k+1 uniformly with respect to the shell thickness for smooth solutions, at least under the assumption that certain geometrical coefficients in the Nagdhi model are replaced by piecewise constants.

How accurate is the streamline diffusion finite element method

Abstract: We investigate the optimal accuracy of the streamline diffusion finite element method applied to convection-dominated problems. For linear/bilinear elements the theoretical order of convergence given in the literature is either O(h 3/2 ) for quasi-uniform meshes or O(h 2 ) for some uniform meshes. The determination of the optimal order in general was an open problem. By studying a special type of meshes, it is shown that the streamline diffusion method may actually converge with any order within this range depending on the characterization of the meshes.

Optimal, quasi-optimal and superlinear band-Toeplitz preconditioners for asymptotically ill-conditioned positive definite Toeplitz systems

Abstract: In this paper we are concerned with the solution of n x n Hermitian Toeplitz systems with nonnegative generating functions f. The preconditioned conjugate gradient (PCG) method with the well-known circulant preconditioners fails in the case where f has zeros. In this paper we consider as preconditioners band-Toeplitz matrices generated by trigonometric polynomials g of fixed degree l. We use different strategies of approximation of f to devise a polynomial g which has some analytical properties of f, is easily computable and is such that the corresponding preconditioned system has a condition number bounded by a constant independent of n. For each strategy we analyze the cost per iteration and the number of iterations required for the convergence within a preassigned accuracy. We obtain different estimates of l for which the total cost of the proposed PCG methods is optimal and the related rates of convergence are superlinear. Finally, for the most economical strategy, we perform various numerical experiments which fully confirm the effectiveness of approximation theory tools in the solution of this kind of linear algebra problems.

On the singular values of Weber modular functions

Abstract: The minimal polynomials of the singular values of the classical Weber modular functions give far simpler defining polynomials for the class fields of imaginary quadratic fields than the minimal polynomials of singular moduli of level 1. We describe computations of these polynomials and give conjectural formulas describing the prime decomposition of their resultants and discriminants, extending the formulas of Gross-Zagier for the level 1 case.

A general framework for high-accuracy parametric interpolation

Abstract: In this paper we establish a general framework for so-called parametric, polynomial, interpolation methods for parametric curves. In contrast to traditional methods, which typically approximate the components of the curve separately, parametric methods utilize geometric information (which depends on all the components) about the curve to generate the interpolant. The general framework suggests a multitude of interpolation methods in all space dimensions, and some of these have been studied by other authors as independent methods of approximation. Since the approximation methods are nonlinear, questions of solvability and stability have to be considered. As a special case of a general result, we prove that four points on a planar curve can be interpolated by a quadratic with fourth-order accuracy, if the points are sufficiently close to a point with nonvanishing curvature. We also find that six points on a planar curve can be interpolated by a cubic, with sixth-order accuracy, provided the points are sufficiently close to a point where the curvature does not have a double zero. In space it turns out that five points sufficiently close to a point with nonvanishing torsion can be interpolated by a cubic, with fifth-order accuracy.

New estimates for Ritz vectors

Abstract: The following estimate for the Rayleigh-Ritz method is proved: |λ-λ||(u, u)| ≤:|Au - λu|| sin <{u; U}, ||u|| = 1. Here A is a bounded self-adjoint operator in a real Hilbert/euclidian space, {λ, u} one of its eigenpairs, U a trial subspace for the Rayleigh-Ritz method, and {λ, u} a Ritz pair. This inequality makes it possible to analyze the fine structure of the error of the Rayleigh-Ritz method, in particular, it shows that |(u, u)| ≤ C ∈ 2 , if an eigenvector u is close to the trial subspace with accuracy ∈ and a Ritz vector u is an ∈ approximation to another eigenvector, with a different eigenvalue. Generalizations of the estimate to the cases of eigenspaces and invariant subspaces are suggested, and estimates of approximation of eigenspaces and invariant subspaces are proved.

A subspace limited memory Quasi-Newton algorithm for large-scale nonlinear bound constrained optimization

Abstract: In this paper we propose a subspace limited memory quasi-Newton method for solving large-scale optimization with simple bounds on the variables. The limited memory quasi-Newton method is used to update the variables with indices outside of the active set, while the projected gradient method is used to update the active variables. The search direction consists of three parts: a subspace quasi-Newton direction, and two subspace gradient and modified gradient directions. Our algorithm can be applied to large-scale problems as there is no need to solve any subproblems. The global convergence of the method is proved and some numerical results are also given.

On the p -Divisibility of Fermat quotients

Abstract: The authors carried out a numerical search for Fermat quotients Q a = (a p-1 - 1)/p vanishing mod p, for 1 ≤ a ≤ p - 1, up top < 10 6 . This article reports on the results and surveys the associated theoretical properties of Q a . The approach of fixing the prime p rather than the base a leads to some aspects of the theory apparently not published before.

Implicitly restarted Arnoldi with purification for the shift-invert transformation

Abstract: The need to determine a few eigenvalues of a large sparse generalised eigenvalue problem AZ = λBx with positive semidefinite B arises in many physical situations, for example, in a stability analysis of the discretised Navier-Stokes equation. A common technique is to apply Arnoldi's method to the shift-invert transformation, but this can suffer from numerical instabilities as is illustrated by a numerical example. In this paper, a new method that avoids instabilities is presented which is based on applying the implicitly restarted Arnoldi method with the B semi-inner product and a purification step. The paper contains a rounding error analysis and ends with brief comments on some extensions.

On some computational problems in finite Abelian groups

Abstract: We present new algorithms for computing orders of elements, discrete logarithms, and structures of finite abelian groups. We estimate the computational complexity and storage requirements, and we explicitly determine the O-constants and Ω-constants. We implemented the algorithms for class groups of imaginary quadratic orders and present a selection of our experimental results. Our algorithms are based on a modification of Shanks' baby-step giant-step strategy, and have the advantage that their computational complexity and storage requirements are relative to the actual order, discrete logarithm, or size of the group, rather than relative to an upper bound on the group order.

The resolution of the Gibbs phenomenon for spherical harmonics

Abstract: Spherical harmonics have been important tools for solving geophysical and astrophysical problems. Methods have been developed to effectively implement spherical harmonic expansion approximations. However, the Gibbs phenomenon was already observed by Weyl for spherical harmonic expansion approximations to functions with discontinuities, causing undesirable oscillations over the entire sphere. Recently, methods for removing the Gibbs phenomenon for one-dimensional discontinuous functions have been successfully developed by Gottlieb and Shu. They proved that the knowledge of the first N expansion coefficients (either Fourier or Gegenbauer) of a piecewise analytic function f(x) is enough to recover an exponentially convergent approximation to the point values of f(x) in any subinterval in which the function is analytic. Here we take a similar approach, proving that knowledge of the first N spherical harmonic coefficients yield an exponentially convergent approximation to a spherical piecewise smooth function f(θ, Φ) in any subinterval [θ 1 , θ 2 ], Φ ∈ [0,2π], where the function is analytic. Thus we entirely overcome the Gibbs phenomenon.

On generalized bisection of n -simplices

Abstract: A generalized procedure of bisection of n-simplices is introduced, where the bisection point can be an (almost) arbitrary point at one of the longest edges. It is shown that nested sequences of simplices generated by successive generalized bisection converge to a singleton, and an exact bound of the convergence speed in terms of diameter reduction is given. For regular simplices, which mark the worst case, the edge lengths of each worst and best simplex generated by successive bisection are given up to depth n. For n = 2 and 3, the sequence of worst case diameters is provided until it is halved.