Showing papers in "Siam Journal on Scientific and Statistical Computing in 1986"
TL;DR: An iterative method for solving linear systems, which has the property of minimizing at every step the norm of the residual vector over a Krylov subspace.
Abstract: We present an iterative method for solving linear systems, which has the property of minimizing at every step the norm of the residual vector over a Krylov subspace. The algorithm is derived from t...
10,907 citations
TL;DR: In this article, a pressure correction method for viscous incompressible flow is presented that is second order accurate in time and space, and a practical application is given for a numerical example.
Abstract: A pressure correction method for (time-dependent) viscous incompressible flow is presented that is second order accurate in time and space. The order of accuracy is proved for a model scheme and demonstrated for a numerical example. A practical application is given.
716 citations
TL;DR: In this paper, the linear inversion (deconvolution) of band-limited reflection seismograms is studied and a new cost functional is proposed that allows robust profile reconstruction in the presence of noise.
Abstract: We present a method for the linear inversion (deconvolution) of band-limited reflection seismograms. A large convolution problem is first broken up into a sequence of smaller problems. Each small problem is then posed as an optimization problem to resolve the ambiguity presented by the band-limited nature of the data. A new cost functional is proposed that allows robust profile reconstruction in the presence of noise. An algorithm for minimizing the cost functional is described. We present numerical experiments which simulate data interpretation using this procedure.
549 citations
TL;DR: Applications of the polar decomposition to factor analysis, aerospace computations and optimisation are outlined; and a new method is derived for computing the square root of a symmetric positive definite matrix.
Abstract: A quadratically convergent Newton method for computing the polar decomposition of a full-rank matrix is presented and analysed. Acceleration parameters are introduced so as to enhance the initial rate of convergence and it is shown how reliable estimates of the optimal parameters may be computed in practice.To add to the known best approximation property of the unitary polar factor, the Hermitian polar factor H of a nonsingular Hermitian matrix A is shown to be a good positive definite approximation to Aand $\frac{1}{2}(A + H)$ is shown to be a best Hermitian positive semi-definite approximation to A. Perturbation bounds for the polar factors are derived.Applications of the polar decomposition to factor analysis, aerospace computations and optimisation are outlined; and a new method is derived for computing the square root of a symmetric positive definite matrix.
446 citations
TL;DR: The major theme of this work is the development of an iterative scheme for the construction of a smooth surface, presented by global basis functions, which approximates only the smooth components of a set of scattered noisy data.
Abstract: In many applications one encounters the problem of approximating surfaces from data given on a set of scattered points in a two-dimensional domain. The global interpolation methods with Duchon's “thin plate splines” and Hardy's multiquadrics are considered to be of high quality; however, their application is limited, due to computational difficulties, to $ \sim 150$ data points. In this work we develop some efficient iterative schemes for computing global approximation surfaces interpolating a given smooth data. The suggested iterative procedures can, in principle, handle any number of data points, according to computer capacity. These procedures are extensions of a previous work by Dyn and Levin on iterative methods for computing thin-plate spline interpolants for data given on a square grid. Here the procedures are improved significantly and generalized to the case of data given in a general configuration.The major theme of this work is the development of an iterative scheme for the construction of a smooth surface, presented by global basis functions, which approximates only the smooth components of a set of scattered noisy data. The novelty in the suggested method is in the construction of an iterative procedure for low-pass filtering based on detailed spectral properties of a preconditioned matrix. The general concepts of this approach can also be used in designing iterative computation procedures for many other problems.The interpolation and smoothing procedures are tested, and the theoretical results are verified, by many numerical experiments.
429 citations
TL;DR: In this article, the authors present the results of an exhaustive search to find optimal full period multipliers for the multiplicative congruential random number generator with prime modulus $2^{31} - 1$.
Abstract: This paper presents the results of an exhaustive search to find optimal full period multipliers for the multiplicative congruential random number generator with prime modulus $2^{31} - 1$. Here a multiplier is said to be optimal if the distance between adjacent parallel hyperplanes on which k-tuples lie does not exceed the minimal achievable distance by more than 25 percent for $k = 2, \cdots ,6$. This criterion is considerably more stringent than prevailing standards of acceptability and leads to a total of only 414 multipliers among the more than 534 million candidate multipliers.Section 1 reviews the basic properties of linear congruential generators and § 2 describes worst case performance measures. These include the maximal distance between adjacent parallel hyperplanes, the minimal number of parallel hyperplanes, the minimal distance between k-tuples, the lattice ratio and the discrepancy. Section 3 presents the five best multipliers and compares their performances with those of three commonly employed multipliers for all measures but the lattice test. Comparisons using packing measures in the space of k-tuples and in the dual space are also made. Section 4 presents the results of applying a battery of statistical tests to the best five to detect local departures from randomness. None were found. The Appendix contains a list of all optimal multipliers.
215 citations
TL;DR: In this article, the FG-algorithm was proposed to find an orthogonal matrix B such that the deviation of the matrix B from diagonality is minimized by π(B^T A_1 B, \cdots,B^t A_k B;n_1, \cdot,A_k ;n_ 1, \CDots,n_k ) = \prod (n = 1}^k [\det (A_1 )]^{n_i } /[\det(A_i )]^{
Abstract: For $K \geqq 1$ positive definite symmetric matrices $A_1 , \cdots ,A_k $ of dimension $p \times p$ we define the function $\phi (A_1 , \cdots ,A_k ;n_1 , \cdots ,n_k ) = \prod _{i = 1}^k [\det (\operatorname{diag} A_1 )]^{n_i } /[\det (A_i )]^{n_i } $, where $n_i $ are positive constants, as a measure of simultaneous deviation of $A_1 , \cdots ,A_k $ from diagonality. We give an iterative algorithm, called the FG-algorithm, to find an orthogonal $p \times p$-matrix B such that $\phi (B^T A_1 B, \cdots ,B^T A_k B;n_1 , \cdots ,n_k )$ is minimum. The matrix B is said to transform $A_1 , \cdots ,A_k $ simultaneously to nearly diagonal form. Conditions for the uniqueness of the solution are given.The FG-algorithm can be used to find the maximum likelihood estimates of common principal components in k groups (Flury (1984)). For $k = 1$, the FG-algorithm computes the characteristic vectors of the single positive definite symmetric matrix $A_i $.
194 citations
TL;DR: In this paper, it was shown that the reacting compressible Navier-Stokes equations have dynamically stable weak detonations which occur in bifurcating wave patterns from strong detonation initial data.
Abstract: Several remarkable theoretical and computational properties of reacting shock waves are both documented and analyzed. In particular, for sufficiently small heat release or large reaction rate, we demonstrate that the reacting compressible Navier–Stokes equations have dynamically stable weak detonations which occur in bifurcating wave patterns from strong detonation initial data. In the reported calculations, an increase in reaction rate by a factor of 5 is sufficient to create the bifurcation from a spiked nearly Z-N-D detonation to the wave pattern with a precursor weak detonation. The numerical schemes used in the calculations are fractional step methods based on the use of a second order Godunov method in the inviscid hydrodynamic sweep; on sufficiently coarse meshes in inviscid calculations, these fractional step schemes exhibit qualitatively similar but purely numerical bifurcating wave patterns with numerical weak detonations. We explain this computational phenomenon theoretically through a new class of nonphysical discrete travelling waves for the difference scheme which are numerical weak detonations. The use of simplified model equations both to predict and analyze the theoretical and numerical phenomena is emphasized.
191 citations
TL;DR: The key idea is a new triangular solver that uses depth-first search and topological ordering to take advantage of sparsity in the right-hand side to factor sparse matrices by Gaussian elimination with partial privoting in time proportional to the number of arithmetic operations.
Abstract: Existing sparse partial pivoting algorithms can spend asymptomatically more time manipulating data structures than doing arithmetic, although they are tuned to be efficient on many large problems. We present an algorithm to factor sparse matrices by Gaussian elimination with partial privoting in time proportional to the number of arithmetic operations. Implementing this algorithm requires only simple data structures and gives a code that is competitive with, and often faster than, existing sparse codes. The key idea is a new triangular solver that uses depth-first search and topological ordering to take advantage of sparsity in the right-hand side.
188 citations
TL;DR: Davidson’s method for computing a few eigenpairs of large sparse symmetric matrices is generalized to a method which offers a powerful way of applying preconditioning techniques developed for solving systems of linear equations to solving eigenvalue problems.
Abstract: This paper analyzes Davidson’s method for computing a few eigenpairs of large sparse symmetric matrices. An explanation is given for why Davidson’s method often performs well but occasionally performs very badly. Davidson’s method is then generalized to a method which offers a powerful way of applying preconditioning techniques developed for solving systems of linear equations to solving eigenvalue problems.
175 citations
TL;DR: With the correct choice of ordering the algorithm can be implemented using systolic array processors (Gentleman, personal communication), and can also be used to compute any CS decomposition of a unitary matrix.
Abstract: An algorithm is described for computing the generalized singular value decomposition of $A(m \times n)$ and $B(p \times n)$. Unitary matrices U, V and Q are developed so that $U^H AQ$ and $V^H BQ$ have as many nonzero parallel rows as possible, and these correspond to the common row space of the two matrices. The algorithm consists of an iterative sequence of cycles where each cycle is made up of the serial application of $2 \times 2$ generalized singular value decompositions. Convergence appears to be at least quadratic. With the correct choice of ordering the algorithm can be implemented using systolic array processors (Gentleman,personal communication). The algorithm can also be used to compute any CS decomposition of a unitary matrix.
TL;DR: In this article, the behavior of numerical ODE methods for the solution of systems of differential equations coupled with algebraic constraints is investigated, showing how to overcome problems of matrix ill-conditioning, and giving convergence tests and error tests which are supported by theory.
Abstract: In this paper we investigate the behavior of numerical ODE methods for the solution of systems of differential equations coupled with algebraic constraints. Systems of this form arise frequently in the modelling of problems from physics and engineering; we study some particular examples from fluid dynamics and constrained mechanical systems. We investigate some of the practical difficulties of implementing variable-stepsize backward differentiation formulas for the solution of these equations, showing how to overcome problems of matrix ill-conditioning, and giving convergence tests and error tests which are supported by theory.
TL;DR: In this article, a numerical solution of the two-fluid incompressible Euler equation is used to study the Rayleigh-Taylor instability, which has the distinguishing feature of being a predominantly Eulerian method in which sharp interfaces are preserved with zero numerical diffusion.
Abstract: A numerical solution of the two-fluid incompressible Euler equation is used to study the Rayleigh–Taylor instability. The solution is based on the method of front tracking, which has the distinguishing feature of being a predominantly Eulerian method in which sharp interfaces are preserved with zero numerical diffusion. In this paper validation of the method is obtained by comparison with existing numerical solutions based on conformal mapping. An initial study of heterogeneity is presented.
TL;DR: The Hamiltonian-Schur decomposition (HSC) decomposition as discussed by the authors is a variant of Schur decompositions for Hamiltonian matrices that arise from single input control systems.
Abstract: This paper presents a variant $QR$ algorithm for calculating a Hamiltonian–Schur decomposition [10]. It is defined for Hamiltonian matrices that arise from single input control systems. Numerical stability and Hamiltonian structure are preserved by using unitary symplectic similarity transformations. Following a finite step reduction to a Hessenberg-like condensed form, a sequence of similarity transformations yields a permuted triangular matrix. As the iteration converges, it deflates into problems of lower dimension. Convergence is accelerated by varying a scalar shift. When the Hamiltonian matrix is real, complex arithmetic can be avoided by using an implicit double shift technique. The Hamiltonian-Schur decomposition yields the same invariant subspace information as a Schur decomposition but requires significantly less work and storage for problems of size greater than about 20.
TL;DR: In this article, a rigorous analysis is presented for the method of modified equations whereby its range of applicability and its shortcomings are delineated, and the theoretical findings are confirmed throughout by computational experiments.
Abstract: A rigorous analysis is presented for the method of modified equations whereby its range of applicability and its shortcomings are delineated. Numerous examples from different areas are presented and the theoretical findings are confirmed throughout by computational experiments.
TL;DR: In this article, the singular value decomposition of a product of two general matrices without explicitly forming the product is computed using plane rotations applied to the two matrices separately, based on an earlier Jacobi-like method due to Kogbetliantz.
Abstract: An algorithm is developed for computing the singular value decomposition of a product of two general matrices without explicitly forming the product. The algorithm is based on an earlier Jacobi-like method due to Kogbetliantz and uses plane rotations applied to the two matrices separately. A triangular variant of the basic algorithm is developed that reduces the amount of work required.
[...]
TL;DR: Multicolor SOR methods can be found that have the same rate of convergence as the natural rowwise SOR method for a wide range of stencils used to discretize partial differential equations and can be efficiently implemented on a wide class of parallel computers.
Abstract: The work of Young in 1950, see Young [1950], [1971], showed that the Red/ Black ordering and the natural rowwise ordering of matrices with Property A, such as those arising from the 5-point discretization of Poisson's equation, lead to SOR iteration matrices with identical eigenvalues. With the advent of parallel computers, multicolor point SOR schemes have been proposed for more complicated stencils on 2-dimensional rectangular grids, see Adams and Ortega [1982] for example, but to our knowledge, no theory has been provided for the rate of convergence of these methods relative to that of the natural rowwise scheme.New results show that certain matrices may be reordered so the resulting multicolor SOR matrix has the same eigenvalues as that for the original ordering. In addition, for a wide range of stencils, we show how to choose multicolor orderings so the multicolor SOR matrices have the same eigenvalues as the natural rowwise SOR matrix. The strategy for obtaining these orderings is based on “data flow” concepts and can be used to reach Young's conclusions above for the 5-point stencil.The importance of these results is threefold. Firstly, a constructive and easy means of finding these multicolorings is a direct consequence of the theory; secondly, multicolor SOR methods can be found that have the same rate of convergence as the natural rowwise SOR method for a wide range of stencils used to discretize partial differential equations; and thirdly, these multicolor SOR methods can be efficiently implemented on a wide class of parallel computers.
TL;DR: In this article, the numerical consistency of equality and inequality invariants known to be obeyed by the true solutions of ODES is discussed, and a general framework is used to justify a particular form of Shampine's technique for both one-step and multistep methods in the case of equality invariants.
Abstract: Techniques for maintaining the numerical consistency of equality and inequality invariants known to be obeyed by the true solutions of ODES are discussed. Shampine [11] has examined the problem and made some recommendations for one-step methods. A general framework is used to justify a.particular form of his technique for both one-step and multistep methods in the case of equality invariants. This framework can also be used for inequality invariants, but is not satisfactory for multistep methods. However, in a large class of inequality invariants, typified by those arising in chemical kinetic problems, it is shown that the invariant can be satisfied automatically by the numerical method.
TL;DR: In this paper, an upwind second-order scheme is proposed for time-dependent, inviscid, compressible fluid flow in one space dimension and variable cross-section, where the basic ingredient of the method is an analytic solution of the corresponding GRP (Generalized Riemann Problem), leading to explicit expressions for the fluxes.
Abstract: An upwind second-order scheme is proposed for time-dependent, inviscid, compressible fluid flow in one space dimension and variable cross-section. The basic ingredient of the method is an analytic solution of the corresponding GRP (Generalized Riemann Problem), leading to explicit expressions for the fluxes. Sonic lines are given a special treatment. The only accessory technique used is a simple monotonicity algorithm. Results of two test problems are shown: (a) the cylindrical converging shock problem, where boundary conditions at $r = 0$ are shown to be natural extensions of the analytic method; (b) a quasi 1-D nozzle problem.
TL;DR: This paper describes data structures and algorithms for the automatic generation of adaptive subgrids, a technique used with adaptive mesh refinement for solving partial differential equations, which generate a nested sequence of finer and finer grids on an underlying coarse grid.
Abstract: This paper describes data structures and algorithms for the automatic generation of adaptive subgrids, a technique used with adaptive mesh refinement for solving partial differential equations. Our algorithms generate a nested sequence of finer and finer grids on an underlying coarse grid. There are two aspects to the data structures. Trees are used to do the grid management for this type of grid structure.Secondly, the automatic grid generation algorithms use data structures with special nearest neighbor properties. Examples of grids from actual adaptive numerical computations are shown.
TL;DR: An iterative method for solving large sparse nonsymmetric linear systems of equations that enhances Manteuflel's adaptive Chebyshev method with a conjugate gradient-like method that replaces the modified power method for computing needed eigenvalue estimates with Arnoldi’s method.
Abstract: We present an iterative method for solving large sparse nonsymmetric linear systems of equations that enhances Manteuflel’s adaptive Chebyshev method with a conjugate gradient-like method. The new method replaces the modified power method for computing needed eigenvalue estimates with Arnoldi’s method, which can be used to simultaneously compute eigenvalues and to improve the approximate solution. Convergence analysis and numerical experiments suggest that the method is more efficient than the original adaptive Chebyshev algorithm.
TL;DR: A more general family of nodal schemes with O(h) convergence for any positive integer k under appropriate smoothness assumptions is presented, first introduced within a straightforward nonconforming finite element framework and then led to special numerical quadrature schemes which can be obtained directly through basic physical principles.
Abstract: Nodal schemes have been originally developed in numerical reactor calculation, especially in the area of neutron diffusion problems. Broadly speaking, they constitute accurate and fast methods sharing many attractive features of the finite element as well as of the finite difference methods. In a previous work, some of the simplest nodal schemes were identified with nonstandard nonconforming finite element schemes exhibiting $O(h)$ convergence in $H^1 $ norm. We present here a more general family of nodal schemes with $O(h^k )$ convergence for any positive integer k under appropriate smoothness assumptions. This new family is first introduced within a straightforward nonconforming finite element framework. Under special numerical quadrature schemes, we are then led to nodal schemes which can be obtained directly through basic physical principles. Finally, dimensionally reduced versions are obtained by transverse integration and stand as strong candidates to practical implementations of the alternating direction type.
TL;DR: A method for computing the smallest eigenvalue of a symmetric positive definite Toeplitz matrix relies solely upon the Levinson–Durbin algorithm.
Abstract: A method for computing the smallest eigenvalue of a symmetric positive definite Toeplitz matrix is given. It relies solely upon the Levinson–Durbin algorithm. The procedure involves a combination of bisection and Newton's method. Good starting values are also shown to be obtainable from the Levinson–Durbin algorithm.
TL;DR: This work introduces a new multigrid continuation method for computing solutions of nonlinear elliptic eigenvalue problems which contain limit points (also called turning points or folds), which produces considerable storage savings over direct continuation methods, as well as better initial coarse grid approximations, and avoids complicated algorithms for determining the parameter on finer grids.
Abstract: We introduce a new multigrid continuation method for computing solutions of nonlinear elliptic eigenvalue problems which contain limit points (also called turning points or folds). Our method combines the frozen tau technique of Brandt with pseudo-arc length continuation and correction of the parameter on the coarsest grid. This produces considerable storage savings over direct continuation methods,as well as better initial coarse grid approximations, and avoids complicated algorithms for determining the parameter on finer grids. We provide numerical results for second, fourth and sixth order approximations to the two-parameter, two-dimensional stationary reaction-diffusion problem:\[ \Delta u + \lambda \exp (u/(1 + \alpha u)) = 0. \]For the higher order interpolations we use bicubic and biquintic splines. The convergence rate is observed to be independent of the occurrence of limit points.
TL;DR: In this article, an importance sampling Monte Carlo method is proposed to approximate each of the weighted average of products of one-dimensional integrals and a prior error bound and a posterior error bound for the ratio are developed to measure the efficiency of the Monte Carlo approach.
Abstract: The computation of Bayes estimators based on mixtures of Dirichlet processes is treated in this article. These estimators may be written as ratios of two multidimensional integrals, each of which may be decomposed into a weighted average of products of one-dimensional integrals. An importance sampling Monte Carlo method is proposed to approximate each of the weighted averages. A prior error bound for each of the Monte Carlo estimators and a posterior error bound for the ratio are developed to measure the efficiency of the Monte Carlo method. Jackknife and random group error estimates are also considered. Two examples are given which illustrate the computation of the Bayes estimators.
TL;DR: PLTMGC as discussed by the authors is a program package for solving nonlinear elliptic systems that have explicit dependence on a scalar parameter, in addition to being able to compute solutions for fixed parameter values, it can be used to solve the linear eigenvalue problem, trace solution branches, locate singular points (simple turning points and bifurcation points).
Abstract: PLTMGC is a program package for solving nonlinear elliptic systems that have explicit dependence on a scalar parameter. In addition to being able to compute solutions for fixed parameter values, it can be used to solve the linear eigenvalue problem, trace solution branches, locate singular points (simpleturning points and bifurcation points) and switch branch at simple bifurcation points. A multi-grid continuation approach is employed in which a continuation procedure is used to follow the solution curve on the coarsest grid and a multi-grid algorithm is used to refine the solution at selected points using an adaptive mesh refinement strategy. Some numerical examples illustrating the performance of the package are given.
TL;DR: An algorithm for producing a triangular mesh in a convex polygon is presented and the expected time complexity is shown to be linear in the number of triangles in the mesh.
Abstract: An algorithm for producing a triangular mesh in a convex polygon is presented. It is used in a method for the finite element triangulation of a complex polygonal region of the plane in which the region is decomposed into convex polygons. The interior vertices of the mesh are chosen to be on a quasi-uniform grid, different mesh spacings are specified for the edges of the polygon, and the mesh is a Delaunay triangulation. The correctness of the algorithm is proved and the expected time complexity is shown to be linear in the number of triangles in the mesh.
TL;DR: In this paper, the stability of the Sherman-Morrison-Woodbury formula was studied and it was shown that if the original matrices, A and B, are well conditioned, then there exists matrices U and V such that the Sherman -Morrison Woodbury formula is stable when applied to $A = B - UV^T $.
Abstract: In this paper, we address the stability of the Sherman–Morrison–Woodbury formula. Our main result states that if the original matrices, A and B, are well conditioned, then there exists matrices U and V such that the Sherman–Morrison–Woodbury formula is stable when applied to $A = B - UV^T $.
TL;DR: In this paper, an algorithm for computing the structure elements associated with the Kronecker canonical form (KCF) of a matrix pencil, where A and B are complex m by n matrices, is presented.
Abstract: An algorithm (RGSVD) for computing the structure elements associated with the Kronecker canonical form (KCF) of a matrix pencil $A - \lambda B$, where A and B are complex m by n matrices, is presented. RGSVD is based on repeated generalized singular value decompositions (or more precisely cosine-sine decompositions of partitioned orthonormal matrices). It extracts the structures of the zero and/or the infinite eigenvalues together with the left (row) or right (column) minimal indices of $A - \lambda B$. By accumulating equivalence transformations, RGSVD also produces pairs of reducing subspaces associated with e.g. the zero structure and the right Kronecker indices.
TL;DR: A parallel method for computing the QR-decomposition of an $n \times n$ matrix is proposed and can be extended to handle an $m \timesn$ matrix $(m \geqq n)$.
Abstract: A parallel method for computing the QR-decomposition of an $n \times n$ matrix is proposed. It requires $O(n^2 )$ processors and $O(n)$ units of time. The method can be extended to handle an $m \times n$ matrix $(m \geqq n)$. The requirements then become $O(n^2 )$ processors and $O(m)$ time.