Past Computational Mathematics and Applications Seminar

Abstract 1 Another Orthogonal Matrix

A householder reflection and a suitable product of Givens rotations are two well known examples of an orthogonal matrix with given first column. We present another way to produce such a matrix and apply it to produce a "fast Givens" method to compute the R factor of A, A = QR. This approach avoids the danger of under/overflow.
(joint work with Eric Barszcz)

Abstract 2 An application of Pfaff's Theorem (on skew-symmetric matrices)

There are no constraints on the eigenvalues of a product of two real symmetric matrices but what about the product of two real skew-symmetric matrices?
(joint work with A Dubrulle)

When studying invariant quantities and stability of discretization schemes for time-dependent differential equations(ODEs), Backward error analysis (BEA) has proven itself an invaluable tool. Although the established results give very accurate estimates, the known results are generally given for "worst case" scenarios. By taking into account the structure of the differential equations themselves further improvements on the estimates can be established, and sharper estimates on invariant quantities and stability can be established. In the talk I will give an overview of BEA, and its applications as it stands emphasizing the shortcoming in the estimates. An alternative strategy is then proposed overcoming these shortcomings, resulting in a tool which when used in connection with results from dynamical systems theory gives a very good insight into the dynamics of discretized differential equations.

The convergence of iterative methods for solving the linear system Ax = b with a Hermitian positive definite matrix A depends on the condition number of A: the smaller the latter the faster the former. Hence the idea to multiply the equation by a matrix T such that the condition number of TA is much smaller than that of A. The above is a common interpretation of the technique known as preconditioning, the matrix T being referred to as the preconditioner for A.
The eigenvalue computation does not seem to benefit from the direct application of such a technique. Indeed, what is the point in replacing the standard eigenvalue problem Ax = λx with the generalized one TAx = λTx that does not appear to be any easier to solve? It is hardly surprising then that modern eigensolvers, such as ARPACK, do not use preconditioning directly. Instead, an option is provided to accelerate the convergence to the sought eigenpairs by applying spectral transformation, which generally requires the user to supply a subroutine that solves the system (A−σI)y = z, and it is entirely up to the user to employ preconditioning if they opt to solve this system iteratively.
In this talk we discuss some alternative views on the preconditioning technique that are more general and more useful in the convergence analysis of iterative methods and that show, in particular, that the direct preconditioning approach does make sense in eigenvalue computation. We review some iterative algorithms that can benefit from the direct preconditioning, present available convergence results and demonstrate both theoretically and numerically that the direct preconditioning approach has advantages over the two-level approach. Finally, we discuss the role that preconditioning can play in the a posteriori error analysis, present some a posteriori error estimates that use preconditioning and compare them with commonly used estimates in terms of the Euclidean norm of residual.

We consider the problem of computing ratings using the results of games (such as chess) played between a set of n players, and show how this problem can be reduced to computing the positive eigenvectors corresponding to the dominant eigenvalues of certain n by n matrices. There is a close connection with the stationary probability distributions of certain Markov chains. In practice, if n is large, then the matrices involved will be sparse, and the power method may be used to solve the eigenvalue problems efficiently.

Many industrial flow problems, expecially in the minerals and process

industries, are very complex, with strong interactions between phases

and components, and with very different length and time scales. This

presentation outlines the algorithms used in the CFX-5 software, and

describes the extension of its coupled solver approach to some

multi-scale industrial problems. including Population Balance modelling

to predict size distributions of a disperse phase. These results will be

illustrated on some practical industrial problems.

Spectral projections enjoy high order convergence for globally smooth functions. However, a single discontinuity introduces O(1) spurious oscillations near the discontinuity and reduces the high order convergence rate to first order, Gibbs' Phenomena. Although a direct expansion of the function in terms of its global moments yields this low order approximation, high resolution information is retained in the global moments. Two techniques for the resolution of the Gibbs' phenomenon are discussed, filtering and reprojection methods. An adaptive filter with optimal joint time-frequency localization is presented, which recovers a function from its N term Fourier projection within the error bound \exp(-Nd(x)), where d(x) is the distance from the point being recovered to the nearest discontinuity. Symmetric filtering, however, must sacrifice accuracy when approaching a discontinuity. To overcome this limitation, Gegenbauer postprocessing was introduced by Gottlieb, Shu, et al, which recovers a function from its N term Fourier projection within the error bound \exp(-N). An extension of Gegenbauer postprocessing with improved convergence and robustness properties is presented, the robust Gibbs complements. Filtering and reprojection methods will be put in a unifying framework, and their properties such as robustness and computational cost contrasted. This research was conducted jointly with Eitan Tadmor and Anne Gelb.

The use of weighted matchings is becoming increasingly standard in the

solution of sparse linear systems. While non-symmetric permutations based on these

matchings have been the state-of-the-art for

several years (especially for direct solvers), approaches for symmetric

matrices have only recently gained attention.

\\

In this talk we discuss results of our work on using weighted matchings in

the preconditioning of symmetric indefinite linear systems, following ideas

introduced by Duff and Gilbert. In order to maintain symmetry,

