SIAM Digital Library
 
 
 

SIAM J. on Scientific Computing

BROWSE VOLUMES

Year Range: 

All Journal content published prior to 1997 is part of LOCUS.

Search Issue | RSS Feeds RSS
Previous Issue

2007

Volume 29, Issue 6, pp. 2241-2704


Efficient Numerical Solution of Cahn–Hilliard–Navier–Stokes Fluids in 2D

David Kay and Richard Welford

SIAM J. Sci. Comput. 29, pp. 2241-2257 (17 pages)

Online Publication Date: October 05, 2007

Full Text: | Download PDF

Show Abstract
A finite element discretization of the variable density Cahn–Hilliard–Navier–Stokes system is presented. This system is then decoupled with the Cahn–Hilliard part solved via the nonlinear multigrid method proposed by the authors in [J. Comput. Phys., 212 (2006), pp. 288–304] and the Navier–Stokes part via a multigrid preconditioned GMRES scheme. This new solver is examined numerically and is shown to be efficient with respect to mesh refinement and robust with respect to the problem parameters.

A Discontinuous Galerkin Moving Mesh Method for Hamilton–Jacobi Equations

J. A. Mackenzie and A. Nicola

SIAM J. Sci. Comput. 29, pp. 2258-2282 (25 pages)

Online Publication Date: October 05, 2007

Full Text: | Download PDF

Show Abstract
In this paper we consider the numerical solution of first-order Hamilton–Jacobi equations using the combination of a discontinuous Galerkin finite element method and an adaptive $r$-refinement (mesh movement) strategy. Particular attention is given to the choice of an appropriate adaptivity criterion when the solution becomes discontinuous. Numerical examples in one and two dimensions are presented to demonstrate the effectiveness of the adaptive procedure.

Simultaneous Higher-Order Optical Flow Estimation and Decomposition

Jing Yuan, Christoph Schörr, and Gabriele Steidl

SIAM J. Sci. Comput. 29, pp. 2283-2304 (22 pages)

Online Publication Date: October 05, 2007

Full Text: | Download PDF

Show Abstract
We study the estimation and decomposition of optical flows from highly nonrigid motions. To this end, recent methods from image decomposition into structural and textural parts are combined with variational optical flow estimation. The approaches we suggest amount to minimizing discrete convex functionals using second-order cone programming. Higher-order regularization is necessary in order to accurately recover important flow structure like vortices, and to incorporate key physical properties such as vanishing divergence. For proper discretization, we apply the finite mimetic difference method, which preserves the identities fulfilled by the continuous differential operators. Numerical examples demonstrate the feasibility of the complex approaches.

A Sparse Discretization for Integral Equation Formulations of High Frequency Scattering Problems

Daan Huybrechs and Stefan Vandewalle

SIAM J. Sci. Comput. 29, pp. 2305-2328 (24 pages) | Cited 1 time

Online Publication Date: October 05, 2007

Full Text: | Download PDF

Show Abstract
We consider two-dimensional scattering problems, formulated as an integral equation defined on the boundary of the scattering obstacle. The oscillatory nature of high-frequency scattering problems necessitates a large number of unknowns in classical boundary element methods. In addition, the corresponding discretization matrix of the integral equation is dense. We formulate a boundary element method with basis functions that incorporate the asymptotic behavior of the solution at high frequencies. The method exhibits the effectiveness of asymptotic methods at high frequencies with only few unknowns, but retains accuracy for lower frequencies. New in our approach is that we combine this hybrid method with very effective quadrature rules for oscillatory integrals. As a result, we obtain a sparse discretization matrix for the oscillatory problem. Moreover, numerical experiments indicate that the accuracy of the solution actually increases with increasing frequency. The sparse discretization applies to problems where the phase of the solution can be predicted a priori, for example in the case of smooth and convex scatterers.

On the Performance of Tensor Methods for Solving Ill-Conditioned Problems

Brett W. Bader and Robert B. Schnabel

SIAM J. Sci. Comput. 29, pp. 2329-2351 (23 pages)

