SIAM Digital Library
 
 
 

SIAM J. on Scientific Computing

Year Range: 

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

Search Issue | RSS Feeds RSS
Previous Issue Next Issue

2007

Volume 29, Issue 4, pp. 1355-1824


Numerical Methods for Computing Nonlinear Eigenpairs: Part I. Iso-Homogeneous Cases

Xudong Yao and Jianxin Zhou

SIAM J. Sci. Comput. 29, pp. 1355-1374 (20 pages) | Cited 4 times

Online Publication Date: June 12, 2007

Full Text: | Download PDF

Show Abstract
With a Rayleigh quotient formulation, a local minimax method is developed to solve a class of (iso-homogeneous) nonlinear eigenpair problems for multiple solutions in Banach spaces following their instability order. The algorithm is implemented to compute (weighted) eigenpairs of the p-Laplacian. Numerical eigenfunctions are illustrated by their graphics. Several interesting phenomena have been observed and are open for further investigation. Mathematical analysis related to convergence and instability order of computed eigenfunctions, etc., is also presented.

Regularity Theory and Superalgebraic Solvers for Wire Antenna Problems

Oscar P. Bruno and Michael C. Haslam

SIAM J. Sci. Comput. 29, pp. 1375-1402 (28 pages) | Cited 3 times

Online Publication Date: June 12, 2007

Full Text: | Download PDF

Show Abstract
We consider the problem of evaluating the current distribution $J(z)$ that is induced on a straight wire antenna by a time-harmonic incident electromagnetic field. The scope of this paper is twofold. One of its main contributions is a regularity proof for a straight wire occupying the interval $[-1,1]$. In particular, for a smooth time-harmonic incident field this theorem implies that $J(z) = I(z)/\sqrt{1-z^2}$, where $I(z)$ is an infinitely differentiable function—the previous state of the art in this regard placed $I$ in the Sobolev space $W^{1,p}$, $p>1$. The second focus of this work is on numerics: we present three superalgebraically convergent algorithms for the solution of wire problems, two based on Hallén's integral equation and one based on the Pocklington integrodifferential equation. Both our proof and our algorithms are based on two main elements: (1) a new decomposition of the kernel of the form $G(z) = F_1(z) \ln\! |z| + F_2(z)$, where $F_1(z)$ and $F_2(z)$ are analytic functions on the real line; and (2) removal of the end-point square root singularities by means of a coordinate transformation. The Hallén- and Pocklington-based algorithms we propose converge superalgebraically: faster than $\mathcal{O}(N^{-m})$ and $\mathcal{O}(M^{-m})$ for any positive integer $m$, where $N$ and $M$ are the numbers of unknowns and the number of integration points required for construction of the discretized operator, respectively. In previous studies, at most the leading-order contribution to the logarithmic singular term was extracted from the kernel and treated analytically, the higher-order singular derivatives were left untreated, and the resulting integration methods for the kernel exhibit $\mathcal{O}(M^{-3})$ convergence at best. A rather comprehensive set of tests we consider shows that, in many cases, to achieve a given accuracy, the numbers $N$ of unknowns required by our codes are up to a factor of five times smaller than those required by the best solvers previously available; the required number $M$ of integration points, in turn, can be several orders of magnitude smaller than those required in previous methods. In particular, four-digit solutions were found in computational times of the order of four seconds and, in most cases, of the order of a fraction of a second on a contemporary personal computer; much higher accuracies result in very small additional computing times.

Stability Results for Scattered Data Interpolation by Trigonometric Polynomials

Stefan Kunis and Daniel Potts

SIAM J. Sci. Comput. 29, pp. 1403-1419 (17 pages) | Cited 4 times

Online Publication Date: June 12, 2007

Full Text: | Download PDF

Show Abstract
A fast and reliable algorithm for the optimal interpolation of scattered data on the torus $\mathbb{T}^d$ by multivariate trigonometric polynomials is presented. The algorithm is based on a variant of the conjugate gradient method in combination with the fast Fourier transforms for nonequispaced nodes. The main result is that under mild assumptions the total complexity for solving the interpolation problem at $M$ arbitrary nodes is of order ${\cal O}(M\log M)$. This result is obtained by the use of localized trigonometric kernels where the localization is chosen in accordance with the spatial dimension $d$. Numerical examples show the efficiency of the new algorithm.

