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 Next Issue

2009

Volume 31, Issue 3, pp. 1603-2417


Two-Layer Shallow Water System: A Relaxation Approach

Rémi Abgrall and Smadar Karni

SIAM J. Sci. Comput. 31, pp. 1603-1627 (25 pages)

Online Publication Date: February 27, 2009

Full Text: | Download PDF

Show Abstract
The two-layer shallow water system is an averaged flow model. It forms a nonconservative system which is only conditionally hyperbolic. The coupling between the layers, due to the hydrostatic pressure assumption, does not provide explicit access to the system eigenstructure, which is inconvenient for Riemann solution based numerical schemes. We consider a relaxation approach which offers greater decoupling and accessible eigenstructure. The stability of the model is discussed. Numerical results are shown for unsteady flows as well as for smooth and nonsmooth steady flows.

Spectral Projected Gradient Method with Inexact Restoration for Minimization with Nonconvex Constraints

M. A. Gomes-Ruggiero, J. M. Martínez, and S. A. Santos

SIAM J. Sci. Comput. 31, pp. 1628-1652 (25 pages) | Cited 2 times

Online Publication Date: February 27, 2009

Full Text: | Download PDF

Show Abstract
This work takes advantage of the spectral projected gradient direction within the inexact restoration framework to address nonlinear optimization problems with nonconvex constraints. The proposed strategy includes a convenient handling of the constraints, together with nonmonotonic features to speed up convergence. The numerical performance is assessed by experiments with hard-spheres problems, pointing out that the inexact restoration framework provides an adequate environment for the extension of the spectral projected gradient method for general nonlinearly constrained optimization.

Instabilities and Inaccuracies in the Integration of Highly Oscillatory Problems

M. P. Calvo and J. M. Sanz-Serna

SIAM J. Sci. Comput. 31, pp. 1653-1677 (25 pages) | Cited 2 times

Online Publication Date: February 27, 2009

Full Text: | Download PDF

Show Abstract
By means of a two-frequency test problem, we analyze the instabilities and inaccuracies that may impair the performance of multiple time steps/split operator integrators in highly oscillatory situations, such as those encountered in molecular dynamics, astrophysics, or partial differential equations describing waves. Considered are the impulse (Verlet-I/r-RESPA) method, the mollified impulse method, and the reversible averaging integrator. The analysis covers errors in positions, momenta, and energy.

Sparse Fourier Transform via Butterfly Algorithm

Lexing Ying

SIAM J. Sci. Comput. 31, pp. 1678-1694 (17 pages) | Cited 4 times

Online Publication Date: February 27, 2009

Full Text: | Download PDF

Show Abstract
This paper introduces a fast algorithm for computing sparse Fourier transforms with spatial and Fourier data supported on curves or surfaces. This problem appears naturally in several important applications of wave scattering, digital signal processing, and reflection seismology. The main idea of the algorithm is that the interaction between a frequency region and a spatial region is approximately low rank if the product of their widths are bounded by the maximum frequency. Based on this property, we can approximate the interaction between these two boxes accurately and compactly using a small number of equivalent sources. The overall structure of the algorithm follows the butterfly algorithm. The computation is further accelerated by exploiting the tensor-product property of the Fourier kernel in two and three dimensions. The proposed algorithm is accurate and has the optimal complexity. We present numerical results in both two and three dimensions.

Numerical Conformal Mapping via a Boundary Integral Equation with the Generalized Neumann Kernel

Mohamed M. S. Nasser

SIAM J. Sci. Comput. 31, pp. 1695-1715 (21 pages)

Online Publication Date: February 27, 2009

Full Text: | Download PDF

Show Abstract
We present a unified boundary integral method for approximating the conformal mappings from any bounded or unbounded multiply connected region $G$ onto the five classical canonical slit domains. The method is based on a uniquely solvable boundary integral equation with the generalized Neumann kernel. Using the proposed method, the approximate mapping functions onto the five canonical slit domains can be computed in a unified way by solving linear systems with a common coefficient matrix. The method can be also used for calculating the conformal mappings of simply and doubly connected regions. The performance of the method is illustrated by several examples for regions with smooth boundaries and with piecewise smooth boundaries.

Computing the Conical Function $P^{\mu}_{-1/2+i\tau}(x)$

Amparo Gil, Javier Segura, and Nico M. Temme

SIAM J. Sci. Comput. 31, pp. 1716-1741 (26 pages)

Online Publication Date: February 27, 2009

Full Text: | Download PDF

Show Abstract
A stable computational scheme for the conical function $P^{\mu}_{-1/2+i\tau}(x)$ for $x>-1$, real $\tau$, and $\mu\le 0$ or $\mu\in\mathbb{N}$ is presented. The scheme combines uniform asymptotic expansions for large $|\mu|$ with the application of the three-term recurrence relation on the $\mu$ index in the direction of decreasing $|\mu|$ when $x>0$. When $x<0$, the conditioning of recursion is the opposite, and conical functions can be computed in the direction of increasing $|\mu|$.