Online Publication Date: October 05, 2007

Full Text: | Download PDF

Show Abstract
This paper investigates the performance of tensor methods for solving small- to large-scale systems of nonlinear equations where the Jacobian matrix at the root is ill-conditioned or singular. This condition occurs on many classes of problems, such as identifying or approaching turning points in path-following problems. The singular case has been studied more than the highly ill-conditioned case, for both Newton and tensor methods. It is known that Newton-based methods do not work well with singular problems because they converge linearly to the solution and, in some cases, with poor accuracy. On the other hand, direct tensor methods have performed well on singular problems and have superlinear convergence on such problems under certain conditions. This behavior originates from the use of a special, restricted form of the second-order term included in the local tensor model that provides information lacking in a (nearly) singular Jacobian. With several implementations available for large-scale problems, tensor methods now are capable of solving larger problems. We compare the performance of tensor methods and Newton-based methods for small- to large-scale problems over a range of conditionings, from well-conditioned to ill-conditioned to singular. Previous studies with tensor methods concerned only the ends of this spectrum. Our results show that tensor methods are increasingly superior to Newton-based methods as the problem grows more ill-conditioned.

Numerical Simulation of Acoustic Streaming on Surface Acoustic Wave-driven Biochips

Daniel Köster

SIAM J. Sci. Comput. 29, pp. 2352-2380 (29 pages) | Cited 1 time

Online Publication Date: October 05, 2007

Full Text: | Download PDF

Show Abstract
A novel type of microfluidic biochip employs acoustic surface waves as a pumping mechanism for fluids. Piezoelectric substrate materials in the chip efficiently convert an electromagnetic input signal into an elastic wave confined to the surface layer of the substrate, hence the surface acoustic wave. The interaction of these waves with an adjoining fluid volume leads to streaming patterns in the fluid or of the motion of the fluid as a whole. In this work we study this acoustic streaming effect. We present a two-scale mathematical model, some theoretical results on the well-posedness of the model problem, as well as numerical simulations based on the finite element method. The model will include the case of free capillary fluid boundaries.

Adaptive Semidiscrete Central-Upwind Schemes for Nonconvex Hyperbolic Conservation Laws

Alexander Kurganov, Guergana Petrova, and Bojan Popov

SIAM J. Sci. Comput. 29, pp. 2381-2401 (21 pages)

Online Publication Date: October 17, 2007

Full Text: | Download PDF

Show Abstract
We discover that the choice of a piecewise polynomial reconstruction is crucial in computing solutions of nonconvex hyperbolic (systems of) conservation laws. Using semidiscrete central-upwind schemes, we illustrate that the obtained numerical approximations may fail to converge to the unique entropy solution or the convergence may be so slow that achieving a proper resolution would require the use of (almost) impractically fine meshes. For example, in the scalar case, all computed solutions seem to converge to solutions that are entropy solutions for some entropy pairs. However, in most applications, one is interested in capturing the unique (Kruzhkov) solution that satisfies the entropy condition for all convex entropies. We present a number of numerical examples that demonstrate the convergence of the solutions, computed with the dissipative second-order minmod reconstruction, to the unique entropy solution. At the same time, more compressive and/or higher-order reconstructions may fail to resolve composite waves, typically present in solutions of nonconvex conservation laws, and thus may fail to recover the Kruzhkov solution. In this paper, we propose a simple and computationally inexpensive adaptive strategy that allows us to simultaneously capture the unique entropy solution and to achieve a high resolution of the computed solution. We use the dissipative minmod reconstruction near the points where convexity changes and utilize a fifth-order weighted essentially nonoscillatory (WENO5) reconstruction in the rest of the computational domain. Our numerical examples (for one- and two-dimensional scalar and systems of conservation laws) demonstrate the robustness, reliability, and nonoscillatory nature of the proposed adaptive method.

Optimized Multiplicative, Additive, and Restricted Additive Schwarz Preconditioning

A. St-Cyr, M. J. Gander, and S. J. Thomas

SIAM J. Sci. Comput. 29, pp. 2402-2425 (24 pages) | Cited 4 times