A Fast Algorithm for the Calculation of the Roots of Special Functions

Andreas Glaser, Xiangtao Liu, and Vladimir Rokhlin

SIAM J. Sci. Comput. 29, pp. 1420-1438 (19 pages) | Cited 9 times

Online Publication Date: June 12, 2007

Full Text: | Download PDF

Show Abstract
We describe a procedure for the determination of the roots of functions satisfying second-order ordinary differential equations, including the classical special functions. The scheme is based on a combination of the Prüfer transform with the classical Taylor series method for the solution of ordinary differential equations and requires $O(1)$ operations for the determination of each root. When the functions in question are classical orthogonal polynomials (Legendre, Hermite, etc.), the techniques presented here also provide tools for the evaluation of the weights for the associated Gaussian quadratures. The performance of the scheme for several classical special functions (prolate spheroidal wave functions, Bessel functions, and Legendre, Hermite, and Laguerre polynomials) is illustrated with numerical examples.

A New Approach to Simulating Flow in Discrete Fracture Networks with an Optimized Mesh

Hussein Mustapha and Kassem Mustapha

SIAM J. Sci. Comput. 29, pp. 1439-1459 (21 pages) | Cited 7 times

Online Publication Date: June 19, 2007

Full Text: | Download PDF

Show Abstract
Natural fractured media are highly unpredictable because of existing complex structures at the fracture and at the network levels. Fractures are by themselves heterogeneous objects of broadly distributed sizes, shapes, and orientations that are interconnected in large correlated networks. With little field data and evidence, numerical modeling can provide important information on the underground hydraulic phenomena. However, it must overcome several barriers. First, the complex network structure produces a structure difficult to mesh. Second, the absence of a priori homogenization scale, along with the double fracture and network heterogeneity levels, requires the calculation of large but finely resolved fracture networks resulting in very large simulation domains. To tackle these two related issues, we reduce the highly complex geometry of the fractures by applying a local transformation that suppresses the cumbersome meshing configurations while keeping the networks fundamental, geological, and geometrical characteristics. We show that the flow properties are marginally affected while the problem complexity (i.e., memory capacity and resolution time) can be divided by orders of magnitude. The goal of this article is to propose a method of resolution which takes into account the geometrical complexity met in the networks and which makes it possible to treat a few thousand fractures. The principal aim of this article is to present a tool to slowly modify the structures of the fracture networks to have a good quality mesh with a marginal loss in precision.

On Global Error Estimation and Control for Initial Value Problems

J. Lang and J. G. Verwer

SIAM J. Sci. Comput. 29, pp. 1460-1475 (16 pages) | Cited 4 times

Online Publication Date: June 19, 2007

Full Text: | Download PDF

Show Abstract
This paper addresses global error estimation and control for initial value problems for ordinary differential equations. The focus lies on a comparison between a novel approach based on the adjoint method combined with a small sample statistical initialization and the classical approach based on the first variational equation. Control is achieved through tolerance proportionality. Both approaches are found to work well and to enable estimation and control in a reliable manner.

Incomplete LU Preconditioning with the Multilevel Fast Multipole Algorithm for Electromagnetic Scattering

Tahİr Malas and Levent Gürel

SIAM J. Sci. Comput. 29, pp. 1476-1494 (19 pages) | Cited 13 times

Online Publication Date: June 21, 2007

Full Text: | Download PDF

Show Abstract
Iterative solution of large-scale scattering problems in computational electromagnetics with the multilevel fast multipole algorithm (MLFMA) requires strong preconditioners, especially for the electric-field integral equation (EFIE) formulation. Incomplete LU (ILU) preconditioners are widely used and available in several solver packages. However, they lack robustness due to potential instability problems. In this study, we consider various ILU-class preconditioners and investigate the parameters that render them safely applicable to common surface integral formulations without increasing the ${\cal O}(n\log n)$ complexity of MLFMA. We conclude that the no-fill ILU(0) preconditioner is an optimal choice for the combined-field integral equation (CFIE). For EFIE, we establish the need to resort to methods depending on drop tolerance and apply pivoting for problems with high condition estimate. We propose a strategy for the selection of the parameters so that the preconditioner can be used as a black-box method. Robustness and efficiency of the employed preconditioners are demonstrated over several test problems.