Central-Upwind Schemes for Two-Layer Shallow Water Equations

Alexander Kurganov and Guergana Petrova

SIAM J. Sci. Comput. 31, pp. 1742-1773 (32 pages)

Online Publication Date: February 27, 2009

Full Text: | Download PDF

Show Abstract
We derive a second-order semidiscrete central-upwind scheme for one- and two-dimensional systems of two-layer shallow water equations. We prove that the presented scheme is well-balanced in the sense that stationary steady-state solutions are exactly preserved by the scheme and positivity preserving; that is, the depth of each fluid layer is guaranteed to be nonnegative. We also propose a new technique for the treatment of the nonconservative products describing the momentum exchange between the layers. The performance of the proposed method is illustrated on a number of numerical examples, in which we successfully capture (quasi) steady-state solutions and propagating interfaces.

Simulation of Stochastic Reaction-Diffusion Processes on Unstructured Meshes

Stefan Engblom, Lars Ferm, Andreas Hellander, and Per Lötstedt

SIAM J. Sci. Comput. 31, pp. 1774-1797 (24 pages) | Cited 9 times

Online Publication Date: February 27, 2009

Full Text: | Download PDF

Show Abstract
We model stochastic chemical systems with diffusion by a reaction-diffusion master equation. On a macroscopic level, the governing equation is a reaction-diffusion equation for the averages of the chemical species. On a mesoscopic level, the master equation for a well stirred chemical system is combined with a discretized Brownian motion in space to obtain the reaction-diffusion master equation. The space is covered in our method by an unstructured mesh, and the diffusion coefficients on the mesoscale are obtained from a finite element discretization of the Laplace operator on the macroscale. The resulting method is a flexible hybrid algorithm in that the diffusion can be handled either on the meso- or on the macroscale level. The accuracy and the efficiency of the method are illustrated in three numerical examples inspired by molecular biology.

An Immersed Interface Method for the Incompressible Navier–Stokes Equations with Discontinuous Viscosity Across the Interface

Zhijun Tan, D.V. Le, K.M. Lim, and B.C. Khoo

SIAM J. Sci. Comput. 31, pp. 1798-1819 (22 pages)

Online Publication Date: February 27, 2009

Full Text: | Download PDF

Show Abstract
We present an immersed interface method for solving the incompressible Navier–Stokes equations with discontinuous viscosity across the interface and singular forces. The method is based on the augmented strategy proposed by Li, Ito, and Lai [Comput. Fluids, 36 (2007), pp. 622–635] to decouple the jump conditions of the fluid variables through the introduction of two augmented variables. In the proposed method, the immersed interface is represented by a number of Lagrangian control points, and the augmented interface variables applicable along the interface are solved by the LU method or the GMRES iterative method. In addition, we also derive a new second order accurate bilinear interpolating scheme for the discontinuous velocity and a new jump condition for the normal derivative of the pressure. The jumps in both pressure and velocity and the jumps in their derivatives are related to the augmented interface variables and/or the forces which are interpolated using cubic splines. For a flexible interface, the forces that the interface exerts on the fluid are computed from the constitutive relation of the flexible interface and are applied to the fluid through the jump conditions. The position of the flexible interface is updated implicitly within each time step. The Navier–Stokes equations are discretized on a staggered Cartesian grid by a second order accurate projection method. The numerical results show that the overall scheme is second order accurate.

Construction of Data-Sparse $\mathcal{H}^2$-Matrices by Hierarchical Compression

Steffen Börm

SIAM J. Sci. Comput. 31, pp. 1820-1839 (20 pages) | Cited 2 times

Online Publication Date: February 27, 2009

Full Text: | Download PDF

Show Abstract
Hierarchical matrices ($\mathcal{H}$-matrices) provide an elegant approach to handling large densely populated matrices: The matrix is split into a hierarchy of blocks, and each block is approximated by a low-rank matrix in factorized form. It has been demonstrated that this representation can be used to treat integral and partial differential equations, solve matrix equations from the field of control theory, and evaluate matrix functions efficiently. $\mathcal{H}^2$-matrices use a refined representation that employs a multilevel structure in order to reduce the storage requirements of hierarchical matrices. It has been shown that $\mathcal{H}^2$-matrices can significantly reduce storage requirements for large problems, in particular when combined with modern error control schemes. Until now, all algorithms for constructing an efficient approximation of a general matrix by an $\mathcal{H}^2$-matrix required a representation of the entire original matrix to be kept in storage; therefore the storage requirements of $\mathcal{H}^2$-matrix algorithms could be far larger than those of the final approximation. This paper presents a new approach that allows us to construct an $\mathcal{H}^2$-matrix without storing the entire original matrix. The central idea is to approximate submatrices and combine them by an efficient new algorithm to form approximations of larger matrices until the entire matrix has been treated. Using this new approach, many $\mathcal{H}$-matrix algorithms can be “refitted” easily to compute results in the more efficient representation. Possible applications include efficient matrix arithmetics for the construction of preconditioners, the approximation of matrix functions or solutions of matrix equations, or efficient compression schemes based on the popular cross approximation algorithms.