Online Publication Date: October 17, 2007

Full Text: | Download PDF

Show Abstract
We demonstrate that a small modification of the multiplicative, additive, and restricted additive Schwarz preconditioner at the algebraic level, motivated by optimized Schwarz methods defined at the continuous level, leads to a significant reduction in the iteration count of the iterative solver. Numerical experiments using finite difference and spectral element discretizations of the positive definite Helmholtz problem and an idealized climate simulation illustrate the effectiveness of the new approach.

On Restart and Error Estimation for Krylov Approximation of $w=f(A)v$

Hillel Tal-Ezer

SIAM J. Sci. Comput. 29, pp. 2426-2441 (16 pages)

Online Publication Date: October 17, 2007

Full Text: | Download PDF

Show Abstract
Krylov algorithms are widely used for the solution of large scale linear systems. A Krylov approach is implemented also for the approximation of the vector $w$ which results from operating with a function of a matrix on a vector $v$. For example, when solving a set of linear ODEs we have $w = f (A)v$, where $f(z) = {\rm exp}(tz)$. The main drawback of the Krylov algorithm lies in the need to store all the vectors spanning the approximation space. When solving linear systems, it is possible to overcome this drawback by restarting. In this paper we show how to apply restarts in the general case of approximating $w=f(A)v$. This new algorithm allows a much more efficient approximation of the exponential operator than the standard Krylov algorithm, and it is especially useful in the case of functions which cannot be factored into a product of functions. Yet another interesting subject is the error estimate. The scheme given here provides error estimates “for free,” a trait which enables us to stop the algorithm whenever the desired accuracy is achieved.

A Numerical Bifurcation Analysis of the Ornstein–Zernike Equation with Hypernetted Chain Closure

R. E. Beardmore, A. T. Peplow, and F. Bresme

SIAM J. Sci. Comput. 29, pp. 2442-2463 (22 pages)

Online Publication Date: October 17, 2007

Full Text: | Download PDF

Show Abstract
We study the codimension-one and -two bifurcations of the Ornstein–Zernike equation with hypernetted chain (HNC) closure with Lennard–Jones intermolecular interaction potential. The main purpose of the paper is to present the results of a numerical study undertaken using a suite of algorithms implemented in MATLAB and based on pseudo arc-length continuation for the codimension-one case and a Newton-GMRES method for the codimension-two case. Through careful consideration of the results of our computations, an argument is formulated which shows that spinodal isothermal solution branches arising in this model cannot be reproduced numerically. Furthermore, we show that the existence of an upper bound on the density that can be realized on a vapor isothermal solution branch, which must be present at a spinodal, causes the existence of at least one fold bifurcation along that vapor branch when density is used as the bifurcation parameter. This provides an explanation for previous inconclusive attempts to compute solutions using Newton–Picard methods that are popular in the physical chemistry literature.

Fast Computation of Fourier Integral Operators

Emmanuel Candès, Laurent Demanet, and Lexing Ying

SIAM J. Sci. Comput. 29, pp. 2464-2493 (30 pages) | Cited 9 times

Online Publication Date: October 17, 2007

Full Text: | Download PDF

Show Abstract
We introduce a general purpose algorithm for rapidly computing certain types of oscillatory integrals which frequently arise in problems connected to wave propagation, general hyperbolic equations, and curvilinear tomography. The problem is to numerically evaluate a so-called Fourier integral operator (FIO) of the form $\int e^{2\pi i \Phi(x,\xi)} a(x,\xi) \hat{f}(\xi) \d\xi$ at points given on a Cartesian grid. Here, $\xi$ is a frequency variable, $\hat f(\xi)$ is the Fourier transform of the input $f$, $a(x,\xi)$ is an amplitude, and $\Phi(x,\xi)$ is a phase function, which is typically as large as $|\xi|$; hence the integral is highly oscillatory. Because a FIO is a dense matrix, a naive matrix vector product with an input given on a Cartesian grid of size $N$ by $N$ would require $O(N^4)$ operations. This paper develops a new numerical algorithm which requires $O(N^{2.5} \log N)$ operations and as low as $O(\sqrt{N})$ in storage space (the constants in front of these estimates are small). It operates by localizing the integral over polar wedges with small angular aperture in the frequency plane. On each wedge, the algorithm factorizes the kernel $e^{2 \pi i \Phi(x,\xi)} a(x,\xi)$ into two components: (1) a diffeomorphism which is handled by means of a nonuniform FFT and (2) a residual factor which is handled by numerical separation of the spatial and frequency variables. The key to the complexity and accuracy estimates is the fact that the separation rank of the residual kernel is provably independent of the problem size. Several numerical examples demonstrate the numerical accuracy and low computational complexity of the proposed methodology. We also discuss the potential of our ideas for various applications such as reflection seismology.

