Past events

Non-selfadjoint singular differential equation eigenproblems arise in a number of contexts, including scattering theory, the study of quantum-mechanical resonances, and hydrodynamic and magnetohydrodynamic stability theory.

\\

\\

It is well known that the spectra of non-selfadjoint operators can be pathologically sensitive to perturbation of the

operator. Wilkinson provides matrix examples in his classical text, while Trefethen has studied the phenomenon extensively through pseudospectra, which he argues are often of more physical relevance than the spectrum itself. E.B. Davies has studied the phenomenon particularly in the context of Sturm-Liouville operators and has shown that the eigenfunctions and associated functions of non-selfadjoint singular Sturm-Liouville operators may not even form a complete set in $L^2$.

\\

\\

In this work we ask the question: under what conditions can one expect the regularization process used for selfadjoint singular Sturm-Liouville operators to be successful for non-selfadjoint operators? The answer turns out to depend in part on the so-called Sims Classification of the problem. For Sims Case I the process is not guaranteed to work, and indeed Davies has very recently described the way in which spurious eigenvalues may be generated and converge to certain curves in the complex plane.

\\

\\

Using the Titchmarsh-Weyl theory we develop a very simple numerical procedure which can be used a-posteriori to distinguish genuine eigenvalues from spurious ones. Numerical results indicate that it is able to detect not only the spurious eigenvalues due to the regularization process, but also spurious eigenvalues due to the numerics on an already-regular problem. We present applications to quantum mechanical resonance calculations and to the Orr-Sommerfeld equation.

\\

\\

This work, in collaboration with B.M. Brown in Cardiff, has recently been generalized to Hamiltonian systems.

Stiff systems of ODEs arise commonly when solving PDEs by spectral methods,

so conventional explicit time-stepping methods require very small time steps.

The stiffness arises predominantly through the linear terms, and these

terms can be handled implicitly or exactly, permitting larger time steps.

This work develops and investigates a class of methods known as

'exponential time differencing'. These methods are shown to have a

number of advantages over the more well-known linearly implicit

methods and integrating factor methods.

The Kestrel interface for submitting optimization problems to the NEOS Server augments the established e-mail, socket, and web interfaces by enabling easy usage of remote solvers from a local modeling environment.

\\

\\

Problem generation, including the run-time detection of syntax errors, occurs on the local machine using any available modeling language facilities. Finding a solution to the problem takes place on a remote machine, with the result returned in the native modeling language format for further processing. A byproduct of the Kestrel interface is the ability to solve multiple problems generated by a modeling language in parallel.

\\

\\

This mechanism is used, for example, in the GAMS/AMPL solver available through the NEOS Server, which internally translates a submitted GAMS problem into AMPL. The resulting AMPL problem is then solved through the NEOS Server via the Kestrel interface. An advantage of this design is that the GAMS to AMPL translator does not need to be collocated with the AMPL solver used, removing restrictions on solver choice and reducing administrative costs.

\\

\\

This talk is joint work with Elizabeth Dolan.

Tridiagonal matrices and three term recurrences and second order equations appear amazingly often, throughout all of mathematics. We won't try to review this subject. Instead we look in two less familiar directions.

\\

\\

Here is a tridiagonal matrix problem that waited surprisingly long for a solution. Forward elimination factors T into LDU, with the pivots in D as usual. Backward elimination, from row n to row 1, factors T into U_D_L_. Parlett asked for a proof that diag(D + D_) = diag(T) + diag(T^-1).^-1. In an excellent paper (Lin Alg Appl 1997) Dhillon and Parlett extended this four-diagonal identity to block tridiagonal matrices, and also applied it to their "Holy Grail" algorithm for the eigenproblem. I would like to make a different connection, to the Kalman filter.

\\

\\

The second topic is a generalization of tridiagonal to "tree-diagonal". Unlike the interval, the tree can branch. In the matrix T, each vertex is connected only to its neighbors (but a branch point has more than two neighbors). The continuous analogue is a second order differential equation on a tree. The "non-jump" conditions at a meeting of N edges are continuity of the potential (N-1 equations) and Kirchhoff's Current Law (1 equation). Several important properties of tridiagonal matrices, including O(N) algorithms, survive on trees.

No seminar this week

Let the thin plate spline radial basis function method be applied to

interpolate values of a smooth function $f(x)$, $x \!\in\! {\cal R}^d$.

It is known that, if the data are the values $f(jh)$, $j \in {\cal Z}^d$,

where $h$ is the spacing between data points and ${\cal Z}^d$ is the