Multivariate Regression and Machine Learning with Sums of Separable Functions

Gregory Beylkin, Jochen Garcke, and Martin J. Mohlenkamp

SIAM J. Sci. Comput. 31, pp. 1840-1857 (18 pages) | Cited 1 time

Online Publication Date: March 13, 2009

Full Text: | Download PDF

Show Abstract
We present an algorithm for learning (or estimating) a function of many variables from scattered data. The function is approximated by a sum of separable functions, following the paradigm of separated representations. The central fitting algorithm is linear in both the number of data points and the number of variables and, thus, is suitable for large data sets in high dimensions. We present numerical evidence for the utility of these representations. In particular, we show that our method outperforms other methods on several benchmark data sets.

Iterative Solution of Piecewise Linear Systems and Applications to Flows in Porous Media

Luigi Brugnano and Vincenzo Casulli

SIAM J. Sci. Comput. 31, pp. 1858-1873 (16 pages) | Cited 1 time

Online Publication Date: March 13, 2009

Full Text: | Download PDF

Show Abstract
The correct numerical modeling of free-surface hydrodynamic problems often requires to have the solution of special linear systems whose coefficient matrix is a piecewise constant function of the solution itself. In doing so, one may fulfill relevant physical constraints. The existence, the uniqueness, and two constructive iterative methods to solve a piecewise linear system of the form $\max[\boldsymbol{l},\min(\boldsymbol{u},\mathbf{x})]+T\mathbf{x}=\mathbf{b}$ are analyzed. The methods are shown to have a finite termination property; i.e., they converge to an exact solution in a finite number of steps and, actually, they converge very quickly, as confirmed by a few numerical tests, which are derived from the mathematical modeling of flows in porous media.

An Asymptotic Fitting Finite Element Method with Exponential Mesh Refinement for Accurate Computation of Corner Eddies in Viscous Flows

Alexander V. Shapeev and Ping Lin

SIAM J. Sci. Comput. 31, pp. 1874-1900 (27 pages)

Online Publication Date: March 13, 2009

Full Text: | Download PDF

Show Abstract
It is well known that any viscous fluid flow near a corner consists of infinite series of eddies with decreasing size and intensity, unless the angle is larger than a certain critical angle [H. K. Moffat, J. Fluid Mech., 18 (1964), pp. 1–18]. The objective of the current work is to simulate such infinite series of eddies occurring in steady flows in domains with corners. The problem is approached by high-order finite element method with exponential mesh refinement near the corners, coupled with analytical asymptotics of the flow near the corners. Such approach allows one to compute position and intensity of the eddies near the corners in addition to the other main features of the flow. The method was tested on the problem of the lid-driven cavity flow as well as on the problem of the backward-facing step flow. The results of computations of the lid-driven cavity problem show that the proposed method computes the central eddy with accuracy comparable to the best of existing methods and is more accurate for computing the corner eddies than the existing methods. The results also indicate that the relative error of finding the eddies' intensity and position decreases uniformly for all the eddies as the mesh is refined (i.e., the relative error in computation of different eddies does not depend on their size).

Simulations of Valveless Pumping in an Open Elastic Tube

Wanho Lee, Eunok Jung, and Sunmi Lee

SIAM J. Sci. Comput. 31, pp. 1901-1925 (25 pages)

Online Publication Date: March 27, 2009

Full Text: | Download PDF

Show Abstract
Mathematical models and numerical simulations of flows driven by pumping without valves (valveless pumping) are presented. This work has been originally motivated by a biomedical objective to explain the complicated valveless blood flow mechanism in the circulation. For instance, the fetal circulation before the formation of valves and blood circulation during cardiopulmonary resuscitation. A mathematical model of valveless pumping either in a closed loop system or in an open system consists of a couple of tubes with different elasticities or radii. In this work, we develop new models which are not necessarily connected by two different materials: only one open elastic tube and an open elastic tube contained in a closed (almost) rigid tank. Although only one soft material is used or the soft and rigid materials are not connected, we have observed the existence of a net flow driven by the periodic compress-and-release action and the important features of valveless pumping that have been observed in earlier models or experiments. We have also shown that there exist mostly one-directional net flows inside the elastic tube if a long and thin tube is considered and the pumping is applied to the short portion near the edge of the tube. This might explain why most earlier reports on the physical experiments could not easily observe the existence of both directional net flows depending on the parameters, such as frequency. Another important result is that the direction and the magnitude of a net flow can be explained by the sign and the amount of power, which is work done on the fluid by the fluid pressure and the elastic wall over one period, respectively. A new feature in this work is that only one elastic material is sufficient to show the existence of a net flow in a valveless pump system. Because of a simple structure in our new model, it is much easier to construct a valveless pump system applicable to real world applications, such as a microelectromechanical system device.

