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

2011

Volume 33, Issue 6, pp. 3087-3561


Bernstein–Bézier Finite Elements of Arbitrary Order and Optimal Assembly Procedures

Mark Ainsworth, Gaelle Andriamaro, and Oleg Davydov

SIAM J. Sci. Comput. 33, pp. 3087-3109 (23 pages)

Online Publication Date: November 01, 2011

Full Text: | Download PDF

Show Abstract
Algorithms are presented that enable the element matrices for the standard finite element space, consisting of continuous piecewise polynomials of degree $n$ on simplicial elements in $\mathbb{R}^d$, to be computed in optimal complexity $\mathcal{O}(n^{2d})$. The algorithms (i) take into account numerical quadrature; (ii) are applicable to nonlinear problems; and (iii) do not rely on precomputed arrays containing values of one-dimensional basis functions at quadrature points (although these can be used if desired). The elements are based on Bernstein polynomials and are the first to achieve optimal complexity for the standard finite element spaces on simplicial elements.

An $L^2$-Projection for the Galerkin-Characteristic Solution of Incompressible Flows

Mofdi El-Amrani and Mohammed Seaid

SIAM J. Sci. Comput. 33, pp. 3110-3131 (22 pages)

Online Publication Date: November 01, 2011

Full Text: | Download PDF

Show Abstract
This work is devoted to the development of an accurate Galerkin-characteristic method for simulation of incompressible viscous flows. The Galerkin-characteristic formulation is obtained by using a semi-Lagrangian discretization of the total derivative. The spatial discretization is performed using the finite element method on unstructured meshes. It can be interpreted as a fractional step technique where the convective part and the generalized Stokes part are treated separately. The semi-Lagrangian method requires high-order interpolating procedures for accuracy. A class of $L^2$-projection methods is developed and two interpolating schemes are considered by tracking the feet of the characteristic lines from the integration nodes. To solve the generalized Stokes problem we used a conjugate gradient algorithm. This method does not require any special correction for the pressure. We illustrate the performance of the proposed approach for an advection-diffusion equation with known analytical solution and for the benchmark problem of flow past a circular cylinder. We also present numerical results for a problem of incompressible buoyancy flow. The Galerkin-characteristic method using the $L^2$-projection has been found to be feasible and satisfactory.

Designing Optimal Spectral Filters for Inverse Problems

Julianne Chung, Matthias Chung, and Dianne P. O'Leary

SIAM J. Sci. Comput. 33, pp. 3132-3152 (21 pages)

Online Publication Date: November 01, 2011

Full Text: | Download PDF

Show Abstract
Spectral filtering suppresses the amplification of errors when computing solutions to ill-posed inverse problems; however, selecting good regularization parameters is often expensive. In many applications, data are available from calibration experiments. In this paper, we describe how to use such data to precompute optimal spectral filters. We formulate the problem in an empirical Bayes risk minimization framework and use efficient methods from stochastic and numerical optimization to compute optimal filters. Our formulation of the optimal filter problem is general enough to use a variety of assessments of goodness of the solution estimate, not just the mean square error. The relationship with the Wiener filter is discussed, and numerical examples from signal and image deconvolution illustrate that our proposed filters perform consistently better than well-established filtering methods. Furthermore, we show how our approach leads to easily computed uncertainty estimates for the pixel values.

Developing Finite Element Methods for Maxwell's Equations in a Cole–Cole Dispersive Medium

Jichun Li, Yunqing Huang, and Yanping Lin

SIAM J. Sci. Comput. 33, pp. 3153-3174 (22 pages)

Online Publication Date: November 01, 2011

Full Text: | Download PDF

Show Abstract
In this paper, we consider the time-dependent Maxwell's equations when Cole–Cole dispersive medium is involved. The Cole–Cole model contains a fractional time derivative term, which couples with the standard Maxwell's equations in free space and creates some challenges in developing and analyzing time-domain finite element methods for solving this model as mentioned in our earlier work [J. Li, J. Sci. Comput., 47 (2001), pp. 1–26]. By adopting some techniques developed for the fractional diffusion equations [V.J. Ervin, N. Heuer, and J.P. Roop, SIAM J. Numer. Anal., 45 (2007), pp. 572–591], [Y. Lin and C. Xu, J. Comput. Phys., 225 (2007), pp. 1533–1552], [F. Liu, P. Zhuang, V. Anh, I. Turner, and K. Burrage, Appl. Math. Comput., 191 (2007), pp. 12–20], we propose two fully discrete mixed finite element schemes for the Cole–Cole model. Numerical stability and optimal error estimates are proved for both schemes. The proposed algorithms are implemented and detailed numerical results are provided to justify our theoretical analysis.

