Showing papers in "Mathematics of Computation in 1996"

Maximally equidistributed combined Tausworthe generators

Abstract: Tausworthe random number generators based on a primitive trinomial allow an easy and fast implementation when their parameters obey certain restrictions. However, such generators, with those restrictions, have bad statistical properties unless we combine them. A generator is called maximally equidistributed if its vectors of successive values have the best possible equidistribution in all dimensions. This paper shows how to find maximally equidistributed combinations in an efficient manner, and gives a list of generators with that property. Such generators have a strong theoretical support and lend themselves to very fast software implementations.

Nonsmooth data error estimates for approximations of an evolution equation with a positive-type memory term

Abstract: We study the numerical approximation of an integro-differential equation which is intermediate between the heat and wave equations. The proposed discretization uses convolution quadrature based on the first- and second-order backward difference methods in time, and piecewise linear finite elements in space. Optimal-order error bounds in terms of the initial data and the inhomogeneity are shown for positive times, without assumptions of spatial regularity of the data.

On error estimates of the projection methods for the Navier-Stokes equations: second-order schemes

Abstract: We present in this paper a rigorous error analysis of several projection schemes for the approximation of the unsteady incompressible Navier-Stokes equations. The error analysis is accomplished by interpreting the respective projection schemes as second-order time discretizations of a perturbed system which approximates the Navier-Stokes equations. Numerical results in agreement with the error analysis are also presented.

Balancing domain decomposition for problems with large jumps in coefficients

Abstract: The Balancing Domain Decomposition algorithm uses in each iteration solution of local problems on the subdomains coupled with a coarse problem that is used to propagate the error globally and to guarantee that the possibly singular local problems are consistent. The abstract theory introduced recently by the first-named author is used to develop condition number bounds for conforming linear elements in two and three dimensions. The bounds are independent of arbitrary coefficient jumps between subdomains and of the number of subdomains, and grow only as the squared logarithm of the mesh size h. Computational experiments for two- and three-dimensional problems confirm the theory.

On restarting the Arnoldi method for large nonsymmetric eigenvalue problems

Abstract: The Arnoldi method computes eigenvalues of large nonsymmetric matrices. Restarting is generally needed to reduce storage requirements and orthogonalization costs. However, restarting slows down the convergence and makes the choice of the new starting vector difficult if several eigenvalues are desired. We analyze several approaches to restarting and show why Sorensen's implicit QR approach is generally far superior to the others. Ritz vectors are combined in precisely the right way for an effective new starting vector. Also, a new method for restarting Arnoldi is presented. It is mathematically equivalent to the Sorensen approach but has additional uses.

Primes in arithmetic progressions

Abstract: Strengthening work of Rosser, Schoenfeld, and McCurley, we establish explicit Chebyshev-type estimates in the prime number theorem for arithmetic progressions, for all moduli k ≤ 72 and other small moduli.

On the optimal stability of the Bernstein basis

Abstract: We show that the Bernstein polynomial basis on a given interval is optimally stable, in the sense that no other nonnegative basis yields systematically smaller condition numbers for the values or roots of arbitrary polynomials on that interval. This result follows from a partial ordering of the set of all nonnegative bases that is induced by nonnegative basis transformations. We further show, by means of some low-degree examples, that the Bernstein form is not uniquely optimal in this respect. However, it is the only optimally stable basis whose elements have no roots on the interior of the chosen interval. These ideas are illustrated by comparing the stability properties of the power, Bernstein, and generalized Ball bases.

The p and hp versions of the finite element method for problems with boundary layers

