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

2006

Volume 27, Issue 6, pp. 1803-2182


Computational Study of Fast Methods for the Eikonal Equation

Pierre A. Gremaud and Christopher M. Kuster

SIAM J. Sci. Comput. 27, pp. 1803-1816 (14 pages) | Cited 15 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
A computational study of the fast marching and the fast sweeping methods for the eikonal equation is given. It is stressed that both algorithms should be considered as "direct" (as opposed to iterative) methods. On realistic grids, fast sweeping is faster than fast marching for problems with simple geometry. For strongly nonuniform problems and/or complex geometry, the situation may be reversed. Finally, fully second order generalizations of methods of this type for problems with obstacles are proposed and implemented.

Composition Methods for Differential Equations with Processing

S. Blanes, F. Casas, and A. Murua

SIAM J. Sci. Comput. 27, pp. 1817-1843 (27 pages) | Cited 3 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
We construct numerical integrators for differential equations up to order 12 obtained by composition of basic integrators. The following cases are considered: (i) composition for a system separable in two solvable parts, (ii) composition using as basic methods a first-order integrator and its adjoint, (iii) composition using second-order symmetric methods, and (iv) composition using fourth-order symmetric methods. Each scheme is implemented with a processor or corrector to improve their efficiency, and this can be done virtually cost-free.

Isoperimetric Partitioning: A New Algorithm for Graph Partitioning

Leo Grady and Eric L. Schwartz

SIAM J. Sci. Comput. 27, pp. 1844-1866 (23 pages) | Cited 3 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
We present a new algorithm for graph partitioning based on optimization of the combinatorial isoperimetric constant. It is shown empirically that this algorithm is competitive with other global partitioning algorithms in terms of partition quality. The isoperimetric algorithm is easy to parallelize, does not require coordinate information, and handles nonplanar graphs, weighted graphs, and families of graphs which are known to cause problems for other methods. Compared to spectral partitioning, the isoperimetric algorithm is faster and more stable. An exact circuit analogy to the algorithm is also developed with a natural random walks interpretation. The isoperimetric algorithm for graph partitioning is implemented in our publicly available Graph Analysis Toolbox [L. Grady, Ph.D. thesis, Boston University, Boston, MA, 2004], [L. Grady and E. L. Schwartz, Technical report TR-03-021, Boston University, Boston, MA, 2003] for MATLAB obtainable at http://eslab.bu.edu/software/graphanalysis/. This toolbox was used to generate all of the results compiled in the tables of this paper.

A Block Preconditioner for High-Order Mixed Finite Element Approximations to the Navier--Stokes Equations

David Kay and Emil Lungu

SIAM J. Sci. Comput. 27, pp. 1867-1880 (14 pages) | Cited 1 time

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
An iterative solver for the linear system arising from a discretization by high-order mixed finite elements for the linearized steady state Navier--Stokes problem is presented. Two mixed finite elements are considered: $\mathbf{P}_{k+1} P_{k}$ and $\mathbf{P}_{k+2} P_{k}$. The iterative solver is based on the GMRES method with right preconditioning where the preconditioner has a block upper triangular form. Computational results show that the convergence rate seems to be independent of the grid size parameter $h$ and the polynomial degree $k$ but is slightly affected by the increasing of the Reynolds number.

Efficient Minimization Methods of Mixed l2-l1 and l1-l1 Norms for Image Restoration

Haoying Fu, Michael K. Ng, Mila Nikolova, and Jesse L. Barlow

SIAM J. Sci. Comput. 27, pp. 1881-1902 (22 pages) | Cited 39 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
Image restoration problems are often solved by finding the minimizer of a suitable objective function. Usually this function consists of a data-fitting term and a regularization term. For the least squares solution, both the data-fitting and the regularization terms are in the $\ell$2 norm. In this paper, we consider the least absolute deviation (LAD) solution and the least mixed norm (LMN) solution. For the LAD solution, both the data-fitting and the regularization terms are in the $\ell$1 norm. For the LMN solution, the regularization term is in the $\ell$1 norm but the data-fitting term is in the $\ell$2 norm. Since images often have nonnegative intensity values, the proposed algorithms provide the option of taking into account the nonnegativity constraint. The LMN and LAD solutions are formulated as the solution to a linear or quadratic programming problem which is solved by interior point methods. At each iteration of the interior point method, a structured linear system must be solved. The preconditioned conjugate gradient method with factorized sparse inverse preconditioners is employed to solve such structured inner systems. Experimental results are used to demonstrate the effectiveness of our approach. We also show the quality of the restored images, using the minimization of mixed $\ell$2-$\ell$1 and $\ell1$-$\ell$1 norms, is better than that using only the $\ell$2 norm.