Windowed Spectral Regularization of Inverse Problems

Julianne Chung, Glenn Easley, and Dianne P. O'Leary

SIAM J. Sci. Comput. 33, pp. 3175-3200 (26 pages)

Online Publication Date: November 01, 2011

Full Text: | Download PDF

Show Abstract
Regularization is used in order to obtain a reasonable estimate of the solution to an ill-posed inverse problem. One common form of regularization is to use a filter to reduce the influence of components corresponding to small singular values, perhaps using a Tikhonov least squares formulation. In this work, we break the problem into subproblems with narrower bands of singular values using spectrally defined windows, and we regularize each subproblem individually. We show how to use standard parameter-choice methods, such as the discrepancy principle and generalized cross-validation, in a windowed regularization framework. A perturbation analysis gives sensitivity estimates. We demonstrate the effectiveness of our algorithms on deblurring images and on the backward heat equation.

A Note on the Convergence of SOR for the PageRank Problem

Chen Greif and David Kurokawa

SIAM J. Sci. Comput. 33, pp. 3201-3209 (9 pages)

Online Publication Date: November 08, 2011

Full Text: | Download PDF

Show Abstract
A curious phenomenon when it comes to solving the linear system formulation of the PageRank problem is that while the convergence rate of Gauss–Seidel shows an improvement over Jacobi by a factor of approximately two, successive overrelaxation (SOR) does not seem to offer a meaningful improvement over Gauss–Seidel. This has been observed experimentally and noted in the literature, but to the best of our knowledge there has been no analytical explanation for this thus far. This convergence behavior is surprising because there are classes of matrices for which Gauss–Seidel is faster than Jacobi by a similar factor of two, and SOR accelerates convergence by an order of magnitude compared to Gauss–Seidel. In this short paper we prove analytically that the PageRank model has the unique feature that there exist PageRank linear systems for which SOR does not converge outside a very narrow interval depending on the damping factor, and that in such situations Gauss–Seidel may be the best choice among the relaxation parameters. Conversely, we show that within that narrow interval, there exists no PageRank problem for which SOR does not converge. Our result may give an analytical justification for the popularity of Gauss–Seidel as a solver for the linear system formulation of PageRank.

Alternating Evolution Schemes for Hyperbolic Conservation Laws

Haseena Saran and Hailiang Liu

SIAM J. Sci. Comput. 33, pp. 3210-3240 (31 pages)

Online Publication Date: November 10, 2011

Full Text: | Download PDF

Show Abstract
The alternating evolution (AE) system of Liu [J. Hyperbolic Differ. Equ., 5 (2008), pp. 421–447], $\partial_t u+ \partial_x f(v) =\frac{1}{\epsilon}(v-u)$, $\partial_t v + \partial_x f(u) =\frac{1}{\epsilon}(u-v)$, serves as a refined description of systems of hyperbolic conservation laws $\partial_t \phi+ \partial_x f(\phi)=0$, $\phi(x, 0)=\phi_0(x)$. The solution of conservation laws is precisely captured when two components take the same initial value as $\phi_0$, or is approached by two components exponentially fast when $\epsilon \downarrow 0$ if two initial states are sufficiently close. This nice property enables us to construct novel shock capturing schemes by sampling the AE system on alternating grids. In this paper we develop a class of local AE schemes by taking advantage of the AE system. Our approach is based on an average of the AE system over a hypercube centered at $x$ with vertices at $x\pm\Delta x$. The numerical scheme is then constructed by sampling the averaged system over alternating grids. Higher order accuracy is achieved by a combination of high order nonoscillatory polynomial reconstruction from the obtained averages and a matching Runge–Kutta solver in time discretization. Local AE schemes are made possible by letting the scale parameter $\epsilon$ reflect the local distribution of nonlinear waves. The AE schemes have the advantage of easier formulation and implementation, and efficient computation of the solution. For the first and second order local AE schemes applied to one-dimensional scalar conservation laws, we prove the numerical stability in the sense of satisfying the maximum principle and total variation diminishing property. The formulation procedure of AE schemes in multiple dimensions is given, followed by both the first and second order AE schemes for two-dimensional conservation laws. Numerical experiments for both scalar conservation laws and compressible Euler equations are presented to demonstrate the high order accuracy and capacity of these AE schemes.