Abstract: We study the uniform approximation of boundary layer functions exp(-x/d) for x ∈ (0,1), d ∈ (0,1], by the p and hp versions of the finite element method. For the p version (with fixed mesh), we prove super-exponential convergence in the range p + 1/2 > e/(2d). We also establish, for this version, an overall convergence rate of O(p -1 √ln p) in the energy norm error which is uniform in d, and show that this rate is sharp (up to the √ln p term) when robust estimates uniform in d ∈ (0,1] are considered. For the p version with variable mesh (i.e., the hp version), we show that exponential convergence, uniform in d ∈ (0,1], is achieved by taking the first element at the boundary layer to be of size O(pd). Numerical experiments for a model elliptic singular perturbation problem show good agreement with our convergence estimates, even when few degrees of freedom are used and when d is as small as, e.g., 10 -8 . They also illustrate the superiority of the hp approach over other methods, including a low-order h version with optimal exponential mesh refinement. The estimates established in this paper are also applicable in the context of corresponding spectral element methods.

Two-level additive Schwarz preconditioners for nonconforming finite element methods

Abstract: Two-level additive Schwarz preconditioners are developed for the nonconforming P1 finite element approximation of scalar second-order symmetric positive definite elliptic boundary value problems, the Morley finite element approximation of the biharmonic equation, and the divergence-free nonconforming P1 finite element approximation of the stationary Stokes equations. The condition numbers of the preconditioned systems are shown to be bounded independent of mesh sizes and the number of subdomains in the case of generous overlap.

An algorithm for matrix extension and wavelet construction

Abstract: This paper gives a practical method of extending an n x r matrix P(z), r ≤ n, with Laurent polynomial entries in one complex variable z, to a square matrix also with Laurent polynomial entries. If P(z) has orthonormal columns when z is restricted to the torus T, it can be extended to a paraunitary matrix. If P(z) has rank r for each z ∈ T, it can be extended to a matrix with nonvanishing determinant on T. The method is easily implemented in the computer. It is applied to the construction of compactly supported wavelets and prewavelets from multiresolutions generated by several univariate scaling functions with an arbitrary dilation parameter.

A continuous space-time finite element method for the wave equation

Abstract: We consider a finite element method for the nonhomogeneous second-order wave equation, which is formulated in terms of continuous approximation functions in both space and time, thereby giving a unified treatment of the spatial and temporal discretizations. Our analysis uses primarily energy arguments, which are quite common for spatial discretizations but not for time. We present a priori nodal (in time) superconvergence error estimates without any special time step restrictions. Our method is based on tensor-product spaces for the full discretization.

Quality local refinement of tetrahedral meshes based on 8-subtetrahedron subdivision

Abstract: Let T be a tetrahedral mesh. We present a 3-D local refinement algorithm for T which is mainly based on an 8-subtetrahedron subdivision procedure, and discuss the quality of refined meshes generated by the algorithm. It is proved that any tetrahedron T ∈ T produces a finite number of classes of similar tetrahedra, independent of the number of refinement levels. Furthermore, η(T i n ) ≥ cη(T), where T ∈ T, c is a positive constant independent of T and the number of refinement levels, T i n is any refined tetrahedron of T, and η is a tetrahedron shape measure. It is also proved that local refinements on tetrahedra can be smoothly extended to their neighbors to maintain a conforming mesh. Experimental results show that the ratio of the number of tetrahedra actually refined to the number of tetrahedra chosen for refinement is bounded above by a small constant.

Explicit bounds for primes in residue classes

Abstract: Let E/K be an abelian extension of number fields, with E ¬= Q. Let Δ and n denote the absolute discriminant and degree of E. Let σ denote an element of the Galois group of E/K. We prove the following theorems, assuming the Extended Riemann Hypothesis: (1) There is a degree-1 prime p of K such that (p/E/K) = σ, satisfying Np ≤ (1+ o(1))(logΔ + 2n) 2 . (2) There is a degree-1 prime p of K such that (p/E/K) generates the same group as σ, satisfying Np ≤ (1 + o(1))(log Δ) 2 . (3) For K = Q, there is a prime p such that (p/E/Q) = σ, satisfying P ≤ (1 + o(1))(log Δ) 2 . In (1) and (2) we can in fact take p to be unramified in K/Q. A special case of this result is the following. (4) If gcd(m,q) = 1, the least prime p? m (mod q) satisfies p ≤ (1 + o(1))(φ(q)log q) 2 . It follows from our proof that (1)-(3) also hold for arbitrary Galois extensions, provided we replace σ by its conjugacy class . Our theorems lead to explicit versions of (1)-(4), including the following: the least prime p? m (mod q) is less than 2(q log q) 2 .