On a Class of Uniformly Accurate IMEX Runge–Kutta Schemes and Applications to Hyperbolic Systems with Relaxation

Sebastiano Boscarino and Giovanni Russo

SIAM J. Sci. Comput. 31, pp. 1926-1945 (20 pages)

Online Publication Date: March 27, 2009

Full Text: | Download PDF

Show Abstract
In this paper we consider hyperbolic systems with relaxation in which the relaxation time $\varepsilon$ may vary from values of order one to very small values. When $\varepsilon$ is very small, the relaxation term becomes very strong and highly stiff, and underresolved numerical schemes may produce spurious results. In such cases it is important to have schemes that work uniformly with respect to $\varepsilon$. IMplicit-EXplicit (IMEX) Runge–Kutta (R-K) schemes have been widely used for the time evolution of hyperbolic partial differential equations but the schemes existing in literature do not exhibit uniform accuracy with respect to the relaxation time. We develop new IMEX R-K schemes for hyperbolic systems with relaxation that present better uniform accuracy than the ones existing in the literature and in particular produce good behavior with high order accuracy in the asymptotic limit, i.e., when $\varepsilon$ is very small. These schemes are obtained by imposing new additional order conditions to guarantee better accuracy over a wide range of the relaxation time. We propose the construction of new third-order IMEX R-K schemes of type CK [S. Boscarino, SIAM J. Numer. Anal., 45 (2008), pp. 1600–1621]. In several test problems, these schemes, with a fixed spatial discretization, exhibit for all range of the relaxation time an almost uniform third-order accuracy.

MultiStage Approaches for Optimal Offline Checkpointing

Philipp Stumm and Andrea Walther

SIAM J. Sci. Comput. 31, pp. 1946-1967 (22 pages) | Cited 2 times

Online Publication Date: March 27, 2009

Full Text: | Download PDF

Show Abstract
The computation of derivatives for optimizing time-dependent flow problems is often based on the integration of the adjoint differential equation. For this purpose, the knowledge of the complete forward solution is required. In the area of flow control, especially for three-dimensional problems, it may be impossible to keep track of the full forward solution due to the lack of storage capacities. One usual method to overcome this problem is checkpointing. Using a checkpointing strategy, only parts of the forward solution are kept in the main memory and additional time steps have to be performed. If one extends this approach such that some checkpoints can be stored on disc too, one can reduce the number of additional time steps. On the other hand, one has to take the access cost to one checkpoint into account. On parallel machines, one may use parallel I/O facilities to store and retrieve checkpoints. In these cases, the access cost of the checkpoints is also no longer negligible. Therefore, in this paper, the write and read counts for each checkpoint in a binomial checkpointing approach are examined. This way, the overall access cost to checkpoints is minimized. Numerical results illustrate the derived checkpointing approaches. They also show that checkpointing techniques may reduce the overall computing time despite the required recalculations.

Accelerating the Multilevel Fast Multipole Algorithm with the Sparse-Approximate-Inverse (SAI) Preconditioning

Tahi̇r Malas and Levent Gürel

SIAM J. Sci. Comput. 31, pp. 1968-1984 (17 pages) | Cited 2 times

Online Publication Date: March 27, 2009

Full Text: | Download PDF

Show Abstract
With the help of the multilevel fast multipole algorithm, integral-equation methods can be used to solve real-life electromagnetics problems both accurately and efficiently. Increasing problem dimensions, on the other hand, necessitate effective parallel preconditioners with low setup costs. In this paper, we consider sparse approximate inverses generated from the sparse near-field part of the dense coefficient matrix. In particular, we analyze pattern selection strategies that can make efficient use of the block structure of the near-field matrix, and we propose a load-balancing method to obtain high scalability during the setup. We also present some implementation details, which reduce the computational cost of the setup phase. In conclusion, for the open-surface problems that are modeled by the electric-field integral equation, we have been able to solve ill-conditioned linear systems involving millions of unknowns with moderate computational requirements. For closed-surface problems that can be modeled by the combined-field integral equation, we reduce the solution times significantly compared to the commonly used block-diagonal preconditioner.

Energy Conserving Explicit Local Time Stepping for Second-Order Wave Equations

Julien Diaz and Marcus J. Grote

SIAM J. Sci. Comput. 31, pp. 1985-2014 (30 pages) | Cited 1 time

Online Publication Date: March 27, 2009

Full Text: | Download PDF

Show Abstract
Locally refined meshes impose severe stability constraints on explicit time-stepping methods for the numerical simulation of time dependent wave phenomena. To overcome that stability restriction, local time-stepping methods are developed, which allow arbitrarily small time steps precisely where small elements in the mesh are located. When combined with a symmetric finite element discretization in space with an essentially diagonal mass matrix, the resulting discrete numerical scheme is explicit, is inherently parallel, and exactly conserves a discrete energy. Starting from the standard second-order “leap-frog” scheme, time-stepping methods of arbitrary order of accuracy are derived. Numerical experiments illustrate the efficiency and usefulness of these methods and validate the theory.