Fast Algorithms for Spherical Harmonic Expansions

Vladimir Rokhlin and Mark Tygert

SIAM J. Sci. Comput. 27, pp. 1903-1928 (26 pages) | Cited 12 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
An algorithm is introduced for the rapid evaluation at appropriately chosen nodes on the two-dimensional sphere $S^2$ in ${\mathbb R}^3$ of functions specified by their spherical harmonic expansions (known as the inverse spherical harmonic transform), and for the evaluation of the coefficients in spherical harmonic expansions of functions specified by their values at appropriately chosen points on $S^2$ (known as the forward spherical harmonic transform). The procedure is numerically stable and requires an amount of CPU time proportional to $N^2 (\log N) \log(1/\epsilon)$, where $N^2$ is the number of nodes in the discretization of $S^2$, and $\epsilon$ is the precision of computations. The performance of the algorithm is illustrated via several numerical examples.

Sensitivity Analysis of ODES/DAES Using the Taylor Series Method

Roberto Barrio

SIAM J. Sci. Comput. 27, pp. 1929-1947 (19 pages) | Cited 9 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
This paper studies the applicability of the Taylor method for the sensibility analysis of ODEs and DAEs. Extended automatic differentiation rules are introduced for the calculus of partial derivatives of Taylor series. The numerical method is implemented using an efficient variable-step variable-order scheme. Finally, some numerical tests are presented showing the benefits of the formulation.

Systematic Derivation of Jump Conditions for the Immersed Interface Method in Three-Dimensional Flow Simulation

Sheng Xu and Z. Jane Wang

SIAM J. Sci. Comput. 27, pp. 1948-1980 (33 pages) | Cited 20 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
In this paper, we systematically derive jump conditions for the immersed interface method [SIAM J. Numer. Anal., 31 (1994), pp. 1019-1044; SIAM J. Sci. Comput., 18 (1997), pp. 709-735] to simulate three-dimensional incompressible viscous flows subject to moving surfaces. The surfaces are represented as singular forces in the Navier--Stokes equations, which give rise to discontinuities of flow quantities. The principal jump conditions across a closed surface of the velocity, the pressure, and their normal derivatives have been derived by Lai and Li [Appl. Math. Lett., 14 (2001), pp. 149-154]. In this paper, we first extend their derivation to generalized surface parametrization. Starting from the principal jump conditions, we then derive the jump conditions of all first-, second-, and third-order spatial derivatives of the velocity and the pressure. We also derive the jump conditions of first- and second-order temporal derivatives of the velocity. Using these jump conditions, the immersed interface method is applicable to the simulation of three-dimensional incompressible viscous flows subject to moving surfaces, where near the surfaces the first- and second-order spatial derivatives of the velocity and the pressure can be discretized with, respectively, third- and second-order accuracy, and the first-order temporal derivatives of the velocity can be discretized with second-order accuracy.

Deterministic Simulation of the Boltzmann--Poisson System in GaAs-Based Semiconductors

María J. Cáceres, José A. Carrillo, and Armando Majorana

SIAM J. Sci. Comput. 27, pp. 1981-2009 (29 pages) | Cited 1 time

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
In this work we present a deterministic numerical scheme solving the one-dimensional Boltzmann--Poisson system for semiconductor materials with a complicated band structure such as GaAs. The model is based on a system of two kinetic equations for the electron population in each of the relevant valleys coupled through the Poisson equation to consider the self-consistent potential. Interaction of electrons with the crystal is taken into account by treating the most important inter- and intravalley scattering mechanisms in GaAs. The numerical scheme follows the main ideas in [J. A. Carrillo, I. M. Gamba, A. Majorana, and C.-W. Shu, J. Comput. Phys., 184 (2003), pp. 498-525] already developed for Si-based devices. The newest numerical features are the intervalley scattering, the singularity of the kernels of impurities, and optical polar phonon scattering. The charge conservation of the scheme is achieved by a suitable discretization of the collision operator within the accuracy of the middle-point integration rule. A detailed study of this numerical scheme is made by showing results on the bulk material, numerically validating the results by comparison to direct simulation Monte Carlo (DSMC) results, and testing our results for diodes and Gunn oscillations.