A Fast Algorithm for Fourier Continuation

Mark Lyon

SIAM J. Sci. Comput. 33, pp. 3241-3260 (20 pages)

Online Publication Date: November 10, 2011

Full Text: | Download PDF

Show Abstract
A new algorithm is presented which provides a fast method for the computation of recently developed Fourier continuations (a particular type of Fourier extension method) that yield superalgebraically convergent Fourier series approximations of nonperiodic functions. Previously, the coefficients of an approximating Fourier series have been obtained by means of a regularized singular value decomposition (SVD)-based least-squares solution to an overdetermined linear system of equations. These SVD methods are effective when the size of the system does not become too large, but they quickly become unwieldy as the number of unknowns in the system grows. We demonstrate a novel decoupling of the least-squares problem which results in two systems of equations, one of which may be solved quickly by means of fast Fourier transforms (FFTs) and another that is demonstrated to be well approximated by a low-rank system. Utilizing randomized algorithms, the low-rank system is reduced to a significantly smaller system of equations. This new system is then efficiently solved with drastically reduced computational cost and memory requirements while still benefiting from the advantages of using a regularized SVD. The computational cost of the new algorithm in on the order of the cost of a single FFT multiplied by a slowly increasing factor that grows only logarithmically with the size of the system.

Fast Nonnegative Matrix Factorization: An Active-Set-Like Method and Comparisons

Jingu Kim and Haesun Park

SIAM J. Sci. Comput. 33, pp. 3261-3281 (21 pages)

Online Publication Date: November 10, 2011

Full Text: | Download PDF

Show Abstract
Nonnegative matrix factorization (NMF) is a dimension reduction method that has been widely used for numerous applications, including text mining, computer vision, pattern discovery, and bioinformatics. A mathematical formulation for NMF appears as a nonconvex optimization problem, and various types of algorithms have been devised to solve the problem. The alternating nonnegative least squares (ANLS) framework is a block coordinate descent approach for solving NMF, which was recently shown to be theoretically sound and empirically efficient. In this paper, we present a novel algorithm for NMF based on the ANLS framework. Our new algorithm builds upon the block principal pivoting method for the nonnegativity-constrained least squares problem that overcomes a limitation of the active set method. We introduce ideas that efficiently extend the block principal pivoting method within the context of NMF computation. Our algorithm inherits the convergence property of the ANLS framework and can easily be extended to other constrained NMF formulations. Extensive computational comparisons using data sets that are from real life applications as well as those artificially generated show that the proposed algorithm provides state-of-the-art performance in terms of computational speed.

A Boundary Integral Method for Computing the Dynamics of an Epitaxial Island

Shuwang Li and Xiaofan Li

SIAM J. Sci. Comput. 33, pp. 3282-3302 (21 pages)

Online Publication Date: November 29, 2011

Full Text: | Download PDF

Show Abstract
In this paper, we present a boundary integral method for computing the quasi-steady evolution of an epitaxial island. The problem consists of an adatom diffusion equation (with desorption) on terrace and a kinetic boundary condition at the step (island boundary). The normal velocity for step motion is determined by a two-sided flux. The integral formulation of the problem involves both double- and single-layer potentials due to the kinetic boundary condition. Numerical tests on a growing/shrinking circular or a slightly perturbed circular island are in excellent agreement with the linear analysis, demonstrating that the method is stable, efficient, and spectrally accurate in space. Nonlinear simulations for the growth of perturbed circular islands show that sharp tips and flat edges will form during growth instead of the usual dense branching morphology seen throughout physical and biological systems driven out of equilibrium. In particular, Bales–Zangwill instability is manifested in the form of wave-like fronts (meandering instability) around the tip regions. The numerical techniques presented here can be applied generally to a class of free/moving boundary problems in physical and biological science.