the weighted matching is symmetrized and the cycle structure of the

resulting matching is used to build reorderings that form small diagonal

blocks from the matched entries.

\\

For the preconditioning we investigated two approaches. One is an

incomplete $LDL^{T}$ preconditioning, that chooses 1x1 or 2x2 diagonal pivots

based on a simple tridiagonal pivoting criterion. The second approach

targets distributed computing, and is based on factorized sparse approximate

inverses, whose existence, in turn, is based on the existence of an $LDL^{T}$

factorization. Results for a number of comprehensive test sets are given,

including comparisons with sparse direct solvers and other preconditioning

approaches.

An interior-point method for solving mathematical programs with

equilibrium constraints (MPECs) is proposed. At each iteration of the

algorithm, a single primal-dual step is computed from each subproblem of

a sequence. Each subproblem is defined as a relaxation of the MPEC with

a nonempty strictly feasible region. In contrast to previous

approaches, the proposed relaxation scheme preserves the nonempty strict

feasibility of each subproblem even in the limit. Local and superlinear

convergence of the algorithm is proved even with a less restrictive

strict complementarity condition than the standard one. Moreover,

mechanisms for inducing global convergence in practice are proposed.

Numerical results on the MacMPEC test problem set demonstrate the

fast-local convergence properties of the algorithm.

The trapezoidal rule for numerical integration is remarkably accurate when

the integrand under consideration is smooth and periodic. In this

situation it is superior to more sophisticated methods like Simpson's rule

and even the Gauss-Legendre rule. In the first part of the talk we

discuss this phenomenon and give a few elementary examples. In the second

part of the talk we discuss the application of this idea to the numerical

evaluation of contour integrals in the complex plane.

Demonstrations involving numerical differentiation, the computation

of special functions, and the inversion of the Laplace transform will be

presented.

Plane Couette flow - the flow between two infinite parallel plates moving in opposite directions -

undergoes a discontinuous transition from laminar flow to turbulence as the Reynolds number is

increased. Due to its simplicity, this flow has long served as one of the canonical examples for understanding shear turbulence and the subcritical transition process typical of channel and pipe flows. Only recently was it discovered in very large aspect ratio experiments that this flow also exhibits remarkable pattern formation near transition. Steady, spatially periodic patterns of distinct regions of turbulent and laminar flow emerges spontaneously from uniform turbulence as the Reynolds number is decreased. The length scale of these patterns is more than an order of magnitude larger than the plate separation. It now appears that turbulent-laminar patterns are inevitable intermediate states on the route from turbulent to laminar flow in many shear flows. I will explain how we have overcome the difficulty of simulating these large scale patterns and show results from studies of three types of patterns: periodic, localized, and intermittent.

Sparse grids yield numerical solutions to PDEs with a

significantly reduced number of degrees of freedom. The relative

benefit increases with the dimensionality of the problem, which makes

multi-factor models for financial derivatives computationally tractable.

An outline of a convergence proof for the so called combination

technique will be given for a finite difference discretisation of the

heat equation, for which sharp error bounds can be shown.

Numerical examples demonstrate that by an adaptive (heuristic)

choice of the subspaces European and American options with up to thirty

(and most likely many more) independent variables can be priced with

high accuracy.

The computation of flows of compressible fluids will be

discussed, exploiting the symmetric form of the equations describing

compressible flow.

A cell is a wonderously complex object. In this talk I will

give an overview of some of the mathematical frameworks that are needed

in order to make progress to understanding the complex dynamics of a

cell. The talk will consist of a directed random walk through discrete

Markov processes, stochastic differential equations, anomalous diffusion

and fractional differential equations.

An essential first step in many problems of numerical analysis and

computer graphics is to cover a region with a reasonably regular mesh.

We describe a short MATLAB code that begins with a "distance function"

to describe the region: $d(x)$ is the distance to the boundary

(with d

The talk will start with an introduction to recent development in Mathematica, with emphasis on numerical computing. This will be followed by a discussion of graph drawing algorithms for the display of relational information, in particular force directed algorithms. The talk will show that by employing multilevel approach and octree data structure, it is possible to achieve fast display of very large relational information, without compromising the quality.

We consider the solution of a linear system of equations using the GMRES iterative method. In some applications, performing inexact matrix-vector products in this method may be interesting, provided that a reasonable convergence of GMRES is achieved. A GMRES algorithm where the matrix vector product is performed inexactly is termed ”inexact GMRES algorithm”. An application of this idea occurs in computational electromagnetics, where the fast multipole method provides approximations of the matrix-vector product within a user-defined precision, and where these inaccurate matrix-vector products are all the more cheaper (in terms of CPU time) as the user-defined precision is low. The key point is then to design a control strategy of the accuracy of the matrix-vector product so that the GMRES converges better in the sense that 1) the inexact method achieves a satisfactory limiting accuracy, 2) within a reasonable number of steps.

\\