Numerical Solution to the Sine-Gordon Equation Defined on the Whole Real Axis

Chunxiong Zheng

SIAM J. Sci. Comput. 29, pp. 2494-2506 (13 pages) | Cited 1 time

Online Publication Date: October 17, 2007

Full Text: | Download PDF

Show Abstract
Numerical simulation of the solution to the sine-Gordon equation on the whole real axis is considered in this paper. Based on nonlinear spectral analysis, exact nonreflecting boundary conditions are derived at two artificially introduced boundary points. Then a numerical scheme of second order is proposed to approximate the solution. In the end, some numerical examples are provided to demonstrate the effectiveness of the proposed scheme.

Implementing Generating Set Search Methods for Linearly Constrained Minimization

Robert Michael Lewis, Anne Shepherd, and Virginia Torczon

SIAM J. Sci. Comput. 29, pp. 2507-2530 (24 pages) | Cited 3 times

Online Publication Date: October 17, 2007

Full Text: | Download PDF

Show Abstract
We discuss an implementation of a derivative-free generating set search method for linearly constrained minimization with no assumption of nondegeneracy placed on the constraints. The convergence guarantees for generating set search methods require that the set of search directions possesses certain geometrical properties that allow it to approximate the feasible region near the current iterate. In the hard case, the calculation of the search directions corresponds to finding the extreme rays of a cone with a degenerate vertex at the origin, a difficult problem. We discuss here how state-of-the-art computational geometry methods make it tractable to solve this problem in connection with generating set search. We also discuss a number of other practical issues of implementation, such as the careful treatment of equality constraints and the desirability of augmenting the set of search directions beyond the theoretically minimal set. We illustrate the behavior of the implementation on several problems from the CUTEr test suite. We have found it to be successful on problems with several hundred variables and linear constraints.

A Dual Optimization Approach to Inverse Quadratic Eigenvalue Problems with Partial Eigenstructure

Zheng-Jian Bai, Delin Chu, and Defeng Sun

SIAM J. Sci. Comput. 29, pp. 2531-2561 (31 pages) | Cited 2 times

Online Publication Date: October 24, 2007

Full Text: | Download PDF

Show Abstract
The inverse quadratic eigenvalue problem (IQEP) arises in the field of structural dynamics. It aims to find three symmetric matrices, known as the mass, the damping, and the stiffness matrices, such that they are closest to the given analytical matrices and satisfy the measured data. The difficulty of this problem lies in the fact that in applications the mass matrix should be positive definite and the stiffness matrix positive semidefinite. Based on an equivalent dual optimization version of the IQEP, we present a quadratically convergent Newton-type method. Our numerical experiments confirm the high efficiency of the proposed method.

Inexact Uniformization Method for Computing Transient Distributions of Markov Chains

Roger B. Sidje, Kevin Burrage, and Shev MacNamara

SIAM J. Sci. Comput. 29, pp. 2562-2580 (19 pages) | Cited 1 time

Online Publication Date: October 24, 2007

Full Text: | Download PDF

