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

2002

Volume 23, Issue 6, pp. 1817-2182


A Weakly Overlapping Domain Decomposition Preconditioner for the Finite Element Solution of Elliptic Partial Differential Equations

Randolph E. Bank, Peter K. Jimack, Sarfraz A. Nadeem, and Sergei V. Nepomnyaschikh

SIAM J. Sci. Comput. 23, pp. 1817-1841 (25 pages) | Cited 6 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
We present a new two-level additive Schwarz domain decomposition preconditioner which is appropriate for use in the parallel finite element solution of elliptic partial differential equations (PDEs). As with most parallel domain decomposition methods each processor may be assigned one or more subdomains, and the preconditioner is such that the processors are able to solve their own subproblem(s) concurrently. The novel feature of the technique proposed here is that it requires just a single layer of overlap in the elements which make up each subdomain at each level of refinement, and it is shown that this amount of overlap is sufficient to yield an optimal preconditioner. Some numerical experiments---posed in both two and three space dimensions---are included to confirm that the condition number when using the new preconditioner is indeed independent of the level of mesh refinement on the test problems considered.

A Trust-Region Approach to the Regularization of Large-Scale Discrete Forms of Ill-Posed Problems

Marielba Rojas and Danny C. Sorensen

SIAM J. Sci. Comput. 23, pp. 1842-1860 (19 pages) | Cited 14 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
We consider large-scale least squares problems where the coefficient matrix comes from the discretization of an operator in an ill-posed problem, and the right-hand side contains noise. Special techniques known as regularization methods are needed to treat these problems in order to control the effect of the noise on the solution. We pose the regularization problem as a quadratically constrained least squares problem. This formulation is equivalent to Tikhonov regularization, and we note that it is also a special case of the trust-region subproblem from optimization. We analyze the trust-region subproblem in the regularization case and we consider the nontrivial extensions of a recently developed method for general large-scale subproblems that will allow us to handle this case. The method relies on matrix-vector products only, has low and fixed storage requirements, and can handle the singularities arising in ill-posed problems. We present numerical results on test problems, on an inverse interpolation problem with field data, and on a model seismic inversion problem with field data.

Hidden Discontinuities and Parametric Sensitivity Calculations

John E. Tolsma and Paul I. Barton

SIAM J. Sci. Comput. 23, pp. 1861-1874 (14 pages) | Cited 11 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
Parametric sensitivity analysis is becoming a very important tool when studying differential equations. However, potential pitfalls exist if this analysis is applied naively. In particular, if the model contains discontinuities, then the application of standard numerical codes will typically result in incorrect sensitivity trajectories. The purpose of this paper is to demonstrate the problems associated with "hidden discontinuities" during sensitivity analysis (i.e., discontinuities not handled explicitly) and to present a code analysis and transformation approach through which differential equations containing discontinuities can be solved efficiently and correctly with minimal effort required by the modeler.

On the Numerical Solution of $(\lambda^2 A + \lambda B + C), x = b$ and Application to Structural Dynamics

V. Simoncini and F. Perotti

SIAM J. Sci. Comput. 23, pp. 1875-1897 (23 pages) | Cited 14 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
In this paper we address the numerical solution of a large linear system depending quadratically on a parameter that varies in a wide range. We analyze a solution method, whose computational cost grows only sublinearly with the number of parameters, that relies on the use of an indefinite inner product. Important implementation aspects are treated in detail. The problem arises in various application areas: we shall report on our experience with cases in structural dynamics.

Residual and Backward Error Bounds in Minimum Residual Krylov Subspace Methods

Christopher C. Paige and Zdenvek Strakos

SIAM J. Sci. Comput. 23, pp. 1898-1923 (26 pages) | Cited 4 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
Minimum residual norm iterative methods for solving linear systems Ax=b can be viewed as, and are often implemented as, sequences of least squares problems involving Krylov subspaces of increasing dimensions. The minimum residual method (MINRES) [C. Paige and M. Saunders, SIAM J. Numer. Anal., 12 (1975), pp. 617--629] and generalized minimum residual method (GMRES) [Y. Saad and M. Schultz, SIAM J. Sci. Statist. Comput., 7 (1986), pp. 856--869] represent typical examples. In [C. Paige and Z. Strakos, Bounds for the least squares distance using scaled total least squares, Numer. Math., to appear] revealing upper and lower bounds on the residual norm of any linear least squares (LS) problem were derived in terms of the total least squares (TLS) correction of the corresponding scaled TLS problem. In this paper theoretical results of [C. Paige and Z. Strakos, Bounds for the least squares distance using scaled total least squares, Numer. Math., to appear] are extended to the GMRES context. The bounds that are developed are important in theory, but they also have fundamental practical implications for the finite precision behavior of the modified Gram--Schmidt implementation of GMRES, and perhaps for other minimum norm methods.

