Showing papers in "SIAM Journal on Scientific Computing in 1999"

Weighted ENO Schemes for Hamilton--Jacobi Equations

Abstract: In this paper, we present a weighted ENO (essentially nonoscillatory) scheme to approximate the viscosity solution of the Hamilton--Jacobi equation: $$ \phi_t + H(x_1,\ldots,x_d,t,\phi,\phi_{x_1},\ldots,\phi_{x_d}) = 0. $$ This weighted ENO scheme is constructed upon and has the same stencil nodes as the third order ENO scheme but can be as high as fifth order accurate in the smooth part of the solution. In addition to the accuracy improvement, numerical comparisons between the two schemes also demonstrate that the weighted ENO scheme is more robust than the ENO scheme.

A Nonlinear Primal-Dual Method for Total Variation-Based Image Restoration

Abstract: We present a new method for solving total variation (TV) minimization problems in image restoration The main idea is to remove some of the singularity caused by the nondifferentiability of the quantity $| abla u|$ in the definition of the TV-norm before we apply a linearization technique such as Newton's method This is accomplished by introducing an additional variable for the flux quantity appearing in the gradient of the objective function, which can be interpreted as the normal vector to the level sets of the image u Our method can be viewed as a primal-dual method as proposed by Conn and Overton [ A Primal-Dual Interior Point Method for Minimizing a Sum of Euclidean Norms, preprint, 1994] and Andersen [Ph D thesis, Odense University, Denmark, 1995] for the minimization of a sum of Euclidean norms In addition to possessing local quadratic convergence, experimental results show that the new method seems to be globally convergent

A Subspace, Interior, and Conjugate Gradient Method for Large-Scale Bound-Constrained Minimization Problems

Abstract: A subspace adaptation of the Coleman--Li trust region and interior method is proposed for solving large-scale bound-constrained minimization problems. This method can be implemented with either sparse Cholesky factorization or conjugate gradient computation. Under reasonable conditions the convergence properties of this subspace trust region method are as strong as those of its full-space version. Computational performance on various large test problems is reported; advantages of our approach are demonstrated. Our experience indicates that our proposed method represents an efficient way to solve large bound-constrained minimization problems.

An Efficient, Interface-Preserving Level Set Redistancing Algorithm and Its Application to Interfacial Incompressible Fluid Flow

Abstract: In Sussman, Smereka, and Osher [ J. Comp. Phys., 94 (1994), pp. 146--159], a numerical scheme was presented for computing incompressible air--water flows using the level set method. Crucial to the above method was a new iteration method for maintaining the level set function as the signed distance from the zero level set. In this paper we implement a "constraint" along with higher order difference schemes in order to make the iteration method more accurate and efficient. Accuracy is measured in terms of the new computed signed distance function and the original level set function having the same zero level set. We apply our redistancing scheme to incompressible flows with noticeably better resolved results at reduced cost. We validate our results with experiment and theory. We show that our "distance level set scheme" with the added constraint competes well with available interface tracking schemes for basic advection of an interface. We perform basic accuracy checks and more stringent tests involving complicated interfacial structures. As with all level set schemes, our method is easy to implement.

A Note on Preconditioning for Indefinite Linear Systems

Abstract: Preconditioners are often conceived as approximate inverses. For nonsingular indefinite matrices of saddle-point (or KKT) form, we show how preconditioners incorporating an exact Schur complement lead to preconditioned matrices with exactly two or exactly three distinct eigenvalues. Thus approximations of the Schur complement lead to preconditioners which can be very effective even though they are in no sense approximate inverses.

Efficient Asymptotic-Preserving (AP) Schemes For Some Multiscale Kinetic Equations

Abstract: Many kinetic models of the Boltzmann equation have a diffusive scaling that leads to the Navier--Stokes type parabolic equations as the small scaling parameter approaches zero. In practical applications, it is desirable to develop a class of numerical schemes that can work uniformly with respect to this relaxation parameter, from the rarefied kinetic regimes to the hydrodynamic diffusive regimes. An earlier approach in [S. Jin, L. Pareschi, and G. Toscani, SIAM J. Numer. Anal., 35 (1998), pp. 2405--2439] reformulates such systems into the common hyperbolic relaxation system by Jin and Xin for hyperbolic conservation laws used to construct the relaxation schemes and then uses a multistep time-splitting method to solve the relaxation system. Here we observe that the combination of the two time-split steps may yield hyperbolic-parabolic systems that are more advantageous, in both stability and efficiency, for numerical computations. We show that such an approach yields a class of asymptotic-preserving (AP) schemes which are suitable for the computation of multiscale kinetic problems. We use the Goldstein--Taylor and Carleman models to illustrate this approach.