General Tooth Boundary Conditions for Equation Free Modeling

A. J. Roberts and I. G. Kevrekidis

SIAM J. Sci. Comput. 29, pp. 1495-1510 (16 pages) | Cited 1 time

Online Publication Date: June 29, 2007

Full Text: | Download PDF

Show Abstract
We are developing a framework for multiscale computation which enables models at a “microscopic” level of description, for example, lattice Boltzmann, Monte Carlo, or molecular dynamics simulators, to perform modeling tasks at “macroscopic” length scales of interest. The plan is to use the microscopic rules restricted to small “patches” of the domain, the “teeth,” using interpolation to bridge the “gaps.” Here we explore general boundary conditions coupling the widely separated “teeth” of the microscopic simulation that achieve high order accuracy over the macroscale. We present the simplest case when the microscopic simulator is the quintessential example of a PDE. We argue that classic high order interpolation of the macroscopic field provides the correct forcing in whatever boundary condition is required by the microsimulator. Such interpolation leads to tooth boundary conditions, which achieve arbitrarily high order consistency. The high order consistency is demonstrated on a class of linear PDEs in two ways: first through the eigenvalues of the scheme for selected numerical problems, and second using the dynamical systems approach of holistic discretization on a general class of linear PDEs. Analytic modeling shows that, for a wide class of microscopic systems, the subgrid fields and the effective macroscopic model are largely independent of the tooth size and the particular tooth boundary conditions. When applied to patches of microscopic simulations these tooth boundary conditions promise efficient macroscale simulation. We expect the same approach will also accurately couple patch simulations in higher spatial dimensions.

Dual Functions for a Parallel Adaptive Method

Randolph E. Bank and Jeffrey S. Ovall

SIAM J. Sci. Comput. 29, pp. 1511-1524 (14 pages)

Online Publication Date: June 29, 2007

Full Text: | Download PDF

Show Abstract
In this paper, we investigate the effects of pollution error on the performance of the parallel adaptive finite element technique proposed by Bank and Holst in 2000 [R. E. Bank and M. Holst, SIAM J. Sci. Comput., 22 (2000), pp. 1411–1443]. In particular, we consider whether the performance of the algorithm as it was originally proposed can be improved through the use of certain dual functions which give an indication of the global influences on subdomain errors.

Boundary Preserving Semianalytic Numerical Algorithms for Stochastic Differential Equations

Esteban Moro and Henri Schurz

SIAM J. Sci. Comput. 29, pp. 1525-1549 (25 pages) | Cited 10 times

Online Publication Date: July 04, 2007

Full Text: | Download PDF

Show Abstract
Construction of splitting-step methods and properties of related nonnegativity and boundary preserving semianalytic numerical algorithms for solving stochastic differential equations (SDEs) of Itô type are discussed. As the crucial assumption, we oppose conditions such that one can decompose the original system of SDEs into subsystems for which one knows either the exact solution or its conditional transition probability. We present convergence proofs for a newly designed splitting-step algorithm and simulation studies for numerous well-known numerical examples ranging from stochastic dynamics occurring in asset pricing in mathematical finance (Cox–Ingersoll–Ross (CIR) and constant elasticity of variance (CEV) models) to measure-valued diffusion and super-Brownian motion (stochastic PDEs (SPDEs)) as met in biology and physics.

A Variational Approach to Video Segmentation for Botanical Data

Aaron Luttman and John Bardsley

SIAM J. Sci. Comput. 29, pp. 1550-1566 (17 pages) | Cited 2 times

Online Publication Date: July 04, 2007

Full Text: | Download PDF

