SIAM Digital Library
 
 
 

SIAM J. on Scientific Computing

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

Search Issue | RSS Feeds RSS
Previous Issue

2005

Volume 26, Issue 6, pp. 1811-2176


Model Problems for the Multigrid Optimization of Systems Governed by Differential Equations

Robert Michael Lewis and Stephen G. Nash

SIAM J. Sci. Comput. 26, pp. 1811-1837 (27 pages) | Cited 18 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
We discuss a multigrid approach to the optimization of systems governed by differential equations. Such optimization problems appear in many applications and are of a different nature than systems of equations. Our approach uses an optimization-based multigrid algorithm in which the multigrid algorithm relies explicitly on nonlinear optimization models as subproblems on coarser grids. Our goal is not to argue for a particular optimization-based multigrid algorithm, but instead to demonstrate how multigrid can be used to accelerate nonlinear programming algorithms. Furthermore, using several model problems we give evidence (both theoretical and numerical) that the optimization setting is well suited to multigrid algorithms. Some of the model problems show that the optimization problem may be more amenable to multigrid than the governing differential equation. In addition, we relate the multigrid approach to more traditional optimization methods as further justification for the use of an optimization-based multigrid algorithm.

Explicit, Time Reversible, Adaptive Step Size Control

Ernst Hairer and Gustaf Söderlind

SIAM J. Sci. Comput. 26, pp. 1838-1851 (14 pages) | Cited 8 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
Adaptive step size control is difficult to combine with geometric numerical integration. As classical step size control is based on "past" information only, time symmetry is destroyed and with it the qualitative properties of the method. In this paper we develop completely explicit, reversible, symmetry-preserving, adaptive step size selection algorithms for geometric numerical integrators such as the Störmer--Verlet method. A new step density controller is proposed and analyzed using backward error analysis and reversible perturbation theory. For integrable reversible systems we show that the resulting adaptive method nearly preserves all action variables and, in particular, the total energy for Hamiltonian systems. It has the same excellent long-term behavior as that obtained when constant steps are used. With variable steps, however, both accuracy and efficiency are greatly improved.

Factorized Banded Inverse Preconditioners for Matrices with Toeplitz Structure

Fu-Rong Lin, Michael K. Ng, and Wai-Ki Ching

SIAM J. Sci. Comput. 26, pp. 1852-1870 (19 pages) | Cited 10 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
In this paper, we study factorized banded inverse preconditioners for matrices with Toeplitz structure. We show that if a Toeplitz matrix T has certain off-diagonal decay property, then the factorized banded inverse preconditioner approximates T-1 accurately, and the spectra of these preconditioned matrices are clustered around 1. In nonlinear image restoration applications, Toeplitz-related systems of the form I + T*DT are required to solve, where D is a positive nonconstant diagonal matrix. We construct inverse preconditionersfor such matrices. Numerical results show that the performance of our proposed preconditioners are superior to that of circulant preconditioners. A two-dimensional nonlinear image restoration example is also presented to demonstrate the effectiveness of the proposed preconditioner.

Acceleration of the Schwarz Method for Elliptic Problems

M. Garbey

SIAM J. Sci. Comput. 26, pp. 1871-1893 (23 pages) | Cited 5 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
In this paper, we present a family of domain decomposition methods based on Aitken-like acceleration of the Schwarz method, which is an iterative procedure with linear rate of convergence. This paper is a generalization of our method first introduced in 2000 that was restricted to Cartesian grids. We consider the finite volume (FV) approximation on general quadrangle meshes introduced by Faille in 1992. We first present the so-called Aitken--Schwarz procedure for a linear differential operator. This solver is a direct solver, but its computational cost in the general case might be prohibitive. Second, we introduce the Steffensen--Schwarz variant, which is an iterative domain decomposition solver that can be applied to linear and nonlinear problems. We show that this iterative solver has reasonable numerical efficiency and can be applied successfully to several classes of linear and nonlinear elliptic problems. However, the salient feature of our method is that our algorithm has high tolerance for slow networks in the context of distributed parallel computing and is attractive to use with computer architectures for which performance is limited by the memory bandwidth rather than the flop performance of the CPU. This is currently the case for most parallel computers using RISC processor architecture, such as Beowulf systems.

An Adaptive Wavelet Collocation Method for Fluid-Structure Interaction at High Reynolds Numbers

Nicholas K. R. Kevlahan and Oleg V. Vasilyev