A Restricted Additive Schwarz Preconditioner for General Sparse Linear Systems

Abstract: We introduce some cheaper and faster variants of the classical additive Schwarz preconditioner (AS) for general sparse linear systems and show, by numerical examples, that the new methods are superior to AS in terms of both iteration counts and CPU time, as well as the communication cost when implemented on distributed memory computers. This is especially true for harder problems such as indefinite complex linear systems and systems of convection-diffusion equations from three-dimensional compressible flows. Both sequential and parallel results are reported.

A Fast Algorithm for Deblurring Models with Neumann Boundary Conditions

Abstract: Blur removal is an important problem in signal and image processing. The blurring matrices obtained by using the zero boundary condition (corresponding to assuming dark background outside the scene) are Toeplitz matrices for one-dimensional problems and block-Toeplitz--Toeplitz-block matrices for two-dimensional cases. They are computationally intensive to invert especially in the block case. If the periodic boundary condition is used, the matrices become (block) circulant and can be diagonalized by discrete Fourier transform matrices. In this paper, we consider the use of the Neumann boundary condition (corresponding to a reflection of the original scene at the boundary). The resulting matrices are (block) Toeplitz-plus-Hankel matrices. We show that for symmetric blurring functions, these blurring matrices can always be diagonalized by discrete cosine transform matrices. Thus the cost of inversion is significantly lower than that of using the zero or periodic boundary conditions. We also show that the use of the Neumann boundary condition provides an easy way of estimating the regularization parameter when the generalized cross-validation is used. When the blurring function is nonsymmetric, we show that the optimal cosine transform preconditioner of the blurring matrix is equal to the blurring matrix generated by the symmetric part of the blurring function. Numerical results are given to illustrate the efficiency of using the Neumann boundary condition.

A Cyclic Low-Rank Smith Method for Large Sparse Lyapunov Equations

Abstract: In this paper we present the cyclic low-rank Smith method, which is an iterative method for the computation of low-rank approximations to the solution of large, sparse, stable Lyapunov equations. It is based on a generalization of the classical Smith method and profits by the usual low-rank property of the right-hand side matrix. The requirements of the method are moderate with respect to both computational cost and memory. Furthermore, we propose a heuristic for determining a set of suboptimal alternating direction implicit (ADI) shift parameters. This heuristic, which is based on a pair of Arnoldi processes, does not require any a priori knowledge on the spectrum of the coefficient matrix of the Lyapunov equation. Numerical experiments show the efficiency of the iterative scheme combined with the heuristic for the ADI parameters.

A Simple Method for Compressible Multifluid Flows

Abstract: A simple second order accurate and fully Eulerian numerical method is presented for the simulation of multifluid compressible flows, governed by the stiffened gas equation of state, in hydrodynamic regime. Our numerical method relies on a second order Godunov-type scheme, with approximate Riemann solver for the resolution of conservation equations, and a set of nonconservative equations. It is valid for all mesh points and allows the resolution of interfaces. This method works for an arbitrary number of interfaces, for breakup and coalescence. It allows very high density ratios (up to 1000). It is able to compute very strong shock waves (pressure ratio up to 10 5). Contrary to all existing schemes (which consider the interface as a discontinuity) the method considers the interface as a numerical diffusion zone as contact discontinuities are computed in compressible single phase flows, but the variables describing the mixture zone are computed consistently with the density, momentum and energy. Several test problems are presented in one, two, and three dimensions. This method allows, for example, the computation of the interaction of a shock wave propagating in a liquid with a gas cylinder, as well as Richtmeyer--Meshkov instabilities, or hypervelocity impact, with realistic initial conditions. We illustrate our method with the Rusanov flux. However, the same principle can be applied to a more general class of schemes.

Long-Time-Step Methods for Oscillatory Differential Equations