Collocation Methods for General Volterra Functional Integral Equations with Vanishing Delays

Hehu Xie, Ran Zhang, and Hermann Brunner

SIAM J. Sci. Comput. 33, pp. 3303-3332 (30 pages)

Online Publication Date: December 01, 2011

Full Text: | Download PDF

Show Abstract
We analyze the existence, uniqueness, and regularity properties of solutions for general Volterra functional integral equations with the delay function $\theta(t)$ vanishing at the initial point of the given interval $[0,T]$ (with $\theta(t)=qt, \; 0<q<1$, representing an important special case). The focus of the paper is then on the the attainable order of convergence, and the question of possible superconvergence, for collocation solutions in certain piecewise polynomial spaces. Numerical experiments complement the theoretical convergence results.

Fourth Order Time-Stepping for Kadomtsev–Petviashvili and Davey–Stewartson Equations

C. Klein and K. Roidot

SIAM J. Sci. Comput. 33, pp. 3333-3356 (24 pages)

Online Publication Date: December 01, 2011

Full Text: | Download PDF

Show Abstract
Purely dispersive partial differential equations such as the Korteweg–de Vries equation, the nonlinear Schrödinger equation, and higher dimensional generalizations thereof can have solutions which develop a zone of rapid modulated oscillations in the region where the corresponding dispersionless equations have shocks or blow-up. To numerically study such phenomena, fourth order time-stepping in combination with spectral methods is beneficial in resolving the steep gradients in the oscillatory region. We compare the performance of several fourth order methods for the Kadomtsev–Petviashvili and the Davey–Stewartson equations, two integrable equations in $2+1$ dimensions: these methods are exponential time-differencing, integrating factors, time-splitting, implicit Runge–Kutta, and Driscoll's composite Runge–Kutta method. The accuracy in the numerical conservation of integrals of motion is discussed.

Adaptive Pattern Research for Block FSAI Preconditioning

Carlo Janna and Massimiliano Ferronato

SIAM J. Sci. Comput. 33, pp. 3357-3380 (24 pages)

Online Publication Date: December 01, 2011

Full Text: | Download PDF

Show Abstract
An adaptive algorithm is presented to generate automatically the nonzero pattern of the block factored sparse approximate inverse (BFSAI) preconditioner. It is demonstrated that in symmetric positive definite (SPD) problems BFSAI minimizes an upper bound to the Kaporin number of the preconditioned matrix. The mathematical structure of this bound suggests an efficient and easily parallelizable strategy for improving the given nonzero pattern of BFSAI, thus providing a novel adaptive BFSAI (ABF) preconditioner. Numerical experiments performed on large sized finite element problems show that ABF coupled with a block incomplete Cholesky (IC) outperforms BFSAI-IC even by a factor of 4, preserving the same preconditioner density and exhibiting an excellent parallelization degree.

Learning to Predict Physical Properties using Sums of Separable Functions

Mayeul d'Avezac, Ryan Botts, Martin J. Mohlenkamp, and Alex Zunger

SIAM J. Sci. Comput. 33, pp. 3381-3401 (21 pages)

Online Publication Date: December 08, 2011

Full Text: | Download PDF

Show Abstract
We present an algorithm for learning the function that maps a material structure to its value on some property, given the value of this function on several structures. We pose this problem as one of learning (regressing) a function of many variables from scattered data. Each structure is first converted to a weighted set of points by a process that removes irrelevant translations and rotations but otherwise retains full information about the structure. Then, incorporating a weighted average for each structure, we construct the multivariate regression function as a sum of separable functions, following the paradigm of separated representations. The algorithm can treat all finite and periodic structures within a common framework, and in particular does not require all structures to lie on a common lattice. We show how the algorithm simplifies when the structures do lie on a common lattice, and we present numerical results for that case.

Simulation, Characterization, and Optimization of Metabolic Models with the High Performance Systems Biology Toolkit