SIAM J. Sci. Comput. 26, pp. 1894-1915 (22 pages) | Cited 24 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
Two mathematical approaches are combined to calculate high Rey\-nolds number incompressible fluid-structure interaction: a wavelet method to dynamically adapt the computational grid to flow intermittency and obstacle motion, and Brinkman penalization to enforce solid boundaries of arbitrary complexity. We also implement a wavelet-based multilevel solver for the Poisson problem for the pressure at each time step. The method is applied to two-dimensional flow around fixed and moving cylinders for Reynolds numbers in the range $3\times 10^1 \le Re \le 10^5$. The compression ratios of up to 1000 are achieved. For the first time it is demonstrated in actual dynamic simulations that the compression scales like $Re^{1/2}$ over five orders of magnitude, while computational complexity scales like $Re$. This represents a significant improvement over the classical complexity estimate of $Re^{9/4}$ for two-dimensional turbulence.

Computation of a Test Statistic in Data Quality Control

Xiao-Wen Chang, Christopher C. Paige, and Christian C. J. M. Tiberius

SIAM J. Sci. Comput. 26, pp. 1916-1931 (16 pages)

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
When processing observational data, statistical testing is an essential instrument for rendering harmless incidental anomalies and disturbances in the measurements. A commonly used test statistic based on the general linear model is the generalized likelihood ratio test statistic. The standard formula given in the literature for this test statistic is not defined if the noise covariance matrix is singular, and is not suitable for computation if any of the matrices involved are ill-conditioned. Based on Paige's generalized linear least squares method [Comm. Statist. B---Simulation Comput., 7 (1978), pp. 437--453], a numerically stable approach is proposed for the computation of the test statistic, as well as for the estimates of the parameter vectors, and reliable representations of the error covariance matrices for these estimates are presented. This approach allows the noise covariance matrix to be singular and can be applied directly to the linear model with linear equality constraints.

Numerical Normal Forms for Codim 2 Bifurcations of Fixed Points with at Most Two Critical Eigenvalues

Yu. A. Kuznetsov and H. G. E. Meijer

SIAM J. Sci. Comput. 26, pp. 1932-1954 (23 pages) | Cited 11 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
In this paper we derive explicit formulas for the normal form coefficients to verify the nondegeneracy of eight codimension two bifurcations of fixed points with one or two critical eigenvalues. These include all strong resonances, as well as the degenerate flip and Neimark--Sacker bifurcations. Applying our results to n-dimensional maps, one avoids the computation of the center manifold, but one can directly evaluate the critical normal form coefficients in the original basis. The formulas remain valid also for n=2 and allow one to avoid the transformation of the linear part of the map into Jordan form.
The developed techniques are tested on two examples: (1) a three-dimensional map appearing in adaptive control; (2) a periodically forced epidemic model. We compute numerically the critical normal form coefficients for several codim 2 bifurcations occurring in these models. To compute the required derivatives of the Poincaré map for the epidemic model, the automatic differentiation package ADOL-C is used.

Accurate Sum and Dot Product

Takeshi Ogita, Siegfried M. Rump, and Shin'ichi Oishi

SIAM J. Sci. Comput. 26, pp. 1955-1988 (34 pages) | Cited 27 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
Algorithms for summation and dot product of floating-point numbers are presented which are fast in terms of measured computing time. We show that the computed results are as accurate as if computed in twice or K-fold working precision, $K\ge 3$. For twice the working precision our algorithms for summation and dot product are some 40% faster than the corresponding XBLAS routines while sharing similar error estimates. Our algorithms are widely applicable because they require only addition, subtraction, and multiplication of floating-point numbers in the same working precision as the given data. Higher precision is unnecessary, algorithms are straight loops without branch, and no access to mantissa or exponent is necessary.

GMRES Convergence Analysis for a Convection-Diffusion Model Problem

J. Liesen and Z. Strakos

SIAM J. Sci. Comput. 26, pp. 1989-2009 (21 pages) | Cited 5 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
When GMRES [Y. Saad and M. H. Schultz, SIAM J. Sci. Statist. Comput.}, 7 (1986), pp. 856--869] is applied to streamline upwind Petrov--Galerkin (SUPG) discretized convection-diffusion problems, it typically exhibits an initial period of slow convergence followed by a faster decrease of the residual norm. Several approaches were made to understand this behavior. However, the existing analyses are solely based on the matrix of the discretized system and they do not take into account any influence of the right-hand side (determined by the boundary conditions and/or source term in the PDE). Therefore they cannot explain the length of the initial period of slow convergence which is right-hand side dependent.
We concentrate on a frequently used model problem with Dirichlet boundary conditions and with a constant velocity field parallel to one of the axes. Instead of the eigendecomposition of the system matrix, which is ill conditioned, we use its orthogonal transformation into a block-diagonal matrix with nonsymmetric tridiagonal Toeplitz blocks and offer an explanation of GMRES convergence. We show how the initial period of slow convergence is related to the boundary conditions and address the question why the convergence in the second stage accelerates.