Numerical Method for 2D Simulation of a Silicon MESFET with a Hydrodynamical Model Based on the Maximum Entropy Principle

A. M. Blokhin and A. S. Ibragimova

SIAM J. Sci. Comput. 31, pp. 2015-2046 (32 pages)

Online Publication Date: April 03, 2009

Full Text: | Download PDF

Show Abstract
We discuss the application of the method of lines for constructing an original numerical algorithm for two-dimensional (2D) simulation of a silicon metal semiconductor field effect transistor (MESFET) with the new hydrodynamical model based on the maximum entropy principle [A. M. Anile and V. Romano, Contin. Mech. Thermodyn., 11 (1999), pp. 307–325, V. Romano, Contin. Mech. Thermodyn., 12 (2000), pp. 31–51]. The substantiation of the proposed algorithm is conducted on the example of the model problem for the Poisson equation. The results of numerical experiments are given for different parameters, characterizing silicon MESFET.

Efficient Schemes for Total Variation Minimization Under Constraints in Image Processing

Pierre Weiss, Laure Blanc-Féraud, and Gilles Aubert

SIAM J. Sci. Comput. 31, pp. 2047-2080 (34 pages) | Cited 10 times

Online Publication Date: April 03, 2009

Full Text: | Download PDF

Show Abstract
This paper presents new fast algorithms to minimize total variation and more generally $l^1$-norms under a general convex constraint. Such problems are standards of image processing. The algorithms are based on a recent advance in convex optimization proposed by Yurii Nesterov. Depending on the regularity of the data fidelity term, we solve either a primal problem or a dual problem. First we show that standard first-order schemes allow one to get solutions of precision $\epsilon$ in $O(\frac{1}{\epsilon^2})$ iterations at worst. We propose a scheme that allows one to obtain a solution of precision $\epsilon$ in $O(\frac{1}{\epsilon})$ iterations for a general convex constraint. For a strongly convex constraint, we solve a dual problem with a scheme that requires $O(\frac{1}{\sqrt{\epsilon}})$ iterations to get a solution of precision $\epsilon$. Finally we perform some numerical experiments which confirm the theoretical results on various problems of image processing.

Fourier Analysis for Multigrid Methods on Triangular Grids

F.J. Gaspar, J.L. Gracia, and F.J. Lisbona

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

Online Publication Date: April 03, 2009

Full Text: | Download PDF

Show Abstract
In this paper a local Fourier analysis technique for multigrid methods on triangular grids is presented. The analysis is based on an expression of the Fourier transform in new coordinate systems, both in space variables and in frequency variables, associated with reciprocal bases. This tool makes it possible to study different components of the multigrid method in a very similar way to that of rectangular grids. Different smoothers for the discrete Laplace operator obtained with linear finite elements are analyzed. A new three-color smoother has been studied and has proven to be the best choice for “near” equilateral triangles. It is also shown that the block-line smoothers are more appropriate for irregular triangles. Numerical test calculations validate the theoretical predictions.

Quality Triangulations with Locally Optimal Steiner Points

Hale Erten and Alper Üngör

SIAM J. Sci. Comput. 31, pp. 2103-2130 (28 pages) | Cited 1 time

Online Publication Date: April 29, 2009

Full Text: | Download PDF

Show Abstract
We propose two novel ideas to improve the performance of Delaunay refinement algorithms which are used for computing quality triangulations. The first idea is an effective use of the Voronoi diagram and unifies previously suggested Steiner point insertion schemes (circumcenter, sink, off-center) together with a new strategy. The second idea is the integration of a new local smoothing strategy into the refinement process. These lead to two new versions of Delaunay refinement, where the second is simply an extension of the first. For a given input domain and a constraint angle $\alpha$, Delaunay refinement algorithms aim to compute triangulations that have all angles at least $\alpha$. The original Delaunay refinement algorithm of Ruppert is proven to terminate with size-optimal quality triangulations for $\alpha\le20.7^\circ$. In practice, the original and the consequent Delaunay refinement algorithms generally work for $\alpha\le34^\circ$ and fail to terminate for larger constraint angles. Our algorithms provide the same theoretical guarantees as the previous Delaunay refinement algorithms. The second of the proposed algorithms generally terminates for constraint angles up to $42^\circ$. Experiments also indicate that our algorithm computes significantly (about by a factor of two) smaller triangulations than the output of the previous Delaunay refinement algorithms. Moreover, the new algorithms are experimentally shown to outperform the previous algorithms even in the existence of additional constraints, such as the maximum area triangle constraint which is commonly used for computing uniform triangulations.

Design and Implementation of Predictors for Additive Semi-Implicit Runge–Kutta Methods

Inmaculada Higueras, José Miguel Mantas, and Teo Roldán

SIAM J. Sci. Comput. 31, pp. 2131-2150 (20 pages)