Fast evaluation of the Gaunt coefficients

Abstract: Addition theorems for vector spherical harmonics require the determination of the Gaunt coefficients that appear in a linearization expansion of the product of two associated Legendre functions. This paper presents an algorithm for the efficient calculation of these coefficients through solving the most appropriate (lower triangular) linear system and derives all relevant recurrence relations needed in the calculation. This algorithm is also applicable to the calculation of the Clebsch-Gordan coefficients that are closely related to the Gaunt coefficients and are frequently encountered in the quantum theory of angular momentum. The new method described in this paper reduces the computing time to ∼ 1%, compared to the existing formulation that is widely used. This new method can be applied to the calculation of both low- and high-degree coefficients, while the existing formulation works well only for low degrees.

An efficient spectral method for ordinary differential equations with rational function coefficients

Abstract: We present some relations that allow the efficient approximate inversion of linear differential operators with rational function coefficients. We employ expansions in terms of a large class of orthogonal polynomial families, including all the classical orthogonal polynomials. These families obey a simple three-term recurrence relation for differentiation, which implies that on an appropriately restricted domain the differentiation operator has a unique banded inverse. The inverse is an integration operator for the family, and it is simply the tridiagonal coefficient matrix for the recurrence. Since in these families convolution operators (i.e. matrix representations of multiplication by a function) are banded for polynomials, we are able to obtain a banded representation for linear differential operators with rational coefficients. This leads to a method of solution of initial or boundary value problems that, besides having an operation count that scales linearly with the order of truncation N, is computationally well conditioned. Among the applications considered is the use of rational maps for the resolution of sharp interior layers.

Anti-Gaussian quadrature formulas

Abstract: An anti-Gaussian quadrature formula is an (n + 1)-point formula of degree 2n - 1 which integrates polynomials of degree up to 2n + 1 with an error equal in magnitude but of opposite sign to that of the n-point Gaussian formula. Its intended application is to estimate the error incurred in Gaussian integration by halving the difference between the results obtained from the two formulas. We show that an anti-Gaussian formula has positive weights, and that its nodes are in the integration interval and are interlaced by those of the corresponding Gaussian formula. Similar results for Gaussian formulas with respect to a positive weight are given, except that for some weight functions, at most two of the nodes may be outside the integration interval. The anti-Gaussian formula has only interior nodes in many cases when the Kronrod extension does not, and is as easy to compute as the (n + 1)-point Gaussian formula.

Some new error estimates for Ritz-Galerkin methods with minimal regularity assumptions

Abstract: New uniform error estimates are established for finite element approximations u h of solutions u of second-order elliptic equations Lu = f using only the regularity assumption ∥u∥ 1 ≤ c∥f∥ −1 . Using an Aubin-Nitsche type duality argument we show for example that, for arbitrary (fixed) e sufficiently small, there exists an h 0 such that for 0

Product integration for Volterra integral equations of the second kind with weakly singular kernels

Abstract: We introduce a new numerical approach for solving Volterra integral equations of the second kind when the kernel contains a mild singularity. We give a convergence result. We also present numerical examples which show the performance and efficiency of our method.

The 4-triangles longest-side partition of triangles and linear refinement algorithms