A Fourth-Order Time-Splitting Laguerre--Hermite Pseudospectral Method for Bose--Einstein Condensates

Weizhu Bao and Jie Shen

SIAM J. Sci. Comput. 26, pp. 2010-2028 (19 pages) | Cited 24 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
A fourth-order time-splitting Laguerre--Hermite pseudospectral method is introduced for Bose--Einstein condensates (BECs) in three dimensions with cylindrical symmetry. The method is explicit, time reversible, and time transverse invariant. It conserves the position density and is spectral accurate in space and fourth-order accurate in time. Moreover, the new method has two other important advantages: (i) it reduces a three-dimensional (3-D) problem with cylindrical symmetry to an effective two-dimensional (2-D) problem; (ii) it solves the problem in the whole space instead of in a truncated artificial computational domain. The method is applied to vector Gross--Pitaevskii equations (VGPEs) for multicomponent BECs. Extensive numerical tests are presented for the one-dimensional (1-D) GPE, the 2-D GPE with radial symmetry, the 3-D GPE with cylindrical symmetry, as well as 3-D VGPEs for two-component BECs, to show the efficiency and accuracy of the new numerical method.

Finite Element Method for Epitaxial Growth with Thermodynamic Boundary Conditions

Eberhard Bänsch, Frank Hausser, and Axel Voigt

SIAM J. Sci. Comput. 26, pp. 2029-2046 (18 pages) | Cited 3 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
We develop an adaptive finite element method for island dynamics in epitaxial growth. We study a step-flow model, which consists of an adatom (adsorbed atom) diffusion equation on terraces of different height; thermodynamic boundary conditions on terrace boundaries including anisotropic line tension; and the normal velocity law for the motion of such boundaries determined by a two-sided flux, together with the one-dimensional anisotropic ``surface' diffusion (edge diffusion) of edge adatoms along the step edges. The problem is solved using independent meshes: a two-dimensional mesh for the adatom diffusion and one-dimensional meshes for the boundary evolution. A penalty method is used to incorporate the boundary conditions. The evolution of the terrace boundaries includes both the weighted/anisotropic mean curvature flow and the weighted/anisotropic edge diffusion. Its governing equation is solved by a semi-implicit front-tracking method using parametric finite elements.

Error Minimization of Multipole Expansion

Shinichiro Ohnuki and Weng Cho Chew

SIAM J. Sci. Comput. 26, pp. 2047-2065 (19 pages) | Cited 1 time

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
In this paper, we focus on the truncation error of the multipole expansion for the fast multipole method and the multilevel fast multipole algorithm. When the buffer size is large enough, the error can be controlled and minimized by using the conventional selection rules. On the other hand, if the buffer size is small, the conventional selection rules no longer hold, and the new approach which we have recently proposed is needed. However, this method is still not sufficient to minimize the error for small buffer cases. We clarify this fact and show that the information about the placement of true worst-case interaction is needed. A novel algorithm to minimize the truncation error is proposed.

A New Distillation Algorithm for Floating-Point Summation

Yong-Kang Zhu, Jun-Hai Yong, and Guo-Qin Zheng

SIAM J. Sci. Comput. 26, pp. 2066-2078 (13 pages) | Cited 5 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
The summation of n floating-point numbers is ubiquitous in numerical computations. We present a new distillation algorithm for floating-point summation which is stable, efficient, and accurate. The algorithm iteratively "distills" the summands without discarding any significant digit until the partial sums cannot change the whole sum. It uses standard floating-point arithmetic and does not rely on the choice of radix or any other specific assumption. Furthermore, the error bound of our algorithm is independent of n and less than 1 ulp.

Two Interface-Type Numerical Methods for Computing Hyperbolic Systems with Geometrical Source Terms Having Concentrations

Shi Jin and Xin Wen