Show Abstract
The uniformization method (also known as randomization) is a numerically stable algorithm for computing transient distributions of a continuous time Markov chain. When the solution is needed after a long run or when the convergence is slow, the uniformization method involves a large number of matrix-vector products. Despite this, the method remains very popular due to its ease of implementation and its reliability in many practical circumstances. Because calculating the matrix-vector product is the most time-consuming part of the method, overall efficiency in solving large-scale problems can be significantly enhanced if the matrix-vector product is made more economical. In this paper, we incorporate a new relaxation strategy into the uniformization method to compute the matrix-vector products only approximately. We analyze the error introduced by these inexact matrix-vector products and discuss strategies for refining the accuracy of the relaxation while reducing the execution cost. Numerical experiments drawn from computer systems and biological systems are given to show that significant computational savings are achieved in practical applications.

A High-Order Solver for the Heat Equation in 1D domains with Moving Boundaries

Shravan K. Veerapaneni and George Biros

SIAM J. Sci. Comput. 29, pp. 2581-2606 (26 pages) | Cited 1 time

Online Publication Date: October 24, 2007

Full Text: | Download PDF

Show Abstract
We describe a fast high-order accurate method for the solution of the heat equation in domains with moving Dirichlet or Neumann boundaries and distributed forces. We assume that the motion of the boundary is prescribed. Our method extends the work of Greengard and Strain [Comm. Pure Appl. Math., XLIII (1990), pp. 949–963]. Our scheme is based on a time-space Chebyshev pseudo-spectral collocation discretization, which is combined with a recursive product quadrature rule to accurately and efficiently approximate convolutions with Green's function for the heat equation. We present numerical results that exhibit up to eighth-order convergence rates. Assuming $N$ time steps and $M$ spatial discretization points, the evaluation of the solution of the heat equation at the same number of points in space-time requires ${\mathcal{O}} ( N M \log M)$ work. Thus, our scheme can be characterized as “fast”; that is, it is work-optimal up to a logarithmic factor.

Implicit Scheme for Hyperbolic Conservation Laws Using Nonoscillatory Reconstruction in Space and Time

Karthikeyan Duraisamy and James D. Baeder

SIAM J. Sci. Comput. 29, pp. 2607-2620 (14 pages) | Cited 1 time

Online Publication Date: October 31, 2007

Full Text: | Download PDF

Show Abstract
The efficiency of high order accurate schemes for the solution of unsteady hyperbolic conservation laws is adversely affected by time-step restrictions that arise from monotonicity requirements. When applied to the solution of problems involving discontinuities, these restrictions render conventional high order implicit time integration schemes impractical. In the present study, a new single step second order implicit scheme that uses nonoscillatory reconstruction in space and time is presented. Both the spatial and temporal limiters are dependent on the evolving solution, and this nonlinearity allows for a circumvention of total variation diminishing bounds. Numerical results on model scalar and vector hyperbolic equations suggest that the scheme holds promise in achieving accurate and unconditionally nonoscillatory solutions.

The Effect of Problem Perturbations on Nonlinear Dynamical Systems and their Reduced-Order Models

Radu Serban, Chris Homescu, and Linda R. Petzold

SIAM J. Sci. Comput. 29, pp. 2621-2643 (23 pages) | Cited 1 time

Online Publication Date: October 31, 2007

Full Text: | Download PDF

Show Abstract
Reduced-order models are used extensively in many areas of science and engineering for simulation, design, and control. Reduction techniques for nonlinear dynamical systems produce models that depend strongly on the nominal set of parameters for which the reduction is carried out. In this paper we address the following two questions: “What is the effect of perturbations in the problem parameters on the output functional of a nonlinear dynamical system?” and “to what extent does the reduced-order model capture this effect?”

Numerical Methods for Two-Parameter Local Bifurcation Analysis of Maps

W. Govaerts, R. Khoshsiar Ghaziani, Yu. A. Kuznetsov, and H. G. E. Meijer

SIAM J. Sci. Comput. 29, pp. 2644-2667 (24 pages) | Cited 1 time

Online Publication Date: November 02, 2007

Full Text: | Download PDF