General Linear Methods for Volterra Integro-differential Equations with Memory

Chengjian Zhang and Stefan Vandewalle

SIAM J. Sci. Comput. 27, pp. 2010-2031 (22 pages) | Cited 7 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
A new class of numerical methods for Volterra integro-differential equations with memory is developed. The methods are based on the combination of general linear methods with compound quadrature rules. Sufficient conditions that guarantee global and asymptotic stability of the solution of the differential equation and its numerical approximation are established. Numerical examples illustrate the convergence and effectiveness of the numerical methods.

Contact Analysis of Cable Networks by Using Second-Order Cone Programming

Yoshihiro Kanno and Makoto Ohsaki

SIAM J. Sci. Comput. 27, pp. 2032-2052 (21 pages) | Cited 4 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
A method based on mathematical programming is proposed for large deformation and contact analysis of cable networks. By explicitly considering these nonsmooth behaviors, we formulate the linear complementarity problems over symmetric cones under some practically acceptable assumptions. We also present the equivalent second-order cone programming (SOCP) problems, which can be regarded as the minimization problem of total potential energy and complementary energy with the subsidiary constraints on the displacements and contact forces, respectively. By solving the presented SOCP problems by using the primal-dual interior-point method, the equilibrium configurations and internal forces of several cable networks are obtained without any assumptions on stress states and contact conditions.

On the Regularizing Power of Multigrid-type Algorithms

Marco Donatelli and Stefano Serra-Capizzano

SIAM J. Sci. Comput. 27, pp. 2053-2076 (24 pages) | Cited 8 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
We consider the deblurring problem of noisy and blurred images in the case of known space invariant point spread functions with four choices of boundary conditions. We combine an algebraic multigrid previously defined ad hoc for structured matrices related to space invariant operators (Toeplitz, circulants, trigonometric matrix algebras, etc.) and the classical geometric multigrid studied in the partial differential equations context. The resulting technique is parameterized in order to have more degrees of freedom: a simple choice of the parameters allows us to devise a quite powerful regularizing method. It defines an iterative regularizing method where the smoother itself has to be an iterative regularizing method (e.g., conjugate gradient, Landweber, conjugate gradient for normal equations, etc.). More precisely, with respect to the smoother, the regularization properties are improved and the total complexity is lower. Furthermore, in several cases, when it is directly applied to the system $A{\bf f}={\bf g}$, the quality of the restored image is comparable with that of all the best known techniques for the normal equations $A^TA{\bf f}=A^T{\bf g}$, but the related convergence is substantially faster. Finally, the associated curves of the relative errors versus the iteration numbers are "flatter" with respect to the smoother (the estimation of the stop iteration is less crucial). Therefore, we can choose multigrid procedures which are much more efficient than classical techniques without losing accuracy in the restored image (as often occurs when using preconditioning). Several numerical experiments show the effectiveness of our proposals.

AMGe---Coarsening Strategies and Application to the Oseen Equations

Markus Wabro

SIAM J. Sci. Comput. 27, pp. 2077-2097 (21 pages) | Cited 6 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
We provide some extensions to the algebraic multigrid method based on element interpolation (AMGe), concerning the agglomeration process, the application to nonconforming elements, and the application to the mixed finiteelement discretization of the Oseen linearized Navier--Stokes equations.
This last point, using AMGe for mixed finite elements, becomes straightforward because of the availability of coarse-level topologies. We show this exemplarily for the Crouzeix--Raviart element (including a stability result).
The numerical results show that it really pays off to take a closer look at the agglomeration strategy. A "wrong" choice can lead to insufficient convergence or even divergence of the overall multigrid method.