Show Abstract
In order to engage in photosynthesis, plant leaves absorb CO$_2$ via the opening of pores in their surfaces called stomata. Water evaporates through open stomata, however, which is a detriment to plant function. Thus a plant will seek a stomatal aperture that balances its need for CO$_2$ with its aversion to H$_2$O loss. In order to visualize a particular leaf's stomatal aperture, an experimentalist injects the leaf with dye so that it fluoresces. Regions with a higher relative intensity that do not correspond to veins in the leaf then correspond to areas in which the stomata are closed and the darker nonvein regions to areas in which the stomata are open. A camera is used to collect the emitted light, and a fluorescence pattern is measured. Images are continually recorded as these patterns change with time, resulting in a video sequence. The goal of this paper is to segment one such video sequence into fluorescing and nonfluorescing regions. After preprocessing the video, we take a variational approach to the segmentation problem and solve the associated three-dimensional evolution equation using a semi-implicit numerical scheme. A fixed-point formulation of the method is presented as are results of the segmentation for actual leaf data. An erratum to this article has been appended at the end of the pdf file.

Boundary Conditions in the Schur Complement Preconditioning of Dissipative Acoustic Equations

Mika Malinen

SIAM J. Sci. Comput. 29, pp. 1567-1592 (26 pages)

Online Publication Date: July 11, 2007

Full Text: | Download PDF

Show Abstract
This paper is concerned with the development of robust and efficient iterative techniques for solving the coupled systems of linear equations which arise from a mixed finite element discretization of dissipative acoustic equations. The system of linearized, time-harmonic acoustic equations we consider is based on the Navier–Stokes, continuity, and energy equations accompanied by suitable equations of state. In this paper, we focus on the design of preconditioners that can be used in combination with the Krylov subspace methods. Our starting point is the Schur complement formulation of the mixed problem. One of the difficulties encountered with the approach we consider is generally associated with the issue of defining boundary conditions for the associated Schur complement problem. To derive suitable boundary conditions for this problem, we first consider an iterative solution method which is based on the direct use of the Schur complement formulation. This development naturally leads us to address the key question of what boundary conditions should be adjoined to the Schur complement preconditioner. We give numerical examples which demonstrate the effectiveness and robustness of the proposed preconditioner. The examples include both two- and three-dimensional problems with boundary conditions of various types.

Additive and Multiplicative Two-Level Spectral Preconditioning for General Linear Systems

B. Carpentieri, L. Giraud, and S. Gratton

SIAM J. Sci. Comput. 29, pp. 1593-1612 (20 pages) | Cited 2 times

Online Publication Date: July 18, 2007

Full Text: | Download PDF

Show Abstract
In this paper we introduce new preconditioning techniques for the solution of general symmetric and unsymmetric linear systems $A x = b$. These approaches borrow some ideas of the multigrid philosophy designed for the solution of linear systems arising from the discretization of elliptic partial differential equations. We attempt to improve the convergence rate of a prescribed preconditioner $M_1$. In a two-grid framework, this preconditioner is viewed as a smoother and the coarse space is spanned by the eigenvectors associated with the smallest eigenvalues of $M_1A$. We derive both additive and multiplicative variants of the resulting iterated two-level preconditioners for unsymmetric linear systems that can also be adapted for Hermitian positive definite problems. We show that these two-level preconditioners shift the smallest eigenvalues to one and tend to better cluster around one those eigenvalues that $M_1$ already succeeded in moving into the neighborhood of one. We illustrate the behavior of our method through extensive numerical experiments on a set of general linear systems. Finally, we show the effectiveness of these approaches on two challenging real applications; the first comes from a nonoverlapping domain decomposition method in semiconductor device modeling, the second from industrial electromagnetism applications.

Multigrid for High-Dimensional Elliptic Partial Differential Equations on Non-equidistant Grids

H. bin Zubair, C. W. Oosterlee, and R. Wienands

SIAM J. Sci. Comput. 29, pp. 1613-1636 (24 pages) | Cited 5 times

Online Publication Date: July 18, 2007

Full Text: | Download PDF

Show Abstract
This work presents techniques, theory, and numbers for multigrid in a general $d$-dimensional setting. The main focus of this paper is the multigrid convergence for high-dimensional partial differential equations on non-equidistant grids such as may be encountered in a sparse-grid solution. As a model problem we have chosen the anisotropic stationary diffusion equation on a rectangular hypercube. We present some techniques for building the general $d$-dimensional adaptations of the multigrid components and propose grid-coarsening strategies to handle anisotropies that are induced due to discretization on a non-equidistant grid. Apart from the practical formulas and techniques, we present—in detail—the smoothing analysis of the point $\omega$-red-black Jacobi method for a general multidimensional case. We show how relaxation parameters may be evaluated efficiently and used for better convergence. This analysis incorporates full and partial doubling and quadrupling coarsening strategies as well as the second- and the fourth-order finite difference operators. Finally we present some results derived from numerical experiments based on the test problem.