Abstract: Considered are numerical integration schemes for nondissipative dynamical systems in which multiple time scales are present. It is assumed that one can do an explicit separation of the RHS "forces" into fast forces and slow forces such that (i) the fast forces contain the high frequency part of the solution, (ii) the fast forces are conservative, and (iii) the reduced problem consisting only of the fast forces can be integrated much more cheaply than the full problem. The fast forces are allowed to have low frequency components. Particular applications for which the schemes are intended include N-body problems (for which most of the forces are slow) and nonlinear wave phenomena (for which the fast forces can be propagated by spectral methods). The assumption of cheap integration of fast forces implies that the overall cost of integration is primarily determined by the step size used to sample the slow forces. A long-time-step method is one in which this step size exceeds half the period of the fastest normal mode present in the full system. An existing method that comes close to qualifying is the "impulse" method, also known as Verlet-I and r-RESPA. It is shown that it might fail, though, for a couple of reasons. First, it suffers a serious loss of accuracy if the step size is near a multiple of the period of a normal mode, and, second, it is unstable if the step size is near a multiple of half the period of a normal mode. Proposed in this paper is a "mollified" impulse method having an error bound that is independent of the frequency of the fast forces. It is also shown to possess superior stability properties. Theoretical results are supplemented by numerical experiments. The method is efficient and reasonably easy to implement.

Wavelet and Fourier Methods for Solving the Sideways Heat Equation

Abstract: We consider an inverse heat conduction problem, the sideways heat equation, which is a model of a problem, where one wants to determine the temperature on both sides of a thick wall, but where one side is inaccessible to measurements. Mathematically it is formulated as a Cauchy problem for the heat equation in a quarter plane, with data given along the line x=1, where the solution is wanted for $0 \leq x < 1$. The problem is ill-posed, in the sense that the solution (if it exists) does not depend continuously on the data. We consider stabilizations based on replacing the time derivative in the heat equation by wavelet-based approximations or a Fourier-based approximation. The resulting problem is an initial value problem for an ordinary differential equation, which can be solved by standard numerical methods, e.g., a Runge--Kutta method. We discuss the numerical implementation of Fourier and wavelet methods for solving the sideways heat equation. Theory predicts that the Fourier method and a method based on Meyer wavelets will give equally good results. Our numerical experiments indicate that also a method based on Daubechies wavelets gives comparable accuracy. As test problems we take model equations with constant and variable coefficients. We also solve a problem from an industrial application with actual measured data.

Preconditioning for the Steady-State Navier--Stokes Equations with Low Viscosity

Abstract: We introduce a preconditioner for the linearized Navier--Stokes equations that is effective when either the discretization mesh size or the viscosity approaches zero. For constant coefficient problems with periodic boundary conditions, we show that the preconditioning yields a system with a single eigenvalue equal to 1, so that performance is independent of both viscosity and mesh size. For other boundary conditions, we demonstrate empirically that convergence depends only mildly on these parameters and we give a partial analysis of this phenomenon. We also show that some expensive subsidiary computations required by the new method can be replaced by inexpensive approximate versions of these tasks based on iteration, with virtually no degradation of performance.

Hybrid Gauss-Trapezoidal Quadrature Rules

Abstract: A new class of quadrature rules for the integration of both regular and singular functions is constructed and analyzed. For each rule the quadrature weights are positive and the class includes rules of arbitrarily high-order convergence. The quadratures result from alterations to the trapezoidal rule, in which a small number of nodes and weights at the ends of the integration interval are replaced. The new nodes and weights are determined so that the asymptotic expansion of the resulting rule, provided by a generalization of the Euler--Maclaurin summation formula, has a prescribed number of vanishing terms. The superior performance of the rules is demonstrated with numerical examples and application to several problems is discussed.

Inexact Preconditioned Conjugate Gradient Method with Inner-Outer Iteration

Abstract: An important variation of preconditioned conjugate gradient algorithms is inexact preconditioner implemented with inner-outer iterations [G. H. Golub and M. L. Overton, Numerical Analysis, Lecture Notes in Math. 912, Springer, Berlin, New York, 1982], where the preconditioner is solved by an inner iteration to a prescribed precision. In this paper, we formulate an inexact preconditioned conjugate gradient algorithm for a symmetric positive definite system and analyze its convergence property. We establish a linear convergence result using a local relation of residual norms. We also analyze the algorithm using a global equation and show that the algorithm may have the superlinear convergence property when the inner iteration is solved to high accuracy. The analysis is in agreement with observed numerical behavior of the algorithm. In particular, it suggests a heuristic choice of the stopping threshold for the inner iteration. Numerical examples are given to show the effectiveness of this choice and to compare the convergence bound.

A Priori Sparsity Patterns for Parallel Sparse Approximate Inverse Preconditioners