set of points in $d$ dimensions with integer coordinates, then the

accuracy of the interpolant is of magnitude $h^{d+2}$. This beautiful

result, due to Buhmann, will be explained briefly. We will also survey

some recent findings of Bejancu on Lagrange functions in two dimensions

when interpolating at the integer points of the half-plane ${\cal Z}^2

\cap \{ x : x_2 \!\geq\! 0 \}$. Most of our attention, however, will

be given to the current research of the author on interpolation in one

dimension at the points $h {\cal Z} \cap [0,1]$, the purpose of the work

being to establish theoretically the apparent deterioration in accuracy

at the ends of the range from ${\cal O} ( h^3 )$ to ${\cal O} ( h^{3/2}

)$ that has been observed in practice. The analysis includes a study of

the Lagrange functions of the semi-infinite grid ${\cal Z} \cap \{ x :

x \!\geq\! 0 \}$ in one dimension.

In this talk we review experiences of using the Harwell Subroutine

Library and other numerical software codes in implementing large scale

solvers for commercial industrial process simulation packages. Such

packages are required to solve problems in an efficient and robust

manner. A core requirement is the solution of sparse systems of linear

equations; various HSL routines have been used and are compared.

Additionally, the requirement for fast small dense matrix solvers is

examined.

The sphere is an important setting for applied mathematics, yet the underlying approximation theory and numerical analysis needed for serious applications (such as, for example, global weather models) is much less developed than, for example, for the cube.

\\

\\

This lecture will apply recent developments in approximation theory on the sphere to two different problems in scientific computing.

\\

\\

First, in geodesy there is often the need to evaluate integrals using data selected from the vast amount collected by orbiting satellites. Sometimes the need is for quadrature rules that integrate exactly all spherical polynomials up to a specified degree $n$ (or equivalently, that integrate exactly all spherical harmonies $Y_{\ell ,k}(\theta ,\phi)$ with $\ell \le n).$ We shall demonstrate (using results of M. Reimer, I. Sloan and R. Womersley in collaboration with

W. Freeden) that excellent quadrature rules of this kind can be obtained from recent results on polynomial interpolation on the sphere, if the interpolation points (and thus the quadrature points) are chosen to be points of a so-called extremal fundamental system.

\\

\\

The second application is to the scattering of sound by smooth three-dimensional objects, and to the inverse problem of finding the shape of a scattering object by observing the pattern of the scattered sound waves. For these problems a methodology has been developed, in joint work with I.G. Graham, M. Ganesh and R. Womersley, by applying recent results on constructive polynomial approximation on the sphere. (The scattering object is treated as a deformed sphere.)

Continuing advances in computing technology provide the power not only to solve

increasingly large and complex process modeling and optimization problems, but also

to address issues concerning the reliability with which such problems can be solved.

For example, in solving process optimization problems, a persistent issue

concerning reliability is whether or not a global, as opposed to local,

optimum has been achieved. In modeling problems, especially with the

use of complex nonlinear models, the issue of whether a solution is unique

is of concern, and if no solution is found numerically, of whether there

actually exists a solution to the posed problem. This presentation

focuses on an approach, based on interval mathematics,

that is capable of dealing with these issues, and which

can provide mathematical and computational guarantees of reliability.

That is, the technique is guaranteed to find all solutions to nonlinear

equation solving problems and to find the global optimum in nonlinear

optimization problems. The methodology is demonstrated using several

examples, drawn primarily from the modeling of phase behavior, the

estimation of parameters in models, and the modeling, using lattice

density-functional theory, of phase transitions in nanoporous materials.

This talk reviews some recent joint work with Michael Eiermann and Olaf

Schneider which introduced a framework for analyzing some popular

techniques for accelerating restarted Krylov subspace methods for

solving linear systems of equations. Such techniques attempt to compensate

for the loss of information due to restarting methods like GMRES, the

memory demands of which are usually too high for it to be applied to

large problems in unmodified form. We summarize the basic strategies which

have been proposed and present both theoretical and numerical comparisons.

Support Vector Machines are a new and very promising approach to

machine learning. They can be applied to a wide range of tasks such as

classification, regression, novelty detection, density estimation,

etc. The approach is motivated by statistical learning theory and the

algorithms have performed well in practice on important applications

such as handwritten character recognition (where they currently give

state-of-the-art performance), bioinformatics and machine vision. The

learning task typically involves optimisation theory (linear, quadratic

and general nonlinear programming, depending on the algorithm used).

In fact, the approach has stimulated new questions in optimisation