Adaptive Reduced-Order Controllers for a Thermal Flow System Using Proper Orthogonal Decomposition

S. S. Ravindran

SIAM J. Sci. Comput. 23, pp. 1924-1942 (19 pages) | Cited 29 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
An adaptive reduced-order controller design is presented for flow control using proper orthogonal decomposition (POD). In reduced-order controller design, the idea is to start with an ensemble of data obtained from numerical simulation of the underlying partial differential equations (PDEs). POD is then used to obtain a reduced set of basis functions which is then used to derive a reduced-order model of the PDEs via Galerkin projection. This reduced-order model allows us to derive a reduced-order controller. However, it is not clear, a priori, what is the best way to obtain an ensemble of data that would give basis functions that represent the influence of the control action on the system. In this paper we explore an adaptive procedure for reduced-order controller design that improves the reduced-order model by successively updating the ensemble of data during the optimization iterations. We illustrate this method on a control problem in thermal flow system modeled by a thermally coupled Navier--Stokes equations. Numerical results are presented for a vorticity regulation problem in fluid flows using boundary temperature as control mechanism. Through our numerical experiments we demonstrate the feasibility and applicability of the adaptive reduced-order controllers.

Implicit Integration of the Time-Dependent Ginzburg--Landau Equations of Superconductivity

D. O. Gunter, H. G. Kaper, and G. K. Leaf

SIAM J. Sci. Comput. 23, pp. 1943-1958 (16 pages) | Cited 2 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
This article is concerned with the integration of the time-dependent Ginzburg--Landau (TDGL) equations of superconductivity. Four algorithms, ranging from fully explicit to nonlinearly implicit, are presented and evaluated for stability, accuracy, and compute time. The benchmark problem for the evaluation is the equilibration of a vortex configuration in a superconductor that is embedded in a thin insulator and subject to an applied magnetic field.

A Domain Decomposition Method for Advection-Diffusion Processes with Application to Blood Solutes

Alfio Quarteroni, Alessandro Veneziani, and Paolo Zunino

SIAM J. Sci. Comput. 23, pp. 1959-1980 (22 pages) | Cited 4 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
In the present paper we consider a heterogeneous model for the dynamics of a blood solute both in the vascular lumen and inside the arterial wall. In the lumen, we consider an advection-diffusion equation, where the convective field is provided by the velocity of blood, which is in turn obtained by solving the Navier--Stokes equations. Inside the arterial wall we consider a pure diffusive dynamics. Since the endothelial layer at the interface between the lumen and the wall acts as a permeable membrane, whose permeability depends on the shear rate exerted by the blood, the solute concentration is discontinuous across this membrane. A possible approach for the numerical study of this kind of problem is inspired by domain decomposition techniques. In particular, we introduce a splitting in the computation and alternate the solution of the advection-diffusion equation in the lumen with that of the diffusion equation in the wall. We set up an efficient iterative method, based on a suitable reformulation of the problem in terms of a Steklov--Poincaré interface equation. This formulation is a nonstandard one because of the concentration discontinuity at the lumen-wall interface and plays a key role in the proof of convergence of our method. In particular, we prove that the convergence rate performed by the proposed method is independent of the finite element space discretization and provides a criterion for the selection of an acceleration parameter.
Several numerical results, referred to as biomedical applications, support our theoretical conclusions and illustrate the efficiency of this algorithm.

Convergence of Nonconvergent IRK Discretizations of Optimal Control Problems with State Inequality Constraints

J. T. Betts, N. Biehn, and S. L. Campbell