Monte Lunacek, Ambarish Nag, David M. Alber, Kenny Gruchalla, Christopher H. Chang, and Peter A. Graf

SIAM J. Sci. Comput. 33, pp. 3402-3424 (23 pages)

Online Publication Date: December 13, 2011

Full Text: | Download PDF

Show Abstract
The High Performance Systems Biology Toolkit (HiPer SBTK) is a collection of simulation and optimization components for metabolic modeling and the means to assemble these components into large parallel processing hierarchies suiting a particular simulation and optimization need. The components come in a variety of different categories: model translation, model simulation, parameter sampling, sensitivity analysis, parameter estimation, and optimization. They can be configured at runtime into hierarchically parallel arrangements to perform nested combinations of simulation characterization tasks with excellent parallel scaling to thousands of processors. We describe the observations that led to the system, the components, and how one can arrange them. We show nearly 90% efficient scaling to over 13,000 processors, and we demonstrate three complex yet typical examples that have run on $\sim$1000 processors and accomplished billions of stiff ordinary differential equation simulations. This work opens the door for the systems biology metabolic modeling community to take effective advantage of large scale high performance computing resources for the first time.

A Bootstrap Algebraic Multilevel Method for Markov Chains

Matthias Bolten, Achi Brandt, James Brannick, Andreas Frommer, Karsten Kahl, and Ira Livshits

SIAM J. Sci. Comput. 33, pp. 3425-3446 (22 pages)

Online Publication Date: December 13, 2011

Full Text: | Download PDF

Show Abstract
This work concerns the development of an algebraic multilevel method for computing state vectors of Markov chains. We present an efficient bootstrap algebraic multigrid (AMG) method for this task. In our proposed approach, we employ a multilevel eigensolver, with interpolation built using ideas based on compatible relaxation, algebraic distances, and least squares fitting of test vectors. Our adaptive variational strategy for computation of the state vector of a given Markov chain is then a combination of this multilevel eigensolver and an associated additive multilevel preconditioned correction process. We show that the bootstrap AMG eigensolver by itself can efficiently compute accurate approximations to the steady state vector. An additional benefit of the bootstrap approach is that it yields an accurate interpolation operator for many other eigenmodes. This in turn allows for the use of the resulting multigrid hierarchy as a preconditioner to accelerate the generalized minimal residual (GMRES) iteration for computing an additive correction equation for the approximation to the steady state vector. Unlike other existing multilevel methods for Markov chains, our method does not employ any special processing of the coarse-level systems to ensure that stochastic properties of the fine-level system are maintained there. The proposed approach is applied to a range of test problems involving nonsymmetric M-matrices arising from stochastic matrices and showing promising results.

Detecting Localization in an Invariant Subspace

Christof Vömel and Beresford N. Parlett

SIAM J. Sci. Comput. 33, pp. 3447-3467 (21 pages)

Online Publication Date: December 13, 2011

Full Text: | Download PDF

Show Abstract
A normalized eigenvector, or, more interestingly, an invariant subspace, is localized if its significant entries are defined by just part(s) of the matrix and negligible elsewhere. This paper presents two new procedures to detect such localization in eigenvectors of a symmetric tridiagonal matrix. The procedures are intended for use before the actual eigenvector computation. If localization is found, one may reduce costs by computing the vectors just from the relevant matrix regions. Practical eigensolvers from numerical libraries such as LAPACK and ScaLAPACK already inspect a given tridiagonal $T$ for off-diagonal entries that are of small magnitude relative to the matrix norm. These so-called splitting points indicate that $T$ breaks into smaller blocks, each one defining a subset of eigenvalues and localized eigenvectors. However, localization can occur even when none of the off-diagonals is particularly small. Our study investigates this more complicated phenomenon in the context of invariant subspaces belonging to isolated eigenvalue clusters.

Algebraic Distance on Graphs

Jie Chen and Ilya Safro

SIAM J. Sci. Comput. 33, pp. 3468-3490 (23 pages)

Online Publication Date: December 13, 2011

Full Text: | Download PDF