theory, principally concerned with the issue of how to handle problems

with a large numbers of variables. In the first part of the talk I will

overview this subject, in the second part I will describe some of the

speaker's contributions to this subject (principally, novelty

detection, query learning and new algorithms) and in the third part I

will outline future directions and new questions stimulated by this

research.

The flow in a cylinder with a rotating endwall has continued to

attract much attention since Vogel (1968) first observed the vortex

breakdown of the central core vortex that forms. Recent experiments

have observed a multiplicity of unsteady states that coexist over a

range of the governing parameters. In spite of numerous numerical and

experimental studies, there continues to be considerable controversy

with fundamental aspects of this flow, particularly with regards to

symmetry breaking. Also, it is not well understood where these

oscillatory states originate from, how they are interrelated, nor how

they are related to the steady, axisymmetric basic state.

\\

\\

In the aspect ratio (height/radius) range 1.6 2.8. An efficient and

accurate numerical scheme is presented for the three-dimensional

Navier-Stokes equations in primitive variables in a cylinder. Using

these code, primary and secondary bifurcations breaking the SO(2)

symmetry are analyzed.

\\

\\

We have located a double Hopf bifurcation, where an axisymmetric limit

cycle and a rotating wave bifurcate simultaneously. This codimension-2

bifurcation is very rich, allowing for several different scenarios. By

a comprehensive two-parameter exploration about this point we have

identified precisely to which scenario this case corresponds. The mode

interaction generates an unstable two-torus modulate rotating wave

solution and gives a wedge-shaped region in parameter space where the

two periodic solutions are both stable.

\\

\\

For aspect ratios around three, experimental observations suggest that

the first mode of instability is a precession of the central vortex

core, whereas recent linear stability analysis suggest a Hopf

bifurcation to a rotating wave at lower rotation rates. This apparent

discrepancy is resolved with the aid of the 3D Navier-Stokes

solver. The primary bifurcation to an m=4 traveling wave, detected by

the linear stability analysis, is located away from the axis, and a

secondary bifurcation to a modulated rotating wave with dominant modes

m=1 and 4, is seen mainly on the axis as a precessing vortex breakdown

bubble. Experiments and the linear stability analysis detected

different aspects of the same flow, that take place in different

spatial locations.

We combine linear algebra techniques with finite element techniques to obtain a reliable stopping criterion for the Conjugate Gradient algorithm. The finite element method approximates the weak form of an elliptic partial differential equation defined within a Hilbert space by a linear system of equations A x = b, where A is a real N by N symmetric and positive definite matrix. The conjugate gradient method is a very effective iterative algorithm for solving such systems. Nevertheless, our experiments provide very good evidence that the usual stopping criterion based on the Euclidean norm of the residual b - Ax can be totally unsatisfactory and frequently misleading. Owing to the close relationship between the conjugate gradient behaviour and the variational properties of finite element methods, we shall first summarize the principal properties of the latter. Then, we will use the recent results of [1,2,3,4]. In particular, using the conjugate gradient, we will compute the information which is necessary to evaluate the energy norm of the difference between the solution of the continuous problem, and the approximate solution obtained when we stop the iterations by our criterion.

Finally, we will present the numerical experiments we performed on a selected ill-conditioned problem.

References

• [1] M. Arioli, E. Noulard, and A. Russo, Vector Stopping Criteria for Iterative Methods: Applications to PDE's, IAN Tech. Rep. N.967, 1995.
• [2] G.H. Golub and G. Meurant, Matrices, moments and quadrature II; how to compute the norm of the error in iterative methods, BIT., 37 (1997), pp.687-705.
• [3] G.H. Golub and Z. Strakos, Estimates in quadratic formulas, Numerical Algorithms, 8, (1994), pp.~241--268.
• [4] G. Meurant, The computation of bounds for the norm of the error in the conjugate gradient algorithm, Numerical Algorithms, 16, (1997), pp.~77--87.

Boundary layers are often studied with no pressure gradient

or with an imposed pressure gradient. Either of these assumptions

can lead to difficulty in obtaining solutions. A major advance in fluid

dynamics last century (1969) was the development of a triple deck

formulation for boundary layers where the pressure is not

specified but emerges through an interaction between

boundary layer and the inviscid outer flow. This has given rise to

new computational problems and computations have in turn

fed ideas back into theoretical developments. In this survey talk

based on my new book, I will look at three problems:

flow past a plate, flow separation and flow in channels

and discuss the interaction between theory and computation

in advancing boundary layer theory.