A New Implementation of Symplectic Runge–Kutta Methods

Robert I. McLachlan

SIAM J. Sci. Comput. 29, pp. 1637-1649 (13 pages) | Cited 2 times

Online Publication Date: August 01, 2007

Full Text: | Download PDF

Show Abstract
We propose a “Newton–Taylor” iteration for solving the implicit equations of symplectic Runge–Kutta methods, using the Jacobian of the vector field and matrix-vector multiplications whose extra cost for certain structured problems is negligible. The structure of Hamiltonian ODEs allows this very simple iteration to be effective. The iteration reduces the number of vector field evaluations almost to that of Newton's method, often only one or two per time step, making symplectic Runge–Kutta methods more efficient even at relatively large time steps.

Sensitivities in Large Eddy Simulation and Improved Estimates of Turbulent Flow Functionals

Mihai Anitescu and William J. Layton

SIAM J. Sci. Comput. 29, pp. 1650-1667 (18 pages) | Cited 3 times

Online Publication Date: August 01, 2007

Full Text: | Download PDF

Show Abstract
We consider a prototypical problem: simulate velocity and pressure in a turbulent flow using large eddy simulation (LES) and then use the results to estimate the force that the underlying turbulent flow exerts on an immersed body. For eddy viscosity-type LES models we develop the appropriate continuous sensitivity equation and show how it can improve the functional estimate by using the sensitivity with respect to the user-selected length scale to incorporate the effects of the underlying unresolved small-scale fluctuations on the functionals sought.

A Corner Cutting Algorithm for Evaluating Rational Bézier Surfaces and the Optimal Stability of the Basis

Jorge Delgado and J. M. Peña

SIAM J. Sci. Comput. 29, pp. 1668-1682 (15 pages) | Cited 2 times

Online Publication Date: August 01, 2007

Full Text: | Download PDF

Show Abstract
The usual method for evaluating rational Bézier surfaces uses the projection operator and the representation provided by the Bernstein basis. We prove the optimal stability of the basis used in this representation. We also propose an alternative method for evaluating rational surfaces through that representation. We show the stability properties of this last method and prove that it has better properties than other known methods from the point of view of avoiding underflow and overflow.

Partitioning Sparse Matrices for Parallel Preconditioned Iterative Methods

Bora Uçar and Cevdet Aykanat

SIAM J. Sci. Comput. 29, pp. 1683-1709 (27 pages) | Cited 7 times

Online Publication Date: August 01, 2007

Full Text: | Download PDF

Show Abstract
This paper addresses the parallelization of the preconditioned iterative methods that use explicit preconditioners such as approximate inverses. Parallelizing a full step of these methods requires the coefficient and preconditioner matrices to be well partitioned. We first show that different methods impose different partitioning requirements for the matrices. Then we develop hypergraph models to meet those requirements. In particular, we develop models that enable us to obtain partitionings on the coefficient and preconditioner matrices simultaneously. Experiments on a set of unsymmetric sparse matrices show that the proposed models yield effective partitioning results. A parallel implementation of the right preconditioned BiCGStab method on a PC cluster verifies that the theoretical gains obtained by the models hold in practice.

Fast Directional Multilevel Algorithms for Oscillatory Kernels

Björn Engquist and Lexing Ying

SIAM J. Sci. Comput. 29, pp. 1710-1737 (28 pages) | Cited 16 times

Online Publication Date: August 10, 2007

Full Text: | Download PDF