Show Abstract
We discuss new and improved algorithms for the bifurcation analysis of fixed points and periodic orbits (cycles) of maps and their implementation in matcont, a MATLAB toolbox for continuation and bifurcation analysis of dynamical systems. This includes the numerical continuation of fixed points of iterates of the map with one control parameter, detecting and locating their bifurcation points (i.e., limit point, period-doubling, and Neimark–Sacker) and their continuation in two control parameters, as well as detection and location of all codimension 2 bifurcation points on the corresponding curves. For all bifurcations of codim 1 and 2, the critical normal form coefficients are computed, both numerically with finite directional differences and using symbolic derivatives of the original map. Using a parameter-dependent center manifold reduction, explicit asymptotics are derived for bifurcation curves of double and quadruple period cycles rooted at codim 2 points of cycles with arbitrary period. These asymptotics are implemented into the software and allow one to switch at codim 2 points to the continuation of the double and quadruple period bifurcations. We provide two examples illustrating the developed techniques: a generalized Hénon map and a juvenile/adult competition model from mathematical biology.

Structured Quadratic Inverse Eigenvalue Problem, I. Serially Linked Systems

Moody T. Chu, Nicoletta Del Buono, and Bo Yu

SIAM J. Sci. Comput. 29, pp. 2668-2685 (18 pages)

Online Publication Date: November 02, 2007

Full Text: | Download PDF

Show Abstract
Quadratic pencils arising from applications are often inherently structured. Factors contributing to the structure include the connectivity of elements within the underlying physical system and the mandatory nonnegativity of physical parameters. For physical feasibility, structural constraints must be respected. Consequently, they impose additional challenges on the inverse eigenvalue problems which intend to construct a structured quadratic pencil from prescribed eigeninformation. Knowledge of whether a structured quadratic inverse eigenvalue problem is solvable is interesting in both theory and applications. However, the issue of solvability is problem dependent and has to be addressed structure by structure. This paper considers one particular structure where the elements of the physical system, if modeled as a mass-spring system, are serially linked. The discussion recasts both undamped and damped problems in a framework of inequality systems that can be adapted for numerical computation. Some open questions are described.

Pressure Schur Complement Preconditioners for the Discrete Oseen Problem

Maxim A. Olshanskii and Yuri V. Vassilevski

SIAM J. Sci. Comput. 29, pp. 2686-2704 (19 pages) | Cited 5 times

Online Publication Date: November 02, 2007

Full Text: | Download PDF

Show Abstract
We consider several preconditioners for the pressure Schur complement of the discrete steady Oseen problem. Two of the preconditioners are well known from the literature and the other is new. Supplemented with an appropriate approximate solve for an auxiliary velocity subproblem, these approaches give rise to a family of the block preconditioners for the matrix of the discrete Oseen system. In the paper we critically review possible advantages and difficulties of using various Schur complement preconditioners. We recall existing eigenvalue bounds for the preconditioned Schur complement and prove such for the newly proposed preconditioner. These bounds hold both for LBB stable and stabilized finite elements. Results of numerical experiments for several model two-dimensional and three-dimensional problems are presented. In the experiments we use LBB stable finite element methods on uniform triangular and tetrahedral meshes. One particular conclusion is that in spite of essential improvement in comparison with “simple” scaled mass-matrix preconditioners in certain cases, none of the considered approaches provides satisfactory convergence rates in the case of small viscosity coefficients and a sufficiently complex (e.g., circulating) advection vector field.
Close

Your Recent History Recent History Help

Recently Viewed

  1. Multigrid for High-Dimensional Elliptic Partial Differential Equations on Non-equidistant Grids
  2. Additive and Multiplicative Two-Level Spectral Preconditioning for General Linear Systems
  3. Semi‐Implicit Covolume Method in 3D Image Segmentation
  4. Central WENO Schemes for Hamilton–Jacobi Equations on Triangular Meshes
  5. Constructing Embedded Lattice Rules for Multivariate Integration
© SIAM By using SIAM Publications Online you agree to abide by the Terms and Conditions of Use. Banner art adapted from a figure by Hinke M. Osinga and Bernd Krauskopf (University of Bristol, UK).
close