Online Publication Date: April 29, 2009

Full Text: | Download PDF

Show Abstract
Space discretization of some time-dependent partial differential equations gives rise to stiff systems of ordinary differential equations. In this case, implicit methods should be used, and therefore, in general, nonlinear systems must be solved. The solutions to these systems are approximated by iterative schemes, and, in order to obtain an efficient code, good initializers should be used. Recently, a parallel code based on some Runge–Kutta and additive Runge–Kutta methods has been constructed, focusing especially on additive semi-implicit Runge–Kutta (ASIRK) methods. The aim of the present paper is to develop efficient initializers for these methods.

Computing with Expansions in Gegenbauer Polynomials

Jens Keiner

SIAM J. Sci. Comput. 31, pp. 2151-2171 (21 pages) | Cited 1 time

Online Publication Date: April 29, 2009

Full Text: | Download PDF

Show Abstract
We develop fast algorithms for computations involving finite expansions in Gegenbauer polynomials. A method is described to convert any finite expansion between different families of Gegenbauer polynomials. For a degree-$n$ expansion the computational cost is $\mathcal{O}(n(\log(1/\varepsilon)+|\alpha-\beta|))$, where $\varepsilon$ is the prescribed accuracy, and $\alpha$ and $\beta$ are the respective Gegenbauer indices. Special cases involving Chebyshev polynomials of first kind are particularly important. In combination with (nonequispaced) discrete cosine transforms, we obtain efficient methods for the evaluation of expansions at prescribed nodes, and for the projection onto a sequence of Gegenbauer polynomials from given function values, respectively.

Multigrid Methods and Sparse-Grid Collocation Techniques for Parabolic Optimal Control Problems with Random Coefficients

A. Borzì and G. von Winckel

SIAM J. Sci. Comput. 31, pp. 2172-2192 (21 pages) | Cited 1 time

Online Publication Date: April 29, 2009

Full Text: | Download PDF

Show Abstract
An efficient computational framework to solve nonlinear parabolic optimal control problems with random coefficients is presented. This framework allows us to investigate the influence of randomness or uncertainty of problem's parameters values on the control provided by the optimal control theory. The proposed framework combines space-time multigrid methods with sparse-grid collocation techniques. Theoretical and numerical results of computation of stochastic optimal control solutions and formulation of mean control functions are presented.

Optimized Schwarz Methods for Maxwell's Equations

V. Dolean, M.J. Gander, and L. Gerardo-Giorda

SIAM J. Sci. Comput. 31, pp. 2193-2213 (21 pages) | Cited 4 times

Online Publication Date: May 07, 2009

Full Text: | Download PDF

Show Abstract
Over the last two decades, classical Schwarz methods have been extended to systems of hyperbolic partial differential equations, using characteristic transmission conditions, and it has been observed that the classical Schwarz method can be convergent even without overlap in certain cases. This is in strong contrast to the behavior of classical Schwarz methods applied to elliptic problems, for which overlap is essential for convergence. More recently, optimized Schwarz methods have been developed for elliptic partial differential equations. These methods use more effective transmission conditions between subdomains than the classical Dirichlet conditions, and optimized Schwarz methods can be used both with and without overlap for elliptic problems. We show here why the classical Schwarz method applied to both the time harmonic and time discretized Maxwell's equations converges without overlap: the method has the same convergence factor as a simple optimized Schwarz method for a scalar elliptic equation. Based on this insight, we develop an entire new hierarchy of optimized overlapping and nonoverlapping Schwarz methods for Maxwell's equations with greatly enhanced performance compared to the classical Schwarz method. We also derive for each algorithm asymptotic formulas for the optimized transmission conditions, which can easily be used in implementations of the algorithms for problems with variable coefficients. We illustrate our findings with numerical experiments.

A Priori Error Analysis and Spring Arithmetic

Gilles Chabert and Luc Jaulin

SIAM J. Sci. Comput. 31, pp. 2214-2230 (17 pages)

Online Publication Date: May 07, 2009

Full Text: | Download PDF

Show Abstract
Error analysis is defined by the following concern: Bounding the output variation of a (nonlinear) function with respect to a given variation of the input variables. This paper investigates this issue in the framework of interval analysis. The classical way of analyzing the error is to linearize the function around the point corresponding to the actual input, but this method is local and not reliable. Both drawbacks can be easily circumvented by a combined use of interval arithmetic and domain splitting. However, because of the underlying linearization, a standard interval algorithm leads to a pessimistic bound, and even simply fails (i.e., returns an infinite error) in case of singularity. We propose an original nonlinear approach where intervals are used in a more sophisticated way through the so-called springs. This new structure allows one to represent an (infinite) set of intervals constrained by their midpoints and their radius. The output error is then calculated with a spring arithmetic in the same way as the image of a function is calculated with interval arithmetic. Our method is illustrated on two examples, including an application of geopositioning.