Show Abstract
This paper introduces a new directional multilevel algorithm for solving $N$-body or $N$-point problems with highly oscillatory kernels. These systems often result from the boundary integral formulations of scattering problems and are difficult due to the oscillatory nature of the kernel and the non-uniformity of the particle distribution. We address the problem by first proving that the interaction between a ball of radius $r$ and a well-separated region has an approximate low rank representation, as long as the well-separated region belongs to a cone with a spanning angle of $O(1/r)$ and is at a distance which is at least $O(r^2)$ away from from the ball. We then propose an efficient and accurate procedure which utilizes random sampling to generate such a separated, low rank representation. Based on the resulting representations, our new algorithm organizes the high frequency far field computation by a multidirectional and multiscale strategy to achieve maximum efficiency. The algorithm performs well on a large group of highly oscillatory kernels. Our algorithm is proved to have $O(N\log N)$ computational complexity for any given accuracy when the points are sampled from a two dimensional surface. We also provide numerical results to demonstrate these properties.

A Fast Shadowing Algorithm for High-Dimensional ODE Systems

Wayne B. Hayes and Kenneth R. Jackson

SIAM J. Sci. Comput. 29, pp. 1738-1758 (21 pages) | Cited 1 time

Online Publication Date: August 22, 2007

Full Text: | Download PDF

Show Abstract
Numerical solutions to chaotic dynamical systems are suspect because the exponential divergence of nearby trajectories causes numerical errors to be exponentially magnified with time. To bolster confidence in the reliability of such numerical solutions, we turn to the study of shadowing: a shadow is an exact trajectory of a chaotic dynamical system that stays close to a numerical trajectory for a long time. Finding shadows of numerical trajectories is very computationally intensive, and until recently it has been infeasible to study shadows of higher-dimensional systems. This paper introduces several optimizations to previous algorithms, which collectively achieve an average speedup of almost two orders of magnitude with no measurable loss in effectiveness. We test the algorithm on systems with up to 180 dimensions. Its application to large gravitational $N$-body integrations has already led to a deeper understanding of the reliability of galaxy and cosmological simulations. The source code (available from the first author) provides to the experienced user a generic interface capable of searching for shadows of any set of ordinary differential equations.

Three-Level BDDC in Three Dimensions

Xuemin Tu

SIAM J. Sci. Comput. 29, pp. 1759-1780 (22 pages) | Cited 11 times

Online Publication Date: September 19, 2007

Full Text: | Download PDF

Show Abstract
Balancing domain decomposition by constraints (BDDC) methods are nonoverlapping iterative substructuring domain decomposition methods for the solution of large sparse linear algebraic systems arising from the discretization of elliptic boundary value problems. Their coarse problems are given in terms of a small number of continuity constraints for each subdomain, which are enforced across the interface. The coarse problem matrix is generated and factored by a direct solver at the beginning of the computation and it can ultimately become a bottleneck if the number of subdomains is very large. In this paper, two three-level BDDC methods are introduced for solving the coarse problem approximately for problems in three dimensions. This is an extension of previous work for the two-dimensional case. Edge constraints are considered in this work since vertex constraints alone, which work well in two dimensions, result in a noncompetitive algorithm in three dimensions. Some new technical tools are then needed in the analysis and this makes the three-dimensional case more complicated. Estimates of the condition numbers are provided for two three-level BDDC methods, and numerical experiments are also discussed.

A Cell-Centered Lagrangian Scheme for Two-Dimensional Compressible Flow Problems

Pierre-Henri Maire, Rémi Abgrall, Jérôme Breil, and Jean Ovadia

SIAM J. Sci. Comput. 29, pp. 1781-1824 (44 pages) | Cited 67 times

Online Publication Date: September 19, 2007

Full Text: | Download PDF

Show Abstract
We present a new Lagrangian cell-centered scheme for two-dimensional compressible flows. The primary variables in this new scheme are cell-centered, i.e., density, momentum, and total energy are defined by their mean values in the cells. The vertex velocities and the numerical fluxes through the cell interfaces are not computed independently, contrary to standard approaches, but are evaluated in a consistent manner due to an original solver located at the nodes. The main new feature of the algorithm is the introduction of four pressures on each edge, two for each node on each side of the edge. This extra degree of freedom allows us to construct a nodal solver which fulfills two properties. First, the conservation of momentum and total energy is ensured. Second, a semidiscrete entropy inequality is provided. In the case of a one-dimensional flow, the solver reduces to the classical Godunov acoustic solver: it can be considered as its two-dimensional generalization. Many numerical tests are presented. They are representative test cases for compressible flows and demonstrate the robustness and the accuracy of this new solver.
Close

close