Abstract: In this paper we study geometrical properties of the iterative 4-triangles longest-side partition of triangles (and of a 3-triangles partition), as well as practical algorithms based on these partitions, used both directly for the triangulation refinement problem, and as a basis for point insertion strategies in Delaunay refinement algorithms. The 4-triangles partition is obtained by joining the midpoint of the longest side with the opposite vertex and the midpoints of the two remaining sides. By means of simple geometrical arguments we show that the iterative partition of obtuse triangles systematically improves the triangles (while they remain obtuse) in the following sense: the sequence of smallest angles monotonically increases while the sequence of largest angles monotonically decreases in an amount (at least) equal to the smallest angle of each iteration. This allows us to improve the known bound on the smallest angle (without making use of previous results), and to obtain a better a priori bound on the number of similarly distinct triangles, as a function of the geometry of the initial triangle. Numerical evidence, showing that the practical behavior of the 4-triangles partition is in complete agreement with this theory, is included. A 4-triangles refinement algorithm is also discussed and illustrated. Furthermore, we show that the time cost of the algorithm is linear independently of the size of the triangulation.

The BKK root count in C n

Abstract: The root count developed by Bernshtein, Kushnirenko and Khovanskii only counts the number of isolated zeros of a polynomial system in the algebraic torus (C * ) n . In this paper, we modify this bound slightly so that it counts the number of isolated zeros in C n . Our bound is, apparently, significantly sharper than the recent root counts found by Rojas and in many cases easier to compute. As a consequence of our result, the Huber-Sturmfels homotopy for finding all the isolated zeros of a polynomial system in (C * ) n can be slightly modified to obtain all the isolated zeros in C n .

A MUSCL method satisfying all the numerical entropy inequalities

Abstract: We consider here second-order finite volume methods for one-dimensional scalar conservation laws. We give a method to determine a slope reconstruction satisfying all the exact numerical entropy inequalities. It avoids inhomogeneous slope limitations and, at least, gives a convergence rate of Δx l/2 . It is obtained by a theory of second-order entropic projections involving values at the nodes of the grid and a variant of error estimates, which also gives new results for the first-order Engquist-Osher scheme.

Analysis and convergence of the MAC scheme. II: Navier-Stokes equations

Abstract: The MAC discretization scheme for the incompressible Navier-Stokes equations is interpreted as a covolume approximation to the equations. Using some results from earlier papers dealing with covolume error estimates for div-curl equation systems, and under certain conditions on the data and the solutions of the Navier-Stokes equations, we obtain first-order error estimates for both the vorticity and the pressure.

On the beta expansion for Salem numbers of degree 6

Abstract: For a given β > 1, the beta transformation T = T β is defined for x ∈ [0,1] by Tx:= βx (mod 1). The number β is said to be a beta number if the orbit {T n (1)} n≥1 is finite, hence eventually periodic. It is known that all Pisot numbers are beta numbers, and it is conjectured that this is true for Salem numbers, but this is known only for Salem numbers of degree 4. Here we consider some computational and heuristic evidence for the conjecture in the case of Salem numbers of degree 6, by considering the set of 11836 such numbers of trace at most 15. Although the orbit is small for the majority of these numbers, there are some examples for which the orbit size is shown to exceed 10 9 and for which the possibility remains that the orbit is infinite. There are also some very large orbits which have been shown to be finite: an example is given for which the preperiod length is 39420662 and the period length is 93218808. This is in contrast to Salem numbers of degree 4 where the orbit size is bounded by 2β + 3. An heuristic probabilistic model is proposed which explains the difference between the degree-4 and degree-6 cases. The model predicts that all Salem numbers of degree 4 and 6 should be beta numbers but that degree-6 Salem numbers can have orbits which are arbitrarily large relative to the size of β. Furthermore, the model predicts that a positive proportion of Salem numbers of any fixed degree ≥ 8 will not be beta numbers. This latter prediction is not tested here.

Interior penalty preconditioners for mixed finite element approximations of elliptic problems

Abstract: It is established that an interior penalty method applied to second-order elliptic problems gives rise to a local operator which is spectrally equivalent to the corresponding nonlocal operator arising from the mixed finite element method. This relation can be utilized in order to construct preconditioners for the discrete mixed system. As an example, a family of additive Schwarz preconditioners for these systems is constructed. Numerical examples which confirm the theoretical results are also presented.

