Comlab

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)

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.

This is a special talk outside the normal Computational Mathematics and Application seminar series. Please note it takes place on a Wednesday.

Current community models in the geosciences employ a variety of numerical methods from finite-difference, finite-volume, finite- or spectral elements, to pseudospectral methods. All have specialized strengths but also serious weaknesses. The first three methods are generally considered low-order and can involve high algorithmic complexity (as in triangular elements or unstructured meshes). Global spectral methods do not practically allow for local mesh refinement and often involve cumbersome algebra. Radial basis functions have the advantage of being spectrally accurate for irregular node layouts in multi-dimensions with extreme algorithmic simplicity, and naturally permit local node refinement on arbitrary domains. We will show test examples ranging from vortex roll-ups, modeling idealized cyclogenesis, to the unsteady nonlinear flows posed by the shallow water equations to 3-D mantle convection in the earth’s interior. The results will be evaluated based on numerical accuracy, stability and computational performance.

An iterative substructuring method with balancing Neumann-Neumann preconditioner is known as an efficient parallel algorithm for the elasticity equations. This method was extended to the Stokes equations by Pavarino and Widlund [2002]. In their extension, Q2/P0-discontinuous elements are used for velocity/pressure and a Schur complement system within "benign space", where incompressibility satisfied, is solved by CG method.

For the construction of the coarse space for the balancing preconditioner, some supplementary solvability conditions are considered. In our algorithm for 3-D computation, P1/P1 elements for velocity/pressure with pressure stabilization are used to save computational cost in the stiffness matrix. We introduce a simple coarse space similar to the one of elasticity equations. Owing to the stability term, solvabilities of local Dirichlet problem for a Schur complement system, of Neumann problem for the preconditioner, and of the coarse space problem are ensured. In our implementation, local Dirichlet and Neumann problems are solved by a 4x4-block modified Cholesky factorization procedure with an envelope method, which leads to fast computation with small memory requirement. Numerical result on parallel efficiency with a shared memory computer will be shown. Direct use of the Stokes solver in an application of Earth's mantle convection problem will be also shown.

For the task of solving PDEs, finite difference (FD) methods are particularly easy to implement. Finite element (FE) methods are more flexible geometrically, but tend to be difficult to make very accurate. Pseudospectral (PS) methods can be seen as a limit of FD methods if one keeps on increasing their order of accuracy. They are extremely effective in many situations, but this strength comes at the price of very severe geometric restrictions. A more standard introduction to PS methods (rather than via FD methods of increasing orders of accuracy) is in terms of expansions in orthogonal functions (such as Fourier, Chebyshev, etc.).

Radial basis functions (RBFs) were first proposed around 1970 as a tool for interpolating scattered data. Since then, both our knowledge about them and their range of applications have grown tremendously. In the context of solving PDEs, we can see the RBF approach as a major generalization of PS methods, abandoning the orthogonality of the basis functions and in return obtaining much improved simplicity and flexibility. Spectral accuracy becomes now easily available also when using completely unstructured meshes, permitting local node refinements in critical areas. A very counterintuitive parameter range (making all the RBFs very flat) turns out to be of special interest. Computational cost and numerical stability were initially seen as serious difficulties, but major progress have recently been made also in these areas.

Quasicontinuum methods are a prototypical class of atomistic-to-continuum coupling methods. For example, we may wish to model a lattice defect (a vacancy or a dislocation) by an atomistic model, but the elastic far field by a continuum model. If the continuum model is consistent with the atomistic model (e.g., the Cauchy--Born model) then the main question is how the interface treatment affects the method.

In this talk I will introduce three of the main ideas how to treat the interface. I will explain their strengths and weaknesses by formulating the simplest possible non-trivial model problem and then simply analyzing the two classical concerns of numerical analysis: consistency and stability.