Technical Reports (1994)
Effective Computation of the Analytic Center of the Solution Set in Linear Programming Using Primal-Dual Interior-Point Methods
María D. González-Lima
The centrality property satisfied by the analytic center of the solution set makes its computation very valuable for some linear programming applications. One such application coming from the economic and management sciences is Data Envelopment Analysis (DEA). In DEA one desires a solution of the underlying linear programming model that is in the relative interior of the solution set and one that is in some sense as far away as possible from the relative small changes in the data. In this work we study the effective computation of the analytic center solution by the use of primal-dual interior-point methods. We present an unified study of existing theoretical results for primal-dual interior-point algorithms as they concern the convergence of the iteration sequence and the convergence of the iteration sequence to the analytic center. These theoretical results are evaluated from the point of view of the practical computation of the analytic center. We propose a primal-dual interior-point algorithm for effectively computing the analytic center of the solution set. The algorithm proposed combines good theoretical and numerical properties and its ability to solve real world problems from the DEA application is demonstrated.
December 1994
Controllability Issues for Flow-Related Models: A Computational Approach
Martin Berggren & Roman Glowinski
We discuss the numerical solution of some controllability problems for time-dependent flow models. The emphasis is on algorithmic aspects, discretization issues, and memory-saving devices. Numerical results are presented for controllability studies involving the viscous Burgers equation, an advection-diffusion equation, and the unsteady Stokes equations.
December 1994
A Spectral Preconditioner for Control Problems Associated with Linear Evolution Equations
Martin Berggren & Roman Glowinski
We introduce a spectral preconditioner for control problems associated with first-order temporary evolution equations involving an elliptic, selfadjoint operator. Condition number estimates are derived, and we describe in detail how to efficiently implement a conjugate gradient algorithm using the preconditioner. Numerical results of a control problem involving the heat equation in two space dimensions show that a very limited spectral information is sufficient to greatly reduce the number of iterations in the conjugate gradient algorithm.
December 1994
[East-West J. Numer. Math. 3 (1995), no. 2, 81-109.]
Trust-Region Interior-Point SQP Algorithms for a Class of Nonlinear Programming Problems
J.E. Dennis, Matthias Heinkenschloss, & Luís N. Vicente
In this paper a family of trust-region interior-point SQP algorithms for the solution of minimization problems with nonlinear equality constraints and simple bounds on some of the variables is described and analyzed. Such nonlinear programs arise e.g. from the discretization of optimal control problems. The algorithms treat states and controls as independent variables. They are designed to take advantage of the structure of the problem. In particular they do not rely on matrix factorizations of the linearized constraints, but use solutions of the linearized state equation and the adjoint equation. They are well suited for large scale problems arising from optimal control problems governed by partial differential equations.
The algorithms keep strict feasibility with respect to the bound constraints by using a primal-dual affine scaling method proposed for a different class of problems by Coleman and Li and they exploit trust-region techniques for equality-constrained optimization. Thus, they allow the computation of the steps using a variety of methods, including many iterative techniques.
Global convergence of these algorithms to a first-order KKT limit point is proved under very mild conditions on the trial steps. Under reasonable, but more stringent conditions on the quadratic model and on the trial steps, the sequence of iterates generated by the algorithms is shown to have a limit point satisfying the second-order necessary KKT conditions. The local rate of convergence to a nondegenerate strict local minimizer is q-quadratic. The results given here include as special cases current results for only equality constraints and for only simple bounds.
Numerical results for the solution of an optimal control problem governed by a nonlinear heat equation are reported.
December 1994 (Revised November 1995, December 1996, and September 1997)
[SIAM J. Control Optim. 36 (1998), no. 5, 1750-1794.]
The Solution of the Metric STRESS and SSTRESS Problems in Multidimensional Scaling Using Newton's Method
A.J. Kearsley, R.A. Tapia, & M.W. Trosset
This paper considers numerical algorithms for finding local minimizers of metric multidimensional scaling problems. The two most common optimality criteria (STRESS and SSTRESS) are considered, the leading algorithms for each are carefully explicated, and a new algorithm is proposed. The new algorithm is based on Newton's method and relies on a parameterization that has not previously been used in multidimensional scaling algorithms. In contrast to previous algorithms, a very pleasant feature of the new algorithm is that it can be used with either the STRESS or the SSTRESS criterion. Numerical results are presented for the metric STRESSS problem. These results are quite satisfying and, among other things, suggest that the well-known SMACOF-I algorithm tends to stop prematurely.
December 1994
[Computational Statistics, 13, no. 3 (1998), 369-396]
On the Construction of Strong Complementarity Slackness for DEA Linear Programming Problems Using a Primal-Dual Interior-Point Method
María D. González-Lima, Richard A. Tapia, & Robert M. Thrall
A novel approach for solving the DEA linear programming problems using a primal-dual interior-point method is presented. The solution found by this method satisfies the Strong Complementarity Slackness Condition (SCSC) and maximizes the product of the positive components among all SCSC solutions. This first property is critical in the use of DEA and the second one contributes significantly to the reliability of the solution.
November 1994
[Ann. Oper. Res. 66 (1996), 139-162.]
Trust-Region Interior-Point Algorithms for Minimization Problems with Simple Bounds
J.E. Dennis & Luís N. Vicente
Two trust-region interior-point algorithms for the solution of minimization problems with simple bounds are analyzed and tested. The algorithms scale the local model in a way similar to Coleman and Li [1]. The first algorithm is more usual in that the trust region and the local quadratic model are consistently scaled. The second algorithm proposed here uses an unscaled trust region. A global convergence result for these algorithms is given and dogleg and conjugate-gradient algorithms to compute trial steps are introduced. Some numerical examples that show the advantages of the second algorithm are presented.
November 1994
[In Applied Mathematics and Parallel Computing, 97-107, Physica, Heidelberg, 1996.]
A Chemical Compositional Reservoir Simulator on Distributed Memory Parallel Computers: Comparative Parallel-UTCHEM Simulation Performance Study (Part I)
M. Rame & M. Delshad
This paper presents the application of distributed memory parallel computers to field scale reservoir simulations using a parallel version of UTCHEM, The University of Texas Chemical Flooding Simulator. The model is a general purpose highly vectorized chemical compositional simulator that can simulate a wide range of displacement processes at both field and laboratory scales. The original simulator was modified to run on both distributed memory parallel machines (Intel iPSC/860 and Delta, Connection Machine 5, Kendall Square 1 and 2, and CRAY T3D) and a cluster of workstations.
A domain decomposition approach has been taken towards parallelization of the code. A portion of the discrete reservoir model is assigned to each processor by a set-up routine that attempts a data layout as even as possible from the load-balance standpoint. Each of these subdomains is extended so that data can be shared between adjacent processors for stencil computation. The added routines that make parallel execution possible are written in a modular fashion that makes the porting to new parallel platforms straight forward.
Results of the distributed memory computing performance of Parallel simulator are presented for field scale applications such as tracer flood and polymer flood. A comparison of the wall-clock times for same problems on a vector supercomputer is also presented.
November 1994
Logically Rectangular Mixed Methods for Darcy Flow on General Geometry
T. Arbogast, P.T. Keenan, M.F. Wheeler, & I. Yotov
We consider an expanded mixed finite element formulation (cell centered finite differences) for Darcy flow with a tensor absolute permeability. The reservoir can be geometrically general with internal features, but the computational domain is rectangular. The method is defined on a curvilinear grid that need not be orthogonal, obtained by mapping the rectangular, computational grid. The original flow problem becomes a similar problem with a modified permeability on the computational grid. Quadrature rules turn the mixed method into a cell-centered finite difference method with a 9 point stencil in 2-D and 19 in 3-D.
As shown by theory and experiment, if the modified permeability on the computational domain is smooth, then the convergence rate is optimal and both pressure and velocity are superconvergent at certain points. If not, Lagrange multiplier pressures can be introduced on boundaries of elements so that optimal convergence is retained. This modification presents only small changes in the solution process; in fact, the same parallel domain decomposition algorithms can be applied with little or no change to the code if the modified permeability is smooth over the subdomains.
This Lagrange multiplier procedure can be used to extend the difference scheme to multi-block domains, and to give a coupling with unstructured grids. In all cases, the mixed formulation is locally conservative. Computational results illustrate the advantage and convergence of this method.
October 1994
Large Time Asymptotics in Contaminant Transport in Porous Media
C.N. Dawson, C.J. van Duijn, & R.E. Grundy
In this paper we derive large time solutions of the partial differential equations modelling contaminant transport in porous media for initial data with bounded support. While the main emphasis is on two space dimensions, for the sake of completeness we give a brief summary of the corresponding results for one space dimension. The philosophy behind the paper is to compare the results of a formal asymptotic analysis of the governing equations as t -> infinity with numerical solutions of the complete initial value problem. the analytic results are obtained using the method of "asymptotic balancing" which identifies the dominant terms in the model equations determining the behaviour of the solution in the large time limit. These are found in terms of time scaled space similarity variables and the procedure results in a reduction of the number of independent variables in the original partial differential equation. This generates what we call a reduced equation the solution of which depends crucially on the value of a parameter appearing in the problem. In some cases the reduced equation can be solved explicitly while in others they have a particularly intractable structure which inhibits any analytic or numerical progress. However we can extract a number of global and local properties of the solution which enables us to form a reasonably complete picture of what the profiles look like. Extensive comparison with numerical solution of the original initial value problem provides convincing confirmation of our analytic solutions. In the final section of the paper, by way of motivation for the work, we give some results concerning the temporal behaviour of certain moments of the two dimensional profiles commonly used to compute physical parameter characteristics for contaminant transport in porous media.
November 1994
[SIAM J. Appl. Math. 56 (1996), no. 4, 965-993.]
A Trust-Region Approach to Nonlinear Systems of Equalities and Inequalities
J.E. Dennis Jr., Mahmoud El-Alem, & Karen Williamson
In this paper, two new trust-region algorithms for the numerical solution of systems of nonlinear equalities and inequalities are introduced. The formulation is free of arbitrary parameters and possesses sufficient smoothness to exploit the robustness of the trust-region approach. The proposed algorithms are one-sided least-squares trust-region algorithms. The first algorithm is a single-model algorithm, and the second one is a multi-model algorithm where the Cauchy point computation is a model selection procedure.
Global convergence analyses for the two algorithms are presented. Our analysis generalizes to nonlinear systems of equalities and inequalities the well-developed theory for nonlinear least-squares problems.
Numerical experiments on the two algorithms are also presented. The performances of the two algorithms are reported. The numerical results validate the effectiveness of our approach.
October 1994
[SIAM J. Optim. 9 (1999), no. 2, 291-315.]
A Robust Choice of the Lagrange Multiplier in the SQP Newton Method
Debora Cores & R.A. Tapia
We study the choice of the Lagrange multipliers in the successive quadratic programming method (SQP) for equality constrained optimization.
It is known that the augmented Lagrangian SQP-Newton method depends on the penalty parameter only through the multiplier in the Hessian matrix of the Lagrangian function. This effectively reduces the augmented Lagrangian SQP-Newton method to the Lagrangian SQP-Newton method where only the multiplier estimate depends on the penalty parameter. In this work, we construct a multiplier estimate that depends strongly on the penalty parameter and derive a choice for the penalty parameter that attempts to make the Hessian matrix, restricted to the tangent space of the constraints, positive definite and well conditioned. We demonstrate that the SQP-Newton method with this choice of Lagrange multipliers is locally and q-quadratically convergent. Considerable numerical experimentation is included and shows that our approach merits further investigation.
October 1994
On the Convergence Theory of Trust-Region-Based Algorithms for Equality-Constrained Optimization
John E. Dennis & Luís N. Vicente
In a recent paper, Dennis, El-Alem, and Maciel proved global convergence to a stationary point for a general trust-region-based algorithm for equality-constrained optimization. This general algorithm is based on appropriate choices of trust-region subproblems and seems particularly suitable for large problems.
This paper shows global convergence to a point satisfying the second-order necessary optimality conditions for the same general trust-region-based algorithm. The results given here can be seen as a generalization of the convergence results for trust-region methods for unconstrained optimization obtained by Moré and Sorensen. The behavior of the trust radius and the local rate of convergence are analyzed. Some interesting facts concerning the trust-region subproblem for the linearized constraints, the quasi-normal component of the step, and the hard case are presented.
It is shown how these results can be applied to a class of discretized optimal control problems.
October 1994 (Revised September 1995)
[SIAM J. Optim. 7 (1997), no. 4, 927-950.]
Stability Analysis of Finite-Difference Schemes for the Viscoelastic Wave Equation
Joakim O. Blanch & William W. Symes
It is difficult to predict stability properties of a finite difference scheme. It has to be investigated through the roots of the Z-transformed and Fourier transformed difference scheme (modal equation). To simultaneously investigate several schemes for the viscoelastic wave equation, it is possible to derive the modal equation with parameterized coefficients. Several conditionally stable schemes were found, where the most efficient is a staggered scheme with a stability condition closely resembling that of an elastic scheme.
October 1994
Algorithms for Bilevel Optimization
Natalia Alexandrov & J.E. Dennis, Jr.
General multilevel nonlinear optimization problems arise in design of complex systems and can be used as a means of regularization for multicriteria optimization problems. Here for clarity in displaying our ideas, we restrict ourselves to general bilevel optimization problems, and we present two solution approaches. Both approaches use a trust-region globalization strategy, and they can be easily extended to handle the general multilevel problem. We make no convexity assumptions, but we do assume that the problem has a nondegenerate feasible set. We consider necessary optimality conditions for the bilevel problem formulations and discuss results that can be extended to obtain multilevel optimization formulations with constraints at each level.
August 1994
A Comparison of Nonlinear Programming Approaches to an Elliptic Inverse Problem and a New Domain Decomposition Approach
J.E. Dennis, Jr. & Robert Michael Lewis
We compare three nonlinear programming approaches to a well-known elliptic inverse problem in three spatial dimensions. Two of these approaches may be viewed as conventional; the third approach is new and is based on a domain decomposition technique for the solution of the governing elliptic equation. We discuss the benefits that may be obtained from treating the governing differential equation in approaches may be viewed as conventional;the third approach is new an inverse problem as equality constraints in the optimization problem. We present numerical results and discuss the relative efficacy of the three approaches.
September 1994
Numerical Simulation and Optimal Shape for Viscous Flow by a Fictitious Domain Method
Roland Glowinski, Anthony J. Kearsley, Tsorng-Whay Pan, & Jacques Periaux
In this article we discuss the fictitious domain solution of the Navier-Stokes equations modelling unsteady incompressible viscous flow. The method is based on a Lagrange multiplier treatment of the boundary conditions to be satisfied and is particularly well suited to the treatment of no-slip boundary conditions. This approach allows the use of structured meshes and fast specialized solvers for problems on complicated geometries. Another interesting feature of the fictitious domain approach is that it allows the solution of optimal shape problems without regriding. The resulting methodology is applied to the solution of flow problems including external incompressible viscous flow modelled by the Navier-Stokes equations and then to an optimal shape problem for Stokes and Navier-Stokes flow.
August 1994
[Internat. J. Numer. Methods Fluids 20 (1995), no. 8-9, 695-711.]
Implementation of the Conjugate Gradient Algorithm in DSO
Susan E. Minkoff
Computation of the inner state parameters in DSO inversion requires solving a large normal matrix system. A combined conjugate gradient and Lanczos iterative technique can be used to both solve the system and approximate some of the spectrum of the normal operator. At each iteration of the conjugate gradient algorithm, a small tridiagonal matrix (of dimension equal to the number of iterations) is created which has extreme eigenvalues approximating those of the original matrix. A second matrix whose columns are the normalized residual vectors from the conjugate gradient algorithm allows the corresponding eigenvectors to be computed as well if desired. Implemented so that it can be applied to different DSO inversion problems, the conjugate gradient code provides the user with tools to analyze the condition of the problem as well as the quality of the inversion results. Storage of the residual (Lanczos) vectors may be costly if large problems are being solved. Numerical experiments indicate that the largest eigenvalue is the best approximation in the spectrum and the smallest is generally the next best.
August 1994
RUF 1.0 User Manual
Philip T. Keenan
(Abstract not available)
August 1994
Derivative-Free Pattern Search Methods for Multidisciplinary Design Problems
J.E. Dennis, Jr. & Virginia J. Torczon
There have been interesting recent developments in methods for solving optimization problems without making use of derivative (sensitivity) information. While calculus based methods that employ derivative information can be extremely efficient and very effective, they are not applicable to all MDO problems, for instance, when the function to be optimized is nondifferentiable, when sensitivity information is not available or is not reliable, or when the function values are inaccurate. In these settings, we have found that the multidirectional search method, a derivative-free method we have developed for solving nonlinear optimization problems, can be used effectively. Our analysis of the multidirectional search algorithm has led us to discover that its algebraic structure and resulting convergence theory can be related to an entire class of derivative-free methods, which we now call pattern search methods, that have been in use for decades.
The goal of this paper is to give an introduction to pattern search methods, to describe the features they share by using coordinate search, one of the earliest and best-known (if not as effective) pattern search methods, as an example, and to review some recent developments that suggest that these methods can be extended to handle problems with a mix of continuous and discrete variables - another situation that can arise in MDO problems. Finally, we will discuss when these methods are an appropriate choice for solving MDO problems.
August 1994
Sub-Surface Flow and Transport Parameter Identification via Inverse Simulation
Martha M. Carey
Minimization of a nonlinear output least squares objective function consisting of both concentration and pressure residuals is used to estimate parameters representing hydraulic conductivities. Details of objective function construction and behavior, numerical optimization, forward simulation of flow and transport, mapping of parameters to computational hydraulic conductivities, computational needs and numerical results are discussed. Inclusion of concentration data creates a smooth and convex objective function from which hydraulic conductivities can be retrieved with accuracy.
July 1994 (Revised August 1994)
Minimization of Large Scale Quadratic Function Subject to an Ellipsoidal Constraint
D.C. Sorensen
An important problem in linear algebra and optimization is the Trust Region Problem: Minimize a quadratic function subject to an ellipsoidal constraint. This basic problem has several important large scale applications including seismic inversion and forcing convergence in optimization methods. Existing methods to solve the trust region problem require matrix factorizations that are not feasible in the large scale setting. This paper presents an algorithm for solving the large scale trust region problem that requires a fixed size limited storage proportional to order of the quadratic and that relies only on matrix-vector products. The algorithm recasts the trust region problem in terms of parameterized eigenvalue problem and adjusts the parameter with a superlinearly convergent iteration to find the optimal solution from the eigenvector of the parameterized problem. Only the smallest eigenvalue and corresponding eigenvector of the parameterized problem needs to be computed.
July 1994
An Inexact Hybrid Algorithm for Nonlinear Systems of Equations
Mohammedi El Hallabi
In this work we define a hybrid algorithm for approximating zeros of the nonlinear systems F(x)=0, where F: R^n -> R is continuously differentiable. We are concerned with the possibility that n may be large and the Jacobian F´(x) sparse or singular. Trust region globalization methods are known to be robust and can be applied successfully to obtain global convergence results under rather weak hypotheses. However,these algorithms can be expensive, especially for large problems, if the trust region radius needs to be reduced quite often before an acceptable step is obtained. Exploiting the convex structure of the local model subproblem, we propose a hybrid algorithm that uses both trust region and linesearch globalization strategies. It solves, once and not accurately, a local model to obtain a search direction and then uses a linesearch technique to obtain an acceptable steplength. We demonstrate, under rather weak hypotheses, that the algorithm is is globally convergent and that the sequence of residuals converges to zero. Moreover, under standard assumptions of Newton's method theory, we prove that the rate of convergence is q-superlinear. Furthermore, q-quadratic convergence can be obtained by requiring sufficient accuracy in the solution of the local model trust region subproblem.
July 1994 (Revised April 1995)
Choosing the Forcing Terms in an Exact Newton Method
Stanley C. Eisenstat & Homer F. Walker
An inexact Newton method is a generalization of Newton's method for solving F(x)=0, F: R^n -> R^n, in which, at the kth iteration, the step sk from the current approximate solution xk is required to satisfy a condition ||F(xk)+F´(xk) sk|| <= etak ||F(xk)|| for a "forcing term" etak in [0,1). In typical applications, the choice of the forcing terms is critical to the efficiency of the method and can affect robustness as well. Promising choices of the forcing terms are given, their local convergence properties are analyzed, and their practical performance is shown on a representative set of test problems.
June 1994 (Revised February 1995)
[SIAM J. Sci. Comput. 17 (1996), no. 1, 16-32.]
Multilevel Algorithms for Nonlinear Optimization
Natalia Alexandrov & John E. Dennis, Jr.
Multidisciplinary design optimization (MDO) gives rise to nonlinear optimization problems characterized by a large number of constraints that naturally occur in blocks. We propose a class of multilevel optimization methods motivated by the structure and number of constraints and by the expense of the derivative computations for MDO. The algorithms are an extension to the nonlinear programming problem of the successful class of local Brown-Brent algorithms for nonlinear equations. Our extensions allow the user to partition constraints into arbitrary blocks to fit the application, and they separately process each block and the objective function, restricted to certain subspaces. The methods use trust regions as a globalization strategy, and they have been shown to be globally convergent under reasonable assumptions. The multilevel algorithms can be applied to all classes of MDO formulations. Multilevel algorithms for solving nonlinear systems of equations are a special case of the multilevel optimization methods. In this case, they can be viewed as a trust-region globalization of the Brown-Brent class.
June 1994
[In Optimal Design and Control (Blacksburg, VA, 1994), 1-22, Birkhäuser, Boston, Boston, MA, 1995.]
High Resolution Upwind-Mixed Finite Element Methods for Advection-Diffusion Equations with Variable Time-Stepping
Clint N. Dawson
Numerical methods for advection-diffusion equations are discussed based on approximating advection using a high-resolution upwind finite difference method, and incorporating diffusion using a mixed finite element method. In this approach, advection is approximated explicitly and diffusion implicitly. We first describe the basic procedure where each advection time-step is followed by a diffusion step. Because the explicit nature of the advective scheme requires a CFL time-step constraint, the basic procedure may be expensive, especially if the CFL constraint is severe. Two alternative time-stepping approaches are presented for improving computational efficiency while preserving accuracy. In the first approach, several advective time-steps are computed before taking a diffusion step. In the second approach, the advective time steps are also allowed to vary spatially. Numerical results for these three procedures for a model problem arising in flow through porous media are given.
June 1994
[Numer. Methods Partial Differential Equations 11 (1995), no. 5, 525-538.]
An Efficient Postprocessor for Velocities from Mixed Methods on Triangular Elements
Philip T. Keenan
Certain finite difference methods on rectangular grids for second order elliptic equations are known to yield superconvergent flux approximations. A class of related finite difference methods have recently been defined for triangular meshes by applying special quadrature rules to an extended version of a mixed finite element method [1]; the flux vector fields from these methods are not superconvergent. This report presents empirical evidence indicating that a simple local postprocessing technique recovers higher order accurate vector velocities at element centers on many meshes of triangular elements, with approximately second order accuracy on three lines meshes.
May 1994
Model Reduction of State Space Systems via an Implicitly Restarted Lanczos Method
E.J. Grimme, D.C. Sorensen, & P. Van Dooren
The nonsymmetric Lanczos method has recently received significant attention as a model reduction technique for large-scale systems. Unfortunately, the Lanczos method may produce an unstable partial realization for a given, stable system. To remedy this situation, inexpensive implicit restarts are developed which can be employed to stabilize the Lanczos generated model.
May 1994
[Numer. Algorithms 12 (1996), no. 1-2, 1-31.]
Computation of Eigenvalues for Starlike Domains
Robert A. Book
In this paper, we present a software tool for the computation of eigenvalues of starlike domains defined by polar boundary functions. We also offer and numerically test a conjecture on the monotonicity of the fundamental eigenvalues of the member of a family of starlike domains.
1994
Some Domain Decomposition and Multigrid Preconditioners for Hybrid Mixed Finite Elements
Lawrence C. Cowsar
Discretization of self-adjoint, linear, second-order, uniformly elliptic partial differential equations by hybrid mixed finite elements leads to large, ill-conditioned saddle-point problems. By eliminating the flux variable, a reduced problem is formed that is symmetric and positive definite but still large and ill-conditioned. Several domain decomposition and multigrid preconditioners are applied to the reduced problem, and bounds on their asymptotic rates of convergence are derived.
Two Schwartz domain decomposition methods are shown to converge at least as fast asymptotically as the same methods applied to conforming linear finite element discretizations. In particular, for both the standard additive overlapping Schwarz methods of Smith, it is proven that the rates of convergence of the methods are uniformly bounded with respect to the mesh size in both two and three dimensions under standard assumptions.
Several multigrid preconditioners are constructed for the reduced problem including a generalization of a method due to Bramble, Pasciak and Xu and an adaptation of methods of Wohlmuth and Hoppe. A common feature of these multigrid methods is the use of conforming finite element spaces on the coarser grids. Uniform convergence rates are proven for most of the methods and numerical results that verify the bounds are reported.
A mixed finite element discretization of a simplified model of sediment transport in a two dimensional periodic channel is also described. The results of two simulations that employ one of the multigrid preconditioners are reported.
April 1994
An Abstract Analysis of Differential Semblance Optimization
Mark S. Gockenbach
Differential Semblance Optimization (DSO) is a novel way of approaching a class of inverse problems arising in exploration seismology. The promising feature of the DSO method is that it replaces a nonsmooth, highly nonconvex cost function (the Output Least-Squares (OLS) objective function) with a smooth cost function that is amenable to standard (local) optimization algorithms.
The OLS problem can be written abstractly as a partially linear least-squares problem with linear constraints. The DSO objective function is derived from the associated quadratic penalty function. It is shown that one way to view the DSO objective function is as a regularization of a function that is dual (in a certain sense) to the OLS objective function. By viewing the DSO problem as a perturbation of this dual problem, this method can be shown to be effective. In particular, it is demonstrated that, under suitable assumptions, the DSO method defines a parameterized path of minimizers converging to the desired solution, and that for certain values of the parameter, standard optimization techniques can be used to find a point on the path.
The predictions of the theory are motivated and illustrated on two simple model problems for seismic velocity inversion, the plane wave detection problem and the "layer-over-half-space" problem. It is shown that the theory presented in this thesis extends the existing theory for the plane wave detection problem.
April 1994
A Nonlinear Mixed Finite Element Method for a Degenerate Parabolic Equation Arising in Flow in Porous Media
Todd Arbogast, Mary F. Wheeler, & Nai-Ying Zhang
We study a model nonlinear, degenerate, advection-diffusion equation having application in petroleum reservoir and groundwater aquifer simulation. The main difficulty is that the true solution is typically lacking in regularity; therefore, we consider the problem from the point of view of optimal approximation. Through time integration, we develop a mixed variational form that respects the known minimal regularity, and then we develop and analyze two versions of a mixed finite element approximation, a simpler semidiscrete (time continuous) version and a fully discrete version. Our error bounds are optimal in the sense that all but one of the bounding terms reduce to standard approximation error. The exceptional term is a a nonstandard approximation error term. We also consider our new formulation for the nondegenerate problem, showing the usual optimal L2-error bounds; moreover, superconvergence is obtained under special circumstances.
April 1994
[SIAM J. Numer. Anal. 33 (1996), no. 4, 1669-1687.]
Problem Formulations and Other Optimization Issues in Multidisciplinary Optimization
J.E. Dennis, Jr. & Robert Michael Lewis
This paper is about multidisciplinary (design) optimization, or MDO, the coupling of two or more analysis disciplines with numerical optimization. The "individual discipline feasible" (IDF) approaches introduced here make use of existing specialized analysis codes, and they introduce significant opportunities for coarse-grained computational parallelism particularly well-suited to heterogeneous computing environments.
April 1994
Gradient Calculation of the Traveltime Cost Function Without Ray-Tracing
Alain Sei & William W. Symes
The inverse problem of tomography is an iterative procedure. It requires the computation of the gradient of the traveltime misfit cost function many times. This calculation is customarily done by ray tracing, the path length of the rays being closely related to the gradient. We propose in this work an alternative method to compute the gradient of the traveltime cost function without ray tracing. We use upwind finite difference schemes to compute the traveltime field by solving the eikonal equation. Then by adjoint state techniques we derive a closed-form expression of the gradient of the traveltime cost function. This approach allows an accurate computation of the gradient as well as the freedom to change the norm on the model space.
April 1994
Modeling of a Constant Q: Methodology and Algorithm for an Efficient and Optimally Inexpensive Viscoelastic Technique
Joakim O. Blanch, Johan O.A. Robertsson, & William W. Symes
Linear anelastic phenomena in wave propagation problems can be well modeled through a viscoelastic mechanical model consisting of standard linear solids. In this paper we present a method for modeling of constant Q as a function of frequency based on an explicit closed formula for calculation of the parameter fields. The proposed method enables substantial savings in computations and memory requirements. Experiments show that the new method also yields higher accuracy in the modeling of Q than e.g., the Padé approximant method (Day and Minster, 1984).
March 1994
Deflation Techniques for an Implicitly Restarted Arnoldi Iteration
R.B. Lehoucq & Danny C. Sorensen
A deflation procedure is introduced that is designed to improve convergence of an implicitly restarted Arnoldi iteration for computing a few eigenvalues of a large matrix. As the iteration progresses the Ritz value approximations of the eigenvalues of A converge at different rates. A numerically stable deflation scheme is introduced that implicitly deflates the converged approximations from the iteration. We present two forms of implicit deflation. The first, a locking operation, decouples converged Ritz values and associated vectors from the active part of the iteration. The second, a purging operation, removes unwanted but converged Ritz pairs. Convergence of the iteration is improved and a reduction in computational effort is also achieved. The deflation strategies make it possible to compute multiple or clustered eigenvalues with a single vector restart method. A Block method is not required. These schemes are analyzed with respect to numerical stability and computational results are presented.
September 1994 (Revised February 1995)
[SIAM J. Matrix Anal. Appl. 17 (1996), no. 4, 789-821.]
A Polylogarithmic Bound for an Iterative Substructuring Method for Spectral Elements in Three Dimensions
L.F. Pavarino & O.B. Widlund
Iterative substructuring methods form an important family of domain decomposition algorithms for elliptic finite element problems. A p-version finite element method based on continuous, piecewise Qp functions is considered for second order elliptic problems in three dimensions; this special method can also be viewed as a conforming spectral element method. An iterative method is designed for which the condition number of the relevant operator grows only in proportion to (1 + log p)2. This bound is independent of jumps in the coefficient of the elliptic problem across the interfaces between the subregions. Numerical results are also reported which support the theory.
March 1994
[SIAM J. Numer. Anal. 33 (1996), no. 4, 1303-1335.]
Triangular Decomposition Methods for Solving Reducible Nonlinear Systems of Equations
J.E. Dennis Jr., José Mario Martínez, & Xiaodong Zhang
This paper generalizes to the nonlinear case a standard way to solve general sparse systems of linear equations. In particular, Duff [1977] has suggested that row and column interchanges be applied to permute the coefficient matrix of a linear system into block lower triangular form. The linear system then can be solved by using the associated block Gauss-Seidel or forward block substitution scheme. This is the approach taken in the Harwell MA28 routine. If more than one matrix with the same sparsity pattern is to be factored, then the reordering need not be redone. In extending this approach to nonlinear problems, it is necessary to assume as in the linear case that component equations can be evaluated separately from equations in other blocks. The algorithm for doing the reordering is very fast, and if the equations can be put into block lower triangular form with each block size much less than the dimension of the system, then a large savings is in order because only the diagonal blocks need to be factored. In the nonlinear variants considered here, there are additional advantages. Not only are only the diagonal blocks of the Jacobian matrix computed and factored, but the off-diagonal partial derivatives need not even exist. Numerical tests and analytic results affirm the intuition that these variants are superior locally to Newton's method. Current work is concerned with globalizing these methods and with variants especially suited to parallel implementations.
February 1994
[SIAM J. Optim. 4 (1994), no. 2, 358-382.]
Parallelizing Sylvester-like Operations on a Distributed Memory Computer
Danny Y. Hu & Danny C. Sorensen
A SIMD scheme for parallelization of the 2-D array operation M(x) = (D×A + B×I + V) x is developed for a distributed memory computer with 2-D mesh processor network. This scheme is carefully designed in order to approach the best possible speed-up by avoiding unnecessary message communication. The numerical results on a Touchstone DELTA machine are presented to demonstrate the effectiveness of the scheme. This scheme can be generalized to one that would deal with the 3-D array operation M3D(x) + (I×I×A + I×B×I + C×I×I + V) x on a distributed memory computer with 3-D torus processor network.
February 1994
The Computation of Elementary Unitary Matrices
Richard Lehoucq
The construction of elementary unitary matrices that transform a complex vector to a real multiple of e1, the first column of the identity matrix, is studied. We survey the two well known forms and present what appears to be relatively unknown third form implemented by LAPACK and subroutine CLARFG.
February 1994 (revised August 1995)
[ACM Trans. Math. Software 22 (1996), no. 4, 393-400.]
On the Characterization of Q-Superlinear Convergence of Quasi-Newton Interior-Point Methods for Nonlinear Programming
H.J. Martinez, Z. Parada, & R.A. Tapia
In this paper we extend the well-known Boggs-Tolle-Wang characterization of Q-superlinear convergence for quasi-Newton methods for equality constrained optimization to quasi-Newton interior-point methods for nonlinear programming. Critical issues in this extension include the choice of the centering parameter, the choice of the steplength parameter, and the determination of the primary variables.
February 1994 (Revised April 1995)
[Bol. Soc. Mat. Mexicana (3) 1 (1995), no. 2, 137--148.]
A Parallel Multiphase Numerical Model for Subsurface Contaminant Transport with Biodegradation Kinetics
Todd Arbogast, Clint N. Dawson, & Mary F. Wheeler
We discuss the formulation of a simulator in three spatial dimensions for two phase groundwater flow and transport with biodegradation kinetics that has been developed at Rice University for massively parallel, distributed memory, message passing machines. The numerical procedures employed are a fully implicit mixed finite element method for flow and a characteristics-mixed method for transport and reactions of dissolved chemical species in groundwater. Domain decomposition solvers have been employed for solving the systems of equations resulting from the discretization of the model. Results from applying this simulator to a bioremediation field problem using a recirculation well in an air-water system are discussed.
February 1994
GMRES on (Nearly) Singular Systems
Peter Brown & Homer Walker
We consider the behavior of the GMRES method for solving a linear system Ax=b when A is singular or nearly so, i.e., ill-conditioned. The (near) singularity of A may or may not affect the performance of GMRES, depending on the nature of the system and the initial approximate solution. For singular A, we give conditions under which the GMRES iterates converge safely to a least-squares solution or to the pseudo-inverse solution. These results also apply to any residual minimizing Krylov subspace method that is mathematically equivalent to GMRES. A practical procedure is outlined for efficiently and reliably detecting singularity or ill-conditioning when it becomes a threat to the performance of GMRES.
January 1994 (August 1995)
[SIAM J. Matrix Anal. Appl. 18 (1997), no. 1, 37-51.]
A Domain Decomposition Method for the Acoustic Wave Equation Allowing for Discontinuous Coefficients and Grid Change
Alain Bamberger, Roland Glowinski, & Quang Huy Tran
A domain decomposition technique is proposed for the computation of the acoustic wave equation, in which the bulk modulus and density fields are allowed to be discontinuous at the interfaces. Inside each subdomain, the method presented coincides with the second order finite difference schemes traditionally used in geophysical modelling. However, the possibility of assigning to each subdomain its own space-step makes numerical simulations much less expensive.
Another interest of the method lies in the fact that its hybrid variational formulation naturally leads to exact equations for gridpoints on the interfaces. Transposing Babuska-Brezzi's formalism on mixed and hybrid finite elements provides a suitable functional framework for this domain decomposition formulation and shows that the inf-sup condition remains the basic requirement for convergence to occur.
January 1994
[SIAM J. Numer. Anal. 34 (1997), no. 2, 603-639.]
Dispersion and Cost Analysis of Some Finite Difference Schemes in One-Parameter Acoustic Wave Modelling
William W. Symes & Quang Huy Tran
A systematic comparison is carried out between some standard finite difference schemes, regarding their costs and dispersion properties. To be more specific, given a precision threshold to be imposed on the velocity error and a finite difference scheme, it is possible to determine a time-step and a grid-spacing in an optimal manner, i.e. so as to minimize the computational cost. Using this optimal cost as a criterion, it becomes easy to single out the most economical scheme for the purpose of a synthetic seismic campaign.
This survey represents the preliminary part of the larger project Marmousie-3D, the purpose of which is to create a synthetic data base corresponding to one-parameter media. In this paper, one's attention is focused on the 2-2m family of schemes. Since both homogeneous and heterogeneous cases are investigated, the study is expected to provide realistic figures for future simulations.
January 1994
[Comput. Geosci. 1 (1997), no. 1, 1-33.]
Logically Rectangular Mixed Methods for Groundwater Flow and Transport on General Geometry
Todd Arbogast, Mary F. Wheeler, & Ivan Yotov
We consider an extended mixed finite element formulation for groundwater flow and transport problems with either a tensor hydraulic conductivity or a tensor dispersion. While the aquifer domain can be geometrically general, in our formulation the computational domain is rectangular. The approximating spaces for the mixed method are defined on a smooth curvilinear grid, obtained by a global mapping of the rectangular, computational grid. The original problem is mapped to the computational domain, giving a similar problem with a modified tensor coefficient. Special quadrature rules are introduced to transform the mixed method into a simple cell-centered finite difference method with a 9 point stencil in 2-D and 19 point stencil in 3-D. The resulting scheme is locally mass conservative. In the case of flow, linear Galerkin procedures give first order accurate velocities, while our method is second order accurate. Both computational and theoretical results are given.
January 1994
Mixed Finite Elements for Elliptic Problems with Tensor Coefficients as Finite Differences
Todd Arbogast, Mary F. Wheeler, & Ivan Yotov
We develop the theory of an expanded mixed finite element approximation of second order elliptic problems containing a tensor coefficient. The mixed method is expanded in the sense that three variables are explicitly approximated, namely, the scalar unknown, its gradient, and its flux (the tensor coefficient times the gradient). The expected optimal order approximations are obtained in the L² and H^{-s}-norms, and superconvergence is obtained between the L²-projection of the scalar variable and its approximation. The scheme is suitable for the case in which the coefficient is a tensor that may have zeros, since it does not need to be inverted.The resulting linear system is a saddle point problem. In the case of the lowest order Raviart-Thomas elements on rectangular parallelepipeds, we approximate this expanded mixed method by incorporating certain quadrature rules. This enables us to write the system as a simple, cell-centered finite difference method, requiring the solution of a sparse, positive semidefinite linear system for the scalar unknown. For a general tensor coefficient, the sparsity pattern for the scalar unknown is a nine point stencil in two dimensions, and 19 points in three dimensions. We show that the optimal rates of convergence are retained; moreover, superconvergence is obtained for the scalar unknown as well as for its gradient and flux at certain discrete points. Computational results illustrate these theoretical results.
January 1994
Two-Grid Methods for Mixed Finite Element Approximations of Nonlinear Parabolic Equations
Clint N. Dawson & Mary F. Wheeler
Mixed finite element approximation of nonlinear parabolic equations is discussed. The equation considered is a prototype of a model which arises in flow through porous media. A two-grid approximation scheme is developed and analyzed for implicit time discretizations. In this approach, the full nonlinear system is solved on a "coarse" grid of size H. The nonlinearities are expanded about the coarse grid solution, and the resulting linear but nonsymmetric system is solved on a "fine" grid of size h. Error estimates are derived which demonstrate that the error is O (h^{k+1} + H^{2(k+1)-d/2} + Deltat), where k is the degree of the approximating space for the primary variable and d is spatial dimension, with k >= 1 for d >= 2. For the RT0 space (k=0) on rectangular domains, we present a modified scheme for treating the coarse grid problem. Here we show that the error is O (h + H^{3 - d/2} + Deltat). The above estimates are useful for determining an appropriate H for the coarse grid problem.
January 1994
[In Domain Decomposition Methods in Scientific and Engineering Computing (University Park, PA, 1993), 191-203, Contemp. Math.,180, Amer. Math. Soc., Providence, RI, 1994.]