A Conservative Discontinuous Galerkin Semi-Implicit Formulation for the Navier–Stokes Equations in Nonhydrostatic Mesoscale Modeling

Marco Restelli and Francis X. Giraldo

SIAM J. Sci. Comput. 31, pp. 2231-2257 (27 pages) | Cited 2 times

Online Publication Date: May 07, 2009

Full Text: | Download PDF

Show Abstract
A discontinuous Galerkin (DG) finite element formulation is proposed for the solution of the compressible Navier–Stokes equations for a vertically stratified fluid, which are of interest in mesoscale nonhydrostatic atmospheric modeling. The resulting scheme naturally ensures conservation of mass, momentum, and energy. A semi-implicit time-integration approach is adopted to improve the efficiency of the scheme with respect to the explicit Runge–Kutta time integration strategies usually employed in the context of DG formulations. A method is also presented to reformulate the resulting linear system as a pseudo-Helmholtz problem. In doing this, we obtain a DG discretization closely related to those proposed for the solution of elliptic problems, and we show how to take advantage of the numerical integration rules (required in all DG methods for the area and flux integrals) to increase the efficiency of the solution algorithm. The resulting numerical formulation is then validated on a collection of classical two-dimensional test cases, including density driven flows and mountain wave simulations. The performance analysis shows that the semi-implicit method is, indeed, superior to explicit methods and that the pseudo-Helmholtz formulation yields further efficiency improvements.

Finite Element Methods on Very Large, Dynamic Tubular Grid Encoded Implicit Surfaces

Oliver Nemitz, Michael Bang Nielsen, Martin Rumpf, and Ross Whitaker

SIAM J. Sci. Comput. 31, pp. 2258-2281 (24 pages) | Cited 1 time

Online Publication Date: May 20, 2009

Full Text: | Download PDF

Show Abstract
The simulation of physical processes on interfaces and a variety of applications in geometry processing and geometric modeling are based on the solution of partial differential equations on curved and evolving surfaces. Frequently, an implicit level set type representation of these surfaces is the most effective and computationally advantageous approach. This paper addresses the computational problem of how to solve partial differential equations on highly resolved level sets with an underlying very high-resolution discrete grid. These high-resolution grids are represented in a very efficient dynamic tubular grid encoding format for a narrow band. A reaction diffusion model on a fixed surface and surface evolution driven by a nonlinear geometric diffusion approach, by isotropic or truly anisotropic curvature motion, are investigated as characteristic model problems. The proposed methods are based on semi-implicit finite element discretizations directly on these narrow bands, require only standard numerical quadrature, and allow for large time steps. To combine large time steps with a very thin and thus storage inexpensive narrow band, suitable transparent boundary conditions on the boundary of the narrow band and a nested iteration scheme in each time step are investigated. This nested iteration scheme enables the discrete interfaces to move in a single time step significantly beyond the domain of the narrow band of the previous time step. Furthermore, algorithmic tools are provided to assemble finite element matrices and to apply matrix vector operators via fast, cache-coherent access to the dynamic tubular grid encoded data structure. The consistency of the presented approach is evaluated, and various numerical examples show its application potential.

An Accelerating Quasi-Monte Carlo Method for Option Pricing Under the Generalized Hyperbolic Lévy Process

Junichi Imai and Ken Seng Tan

SIAM J. Sci. Comput. 31, pp. 2282-2302 (21 pages) | Cited 1 time

Online Publication Date: May 20, 2009

Full Text: | Download PDF

Show Abstract
In this paper, we develop a simple and yet practically efficient algorithm for simulating high-dimensional exotic options. Our method is based on an extension of Imai and Tan's linear transformation method, which is originally proposed in the context of simulating a Gaussian process. By generalizing this method to other stochastic processes and exploiting the numerical inversion method of Hormann and Leydold, this method can be used to enhance quasi-Monte Carlo method in a wide range of applications. We demonstrate the relative efficiency of our proposed simulation technique using exotic option examples including Asian, lookback, barrier, and cliquet options for which the underlying asset price follows an exponential generalized hyperbolic Lévy process. We also illustrate the impact of our proposed method on dimension reduction.

Enhanced Mass Conservation in Least-Squares Methods for Navier–Stokes Equations

J. J. Heys, E. Lee, T. A. Manteuffel, S. F. McCormick, and J. W. Ruge

SIAM J. Sci. Comput. 31, pp. 2303-2321 (19 pages)

Online Publication Date: May 20, 2009

Full Text: | Download PDF

Show Abstract
There are many applications of the least-squares finite element method for the numerical solution of partial differential equations because of a number of benefits that the least-squares method has. However, one of the most well-known drawbacks of the least-squares finite element method is the lack of exact discrete mass conservation, in some contexts, due to the fact that the least-squares method minimizes the continuity equation in $L^2$-norm. In this paper, we explore the reason for the mass loss and provide new approaches to retain the mass even in a severely underresolved grid.

A Fast $\ell$1-TV Algorithm for Image Restoration

Xiaoxia Guo, Fang Li, and Michael K. Ng