Show Abstract
Measuring the connection strength between a pair of vertices in a graph is one of the most important concerns in many graph applications. Simple measures such as edge weights may not be sufficient for capturing the effects associated with short paths of lengths greater than one. In this paper, we consider an iterative process that smooths an associated value for nearby vertices, and we present a measure of the local connection strength (called the algebraic distance; see [D. Ron, I. Safro, and A. Brandt, Multiscale Model. Simul., 9 (2011), pp. 407–423]) based on this process. The proposed measure is attractive in that the process is simple, linear, and easily parallelized. An analysis of the convergence property of the process reveals that the local neighborhoods play an important role in determining the connectivity between vertices. We demonstrate the practical effectiveness of the proposed measure through several combinatorial optimization problems on graphs and hypergraphs.

Multitissue Tetrahedral Image-to-mesh Conversion with Guaranteed Quality and Fidelity

Andrey N. Chernikov and Nikos P. Chrisochoides

SIAM J. Sci. Comput. 33, pp. 3491-3508 (18 pages)

Online Publication Date: December 15, 2011

Full Text: | Download PDF

Show Abstract
We present a novel algorithm for tetrahedral image-to-mesh conversion which allows for guaranteed bounds on the smallest dihedral angle and on the distance between the boundaries of the mesh and the boundaries of the tissues. The algorithm produces a small number of mesh elements that comply with these bounds. We also describe and evaluate our implementation of the proposed algorithm that is compatible in performance with a state-of-the art Delaunay code, but in addition solves the small dihedral angle problem.

A Particle-in-cell Method with Adaptive Phase-space Remapping for Kinetic Plasmas

B. Wang, G. H. Miller, and P. Colella

SIAM J. Sci. Comput. 33, pp. 3509-3537 (29 pages)

Online Publication Date: December 15, 2011

Full Text: | Download PDF

Show Abstract
We present a new accurate and efficient particle-in-cell (PIC) method for computing the dynamics of one-dimensional kinetic plasmas. The method overcomes the numerical noise inherent in particle-based methods by periodically remapping the distribution function on a hierarchy of locally refined grids in phase space. Remapping on phase-space grids also provides an opportunity to integrate a collisional model and an associated grid-based solver. The positivity of the distribution function is enforced by redistributing excess phase-space density in a local neighborhood. We demonstrate the method on a number of standard plasma physics problems. It is shown that remapping significantly reduces the numerical noise and results in a more consistent second-order method than the standard PIC method. An error analysis is presented which is based on prior results of Cottet and Raviart's work [SIAM J. Numer. Anal., 21 (1984), pp. 52–76].

Local POD Plus Galerkin Projection in the Unsteady Lid-Driven Cavity Problem

Filippo Terragni, Eusebio Valero, and José M. Vega

SIAM J. Sci. Comput. 33, pp. 3538-3561 (24 pages) | Cited 1 time

Online Publication Date: December 22, 2011

Full Text: | Download PDF

Show Abstract
A local proper orthogonal decomposition (POD) plus Galerkin projection method is applied to the unsteady lid-driven cavity problem, namely the incompressible fluid flow in a two-dimensional box whose upper wall is moved back and forth at moderately large values of the Reynolds number. Such a method was recently developed for one-dimensional parabolic problems. Its extension to fluid dynamics problems is nontrivial (especially if rough CFD codes are used) and consists of using a computational fluid dynamics (CFD) code and a Galerkin system (GS) in a sequence of interspersed intervals $I_{CFD}$ and $I_{GS}$, respectively. The POD manifold is calculated retaining the most energetic POD modes resulting from the snapshots computed in the $I_{CFD}$ intervals; in fact, the POD manifold is completely calculated in the first $I_{CFD}$ interval but only updated in subsequent $I_{CFD}$ intervals. Intended to mimic industrial solvers, the CFD code contains unphysical terms that are introduced for purely numerical reasons to accelerate runs. The use of such a CFD code poses the question of which equations should be Galerkin projected onto the POD manifold, the exact governing equations or the CFD numerical scheme. Also, Galerkin projection can be made using either the standard $L_2$ inner product or a nonstandard one, based on a limited number of mesh points. After addressing these issues, a method is constructed that is able to accelerate the CFD code by a factor of the order of 5–15, depending on the Reynolds number and the nature (steady, periodic, or quasi-periodic) of the forcing velocity.
Close

close