SIAM J. Sci. Comput. 26, pp. 2079-2101 (23 pages) | Cited 8 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
We propose two simple well-balanced methods for hyperbolic systems with geometrical source terms having concentrations. Physical problems under consideration include the shallow water equations with topography and the quasi-one-dimensional isothermal nozzle flows. These two methods use the numerical fluxes already obtained from the corresponding homogeneous systems in the source terms, and one needs only a black-box (approximate) Riemann solver for the homogeneous system. Compared with our previous method developed in [S. Jin and X. Wen, J. Comput. Math., 22 (2004), pp. 230--249], these methods avoid the Newton iterations in the evaluation of the source term. Numerical experiments demonstrate that both methods give good numerical approximations to the subcritical and supercritical flows. With a transonic fix, both methods also capture with a high resolution the transonic flows over the concentration. These methods are applicable to both unsteady and steady state computations.

Accurate and Efficient Boundary Integral Methods for Electrified Liquid Bridge Problems

Darko Volkov, Demetrios T. Papageorgiou, and Peter G. Petropoulos

SIAM J. Sci. Comput. 26, pp. 2102-2132 (31 pages) | Cited 1 time

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
We derive and implement boundary integral methods for axisymmetric liquid bridge problems in the presence of an axial electric field. The liquid bridge is bounded by solid parallel electrodes placed perpendicular to the axis of symmetry and held at a constant potential difference. The fluid is assumed to be nonconducting and has permittivity different from that of the passive surrounding medium. The problem reduces to the solution of two harmonic problems for the fluid and voltage potential inside the bridge and another harmonic problem for the voltage potential outside the bridge. The shape of the moving interface is determined by the imposition of stress, as well as kinematic and electric field boundary conditions, the former condition accounting for discontinuous electric stresses across the interface. We propose fast and highly accurate boundary integral methods based on fast summations of appropriate series representations of axisymmetric Green's functions in bounded geometries. We implement our method to calculate equilibrium shapes for electrified liquid bridges in the absence and presence of gravity. Such calculations appear in the literature using finite element methods, and our boundary integral approach is a fast and accurate alternative.

Algorithms for Numerical Analysis in High Dimensions

Gregory Beylkin and Martin J. Mohlenkamp

SIAM J. Sci. Comput. 26, pp. 2133-2159 (27 pages) | Cited 64 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
Nearly every numerical analysis algorithm has computational complexity that scales exponentially in the underlying physical dimension. The separated representation, introduced previously, allows many operations to be performed with scaling that is formally linear in the dimension. In this paper we further develop this representation by
(i) discussing the variety of mechanisms that allow it to be surprisingly efficient;
(ii) addressing the issue of conditioning;
(iii) presenting algorithms for solving linear systems within this framework; and
(iv) demonstrating methods for dealing with antisymmetric functions, as arise in the multiparticle Schrödinger equation in quantum mechanics.
Numerical examples are given.

Fast Computation for Large Magnetostatic Systems Adapted for Micromagnetism

Stéphane Labbé

SIAM J. Sci. Comput. 26, pp. 2160-2175 (16 pages) | Cited 3 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
In this paper, an efficient method is developed for computing the magnetostatic field for ferromagnetic materials on large structured meshes. The problem is discretized using a finite volume approximation. The discrete operator is proved to preserve the main properties of the continuous model, and a lower estimate of its lower eigenvalue is given. Using the fact that the discrete operator has a block-Toeplitz structure for cubic meshes in parallelepipedic domains, a fast solving method is built. Based upon the use of fast Fourier transform, this method allows one to reduce the computational cost from n2 to O(n log(n)) but also to reduce the storage to O(n) instead of n2, where n is the number of cells in the mesh.

Erratum: Spartan Gibbs Random Field Models for Geostatistical Applications

Dionissios T. Hristopulos

SIAM J. Sci. Comput. 26, pp. 2176-2176 (1 page)

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
Typographical errors in certain mathematical formulas that appear in the article "Spartan Gibbs Random Field Models for Geostatistical Applications" [SIAM J. Sci. Comput., 24 (2003), pp. 2125--2162] are corrected.
Close

Your Recent History Recent History Help

Recently Viewed

  1. A New Error Estimate of the Fast Gauss Transform
  2. A Preconditioner for the Steady-State Navier--Stokes Equations
  3. Iterative Solution of Cyclically Reduced Systems Arising from Discretization of the Three-Dimensiona...
  4. Nonoscillatory Central Schemes for Multidimensional Hyperbolic Conservation Laws
  5. A Fast Multilevel Algorithm for the Solution of Nonlinear Systems of Conductive-Radiative Heat Trans...
© 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