A Stability Notion for Lattice Boltzmann Equations

M. K. Banda, W. A. Yong, and A. Klar

SIAM J. Sci. Comput. 27, pp. 2098-2111 (14 pages) | Cited 4 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
Lattice Boltzmann equations are usually constructed to satisfy some physical requirements such as Galilean invariance and isotropy, to possess a velocity-independent pressure and no compressible effects, and so on. In this paper, a stability requirement is introduced as a new criterion for the constructions. With this requirement, we derive some relations of parameters for several lattice Boltzmann models. Interestingly, these relations are satisfied by many choices of parameters used in the literature.

A Reordering for the PageRank Problem

Amy N. Langville and Carl D. Meyer

SIAM J. Sci. Comput. 27, pp. 2112-2120 (9 pages) | Cited 7 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
We describe a reordering particularly suited to the PageRank problem, which reduces the computation of the PageRank vector to that of solving a much smaller system and then using forward substitution to get the full solution vector. We compare the theoretical rates of convergence of the original PageRank algorithm to that of the new reordered PageRank algorithm, showing that the new algorithm can do no worse than the original algorithm. We present results of an experimental comparison on five datasets, which demonstrate that the reordered PageRank algorithm can provide a speedup of as much as a factor of 6. We also note potential additional benefits that result from the proposed reordering.

Numerical Stochastic Integration for Quasi-Symplectic Flows

R. Mannella

SIAM J. Sci. Comput. 27, pp. 2121-2139 (19 pages) | Cited 3 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
Algorithms for the numerical integration of Langevin equations obeying detailed balance are introduced. The algorithms are derived fulfilling the requirements that they should become symplectic in the deterministic frictionless limit, and should reproduce the equilibrium distributions to some higher order in the integration time step when the equilibrium distribution exists. Extensions to the case when the system volume or pressure are kept constant are discussed. Comparisons with other integration schemes are carried out both for static and dynamical quantities.

Recycling Subspace Information for Diffuse Optical Tomography

Misha E. Kilmer and Eric de Sturler

SIAM J. Sci. Comput. 27, pp. 2140-2166 (27 pages) | Cited 14 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
We discuss the efficient solution of a long sequence of slowly varying linear systems arising in computations for diffuse optical tomographic imaging.
The reconstruction of three-dimensional absorption and scattering information by matching computed solutions from a parameterized model to measured data leads to a nonlinear least squares problem that we solve using the Gauss--Newton method with a line search. This algorithm requires the solution of a long sequence of linear systems. Each choice of parameters in the nonlinear least squares algorithm results in a different matrix describing the optical properties of the medium. These matrices change slowly from one step to the next, but may change significantly over many steps. For each matrix we must solve a set of linear systems with multiple shifts and multiple right-hand sides.
For this problem, we derive strategies for recycling Krylov subspace information that exploit properties of the application and the nonlinear optimization algorithm to significantly reduce the total number of iterations over all linear systems. Furthermore, we introduce variants of GCRO that exploit symmetry and that allow simultaneous solution of multiple shifted systems using a single Krylov subspace in combination with recycling. Although we focus on a particular application and optimization algorithm, our approach is applicable generally to problems where sequences of linear systems must be solved. This may guide other researchers to exploit the opportunities of tunable solvers.
We provide results for two sets of numerical experiments to demonstrate the effectiveness of the resulting method.

Bounding the Solutions of Parameter Dependent Nonlinear Ordinary Differential Equations

Adam B. Singer and Paul I. Barton

SIAM J. Sci. Comput. 27, pp. 2167-2182 (16 pages) | Cited 18 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
This paper presents two techniques for generating rigorous bounds on the solution of parameter dependent nonlinear ordinary differential equations. The first technique is an extension of differential inequalities that enables the construction of tight time varying state bounds by utilizing prior knowledge of the solution of the differential equations. The second technique provides a method for constructing pointwise in time convex and concave relaxations of the image of the solution of a system of parameter dependent differential equations on subsets of a Euclidean space. Two examples are presented to demonstrate the construction of the bounds. The examples include a brief discussion of a computer program written to automatically generate the bounds.
Close

close