SIAM J. Sci. Comput. 23, pp. 1981-2007 (27 pages) | Cited 8 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
It has been observed that optimization codes are sometimes able to solve inequality state constrained optimal control problems with discretizations which do not converge when used as integrators on the constrained dynamics. Understanding this phenomenon could lead to a more robust design for direct transcription codes as well as better test problems. This paper examines how this phenomenon can occur. The difference between solving index 3 differential algebraic equations (DAEs) using the trapezoid method in the context of direct transcription for optimal control problems and a straightforward implicit Runge--Kutta (IRK) formulation of the same trapezoidal discretization is analyzed. It is shown through numerical experience and theory that the two can differ as much as O(1/h3) in the control. The optimization can use a small sacrifice in the accuracy of the states in the early stages of the trapezoidal method to produce better accuracy in the control, whereas more precise solutions converge to an incorrect solution. Convergence independent of the index of the constraints is then proven for one class of problems. The theoretical results are used to explain computational observations.

Principal Angles between Subspaces in an A-Based Scalar Product: Algorithms and Perturbation Estimates

Andrew V. Knyazev and Merico E. Argentati

SIAM J. Sci. Comput. 23, pp. 2008-2040 (33 pages) | Cited 5 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
Computation of principal angles between subspaces is important in many applications, e.g., in statistics and information retrieval. In statistics, the angles are closely related to measures of dependency and covariance of random variables. When applied to column-spaces of matrices, the principal angles describe canonical correlations of a matrix pair. We highlight that all popular software codes for canonical correlations compute only cosine of principal angles, thus making impossible, because of round-off errors, finding small angles accurately. We review a combination of sine and cosine based algorithms that provide accurate results for all angles. We generalize the method to the computation of principal angles in an A-based scalar product for a symmetric and positive definite matrix A. We provide a comprehensive overview of interesting properties of principal angles. We prove basic perturbation theorems for absolute errors for sine and cosine of principal angles with improved constants. Numerical examples and a detailed description of our code are given.

Fourth Order Chebyshev Methods with Recurrence Relation

Assyr Abdulle

SIAM J. Sci. Comput. 23, pp. 2041-2054 (14 pages) | Cited 29 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
In this paper, a new family of fourth order Chebyshev methods (also called stabilized methods) is constructed. These methods possess nearly optimal stability regions along the negative real axis and a three-term recurrence relation. The stability properties and the high order make them suitable for large stiff problems, often space discretization of parabolic PDEs. A new code ROCK4 is proposed, illustrated at several examples, and compared to existing programs.

Hamilton--Jacobi Equations on a Manifold and Applications to Grid Generation or Refinement

Ph. Hoch and M. Rascle

SIAM J. Sci. Comput. 23, pp. 2055-2073 (19 pages) | Cited 1 time

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
In this paper, we consider Hamilton--Jacobi equations on a manifold, typically on the graph of some previously computed function z(x,y), and we show how the corresponding level set method allows us to generate and/or to refine a mesh in regions where this function z has large derivatives. Such as it is, the method needs to be strongly improved and accelerated, but the principle is awfully natural, and in principle the method is fully automatic. Similar ideas are also useful in image processing, in particular for the active contours method.

Chebyshev Spectral Methods for Radiative Transfer

Arnold D. Kim and Miguel Moscoso

SIAM J. Sci. Comput. 23, pp. 2074-2094 (21 pages) | Cited 13 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
We study the performance of Chebyshev spectral methods for time-dependent radiative transfer equations. Starting with a method for one-dimensional problems in homogeneous media, we show that the modifications needed to consider more general problems such as inhomogeneous media, polarization, and higher dimensions are straightforward. In this method, we approximate the spatial dependence of the intensity by an expansion of Chebyshev polynomials. This yields a coupled system of integro-differential equations for the expansion coefficients that depend on angle and time. Next, we approximate the integral operation on the angle variables using a Gaussian quadrature rule resulting in a coupled system of differential equations with respect to time. Using a second-order finite difference approximation, we discretize the time variable. We solve the resultant system of equations with an efficient algorithm that makes Chebyshev spectral methods competitive with other methods for radiative transfer equations.

Quadratic Convergence for Valuing American Options Using a Penalty Method

P. A. Forsyth and K. R. Vetzal