Abstract: Parallel algorithms for computing sparse approximations to the inverse of a sparse matrix either use a prescribed sparsity pattern for the approximate inverse or attempt to generate a good pattern as part of the algorithm. This paper demonstrates that, for PDE problems, the patterns of powers of sparsified matrices (PSMs) can be used a priori as effective approximate inverse patterns, and that the additional effort of adaptive sparsity pattern calculations may not be required. PSM patterns are related to various other approximate inverse sparsity patterns through matrix graph theory and heuristics concerning the PDE's Green's function. A parallel implementation shows that PSM-patterned approximate inverses are significantly faster to construct than approximate inverses constructed adaptively, while often giving preconditioners of comparable quality.

A Deflated Version of the Conjugate Gradient Algorithm

Abstract: We present a deflated version of the conjugate gradient algorithm for solving linear systems. The new algorithm can be useful in cases when a small number of eigenvalues of the iteration matrix are very close to the origin. It can also be useful when solving linear systems with multiple right-hand sides, since the eigenvalue information gathered from solving one linear system can be recycled for solving the next systems and then updated.

A Discontinuous Galerkin Finite Element Method for Hamilton-Jacobi Equations

Abstract: In this paper, we present a discontinuous Galerkin finite element method for solving the nonlinear Hamilton--Jacobi equations. This method is based on the Runge--Kutta discontinuous Galerkin finite element method for solving conservation laws. The method has the flexibility of treating complicated geometry by using arbitrary triangulation, can achieve high-order accuracy with a local, compact stencil, and is suited for efficient parallel implementation. One- and two-dimensional numerical examples are given to illustrate the capability of the method. At least kth order of accuracy is observed for smooth problems when kth degree polynomials are used, and derivative singularities are resolved well without oscillations, even without limiters.

The Fringe Region Technique and the Fourier Method Used in the Direct Numerical Simulation of Spatially Evolving Viscous Flows

Abstract: To eliminate the problem with artificial boundary conditions and facilitate the use of Fourier methods, the fringe region (or filter, damping layer, absorbing layer, sponge layer) technique has been used in direct simulations of transitional and turbulent boundary layers. Despite the fact that good computational results have been obtained with this technique, it is not fully understood. The analysis in this paper indicates that the primary importance of the fringe region technique is to damp out the deviation associated with large scales in the direction normal to the wall. The lack of boundary conditions is compensated by the knowledge of an exact solution in the fringe region of the computational domain. The upstream influence from the fringe region is small. Numerical experiments verifying the theoretical predictions are presented.

Robustness and Scalability of Algebraic Multigrid

Abstract: Algebraic multigrid (AMG) is currently undergoing a resurgence in popularity, due in part to the dramatic increase in the need to solve physical problems posed on very large, unstructured grids. While AMG has proved its usefulness on various problem types, it is not commonly understood how wide a range of applicability the method has. In this study, we demonstrate that range of applicability, while describing some of the recent advances in AMG technology. Moreover, in light of the imperatives of modern computer environments, we also examine AMG in terms of algorithmic scalability. Finally, we show some of the situations in which standard AMG does not work well and indicate the current directions taken by AMG researchers to alleviate these difficulties.

A Stable Numerical Method for Inverting Shape from Moments

Abstract: We derive a stable technique, based upon matrix pencils, for the reconstruction of (or approximation by) polygonal shapes from moments. We point out that this problem can be considered the dual of 2-D numerical quadrature over polygonal domains. An analysis of the sensitivity of the problem is presented along with some numerical examples illustrating the relevant points. Finally, an application to the problem of gravimetry is explored where the shape of a gravitationally anomalous region is to be recovered from measurements of its exterior gravitational field.

On Updating Problems in Latent Semantic Indexing

Abstract: We develop new SVD-updating algorithms for three types of updating problems arising from latent semantic indexing (LSI) for information retrieval to deal with rapidly changing text document collections. We also provide theoretical justification for using a reduced-dimension representation of the original document collection in the updating process. Numerical experiments using several standard text document collections show that the new algorithms give higher (interpolated) average precisions than the existing algorithms, and the retrieval accuracy is comparable to that obtained using the complete document collection.

Robust Computational Algorithms for Dynamic Interface Tracking in Three Dimensions

Abstract: Front tracking provides sharp resolution of wave fronts through the active tracking of interfaces between distinct materials. A major challenge to this method is to handle changes in the interface topology. We describe two algorithms, implemented in the front tracking code FronTier, to model dynamic changes in three-dimensional interfaces. The two methods can be combined to give a hybrid method that is superior to each individual method. The success of these algorithms is shown by simulations of Rayleigh--Taylor instability, which is an interfacial instability driven by an acceleration directed across a material interface. Our numerical results are validated by comparing the numerical computation of the velocity of a single rising bubble with an analytic model for the bubble velocity.