SIAM J. Sci. Comput. 31, pp. 2322-2341 (20 pages) | Cited 1 time

Online Publication Date: May 28, 2009

Full Text: | Download PDF

Show Abstract
Image restoration problems are often solved by finding the minimizer of a suitable objective function consisting of a data-fitting term and a regularization term. In this paper, we consider the data-fitting term measured in the $\ell$1 norm to handle non-Gaussian additive noise and the regularization term given by the total variation (TV) to restore image edges. We propose a new algorithm for this image restoration problem by making use of new variables to modify the data-fitting term and the TV regularization term. An alternating minimization method based on the new formulation is employed to restore blurred and noisy images. Our experimental results show that the quality of restored images by the proposed method is competitive with those restored by the other tested methods. We also show the convergence of the alternating minimization algorithm and demonstrate that the proposed algorithm is very efficient.

Using Global Interpolation to Evaluate the Biot-Savart Integral for Deformable Elliptical Gaussian Vortex Elements

Rodrigo B. Platte, Louis F. Rossi, and Travis B. Mitchell

SIAM J. Sci. Comput. 31, pp. 2342-2360 (19 pages)

Online Publication Date: May 28, 2009

Full Text: | Download PDF

Show Abstract
This paper introduces a new method for approximating the Biot-Savart integral for elliptical Gaussian functions using high-order interpolation and compares it to an existing method based on small aspect ratio asymptotics. The new evaluation technique uses polynomials to approximate the kernel corresponding to the integral representation of the streamfunction. We determine the polynomial coefficients by interpolating precomputed values from look-up tables over a wide range of aspect ratios. When implemented in a full nonlinear vortex method, we find that the new technique is almost three times faster and, unlike the asymptotic method, provides uniform accuracy over the full range of aspect ratios. As a proof-of-concept for large scale computations, we use the new technique to calculate inviscid axisymmetrization and filamentation of a two-dimensional elliptical fluid vortex. We compare our results with those from a pseudo-spectral computation and from electron vortex experiments, and find good agreement among the three approaches.

Bounding Data with a Piecewise Linear Band

Allan R. Willms

SIAM J. Sci. Comput. 31, pp. 2361-2367 (7 pages) | Cited 1 time

Online Publication Date: May 28, 2009

Full Text: | Download PDF

Show Abstract
This paper presents an algorithm for finding a piecewise linear band which bounds given two-dimensional data. The band is of constant height, and the distance between adjacent slope discontinuities in opposite directions is restricted by two user-specified parameters.

Radially Projected Finite Elements

Amnon J. Meir and Necibe Tuncer

SIAM J. Sci. Comput. 31, pp. 2368-2385 (18 pages)

Online Publication Date: May 28, 2009

Full Text: | Download PDF

Show Abstract
We describe and analyze a new finite element discretization for domains with spheroidal geometry. In particular, we describe how the method can be used to approximate solutions as well as eigenvalues and eigenfunctions of partial differential equations posed on the sphere, ellipsoidal shells, and cylindrical shells. These novel, so-called “radially projected finite elements” are particularly attractive for numerical simulations since the resulting finite element discretization is “logically rectangular” and may be easily implemented or incorporated into existing finite element codes.

The $\alpha$-Method Direct Transcription in Path Constrained Dynamic Optimization

Nigam Chandra Parida and Soumyendu Raha

SIAM J. Sci. Comput. 31, pp. 2386-2417 (32 pages)

Online Publication Date: May 28, 2009

Full Text: | Download PDF

Show Abstract
Numerically discretized dynamic optimization problems having active inequality and equality path constraints that along with the dynamics induce locally high index differential algebraic equations often cause the optimizer to fail in convergence or to produce degraded control solutions. In many applications, regularization of the numerically discretized problem in direct transcription schemes by perturbing the high index path constraints helps the optimizer to converge to useful control solutions. For complex engineering problems with many constraints it is often difficult to find effective nondegenerate perturbations that produce useful solutions in some neighborhood of the correct solution. In this paper we describe a numerical discretization that regularizes the numerically consistent discretized dynamics and does not perturb the path constraints. For all values of the regularization parameter the discretization remains numerically consistent with the dynamics and the path constraints specified in the original problem. The regularization is quantifiable in terms of time step size in the mesh and the regularization parameter. For fully regularized systems the scheme converges linearly in time step size. The method is illustrated with examples.
Close

Your Recent History Recent History Help

Recently Viewed

  1. A Three-Dimensional Mixed Finite-Element Approximation of the Semiconductor Energy-Transport Equatio...
  2. Augmented Mixed Finite Element Methods for the Stationary Stokes Equations
  3. Iterative Algorithms Based on Decoupling of Deblurring and Denoising for Image Restoration
  4. Constructing Sobol Sequences with Better Two-Dimensional Projections
  5. Asymptotic Sampling Distribution for Polynomial Chaos Representation from Data: A Maximum Entropy an...
© 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