SIAM J. Sci. Comput. 23, pp. 2095-2122 (28 pages) | Cited 45 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
The convergence of a penalty method for solving the discrete regularized American option valuation problem is studied. Sufficient conditions are derived which both guarantee convergence of the nonlinear penalty iteration and ensure that the iterates converge monotonically to the solution. These conditions also ensure that the solution of the penalty problem is an approximate solution to the discrete linear complementarity problem. The efficiency and quality of solutions obtained using the implicit penalty method are compared with those produced with the commonly used technique of handling the American constraint explicitly. Convergence rates are studied as the timestep and mesh size tend to zero. It is observed that an implicit treatment of the American constraint does not converge quadratically (as the timestep is reduced) if constant timesteps are used. A timestep selector is suggested which restores quadratic convergence.

A New Method for Micro-Macro Simulations of Viscoelastic Flows

C. Chauviere

SIAM J. Sci. Comput. 23, pp. 2123-2140 (18 pages) | Cited 4 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
In this paper we develop a robust numerical method derived from the Brownian configuration field method [M.A. Hulsen, A.P.G. Van Heel, and B.H.A.A. Van Den Brule, J. Non-Newtonian Fluid Mech., 70 (1997), pp. 79--101] for the simulation of flows of dilute polymeric solutions. The statistical properties of the Wiener stochastic process in the stochastic differential equation describing the evolution of the configuration vector are exploited in order to derive a simple expression for the polymeric contribution to the stress. The method is tested numerically by solving the benchmark problem of the flow of an Oldroyd B fluid past a cylinder in a channel using a spectral element method. Results are presented to demonstratethe advantages of the proposed method.

Multilevel Method for Mixed Eigenproblems

R. Hiptmair and K. Neymeyr

SIAM J. Sci. Comput. 23, pp. 2141-2164 (24 pages) | Cited 5 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
For a Lipschitz-polyhedron $\Omega\subset\mathbb{R}^3$ we consider eigenvalue problems $\operatorname{\bf curl}\alpha\operatorname{\bf curl}\mathbf{u}=\lambda\mathbf{u}$ and $\operatorname{grad}\alpha\operatorname{div}\mathbf{u} = \lambda\mathbf{u}$, $\lambda>0$, set in $\boldsymbol{H}(\operatorname{\bf curl};\Omega)$ and $\boldsymbol{H}(\operatorname{div};\Omega)$. They are discretized by means of the conforming finite elements introduced by Nédélec. The preconditioned inverse iteration in its subspace variant is adapted to these problems. A standard multigrid scheme serves as the preconditioner. The main challenge arises from the large kernels of the operators curl and div. However, thanks to the choice of finite element spaces these kernels have a direct representation through the gradients/rotations of discrete potentials. This makes it possible to use a multigrid iteration in potential space to obtain approximate projections onto the orthogonal complements of the kernels. There is ample evidence that this will lead to an asymptotically optimal method. Numerical experiments confirm the excellent performance of the method even on very fine grids.

A Block Orthogonalization Procedure with Constant Synchronization Requirements

Andreas Stathopoulos and Kesheng Wu

SIAM J. Sci. Comput. 23, pp. 2165-2182 (18 pages) | Cited 5 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
First, we consider the problem of orthonormalizing skinny (long) matrices. We propose an alternative orthonormalization method that computes the orthonormal basis from the right singular vectors of a matrix. Its advantages are that (a) all operations are matrix-matrix multiplications and thus cache efficient, (b) only one synchronization point is required in parallel implementations, and (c) it is typically more stable than classical Gram--Schmidt (GS). Second, we consider the problem of orthonormalizing a block of vectors against a previously orthonormal set of vectors and among itself. We solve this problem by alternating iteratively between a phase of GS and a phase of the new method. We provide error analysis and use it to derive bounds on how accurately the two successive orthonormalization phases should be performed to minimize total work performed. Our experiments confirm the favorable numerical behavior of the new method and its effectiveness on modern parallel computers.
Close

Your Recent History Recent History Help

Recently Viewed

  1. Iterative Solution of Cyclically Reduced Systems Arising from Discretization of the Three-Dimensiona...
  2. Nonoscillatory Central Schemes for Multidimensional Hyperbolic Conservation Laws
  3. A Fast Multilevel Algorithm for the Solution of Nonlinear Systems of Conductive-Radiative Heat Trans...
  4. Parallel Newton--Krylov--Schwarz Algorithms for the Transonic Full Potential Equation
  5. Numerical Solution of the Riemann Problem for Two-Dimensional Gas Dynamics
© 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