In [1], a relaxation strategy is proposed for general systems and validated on a large set of numerical experiments. This work is based on heuristic considerations and proposes a strategy that enables a convergence of the GMRES iterates $x_{k}$ within a relative normwise backward error $\frac{\|b−Ax_{k}\|}{\|A\| \|x_{k}\| + \|b\|}$ less than a prescribed quantity $\eta$ > 0, on a significant number of numerical experiments. Similar strategies have been applied to the solution of device simulation problems using domain decomposition [2] and to the preconditioning of a radiation diffusion problem in [5].

\\

A step toward a theoretical explanation of the observed behaviour of the inexact GMRES is proposed in [3, 4]. In this talk, we show that in spite of this considerable theoretical study, the experimental work of [1] is not fully understood yet. We give an overview of the questions that still remains open both in exact arithmetic and in floating-point arithmetic, and we provide some insights into the solution of some of them. Possible applications of this work for the preconditioned GMRES method, when the matrix-vector product is accurate but the preconditioning operation is approximate, are also investigated, based on [3].

\\

\\

References

\\

[1] A. Bouras and V. Frayss´e. Inexact matrix-vector products in Krylov methods for solving linear systems: a relaxation strategy. SIAM Journal on Matrix Analysis and Applications, 2004. To appear.

\\

[2] A. Bouras, V. Frayss´e, and L. Giraud. A relaxation strategy for inner-outer linear solvers in domain decomposition methods. Technical Report TR/PA/00/17, CERFACS, Toulouse, France, 2000.

\\

[3] V. Simoncini and D. B. Szyld. Theory of inexact Krylov subspace methods and applications to scientific computing. SIAM Journal Scientific Computing, 25:454–477, 2003.

\\

[4] J. van den Eshof and G. L. G. Sleijpen. Inexact Krylov subspace methods for linear systems. SIAM Journal on Matrix Analysis and Applications, February 2004. To appear.

\\

[5] J. S. Warsa, M. Benzi, T. A. Warein, and J. E. Morel. Preconditioning a mixed discontinuous finite element method for radiation diffusion. Numerical Linear Algebra with Applications, 2004. To appear.

In recent years, there has been considerable interest, especially in the context of

fluid-dynamics, in nonconforming finite element methods that are based on discontinuous

piecewise polynomial approximation spaces; such approaches are referred to as discontinuous

Galerkin (DG) methods. The main advantages of these methods lie in their conservation properties, their ability to treat a wide range of problems within the same unified framework, and their great flexibility in the mesh-design. Indeed, DG methods can easily handle non-matching grids and non-uniform, even anisotropic, polynomial approximation degrees. Moreover, orthogonal bases can easily be constructed which lead to diagonal mass matrices; this is particularly advantageous in unsteady problems. Finally, in combination with block-type preconditioners, DG methods can easily be parallelized.

\\

\\

In this talk, we introduce DG discretizations of mixed field and potential-based formulations of

eddy current problems in the time-harmonic regime. For the electric field formulation, the

divergence-free constraint within non-conductive regions is imposed by means of a Lagrange

multiplier. This allows for the correct capturing of edge and corner singularities in polyhedral domains; in contrast, additional Sobolev regularity must be assumed in the DG formulation, and their conforming counterparts, when regularization techniques are employed. In particular, we present a mixed method involving discontinuous $P^\ell-P^\ell$ elements, which includes a normal jump stabilization term, and a non-stabilized variant employing discontinuous $P^\ell-P^{\ell+1}$ elements.The first formulation delivers optimal convergence rates for the vector-valued unknowns in a suitable energy norm, while the second (non-stabilized) formulation is designed to yield optimal convergence rates in both the $L^2$--norm, as well as in a suitable energy norm. For this latter method, we also develop the {\em a posteriori} error estimation of the mixed DG approximation of the Maxwell operator. Indeed, by employing suitable Helmholtz decompositions of the error, together with the conservation properties of the underlying method, computable upper bounds on the error, measured in terms of the energy norm, are derived.

\\

\\

Numerical examples illustrating the performance of the proposed methods will be presented; here,

both conforming and non-conforming (irregular) meshes will be employed. Our theoretical and

numerical results indicate that the proposed DG methods provide promising alternatives to standard conforming schemes based on edge finite elements.

The satisfiability (SAT) problem is a central problem in mathematical

logic, computing theory, and artificial intelligence. We consider

instances of SAT specified by a set of boolean variables and a

propositional formula in conjunctive normal form. Given such an instance,

the SAT problem asks whether there is a truth assignment to the variables

such that the formula is satisfied. It is well known that SAT is in

general NP-complete, although several important special cases can be

solved in polynomial time. Extending the work of de Klerk, Warners and van

Maaren, we present new linearly sized semidefinite programming (SDP)

relaxations arising from a recently introduced paradigm of higher

semidefinite liftings for discrete optimization problems. These

relaxations yield truth assignments satisfying the SAT instance if a

feasible matrix of sufficiently low rank is computed. The sufficient rank

values differ between relaxations and can be viewed as a measure of the

relative strength of each relaxation. The SDP relaxations also have the

ability to prove that a given SAT formula is unsatisfiable. Computational

results on hard instances from the SAT Competition 2003 show that the SDP

approach has the potential to complement existing techniques for SAT.