Incomplete Cholesky Factorizations with Limited Memory

Abstract: We propose an incomplete Cholesky factorization for the solution of large-scale trust region subproblems and positive definite systems of linear equations. This factorization depends on a parameter p that specifies the amount of additional memory (in multiples of n, the dimension of the problem) that is available; there is no need to specify a drop tolerance. Our numerical results show that the number of conjugate gradient iterations and the computing time are reduced dramatically for small values of p. We also show that in contrast with drop tolerance strategies, the new approach is more stable in terms of number of iterations and memory requirements.

Distributed Schur Complement Techniques for General Sparse Linear Systems

Abstract: This paper presents a few preconditioning techniques for solving general sparse linear systems on distributed memory environments. These techniques utilize the Schur complement system for deriving the preconditioning matrix in a number of ways. Two of these preconditioners consist of an approximate solution process for the global system, which exploits approximate LU factorizations for diagonal blocks of the Schur complement. Another preconditioner uses a sparse approximate-inverse technique to obtain certain local approximations of the Schur complement. Comparisons are reported for systems of varying difficulty.

Orderings for Incomplete Factorization Preconditioning of Nonsymmetric Problems

Abstract: Numerical experiments are presented whereby the effect of reorderings on the convergence of preconditioned Krylov subspace methods for the solution of nonsymmetric linear systems is shown. The preconditioners used in this study are different variants of incomplete factorizations. It is shown that certain reorderings for direct methods, such as reverse Cuthill--McKee, can be very beneficial. The benefit can be seen in the reduction of the number of iterations and also in measuring the deviation of the preconditioned operator from the identity.

Solving Hyperbolic PDEs Using Interpolating Wavelets

Abstract: A method is presented for adaptively solving hyperbolic PDEs. The method is based on an interpolating wavelet transform using polynomial interpolation on dyadic grids. The adaptability is performed automatically by thresholding the wavelet coefficients. Operations such as differentiation and multiplication are fast and simple due to the one-to-one correspondence between point values and wavelet coefficients in the interpolating basis. Treatment of boundary conditions is simplified in this sparse point representation (SPR). Numerical examples are presented for one- and two-dimensional problems. It is found that the proposed method outperforms a finite difference method on a uniform grid for certain problems in terms of flops.

Moving Mesh Strategy Based on a Gradient Flow Equation for Two-Dimensional Problems

Abstract: In this paper we introduce a moving mesh method for solving PDEs in two dimensions. It can be viewed as a higher-dimensional generalization of the moving mesh PDE (MMPDE) strategy developed in our previous work for one-dimensional problems [W. Huang, Y. Ren, and R. D. Russell, SIAM J. Numer. Anal., 31 (1994), pp. 709--730]. The MMPDE is derived from a gradient flow equation which arises using a mesh adaptation functional in turn motivated from the theory of harmonic maps. Geometrical interpretations are given for the gradient equation and functional, and basic properties of this MMPDE are discussed. Numerical examples are presented where the method is used both for mesh generation and for solving time-dependent PDEs. The results demonstrate the potential of the mesh movement strategy to concentrate the mesh points so as to adapt to special problem features and to also preserve a suitable level of mesh orthogonality.

Low-Rank Matrix Approximation Using the Lanczos Bidiagonalization Process with Applications

Abstract: Low-rank approximation of large and/or sparse matrices is important in many applications, and the singular value decomposition (SVD) gives the best low-rank approximations with respect to unitarily-invariant norms In this paper we show that good low-rank approximations can be directly obtained from the Lanczos bidiagonalization process applied to the given matrix without computing any SVD We also demonstrate that a so-called one-sided reorthogonalization process can be used to maintain an adequate level of orthogonality among the Lanczos vectors and produce accurate low-rank approximations This technique reduces the computational cost of the Lanczos bidiagonalization process We illustrate the efficiency and applicability of our algorithm using numerical examples from several applications areas

The Regular Fourier Matrices and Nonuniform Fast Fourier Transforms

Abstract: For any triple of positive integers (m,N,q), the matrix F(m,N,q), called the (m,N,q)-regular Fourier matrix, is defined. The regular Fourier matrices F(m,N,q) are then applied to set up new algorithms for nonuniform fast Fourier transforms. Numerical results show that the accuracies obtained by our algorithms are much better than previously reported results with the same computation complexity. The algorithms require $O(N\cdot\log_2N)$ arithmetic operations, where N is the number of data points.