Computing p( x ): the Meissel, Lehmer, Lagarias, Miller, Odlyzko method

Abstract: Let π(x) denote the number of primes ≤ x. Our aim in this paper is to present some refinements of a combinatorial method for computing single values of π(x), initiated by the German astronomer Meissel in 1870, extended and simplified by Lehmer in 1959, and improved in 1985 by Lagarias, Miller and Odlyzko. We show that it is possible to compute π(x) in O(x 2/3 /log 2 x) time and O(x 1/3l log 3 x log log x) space. The algorithm has been implemented and used to compute π(10 18 ).

Efficient algorithms for computing the L 2 -discrepancy

Abstract: The \(L_2\)-discrepancy is a quantitative measure of precision for multivariate quadrature rules. It can be computed explicitly. Previously known algorithms needed \(O(m^2\)) operations, where \(m\) is the number of nodes. In this paper we present algorithms which require \(O(m(log m)^d)\) operations.

Strictly positive definite functions on spheres in Euclidean spaces

Abstract: In this paper we study strictly positive definite functions on the unit sphere of the m-dimensional Euclidean space. Such functions can be used for solving a scattered data interpolation problem on spheres. Since positive definite functions on the sphere were already characterized by Schoenberg some fifty years ago, the issue here is to determine what kind of positive definite functions are actually strictly positive definite. The study of this problem was initiated recently by Xu and Cheney (Proc. Amer. Math. Soc. 116 (1992), 977-981), where certain sufficient conditions were derived. A new approach, which is based on a critical connection between this problem and that of multivariate polynomial interpolation on spheres, is presented here. The relevant interpolation problem is subsequently analyzed by three different complementary methods. The first is based on the de Boor-Ron general least solution for the multivariate polynomial interpolation problem. The second, which is suitable only for m = 2, is based on the connection between bivariate harmonic polynomials and univariate analytic polynomials, and reduces the problem to the structure of the integer zeros of bounded univariate exponentials. Finally, the last method invokes the realization of harmonic polynomials as the polynomial kernel of the Laplacian, thereby exploiting some basic relations between homogeneous ideals and their polynomial kernels.

Efficiency of a posteriori BEM-error estimates for first-kind integral equations on quasi-uniform meshes

Abstract: In the numerical treatment of integral equations of the first kind using boundary element methods (BEM), the author and E. P. Stephan have derived a posteriori error estimates as tools for both reliable computation and self-adaptive mesh refinement. So far, efficiency of those a posteriori error estimates has been indicated by numerical examples in model situations only. This work affirms efficiency by proving the reverse inequality. Based on best approximation, on inverse inequalities and on stability of the discretization, and complementary to our previous work, an abstract approach yields a converse estimate. This estimate proves efficiency of an a posteriori error estimate in the BEM on quasi uniform meshes for Symm's integral equation, for a hypersingular equation, and for a transmission problem.

A priori error estimates for numerical methods for scalar conservation laws: Part I: the general approach

Abstract: We construct a general theory of a priori error estimates for scalar conservation laws by suitably modifying the original Kuznetsov approximation theory. As a first application of this general technique, we show that error estimates for conservation laws can be obtained without having to use explicitly any regularity properties of the approximate solution. Thus, we obtain optimal error estimates for the Engquist-Osher scheme without using the fact (i) that the solution is uniformly bounded, (ii) that the scheme is total variation diminishing, and (iii) that the discrete semigroup associated with the scheme has the L 1 -contraction property, which guarantees an upper bound for the modulus of continuity in time of the approximate solution.

Asymptotic semismoothness probabilities

Abstract: We call an integer semismooth with respect to y and z if each of its prime factors is ≤ y, and all but one are ≤ z. Such numbers are useful in various factoring algorithms, including the quadratic sieve. Let G(α,β) be the asymptotic probability that a random integer n is semismooth with respect to n β and n α . We present new recurrence relations for G and related functions. We then give numerical methods for computing G, tables of G, and estimates for the error incurred by this asymptotic approximation.