%%% ====================================================================
%%% BibTeX-file{
%%%     author          = "Nicholas J. Higham",
%%%     version         = "1.00",
%%%     date            = "October 24, 2005",
%%%     time            = "",
%%%     filename        = "u-manchester-mccm.bib",
%%%     address         = "Department of Mathematics
%%%                        University of Manchester
%%%                        Manchester M13 9PL
%%%                        England",
%%%     telephone       = "+44 (0)161 275 5800",
%%%     FAX             = "+44 (0)161 275 5819",
%%%     checksum        = ""
%%%     email           = "higham at ma.man.ac.uk (Internet)",
%%%     codetable       = "ISO/ASCII",
%%%     keywords        = "bibliography, Manchester, numerical
%%%                        analysis",
%%%     supported       = "yes",
%%%     docstring       = "BibTeX bibliography for the Numerical
%%%                        Analysis Reports of the Manchester Centre for
%%%                        Computational Mathematics
%%%                        (University of Manchester and UMIST)."
%%%     }
%%% ====================================================================

%%% The URL (Uniform Resource Locator) of this file is
%%%
%%%   URL = "http://www.maths.man.ac.uk/~nareports/u-manchester-mccm.bib"
%%%
%%% This is narep.bib, the index for the Numerical Analysis Reports of the
%%% Manchester Centre for Computational Mathematics
%%% (University of Manchester and UMIST).
%%% It is in BibTeX form, with additional fields.
%%% Entries are in reverse chronological order.

@String{inst-MCCM               = "Manchester Centre for
Computational Mathematics"}
%%% Numerical Analysis Reports: ISSN 1360-1725

@String{t-NAREP                 = "Numerical Analysis Report"}

@String{inst-U-MANCHESTER       = "University of Manchester"}
@String{inst-U-MANCHESTER:adr   = "Manchester M13 9PL, England"}

@String{pub-LONGMAN             = "Longman Scientific and Technical"}

@String{pub-OUP                 = "Oxford University Press"}

@String{pub-SV                  = "Spring{\-}er-Ver{\-}lag"}
Germany~/ London, UK~/ etc."}
%%% ====================================================================

@techreport{Guo:2005:SNM,
author = "Chun-Hua Guo and Nicholas J. Higham",
title = "A {Schur--Newton} Method for the Matrix $p$th Root and its Inverse",
institution =  inst-MCCM,
type = t-NAREP,
number = "No. 475",
pages = "17",
month = oct,
year = "2005",
URL =          "http://www.maths.man.ac.uk/~nareports/narep475.pdf",
abstract = "Newton's method for the inverse matrix $p$th root, $A^{-1/p}$, has
the attraction that it involves only matrix multiplication. We
show that if the starting matrix is $c^{-1}I$ for $c\in\R^+$ then
the iteration converges quadratically to $A^{-1/p}$ if the
eigenvalues of $A$ lie in a wedge-shaped convex set containing the
disc $\{\, z: |z-c^p| < c^p\,\}$. We derive an optimal choice of
$c$ for the case where $A$ has real, positive eigenvalues. An
application is described to roots of transition matrices from
Markov models, in which for certain problems the convergence
condition is satisfied with $c=1$. Although the basic Newton
iteration is numerically unstable, a coupled version is stable and
a simple modification of it provides a new coupled iteration for
the matrix $p$th root. For general matrices we develop a hybrid
algorithm that computes a Schur decomposition, takes square roots
of the upper (quasi)triangular factor, and applies the coupled
Newton iteration to a matrix for which fast convergence is
guaranteed. The new algorithm can be used to compute either
$A^{1/p}$ or $A^{-1/p}$, and for large $p$ that are not highly
composite it is more efficient than the method of Smith based
entirely on the Schur decomposition.",
keywords = "Matrix $p$th root, principal $p$th root, matrix logarithm, inverse,
Newton's method, preprocessing, Schur decomposition, numerical
stability, convergence, Markov model, transition matrix"
}

@techreport{Elman:05:IFI,
author = "Howard Elman and Alison Ramage and David Silvester",
title = "IFISS: a {Matlab} Toolbox for Modelling Incompressible Flow",
institution =  inst-MCCM,
type = t-NAREP,
number = "No. 474",
month = oct,
year = 2005,
note = "Submitted to ACM Transactions on Mathematical Software",
URL =          "http://www.maths.man.ac.uk/~nareports/narep474.pdf",
abstract = "IFISS is a graphical Matlab package for the interactive numerical study
of incompressible flow problems. It includes algorithms for discretisation by
mixed finite element methods and a posteriori error estimation of the
computed solutions. The package can also be used as
a computational laboratory for experimenting with state-of-the-art
preconditioned iterative solvers for the discrete linear equation systems that
arise in incompressible flow modelling. A unique feature of the package is its
comprehensive nature; for each
problem addressed, it enables the study of both discretisation and iterative
solution algorithms as well as the interaction between the two and the
resulting effect on overall efficiency."
}

@techreport{Tisseur:05:SCN,
author = "Fran\c{c}oise Tisseur and Stef Graillat",
title = "Structured Condition Numbers and Backward Errors
in Scalar Product Spaces",
institution =  inst-MCCM,
type = t-NAREP,
number = "No. 473",
pages = "21",
month = sep,
year = "2005",
URL =          "http://www.maths.man.ac.uk/~nareports/narep473.pdf",
abstract = "We investigate the effect of structure-preserving perturbations on the
solution to a linear system, matrix inversion, and distance to
singularity. Particular attention is paid to linear and nonlinear
structures that form Lie algebras, Jordan algebras and automorphism
groups of a scalar product. These include complex symmetric,
pseudo-symmetric, persymmetric, skew-symmetric, Hamiltonian, unitary,
complex orthogonal and symplectic matrices. We show that under  mild
assumptions on the scalar product, there is little or no difference
between  structured and unstructured condition numbers and distance to
singularity for matrices in Lie and Jordan algebras. Hence, for these
classes of matrices, the usual unstructured perturbation analysis is
sufficient. We show this is not true in general for structures in
automorphism groups. Bounds and computable expressions for the
structured condition numbers for a linear system and matrix inversion
are derived for these nonlinear structures.
Structured backward errors for the approximate solution of linear
systems are also considered. Conditions are given for the structured
backward error to be finite. We prove that for Lie and Jordan
algebras, whenever the structured backward error is finite, it is
within a small factor of the unstructured one. The same conclusion
holds for orthogonal and unitary structures but cannot easily be
extended to other matrix groups. For some nonlinear structures such as
complex orthogonal, we show that the structured backward error can be
rewritten in terms of a challenging Procrustes problem.
This work extends and unifies earlier analyses.",
keywords =  "structured matrices, normwise structured perturbations, structured
linear systems, condition number, backward error, distance to
singularity, Lie algebra, Jordan algebra, automorphism group, scalar
product, bilinear form, sesquilinear form, orthosymmetric."
}

@techreport{MCCM:2005:AR4,
author =       "Manchester Centre for Computational Mathematics",
title =        "Annual Report: {January--December} 2004",
type =         t-NAREP,
number =       "No. 472",
institution =  inst-MCCM,
pages =        "25",
month =        sep,
year =         "2005",
URL =          "http://www.maths.man.ac.uk/~nareports/narep472.pdf"
}

@techreport{Bhatia:2005:IDM,
author = "Rajendra Bhatia",
title = "Infinitely Divisible Matrices",
institution =  inst-MCCM,
type = t-NAREP,
number = "No. 471",
pages = "36",
month = aug,
year = "2005",
URL =          "http://www.maths.man.ac.uk/~nareports/narep471.pdf"
}

@techreport{Lancaster:2005:NWS,
author = "Peter Lancaster and Panayiotis Psarrakos",
title = "A Note on Weak and Strong Linearizations of
Regular Matrix Polynomials",
institution =  inst-MCCM,
type = t-NAREP,
number = "No. 470",
pages = "6",
month = jun,
year = "2005",
URL =          "http://www.maths.man.ac.uk/~nareports/narep470.pdf"
}

@techreport{Berhanu:05:PBH,
author = "Michael Berhanu",
title = "Perturbation Bounds for Hyperbolic Matrix Factorizations",
institution =  inst-MCCM,
type = t-NAREP,
number = "No. 469",
pages = "24",
month = jun,
year = "2005",
URL =          "http://www.maths.man.ac.uk/~nareports/narep469.pdf",
abstract = "Several matrix factorizations depend on orthogonal factors,
matrices that preserve the Euclidean scalar product. Some of these
factorizations can be extended and generalized to
$(J,\Jt)$-orthogonal factors, that is, matrices that satisfy
$H^TJH=\Jt$, where $J$ and $\Jt$ are diagonal with diagonal
elements $\pm1$. The purpose of this work is to analyze the
perturbation of matrix factorizations that have a
$(J,\Jt)$-orthogonal or orthogonal factor and to give first order
perturbation bounds. For each factorization analyzed, we give the
sharpest possible first order bound,  which yields a condition
number. The cost of computing these condition numbers is high. It
is usually equivalent to computing the $2$-norm of an $n^2\times n^2$ matrix for a problem of size $n$. Thus, we also propose less
sharp bounds that are less expensive to compute.",
keywords = "Condition number, hyperbolic matrix factorization, manifold"
}

@techreport{Higham:05:IPL,
author = "Nicholas J. Higham",
title = "An Interview with {Peter Lancaster}",
institution =  inst-MCCM,
type = t-NAREP,
number = "No. 468",
pages = "10",
month = jun,
year = "2005",
URL =          "http://www.maths.man.ac.uk/~nareports/narep468.pdf",
}

@techreport{Tisseur:05:SEC,
author = "Michael Karow and Daniel Kressner and Fran\c{c}oise Tisseur",
title = "Structured Eigenvalue Condition Numbers",
institution =  inst-MCCM,
type = t-NAREP,
number = "No. 467",
pages = "16",
month = apr,
year = "2005",
URL =          "http://www.maths.man.ac.uk/~nareports/narep467.pdf",
abstract = "This paper investigates the effect of structure-preserving
perturbations on the eigenvalues of linearly and nonlinearly
structured eigenvalue problems. Particular attention is paid to
structures that form Jordan algebras, Lie algebras, and automorphism
groups of a scalar product. Bounds and computable expressions for
structured eigenvalue condition numbers are derived for these classes
of matrices, which include complex symmetric, pseudo symmetric,
persymmetric, skew-symmetric, Hamiltonian, symplectic, and orthogonal
matrices. In particular we show that under mild assumptions on the
scalar product, the structured and unstructured eigenvalue condition
numbers are equal for structures in Jordan algebras. For Lie
algebras, the effect on the condition number of incorporating
structure varies greatly with the structure. We identify  Lie
algebras for which structure does not affect the eigenvalue condition
number.",
keywords = "Structured eigenvalue problem, condition number, Jordan algebra,
Lie algebra, automorphism group, symplectic, perplectic,
pseudo-orthogonal, pseudo-unitary, complex symmetric, persymmetric,
perskew-symmetric, Hamiltonian, skew-Hamiltonian, structure
preservation."
}

@techreport{Mackey:2005:PPE,
author = "D. Steven Mackey and Niloufer Mackey and Christian Mehl and
Volker Mehrmann",
title = "Palindromic Polynomial Eigenvalue Problems: {Good} Vibrations from
Good Linearizations",
institution =  inst-MCCM,
type =         t-NAREP,
number = "No. 466",
pages = 26,
month = apr,
year = "2005",
URL =          "http://www.maths.man.ac.uk/~nareports/narep466.pdf",
abstract = "Palindromic polynomial eigenvalue problems and related classes of
structured eigenvalue problems are considered. These structures
generalize the concepts of symplectic and Hamiltonian matrices to
matrix polynomials. We discuss several applications where these
matrix polynomials arise, and show how linearizations can be
derived that re ect the structure of all these structured matrix
polynomials and therefore preserve symmetries in the spectrum.",
keywords = "nonlinear eigenvalue problem, palindromic matrix polynomial, even
matrix polynomial, odd matrix polynomial, Cayley transformation,
linearization, Hamiltonian matrix, symplectic matrix"
}

@techreport{Higham:2005:CLM,
author = "Nicholas J. Higham and D. Steven Mackey and Fran\c{c}oise Tisseur",
title = "The Conditioning of Linearizations of Matrix Polynomials",
institution =  inst-MCCM,
type = t-NAREP,
number = "No. 465",
pages = "25",
month = apr,
year = "2005",
URL =          "http://www.maths.man.ac.uk/~nareports/narep465.pdf",
abstract = "The standard way of solving the polynomial eigenvalue problem of
degree $m$ in $n\times n$ matrices is to linearize'' to a pencil
in $mn\times mn$ matrices and solve the generalized eigenvalue
problem. For a given polynomial, $P$, infinitely many
linearizations exist and they can have widely varying eigenvalue
condition numbers. We investigate the conditioning of
linearizations from a vector space $\mathbb{DL}(P)$ of pencils
recently identified and studied by Mackey, Mackey, Mehl, and
Mehrmann. We look for the best conditioned linearization and
compare the conditioning with that of the original polynomial. Two
particular pencils are shown always to be almost optimal over
linearizations in $\mathbb{DL}(P)$ for eigenvalues of modulus
greater than or less than 1, respectively, provided that the
problem is not too badly scaled and that the pencils are
linearizations. Moreover, under this scaling assumption, these
pencils are shown to be about as well conditioned as the original
polynomial. For quadratic eigenvalue problems that are not too
heavily damped, a simple scaling is shown to convert the problem
to one that is well scaled. We also analyze the eigenvalue
conditioning of the widely used first and second companion
linearizations. The conditioning of the first companion
linearization relative to that of $P$ is shown to depend on the
coefficient matrix norms, the eigenvalue, and the left \ev s of
the linearization and of $P$. The companion form is found to be
potentially much more ill conditioned than $P$, but if the 2-norms
of the coefficient matrices are all approximately 1 then the
companion form and $P$ are guaranteed to have similar condition
numbers. Analogous results hold for the second companion form. Our
results are phrased in terms of both the standard relative
condition number and the condition number of Dedieu and Tisseur
for the problem in homogeneous form, this latter condition number
having the advantage of applying to zero and infinite
eigenvalues.",
keywords = "matrix polynomial, matrix pencil, linearization, companion form,
condition number, homogeneous form, quadratic eigenvalue problem,
vector space, scaling"
}

@techreport{Mackey:2005:VSL,
author = "D. Steven Mackey and Niloufer Mackey and Christian Mehl and
Volker Mehrmann",
title = "Vector Spaces of Linearizations for Matrix Polynomials",
institution =  inst-MCCM,
type =         t-NAREP,
number = "No. 464",
pages = "32",
month = apr,
year = "2005",
URL =          "http://www.maths.man.ac.uk/~nareports/narep464.pdf",
abstract = "The classical approach to investigating polynomial eigenvalue
problems is linearization, where the polynomial is converted into
a larger matrix pencil with the same eigenvalues. For any
polynomial there are innitely many linearizations with widely
varying properties, but in practice the companion forms are
typically used. However, these companion forms are not always
entirely satisfactory, and linearizations with special properties
may sometimes be required. In this paper we develop a systematic
approach to generating large classes of linearizations for matrix
polynomials. Given a polynomial P, we show how to simply construct
two vector spaces of pencils that generalize the companion forms
of P, and prove that almost all of these pencils are
linearizations for P. Eigenvectors of these pencils are shown to
be closely related to those of P. A distinguished subspace is then
isolated, and the special properties of these pencils are
investigated. These spaces of pencils provide a convenient arena
in which to look for structured linearizations of structured
polynomials, as well as to try to optimize the conditioning of
linearizations, issues to be addressed in further work.",
keywords = "matrix polynomial, matrix pencil, linearization, strong
linearization, shifted sum, companion form"
}

@techreport{Lancaster:2005:CFS,
author = "Peter Lancaster and Leiba Rodman",
title = "Canonical Forms for Symmetric/Skew-symmetric Real Matrix Pairs
under Strict Equivalence and Congruence",
institution =  inst-MCCM,
type = t-NAREP,
number = "No. 463",
pages = "77",
month = mar,
year = "2005",
URL =          "http://www.maths.man.ac.uk/~nareports/narep463.pdf",
abstract = "A systematic development is made of the simultaneous reduction of
pairs of quadratic forms over the reals, one of which is
skew-symmetric and the other is either symmetric or
skew-symmetric. These reductions are by strict equivalence and by
congruence, over the reals or over the complex numbers, and
essentially complete proofs are presented.  Some closely related
results which can be derived from the canonical forms of pairs of
symmetric/skew-symmetric real forms are also included. They
concern simultaneously neutral subspaces, Hamiltonian and
skew-Hamiltonian matrices, and canonical structures of real
skew-symmetric indefinite inner product, and real matrices which
are skew-adjoint in a regular symmetric indefinite inner product.
The paper is largely expository, and continues the comprehensive
account of the reduction of pairs of matrices started in a
preceding paper.",
keywords = ""
}

@techreport{Lancaster:2005:ISP,
author = "Peter Lancaster",
title = "Inverse Spectral Problems for Semi-simple Damped Vibrating Systems",
institution =  inst-MCCM,
type = t-NAREP,
number = "No. 462",
pages = "24",
month = mar,
year = "2005",
URL =          "http://www.maths.man.ac.uk/~nareports/narep462.pdf",
abstract = "Computational schemes are investigated for the solution of inverse
spectral problems for $n\times n$ real systems of the form
$L(\la)=M\la^2 +D\la +K$. Thus, admissible sets of data concerning
systems of eigenvalues and eigenvectors are examined and
procedures for generating associated (isospectral) families of
systems are developed. The analysis includes symmetric systems,
systems with mixed real/non-real spectrum, systems with  positive
definite coefficients, and hyperbolic systems (with real
spectrum). A one-to-one correspondence between Jordan pairs and
structure preserving similarities is clarified. An examination of
complex symmetric matrices is included.",
keywords = ""
}

@techreport{Hargreaves:2005:EAM,
author = "Gareth I. Hargreaves and Nicholas J. Higham",
title = "Efficient Algorithms for the Matrix Cosine and Sine",
institution =  inst-MCCM,
type = t-NAREP,
number = "461",
pages = "17",
month = feb,
year = "2005",
URL =          "http://www.maths.man.ac.uk/~nareports/narep461.pdf",
abstract = "Several improvements are made to an algorithm of Higham and Smith
for computing the matrix cosine. The original algorithm scales the
matrix by a power of 2 to bring the $\infty$-norm to 1 or less,
evaluates the [8/8] Pad\'e approximant, then uses the double-angle
formula $\cos(2A) = 2\cos^2A - I$ to recover the cosine of the
original matrix. The first improvement is to phrase truncation
error bounds in terms of $\norm{A^2}^{1/2}$ instead of the (no
smaller and potentially much larger quantity) $\norm{A}$. The
second is to choose the degree of the Pad\'e approximant to
minimize the computational cost subject to achieving a desired
truncation error. A third improvement is to use an absolute,
rather than relative, error criterion in the choice of Pad\'e
approximant; this allows the use of higher degree approximants
without worsening an a priori error bound. Our theory and
experiments show that each of these modifications brings a
reduction in computational cost. Moreover, because the
modifications tend to reduce the number of double-angle steps they
usually result in a more accurate computed cosine in floating
point arithmetic. We also derive an algorithm for computing both
$\cos(A)$ and $\sin(A)$, by adapting the ideas developed for the
cosine and intertwining the cosine and sine double angle
recurrences.",
keywords = "matrix function, matrix cosine, matrix sine, matrix exponential,
double-angle formula, rounding error analysis, Schur--Parlett
method, MATLAB "
}

@techreport{Lancaster:2005:CFH,
author = "Peter Lancaster and Leiba Rodman",
title = "Canonical Forms for Hermitian Matrix Pairs
under Strict Equivalence and Congruence",
institution =  inst-MCCM,
type = t-NAREP,
number = "No. 460",
pages = "45",
month = jan,
year = "2005",
URL =          "http://www.maths.man.ac.uk/~nareports/narep460.pdf",
abstract = "After a brief historical review and an account of the canonical
forms attributed to Jordan and Kronecker, a systematic development
over the complex numbers and over the reals. These reductions are
by strict equivalence and by congruence and essentially complete
proofs are presented.  Some closely related results which can be
derived from the canonical forms are also included. They concern
simultaneous diagonalization, a new criterion for the existence of
positive definite linear combinations of a pair of hermitian
matrices, and the canonical structures of matrices which are
selfadjoint in an indefinite inner product. ",
keywords = "",
note = "To appear in SIAM Review."
}

@techreport{Lancaster:2005:IVS2,
author = "Uwe Prells and Peter Lancaster",
title = "Isospectral Vibrating Systems. {Part 2}: Structure Preserving
Transformations",
institution =  inst-MCCM,
type = t-NAREP,
number = "No. 459",
pages = "23",
month = jan,
year = "2005",
URL =          "http://www.maths.man.ac.uk/~nareports/narep459.pdf",
abstract = "The study of inverse problems for $n\times n$ systems of the form
$L(\la):=M\la^2 +D\la +K$ is continued. In this paper it is
assumed that one vibrating system is specified and the objective
is to generate isospectral families of systems, i.e. systems which
reproduce precisely the eigenvalues of the given system together
with their multiplicities. Two central ideas are developed and
used, namely, standard triples of matrices, and structure
preserving transformations.",
keywords = ""
}

@techreport{Lancaster:2005:IVS1,
author = "Peter Lancaster",
title = "Isospectral Vibrating Systems. {Part 1}: The Spectral Method",
institution =  inst-MCCM,
type = t-NAREP,
number = "No. 458",
pages = "17",
month = jan,
year = "2005",
URL =          "http://www.maths.man.ac.uk/~nareports/narep458.pdf",
abstract = "A study is made of inverse problems for $n\times n$ systems of the
form $L(\la)=M\la^2 +D\la +K$.  This paper concerns the
determination of systems in an equivalence class defined by a
fixed $2n\times 2n$ admissible Jordan matrix, i.e. a class of
isospectral systems. Constructive methods are obtained for complex
or real systems with no symmetry constraints. It is also shown how
isospectral families of complex hermitian matrices can be formed.
The case of real symmetric matrices is more difficult. Some
partial solutions are obtained but, in this case, the theory
remains incomplete. Examples are given.",
keywords = "",
note = "To appear in Linear Algebra Appl."
}

@techreport{Mathias:2004:DGE,
author = "Roy Mathias and Chi-Kwong Li",
title = "The Definite Generalized Eigenvalue Problem:
{A} New Perturbation Theory",
institution =  inst-MCCM,
type = t-NAREP,
number = "No. 457",
pages = 51,
month = oct,
year = 2004,
URL =     "http://www.maths.man.ac.uk/~nareports/narep457.pdf",
abstract = "Let $(A,B)$ be a definite pair of $n \times n$ Hermitian matrices.
That is, $|x^*Ax| + |x^*Bx| \neq 0$ for all non-zero vectors
$x \in \IC^n$. It is possible to find an $n \times n$ non-singular
matrix $X$ with unit columns such that
$$X^{*} (A + iB) X = {\rm diag}(\alpha_1 + i \beta_1, \dots, \alpha_n + i \beta_n)$$
where $\alpha_j$ and $\beta_j$ are real numbers.  We call the pairs
$(\alpha_j, \beta_j)$ {\it normalized generalized eigenvalues} of  the
definite pair $(A, B)$.
These pairs have not been studied previously.  We rework the perturbation
theory for  the eigenvalues and eigenvectors  of the definite
generalized eigenvalue problem $\beta Ax = \alpha Bx$ in terms of
these normalized generalized eigenvalues and show that they play
a crucial rule in obtaining the best possible perturbation  bounds.
In particular, in existing perturbation bounds, one can replace most
instances of the Crawford number
$c(A,B) = \min \{ |x^{*} (A + iB)x|: x \in \IC^n, x^{*}x = 1 \}$
with the larger quantity
$d_{\min} = \min\{ |\alpha_j + i\beta_j|: j = 1, \ldots, n \}.$
This results in bounds that can be stronger by an arbitrarily large
factor. We also give a new measure of the separation of the $j$th
eigenvalue from the $k$th:
$|(\alpha_j + i \beta_j) \sin (\arg(\alpha_j + i \beta_j) - \arg(\alpha_k + i \beta_k))| .$
This asymmetric measure is entirely new, and again
results in bounds that can be arbitrarily stronger  than
the existing bounds. We show that all but one of our bounds are attainable.
We also show that the Crawford number is the infimum of the distance
from a definite pencil, {\it a fortiori} diagonalizable, to a
non-diagonalizable pair."
}

@techreport{Shardlow:2004:GED,
author = "Shardlow, Tony and Yan, Yubin",
title = "Geometric Ergodicity for Dissipative Particle Dynamics",
institution =  inst-MCCM,
type = t-NAREP,
number = "No. 456",
pages = 29,
month = oct,
year = 2004,
URL =          "http://www.maths.man.ac.uk/~nareports/narep456.pdf",
abstract = "Dissipative particle dynamics is a system of stochastic
differential equations modelling complex fluid flows. We consider
the problem of $N$ particles evolving on the one dimensional
periodic domain of length $L$ and, if the density of particles is
large, prove geometric convergence to a unique invariant measure.
The proof uses minorization and drift arguments, but allows
elements of the drift and diffusion matrix to have compact support
where hypoellipticity arguments are not directly available."
}

@techreport{Shardlow:2004:MES,
author = "Shardlow, Tony",
title = "Modified Equations for Stochastic Differential Equations",
institution =  inst-MCCM,
type = t-NAREP,
number = "No. 455",
pages = 12,
month = oct,
year = 2004,
URL =          "http://www.maths.man.ac.uk/~nareports/narep455.pdf",
abstract = "We describe a backward error analysis for stochastic differential
equations with respect to weak convergence. Modified equations are
provided for forward and backward Euler approximations to Ito SDEs
with additive noise, and extensions to other types of equations
and approximation discussed."
}

@techreport{Bini:2004:AMP,
author = "Dario A. Bini and Nicholas J. Higham and Beatrice Meini",
title = "Algorithms for the Matrix {$p$}th Root",
institution =  inst-MCCM,
type = t-NAREP,
number = "No. 454",
pages = 26,
month = aug,
year = 2004,
URL =          "http://www.maths.man.ac.uk/~nareports/narep454.pdf",
abstract = "New theoretical results are presented about the principal matrix
$p$th root. In particular, we show that the $p$th root is related
to the matrix sign function and to the Wiener-Hopf factorization,
and that it can be expressed as an integral over the unit circle.
These results are used in the design and analysis of several new
algorithms for the numerical computation of the $p$th root. We
also analyze the convergence and numerical stability properties of
Newton's method for the inverse $p$th root. Preliminary
computational experiments are presented to compare the methods.",
keywords = "Matrix $p$th root, matrix sign function, Wiener-Hopf
factorization, Newton's method, Graeffe iteration, cyclic
reduction, Laurent polynomial"
}

@techreport{Powell:2004:PFP,
author = "Catherine Powell",
title = "Parameter-free $\mathbf{H(div)}$ preconditioning for mixed
finite element formulation of diffusion problems",
institution =  inst-MCCM,
type = t-NAREP,
number = "No. 453",
month = aug,
year = 2004,
URL =  "http://www.maths.man.ac.uk/~nareports/narep453.ps.gz",
abstract = "In this paper, we implement the lowest-order Raviart-Thomas scheme
for generalised diffusion problems and analyse the robustness of
so-called $H(div)$ preconditioning for the resulting saddle-point
systems. Properties of the exact version of the suggested
block-diagonal preconditioner were discussed in some previous
work,  Optimal preconditioning for Raviart-Thomas mixed
formulation of second-order elliptic problems', SIAM J. Matrix
Anal. Appl., 25, 3(2003). Here, new bounds  are established for
the eigenvalue spectrum of the  preconditioned system matrix when
a single V-cycle of any available multigrid scheme is used to
approximately invert the matrix corresponding to the discrete
weighted $H(div)$ operator."
}

@techreport{Higham:2004:SSM,
author = "Nicholas J. Higham",
title = "The Scaling and Squaring Method for the Matrix Exponential
Revisited",
institution =  inst-MCCM,
type = t-NAREP,
number = "No. 452",
pages = "15",
month = jul,
year = "2004",
URL =          "http://www.maths.man.ac.uk/~nareports/narep452.pdf",
abstract = "The scaling and squaring method is the most widely used method for
computing the matrix exponential, not least because it is the
method implemented in MATLAB's \texttt{expm} function. The method
scales the matrix by a power of 2 to reduce the norm to order 1,
computes a Pad\'e approximant to the matrix exponential, then
repeatedly squares to undo the effect of the scaling. We give a
new backward error analysis of the method (in exact arithmetic)
that employs sharp bounds for the truncation errors and leads to
an implementation of essentially optimal efficiency. We also give
new rounding error analysis that shows the computed Pad\'e
approximant of the scaled matrix to be highly accurate. For IEEE
double precision arithmetic the best choice of degree of Pad\'e
approximant turns out to be 13, rather than the 6 or 8 used by
previous authors. Our implementation of the scaling and squaring
method always requires at least two fewer matrix multiplications
than \texttt{expm} when the matrix norm exceeds 1, which can
amount to a 37\% saving in the number of multiplications, and it
is typically more accurate owing to the fewer required squarings.
We also investigate a different scaling and squaring algorithm
proposed by Najfeld and Havel that employs a Pad\'e approximation
to the function $x \coth(x)$. This method is found to be
essentially a variation of the standard one with weaker supporting
error analysis.",
keywords = "matrix function, matrix exponential, Pad\'e approximation, matrix
polynomial evaluation, scaling and squaring method, MATLAB, expm,
backward error analysis, performance profile"
}

@techreport{Roberts:2004:SVO,
author = "Leonid E. Shaikhet and Jason A. Roberts",
title = "Stochastic {Volterra} Integro-Differential Equations:
Stability and Numerical Methods",
institution = inst-MCCM,
type = t-NAREP,
number = "No. 450",
pages = 38,
month = jun,
year = 2004,
url = "http://www.maths.man.ac.uk/~nareports/narep450.pdf",
abstract = "We consider the reliability of some numerical methods in
preserving the stability properties of the linear stochastic
functional differential equation $$\dot{x}(t)=\alpha x(t)+\beta\int^t_0x(s)ds+\sigma x(t-\tau)\dot{W}(t),$$ where
$\alpha, \beta, \sigma, \tau \geq 0$ are real constants, and
$W(t)$ is a standard Wiener process. We adopt the shorthand
notation of $\dot{x}(t)$ to represent the differential $dx(t)$
etc. Our choice of test equation is a stochastic perturbation of
the classical deterministic Brunner \& Lambert test equation for
$\sigma=0$ and so our investigation may be thought of as an
extension of their work. Sufficient conditions for the asymptotic
mean square stability of solutions to both the differential
equation and discrete analogues are derived using the general
method of Lyapunov functionals construction proposed by
Kolmanovskii \& Shaikhet which has previously been successfully
employed for deterministic and stochastic differential and
difference equations with delay. The areas of the regions of
asymptotic stability for each $\theta$ method, indicated by the
sufficient conditions for the discrete system, are shown to be
equal and we show that an upper bound can be put on the time-step
parameter for the numerical method fo which the system is
asymptotically mean-square stable. We illustrate our results by
means of numerical experiments and various stability diagrams. We
examine the extent to which the continuous system can tolerate
stochastic perturbations before losing its stability properties
and we illustrate how one may accurately choose a numerical method
to preserve the stability properties of the original problem in
the numerical solution. Our numerical experiments also indicate
that the quality of the sufficient conditions is very high."
}

@techreport{MCCM:2004:AR3,
author =       "Manchester Centre for Computational Mathematics",
title =        "Annual Report: {January--December} 2003",
type =         t-NAREP,
number =       "No. 449",
institution =  inst-MCCM,
pages =        "25",
month =        may,
year =         "2004",
URL =          "http://www.maths.man.ac.uk/~nareports/narep449.pdf",
}

@techreport{Baker:2004:HAL,
author = "Christopher T. H. Baker and Evelyn Buckwar",
title = "A note on {Halanay}-type stability results for the
{$\theta$}--{Maruyama} method for stochastic delay differential
equations",
institution = inst-MCCM,
type = t-NAREP,
number = "No. 448",
pages = 13,
month = apr,
year = 2004,
url = "http://www.maths.man.ac.uk/~nareports/narep448.pdf",
abstract = "Using an approach that has its origins in work of Halanay, we
consider stability in  mean square of numerical solutions obtained
from the {\em $\theta$--Maruyama discretization} of a test
stochastic delay differential equation $dX(t) = \{f(t) -\alpha X(t) + \beta X(t-\tau)\} {\rd}t + \{g(t) + \eta \ X(t) + \mu X(t-\tau) \} \ {\rd}W(t)$ interpreted in the It\^o sense, where
$W(t)$ denotes a Wiener process (white noise). Our objective in
this report is limited to  demonstrating that we may use
techniques advanced in a recent report by Baker and Buckwar to
obtain criteria for asymptotic and exponential stability, in mean
square, for the solutions of a recurrence ${\widetilde X}_{n+1} - {\widetilde X}_n = \theta\{f_{n+1} -\A h {\widetilde X}_{n+1} + \B h {\widetilde X}_{n+1-N}\} + (1-\theta)\{f_n -\A h{\widetilde X}_{n} + \B h {\widetilde X}_{n-N}\}$ $+ (g_n + \eta \sqrt{h} \, {\widetilde X}_n + \mu \sqrt{h} \, {\widetilde X}_{n-N_{}})\ \xi_n$ where $\xi_n \in N(0,1)$ ($\xi_n$ is normally distributed
with zero mean and variance unity).  This is the $\theta$-Maruyama
formula; recognizable cases arise from taking $\theta = 0$,
$\theta = \half$ or $\theta = 1$; for convenience we assume $\theta \in [0,1]$."
}

@techreport{Hargreaves:2004:CCN,
author = "Gareth I. Hargreaves",
title = "Computing the Condition Number of Tridiagonal and
Diagonal-Plus-Semiseparable Matrices in Linear Time",
institution = inst-MCCM,
type = t-NAREP,
number = "No. 447",
pages = 19,
month = apr,
year = 2004,
url = "http://www.maths.man.ac.uk/~nareports/narep447.pdf",
abstract = "For an $n \times n$ tridiagonal matrix we exploit the structure of
its  QR factorization to devise two new algorithms for computing the 1-norm
condition number in $O(n)$ operations.  The algorithms avoid
underflow and overflow, and are simpler than existing algorithms since
tests are not required for degenerate cases.  An error analysis of the
first algorithm is given, while the second algorithm is shown to be
competitive in speed with existing algorithms.  We then turn our
attention to an $n \times n$ diagonal-plus-semiseparable matrix, $A$,
for which several algorithms have recently been developed to solve
$Ax=b$ in $O(n)$ operations.  We again exploit the QR factorization of
the matrix to present an algorithm that computes the 1-norm condition
number in $O(n)$ operations. ",
}

@techreport{Higham:2004:FPM,
author = "Nicholas J. Higham and D. Steven Mackey and Niloufer Mackey and
Fran\c{c}oise Tisseur",
title = "Functions Preserving Matrix Groups
and Iterations for the Matrix Square Root",
institution =  inst-MCCM,
type = t-NAREP,
number = "446",
pages = 29,
month = mar,
year = "2004",
URL =          "http://www.maths.man.ac.uk/~nareports/narep446.pdf",
abstract = "For any matrix automorphism group $\G$ associated with a bilinear
or sesquilinear form, Mackey, Mackey, and Tisseur have recently
shown that the matrix sign decomposition factors of $A\in\G$ also
lie in $\G$; moreover, the polar factors of $A$ lie in $\G$ if the
matrix of the underlying form is unitary. Groups satisfying the
latter condition include the complex orthogonal, real and complex
symplectic, and pseudo-orthogonal groups. This work is concerned
with exploiting the structure of $\G$ when computing the polar and
matrix sign decompositions of matrices in $\G$. We give sufficient
conditions for a matrix iteration to preserve the group structure
and show that a family of globally convergent rational
Pad\'e-based iterations of Kenney and Laub satisfy these
conditions. The well-known scaled Newton iteration for computing
the unitary polar factor does not preserve group structure, but we
show that the approach of the iterates to the group is precisely
tethered to the approach to unitarity, and that this forces a
different and exploitable structure in the iterates. A similar
relation holds for the Newton iteration for the matrix sign
function. We also prove that the number of iterations needed for
convergence of the structure-preserving methods can be precisely
predicted by running an associated scalar iteration. Numerical
experiments are given to compare the cubically and quintically
converging iterations with Newton's method and to test stopping
criteria. The overall conclusion is that the structure-preserving
iterations and the scaled Newton iteration are all of practical
interest, and which iteration is to be preferred is
problem-dependent.",
keywords = "automorphism group, bilinear form, sesquilinear form, scalar
product, adjoint, Fr\'echet derivative, stability analysis,
perplectic matrix, pseudo-orthogonal matrix, generalized polar
decomposition, matrix sign function, matrix $p$th root, matrix
square root, structure preservation, matrix iteration, Newton
iteration"
}

@techreport{Lancaster:2004:OPM,
author = "Peter Lancaster and Panayiotis Psarrakos",
title = "On the Pseudospectra of Matrix Polynomials",
institution =  inst-MCCM,
type = t-NAREP,
number = "No. 445",
pages = "17",
month = feb,
year = "2004",
URL =          "http://www.maths.man.ac.uk/~nareports/narep445.pdf",
abstract = "The pseudospectra of a matrix polynomial $P(\la)$ are sets of
complex numbers that are eigenvalues of matrix polynomials which
are near to $P(\la)$, i.e., their coefficients are within some
fixed magnitude of the coefficients of $P(\la)$. Pseudospectra
provide important insights into the sensitivity of eigenvalues
under perturbations, and have several applications. First,
qualitative properties concerning boundedness and connected
components of pseudospectra are obtained. Then an accurate
continuation algorithm for the numerical determination of the
boundary of pseudospectra of matrix polynomials is devised and
illustrated. This algorithm is based on a prediction-correction
scheme.",
keywords = "matrix polynomial, eigenvalue, singular value, perturbation,
$\eps$-pseudospectum, boundary, stability"
}

@techreport{Baker:2004:INP,
author = "Christopher T.H.Baker and Eugene I. Parmuzin",
title = "Identification of the initial function for delay differential
equations: Part III: Nonlinear DDEs \& computational results.",
institution = inst-MCCM,
type = t-NAREP,
number = "No. 444",
pages = 27,
month = mar,
year = 2004,
url = "http://www.maths.man.ac.uk/~nareports/narep444.pdf",
abstract = "In this report we consider the data assimilation problem'' for
nonlinear delay differential equations. This problem consists of
finding an initial function that gives rise to a solution of a
given nonlinear delay differential equation, which is a close fit
to observed data. A r\^ole for adjoint equations and fundamental
solutions in the nonlinear case is established. The
Pseudo-Newton'' method to solve the data assimilation
problem'' for nonlinear DDE is presented. This is an extension of
results obtain in Nareport No.431 for linear delay differential
equations."
}

@techreport{Baker:2004:IDP,
author = "Christopher T. H. Baker and Eugene I. Parmuzin",
title = "Identification of the initial function for delay differential
equations: Part II: A discrete analogue of the data assimilation
problem \& computational results.",
institution = inst-MCCM,
type = t-NAREP,
number = "No. 443",
pages = 29,
month = mar,
year = 2004,
url = "http://www.maths.man.ac.uk/~nareports/narep443.pdf",
abstract = "In this report we consider a difference analogue of the data
assimilation problem'' for delay differential equations (DDEs)
presented in the first part (Numerical Analysis Report No. 431).
The problem consists of finding an initial function that gives rise
to a solution of a discretized DDE, which is a close fit to
observed data. A r\^ole for adjoint equations and fundamental
solutions in the discrete case is established, and related discrete
integral equations (summation equations) are obtained. The final
part of  the report consists of a report of numerical experiments
that  demonstrate the performance of the algorithm."
}

@techreport{Lucas:2004:LSC,
author = "Craig Lucas",
title = "{LAPACK}-Style Codes for Level 2 and 3
Pivoted {Cholesky} Factorizations",
institution =  inst-MCCM,
type = t-NAREP,
number = "No. 442",
pages = "31",
month = feb,
year = "2004",
URL =          "http://www.maths.man.ac.uk/~nareports/narep442.pdf",
abstract = "Fortran 77 codes exist in LAPACK for computing the Cholesky
factorization of a symmetric positive definite matrix using
Level 2 and 3 BLAS. In LINPACK there is a Level 1 BLAS
routine for computing the Cholesky factorization with complete
pivoting of a symmetric positive semi-definite matrix. We present
two  new algorithms and Fortran 77 LAPACK-style codes for
computing this pivoted factorization: one using Level 2 BLAS
and one using Level 3 BLAS. We show that on modern machines the
new codes can be many times faster than the LINPACK code. Also,
with a new  stopping criterion they provide more reliable rank
detection and can have a smaller normwise backward error.",
keywords = "pivoted Cholesky factorization, semi-definite matrices,
LAPACK, Level 3 BLAS."
}

@techreport{Davies:2004:SCM,
author = "Philip I. Davies",
title = "Structured Conditioning of Matrix Functions",
institution =  inst-MCCM,
type = t-NAREP,
number = "No. 441",
pages = "33",
month = jan,
year = "2004",
URL =          "http://www.maths.man.ac.uk/~nareports/narep441.pdf",
abstract = "The existing theory of conditioning for matrix functions
$f(X) \colon \C^{\nbyn} \to \C^{\nbyn}$ does not cater for structure
in the matrix $X$. An extension of this theory is presented in which
when $X$ has structure, all perturbations of $X$ are required to
have the same structure. We consider two classes of structured
matrices, those comprising the Jordan algebra $\J$ and the Lie
algebra $\L$ associated with a nondegenerate bilinear or sesquilinear
form on $\R^n$ or $\C^n$. Examples of such classes are the symmetric,
skew-symmetric, Hamiltonian and skew-Hamiltonian matrices. Structured
condition numbers are defined for these two classes. Under certain
conditions on the underlying scalar product, explicit representations
are given for the structured condition numbers. Comparisons between
the unstructured and structured condition numbers are then made.
We show that when the underlying scalar product is a sesquilinear
form, there is no difference between the values of the two condition
numbers for (i) all functions of $X \in \J$,  and (ii) odd and even
functions of $X \in \L$. When the underlying scalar product is a
bilinear form then equality is not guaranteed in all these cases.
Where equality is not guaranteed, bounds are obtained for the ratio
of the unstructured and structured condition numbers.",
keywords = "Matrix functions, Fr\'echet derivative, Condition numbers, Bilinear
forms, Sesquilinear forms, Structured Matrices, Jordan algebra,
Symmetric matrices, Lie algebra, Skew-symmetric matrices"
}

@techreport{Higham:2003:NSB,
author = "Nicholas J. Higham",
title = "The Numerical Stability of Barycentric {Lagrange} Interpolation",
institution =  inst-MCCM,
type = t-NAREP,
number = "No. 440",
pages = "10",
month = dec,
year = "2003",
URL =          "http://www.maths.man.ac.uk/~nareports/narep440.pdf",
abstract = "The Lagrange representation of the interpolating polynomial can be
rewritten in two more computationally attractive forms: a modified
Lagrange form and a barycentric form. We give an error analysis of
the evaluation of the interpolating polynomial using these two
forms. The modified Lagrange formula is shown to be backward
stable. The barycentric formula has a less favourable error
analysis, but is forward stable for any set of interpolating
points with a small Lebesgue constant. Therefore the barycentric
formula can be significantly less accurate than the modified
Lagrange formula only for a poor choice of interpolating points.
This analysis provides further weight to the argument of Berrut
and Trefethen that barycentric Lagrange interpolation should be
the polynomial interpolation method of choice.",
keywords = "polynomial interpolation, Lagrange interpolation, barycentric
formula, rounding error analysis, backward error, forward error,
Lebesgue constant"
}

@techreport{Buckwar:2003:WAS,
author = "Buckwar, Evelyn and Shardlow, Tony",
title = "Weak Approximation of Stochastic Delay Differential Equations",
institution =  inst-MCCM,
type = t-NAREP,
number = "No. 439",
month = dec,
year = "2003",
URL = "http://www.maths.man.ac.uk/~nareports/narep439.pdf"
}

@techreport{Shardlow:2003:NWE,
author = "Shardlow, Tony",
title = "Nucleation of Waves in Excitable Media by Noise",
institution =  inst-MCCM,
type = t-NAREP,
number = "No. 438",
month = dec,
year = "2003",
URL = "http://www.maths.man.ac.uk/~nareports/narep438.pdf"
}

@techreport{Shardlow:2003:NSS,
author = "Shardlow, Tony",
title = "Numerical Simulation of Stochastic {PDE}s for Excitable Media",
institution =  inst-MCCM,
type = t-NAREP,
number = "No. 437",
month = dec,
year = "2003",
URL = "http://www.maths.man.ac.uk/~nareports/narep437.pdf"
}

@techreport{Davies:2003:CFM,
author = "Philip I. Davies and Nicholas J. Higham",
title = "Computing {$f(A)b$} for Matrix Functions {$f$}",
institution =  inst-MCCM,
type = t-NAREP,
number = "No. 436",
pages = "10",
month = nov,
year = "2003",
URL =          "http://www.maths.man.ac.uk/~nareports/narep436.pdf",
abstract = "For matrix functions $f$ we investigate how to compute a
matrix-vector product $f(A)b$ without explicitly
computing $f(A)$. A general method is described that
applies quadrature to the matrix version of the Cauchy
integral theorem. Methods specific to the logarithm,
based on quadrature, and fractional matrix powers, based
on solution of an ordinary differential equation initial
value problem, are also described."
}

@techreport{Ford:2003:SGB,
author = "Judith M. Ford and Neville J. Ford and John Wheeler",
title = "Simulation of Grain Boundary Creep: Analysis of Some New Numerical
Techniques",
institution =  inst-MCCM,
type = t-NAREP,
number = "No. 435",
month = oct,
year = "2003",
URL = "http://www.maths.man.ac.uk/~nareports/narep435.pdf"
}

@techreport{Ford:2003:SLS,
author = "Judith M. Ford and Eugene E. Tyrtyshnikov",
title = "Solving Linear Systems using Wavelet Compression Combined with
{Kronecker} Product Approximation",
institution =  inst-MCCM,
type = t-NAREP,
number = "No. 434",
month = oct,
pages = "11",
year = "2003",
URL = "http://www.maths.man.ac.uk/~nareports/narep434.pdf",
note = "Submitted to special edition of Numerical Algorithms, Proceedings of
Quatrième séminaire sur l'algorithmique numérique appliquée aux
problèmes industriels, Calais, 2003"
}

@techreport{Alam:2003:OCN,
author = "Rafikul Alam",
title = "On the Construction of Nearest Defective Matrices to a Normal
Matrix",
institution =  inst-MCCM,
type = t-NAREP,
number = "No. 433",
pages = "6",
month = sep,
year = "2003",
URL =          "http://www.maths.man.ac.uk/~nareports/narep433.pdf",
abstract = "Let $A$ be an $n$-by-$n$ normal matrix with $n$ distinct
eigenvalues. We describe a simple procedure to construct a
defective matrix whose distance from $A,$ among all the defective
matrices, is the shortest. Our construction requires only an
appropriate pair of eigenvalues of $A$ and their corresponding
eigenvectors. We also show that the nearest defective matrices
influence the evolution of spectral projections of $A$.",
keywords = "Defective matrix, multiple eigenvalue,
pseudospectra, eigenvector, spectral projection"
}

@techreport{Mackey:2003:SFS,
author = "D. Steven Mackey and Niloufer Mackey and Fran\c{c}oise Tisseur",
title =  "Structured Factorizations in Scalar Product Spaces",
type =         t-NAREP,
number =       "No. 432",
institution =  inst-MCCM,
pages =        "",
month =        "",
year = "2003",
URL =          "http://www.maths.man.ac.uk/~nareports/narep432.pdf",
note = "In preparation."
}

@techreport{Baker:2004:ICP,
author = "Christopher T. H. Baker and Eugene I. Parmuzin",
title = "Identification of the initial function for delay differential
equations: Part I: The continuous problem \& an integral equation
analysis.",
institution = inst-MCCM,
type = t-NAREP,
number = "No. 431",
pages = 25,
month = jan,
year = 2004,
url = "http://www.maths.man.ac.uk/~nareports/narep431.pdf",
abstract = "In this report, which is the first of three related reports, we
consider the data assimilation problem'' for delay differential
equations. This problem consists of finding an initial function
that gives rise to a solution which is a close fit to observed
data. A r\^ole for adjoint equations and fundamental solutions is
established, and related integral equations are obtained."
}

@techreport{Baker:2004:GTV,
author = "Christopher T. H. Baker and Eugene I. Parmuzin",
title = "A Guided Tour of Variation of Parameters Formulae for
Continuous and Discretized {DDE}s.",
institution = inst-MCCM,
type = t-NAREP,
number = "No. 430",
pages = 21,
month = jan,
year = 2004,
url = "http://www.maths.man.ac.uk/~nareports/narep430.pdf",
abstract = "In this report we discuss variation of parameters formulae ($i$)
for systems of linear delay differential equations (DDEs), $$\frac{\dss d y(t)}{\dss dt} -A(t) y(t) -B(t) y(t-\tau) = f(t), \; \mbox{for}\ t\in [0,T],$$ subject to an initial condition
$y(t) =\varphi(t) \quad \mbox{for} \ t \in [-\tau,0]$,  and ($ii$)
for the corresponding discretized equations that result from the
application of Euler methods. The explicit Euler method for the
above DDE, using a step $h = \tau / N$ with $T = Kh$ and $N,K \in \Nset$, yields the equations $${\widetilde{y}}_{n+1} - {\widetilde{y}}_{n} -A_n {\widetilde{y}}_n -B_n {\widetilde{y}}_{n-N} =f_n, \quad \mbox{for}\ n \in \{0,1, \cdots, K-1 \},$$ subject to the initial condition ${\widetilde{y}}_n = \varphi_n \; \mbox{for} \ -N \leq n \leq 0,$ where $B_n = B(nh)$,
$A_n = A(nh)$ $f_n = f(nh)$ and $\varphi_n = \varphi(nh)$. Of
interest is the form of the solution to the DDE or its discretized
form, in a representation that depends upon the inhomogeneous term
$f$ or the initial function $\varphi$. Various ways to derive
variation of parameters formulae will be pursued. In the
continuous case, the fundamental solution can be obtained as a
solution of a DDE, a related integral equation, by seeking a
generalized Green's function or by considering adjoint
equations.(A number of results related to the continuous case are
presented in the literature and these are recalled in an
appendix.)  We look for analogous results in the discretized case,
drawing on parallels in linear algebra."
}

@techreport{MCCM:2003:AR2,
author =       "Manchester Centre for Computational Mathematics",
title =        "Annual Report: {January--December} 2002",
type =         t-NAREP,
number =       "No. 429",
institution =  inst-MCCM,
pages =        "19",
month =        jun,
year =         "2003",
URL =          "http://www.maths.man.ac.uk/~nareports/narep429.pdf",
}

@techreport{Tisseur:2003:EAM,
author = "Dario A. Bini and Luca Gemignani and Fran\c{c}oise Tisseur",
title = "The {Ehrlich-Aberth} Method for the Nonsymmetric
Tridiagonal Eigenvalue Problem",
institution =  inst-MCCM,
type = t-NAREP,
number = "No. 428",
pages = "23",
month = jun,
year = "2003",
URL =  "ftp://ftp.ma.man.ac.uk/pub/narep/narep428.pdf",
abstract = "An algorithm based on the Ehrlich-Aberth iteration is presented
for the computation of the zeros of $p(\l)=\det(T-\l I)$, where
$T$ is an irreducible tridiagonal matrix. The algorithm requires
the evaluation of $p(\l)/p'(\l)=-1/\trace(T-\l I)^{-1}$, which is
done here by exploiting the QR factorization of $T-\l I$ and the
semiseparable structure of $(T-\l I)^{-1}$. Two choices of the
initial approximations are considered; the most effective relies
on a divide-and-conquer strategy, and some results motivating this
strategy are given. A Fortran 95 module implementing the algorithm
is provided and numerical experiments that confirm the
effectiveness and the robustness of the approach are presented. In
particular, comparisons with the LAPACK subroutine \texttt{dhseqr}
show that our algorithm is faster for large dimensions.",
keywords = "nonsymmetric eigenvalue problem,
symmetric indefinite generalized eigenvalue problem,
tridiagonal matrix, root finder, QR decomposition,
divide and conquer."
}

@techreport{Mackey:SPA:2003,
author = "D. Steven Mackey and Niloufer Mackey and Daniel M. Dunlavy",
title = "Structure Preserving Algorithms for Perplectic Eigenproblems",
institution =  inst-MCCM,
type = t-NAREP,
number = "No. 427",
pages = "30",
month = may,
year = "2003",
URL =          "http://www.maths.man.ac.uk/~nareports/narep427.pdf"
}

@techreport{Higham:2003:CPD,
author = "Nicholas J. Higham and D. Steven Mackey and Niloufer Mackey and
Fran\c{c}oise Tisseur",
title = "Computing the Polar Decomposition and the Matrix Sign Decomposition
in Matrix Groups",
institution =  inst-MCCM,
type = t-NAREP,
number = "No. 426",
pages = "15",
month = apr,
year = "2003",
URL =          "http://www.maths.man.ac.uk/~nareports/narep426.pdf",
abstract = "For any matrix automorphism group $\G$ associated with a
bilinear or sesquilinear form, Mackey, Mackey, and Tisseur have
recently shown that the matrix sign decomposition factors of
$A\in\G$ also lie in $\G$; moreover, the polar factors of $A$
lie in $\G$ if the matrix of the underlying form is unitary.
Groups satisfying the latter condition include the complex
orthogonal, real and complex symplectic, and pseudo-orthogonal
groups. This work is concerned with exploiting the structure of
$\G$ when computing the polar and matrix sign decompositions of
matrices in $\G$. We give sufficient conditions for a matrix
iteration to preserve the group structure and show that a family
of globally convergent rational Pad\'e-based iterations of
Kenney and Laub satisfy these conditions. The well-known scaled
Newton iteration for computing the unitary polar factor does not
preserve group structure, but we show that the approach of the
iterates to the group is precisely tethered to the approach to
unitarity, and that this forces a different and exploitable
structure in the iterates. A similar relation holds for the
Newton iteration for the matrix sign function. We also prove
that the number of iterations needed for convergence of the
structure-preserving methods can be precisely predicted by
running an associated scalar iteration. Numerical experiments
are given to compare the cubically and quintically converging
iterations with Newton's method and to test stopping criteria.
The overall conclusion is that the structure-preserving
iterations and the scaled Newton iteration are all of practical
interest, and which iteration is to be preferred is
problem-dependent.",
keywords = "automorphism group, bilinear form, sesquilinear form, adjoint,
complex orthogonal matrix, symplectic matrix, perplectic matrix,
pseudo-orthogonal matrix, polar decomposition, matrix sign
decomposition, structure preservation, matrix iteration, Newton
iteration, convergence tests"
}

@techreport{Baker:2003:MDC,
author = "Christopher T. H. Baker and Gennadii A. Bocharov  and Christopher
A. H. Paul and Fathalla A. Rihan",
title = "Models with Delays for Cell Population Dynamics:
Identification, Selection and Analysis. Part I: Computational
Modelling with Functional Differential Equations:
Identification, Selection, and Sensitivity",
institution =  inst-MCCM,
type = t-NAREP,
number = "No. 425",
pages = "28",
month = feb,
year = "2003",
URL =          "http://www.maths.man.ac.uk/~nareports/narep425.pdf",
abstract = "This is the first in a series of reports in which we discuss aspects
of the identification of, discrimination between, and sensitivity of,
mathematical models, and we illustrate with some examples.
Mathematical models based upon certain types of differential
equations, functional differential equations, or systems of
such equations, are often employed to represent the dynamics
of natural, in particular biological, phenomena.  We present
some of the principles underlying the choice of a methodology
for the computational development of quantitatively consistent
models, using scientifically meaningful parameters, based on
observational data.
Exploiting the link between a (weighted) least-squares fit,
information theory, and maximum likelihood, we can select an
appropriate objective function which is minimized over a set
of parameters. To compare the best fit models from amongst a
hierarchy of possible models we can employ indicators that can
be regarded as taking parsimony into account. A sensitivity
analysis provides feedback on the covariances of the
parameters in relation to the best fit and on whether we need
to include all the parameters in the model.  We attempt to
formulate a methodology for assessing the adequacy of
mathematical models, founded on this approach.
Much of our thinking is influenced by discussions that are
dispersed amongst the existing literature. However, we believe
that the value of our material is more than the sum of that of
its separate parts. Additionally, in a number of places we
extend or supplement or reformulate existing results. The
models that we consider generate infinite-dimensional
dynamical systems.",
keywords = " Cell growth, model selection, neutral delay
differential equation, parameter estimation, first \& second
order sensitivity analysis, nonlinear bias, objective
functions, methodologies."
}

@techreport{Higham:2003:SCS,
author = "Nicholas J. Higham and Mihail Konstantinov and Volker Mehrmann
and Petko Petkov",
title = "Sensitivity of Computational Control Problems",
institution =  inst-MCCM,
type = t-NAREP,
number = "No. 424",
pages = "32",
month = mar,
year = "2003",
URL =          "http://www.maths.man.ac.uk/~nareports/narep424.pdf",
abstract = "It is well-known that many factors contribute to the accurate and
efficient numerical solution of mathematical problems such as
those arising in computational control system design. In simple
terms these are the arithmetic of the machine on which the
calculations are carried out, sensitivity (or conditioning) of the
mathematical model to small changes of the data and the numerical
stability of the algorithms. It happens quite often that these
concepts are confused. We define these concepts and demonstrate
some of the subtleties that often lead to confusion. In particular
we demonstrate with several examples what may happen when a
problem is modularized, i.e., split into subproblems for which
computational modules are available. For three classical problems
in computational control, pole placement, linear quadratic control
and optimal $H_\infty$ control, we then discuss the conditioning
of the problems and point out sources of difficulties. We give
some ill-conditioned examples for which even numerically stable
methods fail. We also stress the need for condition and error
estimators that supplement the numerical algorithm and inform the
user about potential or actual difficulties, and we explain what
can be done to avoid these difficulties.",
keywords = "numerical stability, machine
arithmetic, pole placement, linear quadratic control,
algebraic Riccati equation, $H_\infty$ control"
}

@techreport{Mackey:2003:ODS,
author = "D. Steven Mackey and Niloufer Mackey",
title = "On the Determinant of Symplectic Matrices",
institution =  inst-MCCM,
type = t-NAREP,
number = "No. 422",
pages = "12",
month = feb,
year = "2003",
URL =  "http://www.maths.man.ac.uk/~nareports/narep422.pdf",
abstract = "A collection of new and old proofs showing that the determinant of
any symplectic matrix is $+1$ is presented. Structured
factorizations of symplectic matrices play a key role in several
arguments. A constructive derivation of the symplectic analogue of
the Cartan-Dieudonn\'{e} theorem is one of the new proofs in this
essay.",
keywords = "symplectic, determinant, bilinear form, skew-symmetric,
structure-preserving factorizations, symmetries, transvections,
$\G$-reflector, Pfaffian."
}

@techreport{Tisseur:04:SFSPS,
author = "D. Steven Mackey and Niloufer Mackey and Fran\c{c}oise Tisseur",
title = "Structured Factorizations in Scalar Product Spaces",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 421",
pages = "30",
month = Nov,
year = "2004",
URL =  "http://www.maths.man.ac.uk/~nareports/narep421.pdf",
abstract = "
Let $A$ belong to an automorphism group, Lie algebra or Jordan
algebra of a scalar product. When $A$ is factored, to what extent
do the factors inherit structure from $A$? We answer this question
for the principal matrix square root, the matrix sign
decomposition, and the polar decomposition. For general $A$, we
give a simple derivation and characterization of a particular
generalized polar decomposition, and we relate it to other such
decompositions in the literature. Finally, we study
eigendecompositions and structured singular value decompositions,
considering in particular the structure in eigenvalues,
eigenvectors and singular values that persists across a wide range
of scalar products.
A key feature of our analysis is the identification of two
particular classes of scalar products, termed unitary and
orthosymmetric, which serve to unify assumptions for the existence
of %structured factors and structured factorizations. A variety of
different characterizations of these scalar product classes is
given.",
keywords =  "automorphism group, Lie group, Lie algebra, Jordan algebra,
bilinear form, sesquilinear form, scalar product, indefinite inner
Hamiltonian, pseudo-orthogonal, polar decomposition, matrix sign
function, matrix square root, generalized polar decomposition,
eigenvalues, eigenvectors, singular values, structure
preservation."
}

@techreport{Tisseur:2003:GRS,
author = "D. Steven Mackey and Niloufer Mackey
and Fran\c{c}oise Tisseur",
title = "$\mathbb{G}$-Reflectors in Scalar Product Spaces",
institution =  inst-MCCM,
type = t-NAREP,
number = "No. 420",
pages = "24",
month = feb,
year = "2003",
URL =  "http://www.maths.man.ac.uk/~nareports/narep420.ps.gz",
abstract = " We characterize the analogues of Householder reflectors in matrix
groups associated with scalar products. Examples of such groups
include the symplectic and pseudo-unitary groups. We also
precisely delimit the mapping capabilities of these Householder
analogues: given a matrix group $\G$ and vectors $x, y$, we give
necessary and sufficient conditions for the existence of a
Householder-like analogue $G \in \G$ such that $Gx = y$. When $G$
exists, we show how it can be constructed from $x$ and $y$.",
keywords =  "scalar product, bilinear, sesquilinear, isotropic, Householder
reflection, hyperbolic transformation, symmetries, transvections,
symplectic, pseudo-unitary, structure-preserving."
}

@techreport{Tisseur:2003:STS,
author = "D. Steven Mackey and Niloufer Mackey and Fran\c{c}oise Tisseur",
title = "Structured Tools for Structured Matrices",
institution =  inst-MCCM,
type = t-NAREP,
number = "No. 419",
pages = "36",
month = feb,
year = "2003",
URL =  "http://www.maths.man.ac.uk/~nareports/narep419.pdf",
abstract = "We present an extensive and unified collection of
structure-preserving transformations, organized for easy
reference. The structures we work with arise in the context of a
non-degenerate bilinear or sesquilinear form on $\R^n$ or $\C^n$.
We construct a variety of transformations belonging to the
automorphism groups of these forms, that imitate the action of
Givens rotations, Householder reflectors, and Gauss
transformations. We also describe transformations for performing
structured scaling actions. The matrix groups considered in this
paper are the complex orthogonal, real, complex and conjugate
symplectic, real perplectic, real and complex pseudo-orthogonal,
and pseudo-unitary groups. In addition to deriving new
transformations, this paper collects and unifies existing
structure-preserving tools.",
keywords = "structured matrices, matrix groups, Givens rotations, Householder
reflections, complex orthogonal, symplectic, complex symplectic,
conjugate symplectic, real perplectic, pseudo-orthogonal, complex
pseudo-orthogonal, pseudo-unitary, scalar product, bilinear form,
sesquilinear form, automorphism groups, Jordan algebra, Lie
algebra."
}

@techreport{Ford:2003:CKP,
author = "Judith M. Ford and Eugene E. Tyrtyshnikov",
title = "Combining {Kronecker} Product Approximation with Discrete Wavelet
Transforms to Solve Dense, Function-Related Linear Systems",
institution = inst-MCCM,
type =    t-NAREP,
number = "No. 418",
month = jan,
year = 2003,
URL = "http://www.maths.man.ac.uk/~nareports/narep418.pdf"
}

@techreport{Baker:2003:NDE,
author = "Christopher T. H. Baker and Christopher A. H. Paul",
title = "Piecewise Continuous Solutions of Neutral Delay
Differential Equations",
institution = inst-MCCM,
type = t-NAREP,
number = "No. 417",
pages = 35,
month = Dec,
year = 2003,
url =  "http://www.maths.man.ac.uk/~nareports/narep417.pdf",
abstract = "Examples of neutral delay differential equations (NDDEs) are $\frac{d}{dt} u(t) = F\big(t,u(t),u(t-\tau),\frac{d}{dt} u(t-\tau)\big) \quad (t \ge t_0),$ for $\tau > 0$, and (a form
frequently called Hale's form'') $\frac{d}{dt} \{u(t)-g\big(t,u(t),u(t-\tau)\big)\} = f\big(t,u(t),u(t-\tau)\big) \quad (t \ge t_0)$ which is an NDDE when $g$ does not vanish
identically.  A particular solution of such an equation is defined
by initial data on $[-\tau,0]$, and we observe in this report
that, when the problem is suitably interpreted, the solutions of
an NDDE need not be continuous. This fact, largely overlooked in
the numerical analysis literature (many of the texts on the
analysis of NDDEs admit a function as a solution only if it is
continuous), presents some difficulties when considering numerical
methods for approximating solutions. Our discussion of NDDEs
points to the need to be precise in the definition of the term
solution'', and to consider questions of \emph{existence},
\emph{uniqueness} and \emph{smoothness} of any solutions of NDDEs.
Such questions are a precursor to addressing the problem of
constructing numerical methods. Numerical codes for solving
differential equations are generally based on the assumption that
a solution is continuous, and our discussion will indicate
measures that may be taken if one seeks a discontinuous solution
of an NDDE. First we show that the difficulties can sometimes be
alleviated by \emph{perturbing the initial function} to yield a
continuously differentiable solution. Alternatively,
\emph{tracking discontinuities} and converting the problem to a
\emph{delay differential algebraic equation} (DDAE) -- followed
either by direct solution of the DDAE or by treatment after
regularizing the DDAE to a singularly perturbed problem -- are
discussed in the Appendices. The investigation of the simple NDDE
that opens our discussion is followed by a number of results for
general NDDEs."
}

@techreport{Hargreaves:2002:IAM,
author = "Gareth I. Hargreaves",
title = "Interval Analysis in {MATLAB}",
institution = inst-MCCM,
type =    t-NAREP,
number = "No. 416",
pages = 49,
month = dec,
year = 2002,
url = "http://www.maths.man.ac.uk/~nareports/narep416.pdf",
abstract = "The introduction of fast and efficient software for interval
arithmetic, such as the MATLAB toolbox INTLAB, has resulted in the
increased popularity of the use of interval analysis.  We give an
introduction to interval arithmetic and explain how it is
implemented in the toolbox INTLAB.  A tutorial is provided for
those who wish to learn how to use INTLAB.  We then focus on the
interval versions of some important problems in numerical
analysis.  A variety of techniques for solving interval linear
systems of equations are discussed, and these are then tested to
compare timings and accuracy.  We consider univariate and
multivariate interval nonlinear systems and describe algorithms
that enclose all the roots. Finally, we give an application of
interval analysis.  Interval arithmetic is used to take account of
rounding errors in the computation of Viswanath's constant, the
rate at which a random Fibonacci sequence increases.",
keywords = "interval analysis, interval arithmetic, INTLAB, MATLAB, tutorial,
interval linear equations, interval nonlinear equations,
Viswanath's constant"
}

@techreport{Powell:2002:BBP,
author = "Catherine Powell and David Silvester",
title = "Black-Box Preconditioning for Mixed Formulation of Self-Adjoint
Elliptic {PDE}s",
institution =  inst-MCCM,
type = t-NAREP,
number = "No. 415",
month = nov,
year = 2002,
note = "Submitted to the conference proceedings of CISC 2002 for publication
in Springer Lecture Notes in Computer Science and Engineering",
URL =  "http://www.maths.man.ac.uk/~nareports/narep415.ps.gz",
abstract = "Mixed finite element approximation of self-adjoint elliptic PDEs
leads to symmetric indefinite linear systems of equations.
Preconditioning strategies commonly focus on reduced symmetric
positive definite systems and require nested iteration. This
defficiency is avoided if preconditioned \textsc{minres} is
applied to the full indefinite system.  We outline such a
preconditioning strategy, the key building block for which is a
fast solver for a scalar diffusion operator based on black-box
algebraic multigrid. Numerical results are presented for the
Stokes equations arising in incompressible flow modelling and a
variable diffusion equation that arises in modelling potential
flow. We prove that the eigenvalues of the preconditioned matrices
are contained in intervals that are bounded independently of the
discretisation parameter and the PDE coefficients."
}

@techreport{Agyingi:2003:VOP,
author = "Ephraim Agyingi and Christopher T. H. Baker",
title = "Variation of Parameters Formulae for {Volterra} Integral Equations",
institution = inst-MCCM,
type = t-NAREP,
number = "No. 414",
pages = 20,
month = Nov,
year = 2003,
url = "http://www.maths.man.ac.uk/~nareports/narep414.pdf",
abstract = "In this report, the authors obtain variation of parameters
formulae for nonlinear Volterra integral equations and
integro-differential equations. These are related to a variation
of parameters formula for ordinary differential equations that is
due to V.M.~Alekseev.  The authors' results are obtained through
the use of an embedding technique. Some doubts have been cast on
previous results in the literature. Amongst other things, the
authors show that a result on which much of the numerical
analysis for Volterra equations has been founded in the past can
be obtained by a new approach. "
}

@techreport{Bojanczyk:2002:CIL,
author = "Adam Bojanczyk and Nicholas J. Higham and Harikrishna Patel",
title = "The Constrained Indefinite Least Squares Problem: {Theory} and
Algorithms",
institution =  inst-MCCM,
type = t-NAREP,
number = "No. 413",
pages = "11",
month = oct,
year = "2002",
URL =          "http://www.maths.man.ac.uk/~nareports/narep413.pdf",
abstract = "We present theory and algorithms for the equality constrained
indefinite least squares problem, which requires minimization of
an indefinite quadratic form subject to a linear equality
constraint. A generalized hyperbolic QR factorization is
introduced and used in the derivation of perturbation bounds and
to construct a numerical method. An alternative method is obtained
by employing a generalized QR factorization in combination with a
Cholesky factorization. Rounding error analysis is given to show
that both methods have satisfactory numerical stability properties
and numerical experiments are given for illustration. This work
builds on recent work on the unconstrained indefinite least
squares problem by Chandrasekaran, Gu and Sayed and by the
present authors.",
keywords = "equality constrained indefinite least squares problem,
$J$-orthogonal matrix, hyperbolic rotation, hyperbolic QR
factorization, generalized hyperbolic QR factorization, rounding
error analysis, forward stability, perturbation theory, Cholesky
factorization"
}

@techreport{Ford:2002:IDP,
author = "Judith M. Ford",
title = "An Improved {DWT}-Based Preconditioner for Dense Matrix Problems",
institution = inst-MCCM,
type = t-NAREP,
number = "No. 412",
month = sep,
year = 2002,
URL = "http://www.maths.man.ac.uk/~nareports/narep412.ps.gz",
abstract = "We present a new preconditioning method for a class of dense
matrices, based on sparse approximation in a discrete wavelet
basis. We first prove theoretical results that enable us to reduce
the cost of an existing wavelet-based preconditioner and then
incorporate those ideas into the design of our new preconditioner.
We demonstrate the effectiveness of our methods with numerical
examples drawn from physical problems including the
two-dimensional Helmholtz equation for which our preconditioner
gives savings of about 40\% compared with existing discrete
wavelet preconditioners and of up to 87\% compared with direct
solution."
}

@techreport{Higham:2002:CMC,
author = "Nicholas J. Higham and Matthew I. Smith",
title = "Computing the Matrix Cosine",
institution =  inst-MCCM,
type = t-NAREP,
number = "No. 411",
pages = 15,
month = sep,
year = 2002,
URL =          "http://www.maths.man.ac.uk/~nareports/narep411.pdf",
abstract = "An algorithm is developed for computing the matrix cosine,
building on a proposal of Serbin and Blalock. The algorithm scales
the matrix by a power of 2 to make the $\infty$-norm less than or
equal to 1, evaluates a Pad\'e approximant, and then uses the
double angle formula $\cos(2A) = 2\cos(A)^2-I$ to recover the
cosine of the original matrix. In addition, argument reduction and
balancing is used initially to decrease the norm. We give
truncation and rounding error analyses to show that an [8,8]
Pad\'e approximant produces the cosine of the scaled matrix
correct to machine accuracy in IEEE double precision arithmetic,
and we show that this Pad\'e approximant can be more efficiently
evaluated than a corresponding Taylor series approximation. We
also provide error analysis to bound the propagation of errors in
the double angle recurrence. Numerical experiments show that our
algorithm is competitive in accuracy with the Schur--Parlett
method of Davies and Higham, which is designed for general matrix
functions, and it is substantially less expensive than that method
for matrices of $\infty$-norm of order 1. The dominant
computational kernels in the algorithm are matrix multiplication
and solution of a linear system with multiple right-hand sides, so
the algorithm is well suited to modern computer architectures.",
keywords = "matrix function, matrix cosine, matrix exponential, Taylor series,
Pad\'e approximation, double angle formula, rounding error
analysis, Schur--Parlett method, MATLAB"
}

@techreport{Higham:2002:MCT,
author = "Nicholas J. Higham",
title = "The {Matrix Computation Toolbox} for {MATLAB} (Version 1.0)",
institution =  inst-MCCM,
type = t-NAREP,
number = "No. 410",
pages = "19",
month = aug,
year = "2002",
URL =          "http://www.maths.man.ac.uk/~nareports/narep410.pdf",
abstract = "The Matrix Computation Toolbox is a collection of MATLAB M-files
containing functions for constructing test matrices,
computing matrix factorizations, visualizing matrices,
and carrying out direct search optimization.",
keywords = "test matrix, MATLAB, matrix factorizations, Gaussian elimination,
complete pivoting, rook pivoting, symmetric indefinite matrix,
Gram-Schmidt orthogonalization, polar decomposition, matrix sign
function, complete orthogonal decomposition,
generalized QR factorization, Gauss-Jordan elimination,
Gershgorin disks, pseudospectrum, field of values, visualization,
direct search optimization, Strassen's method, vec-permutation
matrix, Vandermonde matrix"
}

@techreport{Tisseur:2002:TDRSIP,
author = "Fran\c{c}oise Tisseur",
title = "Tridiagonal-Diagonal Reduction of Symmetric Indefinite Pairs",
institution =  inst-MCCM,
type = t-NAREP,
number = "No. 409",
pages = "20",
month = sep,
year = "2002",
URL =          "http://www.maths.man.ac.uk/~nareports/narep409.pdf",
abstract = "We consider the reduction of a symmetric indefinite
matrix pair $(A,B)$, with $B$ nonsingular, to
tridiagonal-diagonal form by congruence transformations.
This is an important reduction in solving polynomial
eigenvalue problems with symmetric coefficient matrices
and in frequency response computations.  First the pair
is reduced to symmetric-diagonal form.  Then we describe
three methods reducing the symmetric-diagonal pair to
tridiagonal-diagonal form.  Two of them employ more
stable versions of Brebner and Grad's pseudosymmetric
Givens and pseudosymmetric Householder reductions, while
the third is new and based on a combination of
Householder reflectors and hyperbolic rotations.  We
prove an optimality condition for the transformations
used in the third reduction.  We present numerical
experiments that compare the different approaches and
show improvements over Brebner and Grad's reductions.",
keywords = "symmetric indefinite generalized eigenvalue problem,
tridiagonalization, hyperbolic rotation, unified
rotation, hyperbolic Householder reflector."
}

@techreport{Higham:2002:JOM,
author = "Nicholas J. Higham",
title = "{$J$}-Orthogonal Matrices: {Properties} and Generation",
institution =  inst-MCCM,
type = t-NAREP,
number = "No. 408",
pages = "15",
month = sep,
year = 2002,
URL =          "http://www.maths.man.ac.uk/~nareports/narep408.pdf",
abstract = "A real, square matrix $Q$ is $J$-orthogonal if $Q^TJQ = J$, where
the signature matrix $J = \diag(\pm 1)$. $J$-orthogonal matrices
arise in the analysis and numerical solution of various matrix
problems involving indefinite inner products, including, in
particular, the downdating of Cholesky factorizations. We present
techniques and tools useful in the analysis, application and
construction of these matrices, giving a self-contained treatment
that provides new insights. First, we define and explore the
properties of the exchange operator, which maps $J$-orthogonal
matrices to orthogonal matrices and vice-versa. Then we show how
the exchange operator can be used to obtain a hyperbolic CS
decomposition of a $J$-orthogonal matrix directly from the usual
CS decomposition of an orthogonal matrix. We employ the
decomposition to derive an algorithm for constructing random
$J$-orthogonal matrices with specified norm and condition number.
We also give a short proof of the fact that $J$-orthogonal
matrices are optimally scaled under two-sided diagonal scalings.
We introduce the indefinite polar decomposition and investigate
two iterations for computing the $J$-orthogonal polar factor: a
Newton iteration involving only matrix inversion and a Schulz
iteration involving only matrix multiplication. We show that these
iterations can be used to $J$-orthogonalize a matrix that is not
too far from being $J$-orthogonal.",
keywords = "$J$-orthogonal matrix, exchange operator, gyration operator,
sweep operator, principal pivot transform, hyperbolic CS
decomposition, two-sided scaling, indefinite least squares
problem, hyperbolic QR factorization, Newton's method, Schulz
iteration"
}

@techreport{Tisseur:2002:STTSM,
author = "Seamus. D. Garvey and Fran\c{c}oise Tisseur
and Mike I. Friswell and John E. T. Penny and
Uwe Prells",
title = "Simultaneous Tridiagonalization of Two Symmetric Matrices",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 407",
pages = "18",
month = aug,
year = "2002",
URL =          "http://www.maths.man.ac.uk/~nareports/narep407.pdf",
abstract = "We show how to simultaneously reduce a pair of symmetric
matrices to tridiagonal form by congruence
transformations.  No assumptions are made on the
nonsingularity or definiteness of the two matrices.  The
reduction follows a strategy similar to the one used for
the tridiagonalization of a single symmetric matrix via
Householder reflectors.  Two algorithms are proposed,
one using non-orthogonal rank-one modifications of the
identity matrix and the other, more costly but more
stable, using a combination of Householder reflectors
and non-orthogonal rank-one modifications of the
identity matrix with minimal condition numbers.  Each of
these tridiagonalization processes requires $O(n^3)$
arithmetic operations and respects the symmetry of the
problem.  We illustrate and compare the two algorithms
with some numerical experiments.",
keywords = "Symmetric matrices,
generalized eigenvalue problem,
tridiagonalization,
}

@techreport{Mihajlovic:2002:EPS,
author = "Milan Mihajlovic and David Silvester",
title = "Efficient Parallel Solvers for the Biharmonic Equation",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 406",
month = jul,
year = 2002,
note = "Submitted to Parallel Computing",
URL =          "http://www.maths.man.ac.uk/~nareports/narep406.ps.gz",
abstract = "We examine the convergence characteristics and performance of
parallelised Krylov subspace solvers applied to the linear algebraic
systems that arise from low-order mixed finite element approximation of
the biharmonic problem. Our strategy results in preconditioned systems
that have nearly optimal eigenvalue distribution, which consists of a
tightly clustered set together with a small number of outliers. We
implement the preconditioner operator in a black-box'' fashion using
publicly available parallelised sparse direct solvers and multigrid
solvers for the discrete Dirichlet Laplacian. We present convergence and
timing results that demonstrate the efficiency and scalability of our
strategy when implemented on contemporary computer architectures."
}

@techreport{Davies:2002:USV,
author = "Philip I. Davies and Matthew I. Smith",
title = "Updating the Singular Value Decomposition",
institution =  inst-MCCM,
type = t-NAREP,
number = "No. 405",
pages = "30",
month = aug,
year = "2002",
URL =          "http://www.maths.man.ac.uk/~nareports/narep405.ps.gz",
abstract = "The spectral decomposition of a symmetric matrix $A$
with small off-diagonal and distinct diagonal elements
can be approximated using a direct scheme of R. Davies and Modi
[\emph{Linear Algebra and Appl.,} 77:61-74, 1986].
In this paper a generalization of this method
for computing the singular value decomposition of
close-to-diagonal $A\in\R^{m \times n}$ is presented.
When $A$ has repeated or close'' singular values it is
possible to apply the direct method to split the problem
in two with one part containing the well-separated singular
values and one requiring the computation of the close''
singular values."
}

@techreport{Davies:2002:SPA,
author = "Philip I. Davies and Nicholas J. Higham",
title = "A {Schur--Parlett} Algorithm for Computing Matrix Functions",
institution =  inst-MCCM,
type = t-NAREP,
number = "No. 404",
pages = "22",
month = jul,
year = "2002",
URL =          "http://www.maths.man.ac.uk/~nareports/narep404.pdf",
abstract = "An algorithm for computing matrix functions is presented. It
employs a Schur decomposition with reordering and blocking
followed by the block form of a recurrence of Parlett, with
functions of the nontrivial diagonal blocks evaluated via a Taylor
series. A parameter is used to balance the conflicting
requirements of producing small diagonal blocks and keeping the
separations of the blocks large. The algorithm is intended
primarily for functions having a Taylor series with an infinite
functions, such as the logarithm. Novel features introduced here
include a convergence test that avoids premature termination of
the Taylor series evaluation and an algorithm for reordering and
blocking the Schur form. Numerical experiments show that the
algorithm is competitive with existing special-purpose algorithms
for the matrix exponential, logarithm, and cosine. Nevertheless,
the algorithm can be numerically unstable with the default choice
of its blocking parameter (or in certain cases for all choices),
and we explain why determining the optimal parameter appears to be
a very difficult problem. A MATLAB implementation is available
that is much more reliable than the existing MATLAB function
\texttt{funm}."
}

@techreport{Song:PTD:2002,
author = "Y. Song and C. T. H. Baker",
title = "Perturbation Theory for Discrete {Volterra} Equations",
institution =  inst-MCCM,
type = t-NAREP,
number = "No. 403",
pages = "18",
month = jun,
year = "2002",
URL =          "http://www.maths.man.ac.uk/~nareports/narep403.ps.gz",
abstract = "Fixed  point theory is used to investigate nonlinear discrete
Volterra equations that are perturbed versions of linear equations.
Sufficient conditions are established ($i$) to ensure that stability
(in a sense that is defined) of the solutions of the linear equation
implies a corresponding stability of the zero solution of the
nonlinear equation and ($ii$) to ensure the existence of
asymptotically periodic solutions",
keywords = "Fixed point theorem, perturbation, discrete Volterra equations,
periodic, asymptotically periodic."
}

@techreport{Song:SDV:2002,
author = "Y. Song and C. T. H. Baker",
title = "Stability in Discrete {Volterra} Equations",
institution =  inst-MCCM,
type = t-NAREP,
number = "No. 402",
pages = "19",
month = jun,
year = "2002",
URL =         "ftp://ftp.ma.man.ac.uk/pub/narep/narep402.ps.gz",
abstract = "A linearized stability theory for discrete Volterra
equations of the type $x(n)=f(n)+\sum^n_{j=0}K\big(n, j, x(n)\big), \mbox{(n\ge 0)}$ is developed: Various types of stability are defined and
studied by using the fundamental and resolvent matrices that
are associated with linear equations of the form $x(n)=f(n) + \sum^n_{j=0}B(n, j)x(j), \mbox{(n\ge 0)},$ and several
necessary and sufficient conditions for stability are obtained
for solutions of this linear equation. In conclusion,  a result
concerning the stability of a solution of the nonlinear
equation is provided under certain assumptions.
The results parallel corresponding results for classical
Volterra integral equations of the second kind",
keywords = "Stability,  asymptotic stability,  discrete Volterra equations,
resolvent matrix,  fundamental matrix."
}

@techreport{MCCM:2002:AR1,
author =       "Manchester Centre for Computational Mathematics",
title =        "Annual Report: {January--December} 2001",
type =         t-NAREP,
number =       "No. 401",
institution =  inst-MCCM,
pages =        "21",
month =        may,
year =         "2002",
URL =          "http://www.maths.man.ac.uk/~nareports/narep401.pdf",
}

@techreport{Baker:2002:PSD,
author = "C. T. H. Baker and Y. Song",
title = "Periodic Solutions of Discrete {Volterra} Equations",
institution =  inst-MCCM,
type = t-NAREP,
number = "No. 400",
pages = "21",
month = may,
year = "2002",
URL =          "ftp://ftp.ma.man.ac.uk/pub/narep/narep400.ps.gz",
abstract =
"In this paper we investigate periodic solutions of linear
and nonlinear discrete Volterra equations of convolution,
or non-convolution type with unbounded memory. For linear
discrete Volterra equations,  we establish  Fredholm's
alternative theorem for convolution equations and we show
that, under certain conditions, only one particular bounded
initial function can generate a  unique periodic solution
and this unique periodic solution attracts all other
solutions with bounded initial function.  In other words,
we show that all solutions of linear discrete Volterra
equations are asymptotically periodic under certain conditions.",
keywords =
"Periodic, asymptotically periodic, discrete Volterra
equations, resolvent matrices, Fredholm's alternative."
}

@techreport{Powell:2002:OPR,
author = "Catherine Powell and David Silvester",
title = "Optimal Preconditioning for {Raviart-Thomas} Mixed Formulation
of Second-Order Elliptic Problems",
institution =  inst-MCCM,
type = t-NAREP,
number = "No. 399",
month = mar,
year = 2002,
note = "Submitted to the SIAM Journal on Matrix Analysis and Applications",
URL =          "http://www.maths.man.ac.uk/~nareports/narep399.ps.gz",
abstract = "We evaluate two preconditioning strategies for the indefinite
linear system obtained from Raviart-Thomas mixed finite element
formulation of a second-order elliptic problem with variable
coefficients. In contrast to other approaches, our emphasis is on
linear algebra; we establish the optimality of the
preconditioning using basic properties of the finite element
matrices. The underlying saddle-point problem is well-posed in
two function spaces, $H(div)\times L^{2}$ and $L^{2} \times H^{1}$, leading to the possibility of two distinct types of
preconditioner. For homogeneous Dirichlet boundary conditions,
the discrete problems are identical. This motivates the use of
Raviart-Thomas approximation in both frameworks yielding a
non-conforming method in the second case. We construct two
block-diagonal preconditioners and establish a new uniform
eigenvalue bound. Trials of preconditioned MINRES illustrate that
both preconditioners are robust with respect to the mesh
parameter, and with respect to jumps and anisotropies in the
diffusion coefficients.  "
}

@techreport{Baker:2002:DVO,
author = "C. T. H. Baker and Y. Song",
title = "Discrete Volterra Operators, Fixed Point Theorems \& their
Application",
institution =  inst-MCCM,
type = t-NAREP,
number = "No. 398",
pages = "23",
month = mar,
year = "2002",
URL =          "ftp://ftp.ma.man.ac.uk/pub/narep/narep398.ps.gz",
abstract =

"The existence of solutions of nonlinear discrete Volterra
equations is established.
We define discrete Volterra operators on normed spaces of
infinite sequences of finite-dimensional vectors.  and present
some of their basic properties (continuity, boundedness, and
representation). The treatment relies upon the use of
coordinate functions, and the existence results are obtained
using fixed point theorems for discrete Volterra operators on
infinite-dimensional spaces based on fixed point theorems of
Schauder, Rothe, and Altman, and Banach's contraction mapping
theorem, for finite-dimensional spaces.",
keywords =
"Fixed point theorem, discrete Volterra equations, discrete
Volterra operators."
}

@techreport{Bojanczyk:2002:SIL,
author = "Adam Bojanczyk and Nicholas J. Higham and Harikrishna Patel",
title = "Solving the Indefinite Least Squares Problem by Hyperbolic {QR}
Factorization",
institution =  inst-MCCM,
type = t-NAREP,
number = "No. 397",
pages = "19",
month = jan,
year = "2002",
URL =          "ftp://ftp.ma.man.ac.uk/pub/narep/narep397.pdf",
abstract = "The indefinite least squares (ILS) problem is to minimize a
certain type of indefinite quadratic form. We develop
perturbation theory for the problem and identify a condition
number. We describe and analyze a method for solving the ILS
problem based on hyperbolic QR factorization. This method is
faster than one recently proposed by Chandrasekaran, Gu and Sayed
that employs both QR and Cholesky factorizations. We give a
rounding error analysis of the new method and use the
perturbation theory to show that under a reasonable assumption
the method is forward stable. Our analysis is quite general and
sheds some light on the stability properties of hyperbolic
transformations. Our numerical experiments suggest that the new
method is in practice just as accurate as the method of
Chandrasekaran, Gu and Sayed.",
keywords = "indefinite least squares problem, downdating, hyperbolic
rotation, hyperbolic QR factorization, rounding error analysis,
forward stability, perturbation theory, condition number"
}

@techreport{Baker:2002:BSD,
author = "C. T. H. Baker and J. M. Ford and N. J. Ford",
title = "Bifurcations in Stochastic Delay Differential Equations---{I}:
{Numerical} Investigation ",
institution =  inst-MCCM,
type = t-NAREP,
number = "No. 396",
pages = "18",
month = mar,
year = "2002",
URL =          "ftp://ftp.ma.man.ac.uk/pub/narep/narep396.ps.gz",
abstract =
"We consider {\em stochastic delay differential equations} of
the form $dY(t)=F(t,Y(t),Y(t-\tau))dt + G(t,Y(t),Y(t-\tau)) dW(t) \quad (t \geq t_0),$ interpreted in the It\^{o} sense
where $W(t)$ is a standard {\em Wiener process}, with $Y(t) = \Phi(t)$ for $t \in [-\tau,0]$ (where $\tau > 0$ is the
constant lag, or time-lag').  We are interested in {\em
bifurcations} (that is, changes in the qualitative behaviour
of the solutions to these equations) and we draw on insights
from the deterministic counterpart $y'(t)=f(t,y(t),y(t-\tau)), \quad \mbox{ for t \geq t_0},$
with $y(t) = \phi(t)$ for $t \in [t_0-\tau,t_0]$ and constant
lag $\tau$.  For these deterministic problems we can exploit
known theory and numerical results that enable us to predict
where changes occur in the behaviour of the solutions (exact
and approximate) of the equation.
For those unfamiliar with the area, rather diverse components
of the mathematical background are necessary to understand the
questions of interest.  In this introductory report we first
review some deterministic results and some basic elements of
the stochastic analysis that ($i$) suggests lines of
investigation for the stochastic case and ($ii$) are expected
to facilitate the theoretical investigation of the stochastic
problem. We then give numerical evidence that the presence of
noise, corresponding to the change from the deterministic to
the stochastic form, may change the qualitative behaviour of
the solutions in an interesting way.",
keywords =
"Bifurcation and stability, stochastic and deterministic delay
equations with fixed lag, Euler-type methods with fixed step,
numerical experiments."
}

@techreport{Tisseur:2001:CBE,
author = "Fran\c{c}oise Tisseur",
title = "A Chart of Backward Errors and Condition Numbers for
Singly and Doubly Structured Eigenvalue Problems",
institution =  inst-MCCM,
type = t-NAREP,
number = "No. 395",
pages = "23",
month = dec,
year = "2001",
URL =          "ftp://ftp.ma.man.ac.uk/pub/narep/narep395.ps.gz",
abstract = "We present a chart of structured backward errors for
approximate eigenpairs of singly and doubly structured
eigenvalue problems.  We aim to give, wherever possible,
formulae that are inexpensive to compute so that they
can be used routinely in practice.  We identify a number
of problems for which the structured backward error is
within a factor $\sqrt{2}$ of the unstructured backward
error. This paper collects, unifies and extends existing
work on this subject.",
}

@techreport{Tisseur:2001:IGT,
author = "Jean-Pierre Dedieu and Myong-Hi Kim and Michael Shub
and Fran\c{c}oise Tisseur",
title = "Implicit Gamma Theorems {(I)}: {Pseudoroots} and Pseudospectra",
institution =  inst-MCCM,
type = t-NAREP,
number = "No. 394",
pages = "25",
month = Nov,
year = "2001",
URL =          "ftp://ftp.ma.man.ac.uk/pub/narep/narep394.ps.gz",
abstract = "Let $g: \E \rightarrow \F$ be an analytic function
between two Hilbert spaces $\E$ and $\F$. We study the
set $g(B(x,\eps ))\subset \F$, the image under $g$ of
the closed ball about $x\in\E$ with radius $\eps$. When
$g(x)$ expresses the solution of an equation depending
on $x$ then the elements of $g(B(x,\eps ))$ are
$\eps$-pseudosolutions . Our aim is to investigate the
size of the set $g(B(x,\eps ))$. We derive upper and
lower bounds of the following form $$g(x) + Dg(x)\left ( B(0, c_1 \eps ) \right ) \subseteq g(B(x,\eps )) \subseteq g(x) + Dg(x)\left ( B(0, c_2 \eps ) \right ),$$ where $Dg(x)$ denotes the derivative of $g$ at
$x$. We consider both the case where $g$ is given
explicitly and the case where $g$ is given
implicitly. We apply our results to the implicit
function associated with the evaluation map, namely the
solution map, and to the polynomial eigenvalue
problem. Our results are stated in terms of an invariant
$\gamma$ which has been extensively used by various
authors in the study of Newton's method. The main tool
used here is an implicit $\gamma$ theorem, which
estimates the $\gamma$ of an implicit function in terms
of the $\gamma$ of the function defining it.",
}

@techreport{Cheng:2001:ILBA,
author = "Sheung Hun Cheng and Nicholas J. Higham",
title = "Implementation for {LAPACK} of a Block Algorithm for Matrix 1-Norm
Estimation",
institution =  inst-MCCM,
type = t-NAREP,
number = "No. 393",
pages = "19",
month = aug,
year = "2001",
note = "LAPACK Working Note 152",
URL =          "ftp://ftp.ma.man.ac.uk/pub/narep/narep393.ps.gz",
abstract = "We describe double precision and complex*16 Fortran~77
implementations, in LAPACK style, of a block matrix 1-norm
estimator of Higham and Tisseur. This estimator differs from that
underlying the existing LAPACK code, \texttt{xLACON}, in that it
iterates with a matrix with $t$ columns, where $t\ge 1$ is a
parameter, rather than with a vector, and so the basic
computational kernel is level 3 BLAS operations. Our experiments
with random matrices on a Sun SPARCStation Ultra-5 show that with
$t=2$ or $4$ the new code offers better estimates than
\texttt{xLACON} with a similar execution time. Moreover, with
$t>2$, estimates exact over $95\%$ and $75\%$ of the time are
achieved for the real and complex version respectively, with
execution time growing much slower than $t$. We recommend this
new code be included as an auxiliary routine in LAPACK to
complement the existing LAPACK routine \texttt{xLACON}, upon
which the various drivers should still be based for compatibility
reasons."
}

@techreport{Smith:2001:ASAC,
author =       "Matthew I. Smith",
title =        "A {Schur} Algorithm for Computing Matrix $p$th Roots",
institution =   inst-MCCM,
type =          t-NAREP,
number =       "No. 392",
month =         jul,
year =         "2001",
pages =        "18",
URL =          "ftp://ftp.ma.man.ac.uk/pub/narep/narep392.ps.gz",
abstract =     "Any nonsingular matrix has $p$th roots,
and they can be classified into two groups, those that are
functions of $A$ and those that are not. Matrix $p$th roots can be
computed by a specialized version of Newton's method, but we show
that this iteration has poor convergence and stability properties
in general. We present a Schur algorithm for computing a matrix
$p$th root that generalizes methods of Bj{\"{o}}rck and Hammarling
[\emph{Linear Algebra and Appl.,} 52/53:127-140, 1983] and Higham
[\emph{Linear Algebra and Appl.,} 88/89:405-430, 1987] for the
square root. Our algorithm forms a Schur decomposition of $A$ and
computes a $p$th root of the (quasi)triangular factor by a
recursion. The backward error associated with our Schur method is
examined, and we find that the method has excellent numerical
stability.",
keywords = "matrix $p$th root, Schur algorithm, Newton's method,
commutativity, stability, real arithmetic, rounding error
analysis"
}

@techreport{MCCM:2001:AR0,
author =       "Manchester Centre for Computational Mathematics",
title =        "Annual Report: {January--December} 2000",
type = t-NAREP,
number =       "No. 391",
institution =  inst-MCCM,
pages =        "23",
month =        mar,
year =         "2001",
URL =          "http://www.maths.man.ac.uk/~nareports/narep391.pdf",
}

@techreport{FordJ:2001:PFWT,
author = "J. M. Ford and K. Chen and N. J. Ford",
title = "Parallel Implementation of Fast Wavelet Transforms",
institution =  inst-MCCM,
type = t-NAREP,
number = "No. 389",
pages = "15",
month = feb,
year = "2001",
URL =          "ftp://ftp.ma.man.ac.uk/pub/narep/narep389.ps.gz",
abstract = "We look at the use of parallel techniques for fast wavelet transforms, focusing
particularly on how to make best use of any given number of available  processors.
Efficient algorithms for parallel  FWT have been established in the case where the
lengths of the vectors to be transformed and the number of processors are both powers of 2.
We show that a na{\"}{\i}ve extension to arbitrary vector length and number of processors
may give poor performance. We present new algorithms for the distribution of matrix
elements between processors that overcome these difficulties. We have implemented our
algorithms on a multi-processor Sun Sparc machine, and we present analytical and
experimental results demonstrating significant improvements in performance of FWT
calculation.",
keywords = "parallel, FWT, implementation, MPI"
}

@techreport{FordJ:2001:SPFW,
author = "J. M. Ford and K. Chen and N. J. Ford",
title = "Small-Scale Parallel Implementation of Fast Wavelet Transforms",
institution =  inst-MCCM,
type = t-NAREP,
number = "No. 388",
pages = "7",
month = feb,
year = "2001",
URL =          "ftp://ftp.ma.man.ac.uk/pub/narep/narep388.ps.gz",
abstract = "We look at the use of parallel techniques for Fast Wavelet Transforms, focusing
particularly on how to make best use of a very limited (and, perhaps, variable)
number of available  processors. Efficient algorithms for parallel  FWT have been
established in the case where the lengths of the vectors to be transformed and the
number of processors are both powers of 2. We show that a na{\"}{\i}ve extension to
arbitrary vector length and number of processors may give poor performance or, in
some cases, be impossible. We present new algorithms for the distribution of matrix
elements between processors that overcome these difficulties. We have implemented our
algorithms on a local area network of inexpensive PCs using MPI, and present
experimental results demonstrating significant improvements in performance of FWT
calculation.",
keywords = "PC, parallel, FWT, MPI"
}

@techreport{Ford:2001:NAFDE,
author = "N. J. Ford and A. C. Simpson",
title = "Numerical and Analytical Treatment of Differential Equations of
Fractional Order",
institution =  inst-MCCM,
type = t-NAREP,
number = "No. 387",
pages = "7",
month = feb,
year = "2001",
URL =          "ftp://ftp.ma.man.ac.uk/pub/narep/narep387.ps.gz",
abstract = "We begin with a brief survey of differential equations of fractional order and
their applications. We focus on key issues in the exact and numerical solution
of the equations. We summarise our recent results: We describe the importance
of the appropriate choice of	initial conditions. We show how the behaviour of
solutions changes with small changes in the fractional order. We discuss the use
of {\em systems} of fractional order equations as a tool to solve equations of
higher order. We give results that show how the computational effort required to
solve an equation of fractional order can be reduced, either through careful
application of the  {\em short memory}  principle or through the use of schemes
of varying step-length.",
keywords = "short memory, systems, fractional order"
}

@techreport{Ford:2001:AFDG1,
author = "N. J. Ford and A. C. Simpson",
title = "The Approximate Solution of Fractional Differential Equations of Order
Greater Than 1",
institution =  inst-MCCM,
type = t-NAREP,
number = "No. 386",
pages = "6",
month = feb,
year = "2001",
URL =          "ftp://ftp.ma.man.ac.uk/pub/narep/narep386.ps.gz",
abstract = "Differential equations of fractional order arise, for example, in the generalised
viscoelastic constitutive equations, the generalised Basset force, and
in fractional Brownian motion. Numerical methods for solution of
equations of fixed order $\alpha$ are well established, particularly
for $0<\alpha<1$, however many of the analytical results depend on the
assumption that the equation has a particularly simple form. Recent
work (by Diethelm and Ford and others) has considered approaches for
solving more complicated and multi-term equations. In this paper we
present an alternative approach to numerical solution of order greater
than one and we conclude the paper by comparing the numerical results
for our proposed method with known numerical and analytical solutions
for prototype nonlinear problems.",
keywords = "systems, fractional order"
}

@techreport{Ford:2001:SVA,
author = "N. J. Ford and A. C. Simpson",
title = "The Numerical Solution of Fractional Differential Equations: Speed
Versus Accuracy",
institution =  inst-MCCM,
type = t-NAREP,
number = "No. 385",
pages = "12",
month = feb,
year = "2001",
URL =          "ftp://ftp.ma.man.ac.uk/pub/narep/narep385.ps.gz",
abstract = "This paper is concerned with the development of efficient algorithms for the
approximate solution of fractional differential equations of the
form $D^\alpha y(t)=f(t,y(t)), \mathbb{R}^+-\mathbb{N}$. We
briefly review standard numerical techniques for the solution of
($\dagger$) and we consider how the computational cost may be
reduced by taking into account the structure of the calculations
to be undertaken. We analyse the {\it fixed memory principle} and
present an alternative {\it nested mesh} variant that gives a good
approximation to the true solution at reasonable computational
cost. We conclude with some numerical examples.",
keywords = "short memory, finite memory, logarithmic memory"
}

@techreport{Edwards:2001:BSDE,
author = "J. T. Edwards and N. J. Ford",
title = "Boundedness and Stability of Difference Equations",
institution =  inst-MCCM,
type = t-NAREP,
number = "No. 384",
pages = "14",
month = feb,
year = "2001",
URL =          "ftp://ftp.ma.man.ac.uk/pub/narep/narep384.ps.gz",
abstract = "This paper is concerned with the qualitative behaviour of solutions to
difference equations.	We focus on boundedness and stability of solutions
and we present a unified theory that applies both to autonomous and
non-autonomous equations and to nonlinear equations as well as linear
equations. Our presentation brings together new, established, and
hard-to-find results from the literature and provides a theory that is both
memorable and easy to apply. We show how the theoretical results given here
relate to some of those in the established literature and by means of simple
examples we indicate how the use of Lipschitz constants in this way can provide
useful insights into the qualitative behaviour of solutions to some nonlinear
problems including those arising in numerical analysis.",
keywords = "boundedness, stability, Lipschitz"
}

@techreport{Ford:2001:SHQ,
author = "N. J. Ford and P. D. Crofts and R. H. T. Edwards",
title = "Sensitivity of Hospital Clinic Queues to Patient Non-Attendance",
institution =  inst-MCCM,
type = t-NAREP,
number = "383",
month = feb,
pages = "6",
year = "2001",
URL =          "ftp://ftp.ma.man.ac.uk/pub/narep/narep383.ps.gz",
abstract = "The recently introduced Patients' Charter offers patients in the United
Kingdom a guaranteed minimum quality of service when they require hospital
treatment. For example, the following guarantees are offered: 1. No
patient should be required to spend more than one year	awaiting an offer
of an appointment with a Consultant 2. Patients should not be expected to
wait in clinic for more than 30 minutes without seeing a doctor. Health
economists are active in devising strategies and systems to deal with these
enhanced patient expectations. In this paper, we describe the work of an
interdisciplinary team, funded by the UK National Health Service Research
and Development which investigated computer simulation of clinic queues to
provide further insight into the problem and to produce computer software
that can assist in clinic optimisation.",
keywords = "optimisation, modelling, queues"
}

@techreport{Ford:2001:NDSS,
author = "Neville J Ford and Sjoerd M Verduyn Lunel",
title = "Numerical approximation of delay differential equations with small
solutions",
institution =  inst-MCCM,
type = t-NAREP,
number = "382",
month = feb,
pages = "6",
year = "2001",
URL =          "ftp://ftp.ma.man.ac.uk/pub/narep/narep382.ps.gz",
abstract = "We are concerned with the detection of small solutions of delay differential equations.
We consider, as prototype, an equation of the form $\dot{y}(t)=b(t)y(t-\tau)$ for
which conditions on the function $b$ that govern the existence (or otherwise) of small
solutions are well-known. We apply simple numerical methods to the problem and
investigate the characteristic values of the solution operator. We show how these
characteristic roots may be used to detect the existence of small solutions.",
keywords = "small solutions, floquet theory, characteristic values"
}

@techreport{Ford:2001:CSS,
author = "Neville J. Ford and Sjoerd M. Verduyn Lunel",
title = "Characterising Small Solutions in Delay Differential Equations
Through Numerical Approximations",
institution =  inst-MCCM,
type = t-NAREP,
number = "381",
month = feb,
pages = "15",
year = "2001",
URL =          "ftp://ftp.ma.man.ac.uk/pub/narep/narep381.ps.gz",
abstract = "The existence of small solutions, i.e., solutions that decay
faster than any exponential, for linear time-dependent
differential delay equations with bounded coefficients depends on
specific properties of the coefficients. Although small solutions
do not occur in the finite dimensional approximations of the
differential delay equation we show that the existence of small
solutions for differential delay equations can be predicted from
the behaviour of the spectrum of the finite dimensional
approximations.",
keywords = "small solutions, floquet theory"
}

@techreport{Edwards:2001:IDBP,
author = "J. T. Edwards, J. A. Roberts and N. J. Ford",
title = "The numerical Solution of an Integro-Differential Equation Close to
Bifurcation Points",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 380",
pages = "6",
month = feb,
year = "2001",
URL =          "ftp://ftp.ma.man.ac.uk/pub/narep/narep380.ps.gz",
abstract = "In this paper we  investigate the qualitative behaviour of numerical
approximations to a Volterra integro-differential equation. We consider (as
prototype) a linear problem with fading memory kernel of the form
$$y'(t)=-\int^t_0e^{-\lambda(t-s)}y(s)ds, \qquad \qquad y(0)=1$$ and we consider
the performance of simple numerical schemes applied to solve the equation. We
are concerned with the preservation (or otherwise) of qualitative properties of
the analytical solution in the numerical approximation.   We outline the known
stability behaviour and derive the values of $\lambda$ at which the true solution
bifurcates. We give the corresponding analysis for the discrete schemes and we show
that, as the step size of the numerical scheme decreases, the bifurcation points
tend towards those of the continuous scheme. We illustrate our results with some
numerical examples.",
keywords = "bifurcations, numerical approximations, integro-differential equations"
}

@techreport{Ford:2001:FESO,
author = "N. J. Ford and K. Diethelm",
title = "The numerical Solution of Linear and Non-Linear Fractional
Differential Equations Involving Fractional Derivatives of
several orders",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 379",
month = feb,
pages = 14,
year = 2001,
URL =          "http://www.maths.man.ac.uk/~nareports/narep379.ps.gz",
abstract = "We consider fractional differential equations of the very general
form $y^{(\alpha)}(t)= f(t,y(t),y^{(\beta_1)}(t),y^{(\beta_2)}(t), \ldots,y^{(\beta_n)}(t))$ with $\alpha>\beta_1>\beta_2>\ldots>\beta_n$
and $\alpha-\beta_1\le1$,$\beta_j-\beta_{j+1}\le1$, $0\le\beta_n\le1$,
and where the derivatives are understood in Caputo's sense.For equations of
this type, we discuss the questions of existence and uniqueness of solutions,
and we investigate how the solutions depend on the given data. We propose
convergent and stable numerical methods based on a nearly equivalent system
of fractional differential equations of order not exceeding $\beta_n$. We
give particular emphasis to the practically important linear case and also
present results for nonlinear problems.",
keywords = "milti-term equation, fractional order, fractional Adams method"
}

@techreport{Ford:2001:NBT,
author = "N. J. Ford and K. Diethelm",
title = "Numerical Solution of the {Bagley Torvik} Equation",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 378",
month = feb,
pages = 12,
year = 2001,
note = "Submitted to BIT",
URL = "http://www.maths.man.ac.uk/~nareports/narep378.ps.gz",
abstract = "We consider the numerical solution of the Bagley-Torvik
equation $$A y''(t) + B D_*^{3/2}y(t) + C y(t) = f(t)$$
as a prototype fractional differential equation with two
derivatives. Approximate solutions have recently been proposed in
the book and papers of Podlubny in which the solution obtained with
approximate methods is compared to the exact solution. In this paper
we consider the reformulation of the Bagley-Torvik equation as a
system of fractional differential equations of order $1/2$. This
allows us to propose numerical methods for its solution which are
consistent and stable and have arbitrarily high order.",
keywords = "Bagley Torvik, multi-term, fractional order"
}

@techreport{Ford:2001:AFDE,
author = "N. J. Ford and K. Diethelm",
title = "Analysis of Fractional Differential Equations",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 377",
month = feb,
pages = 17,
year = 2001,
note = "To appear in Journal of Mathematical Analysis and Applications",
URL = "http://www.maths.man.ac.uk/~nareports/narep377.ps.gz",
abstract = "We discuss existence, uniqueness and {\it structural stability}
of solutions of nonlinear differential equations of fractional
order. The differential operators are taken in the
Riemann-Liouville sense and the initial conditions are specified
according to Caputo's suggestion, thus allowing for interpretation
in a physically meaningful way. We investigate in particular the
dependence of the solution on the order of the differential
equation and on the initial condition and we relate our results
to the selection of appropriate numerical schemes for the solution
of fractional differential equations.",
keywords = "fractional order, stability, existence, uniqueness, numerical schemes"
}

@techreport{Bolton:2001:MCLS,
author = "P. Bolton and J. Stratakis and R. W. Thatcher",
title = "Mass Conservation in Least Squares Methods for {Stokes} Flow",
institution =  inst-MCCM,
type = t-NAREP,
number = "No. 376",
pages = "22",
month = jun,
year = "2001",
URL =          "ftp://ftp.ma.man.ac.uk/pub/narep/narep376.pdf"
}

@techreport{Basebi:2001:ASMM,
author =      "Thebe Basebi and Ruth M. Thomas",
title =       "A Study of Moving Mesh Methods Applied to a Thin Flame
Propagating in a Detonator Delay Element",
type = t-NAREP,
number =      "No. 375",
institution =  inst-MCCM,
month =        oct,
year =        "2001",
URL =         "http://www.maths.man.ac.uk/~nareports/narep375.ps.gz",
abstract =    "We study the application of moving mesh methods to a
one-dimensional (time dependent) detonator delay element
problem. We consider moving mesh methods based on the
equidistribution principle derived by
Huang $et$ $al.$ (1994). Adaptive mesh methods have been widely
used recently to solve time dependent partial differential
equations having large solution gradients. Significant
improvements in accuracy and efficiency are achieved by
adapting the nodes (mesh points) so that they are
concentrated about areas of large solution variations.
Each system of equations for the moving mesh methods is
solved in conjunction with the detonator problem. In this
paper, the system of ordinary differential equations that
results (after discretising in space) is solved using the
double precision version of the stiff ordinary differential
equation solver DASSL. The numerical results clearly
demonstrate that the moving mesh methods are capable of
tracking the deflagration wave as it travels down the
detonator delay element more accurately and more efficiently
than a fixed mesh method.",
keywords =    "moving mesh methods, detonator delay element, time dependent
partial differential equations, scaling, monitor functions."
}

@techreport{Cheng:2001:PIB,
author = "Sheung Hun Cheng and Nicholas J. Higham",
title = "Parallel Implementation of a Block Algorithm for Matrix 1-Norm
Estimation",
institution =  inst-MCCM,
type = t-NAREP,
number = "No. 374",
pages = "9",
month = feb,
year = "2001",
URL =          "ftp://ftp.ma.man.ac.uk/pub/narep/narep374.pdf",
abstract = "We describe a parallel Fortran~77 implementation, in ScaLAPACK
style, of a block matrix 1-norm estimator of Higham and Tisseur.
This estimator differs from that underlying the existing
ScaLAPACK code, \texttt{PxLACON}, in that it iterates with a
matrix with $t$ columns, where $t\ge 1$ is a parameter, rather
than with a vector, and so the basic computational kernel is
level 3 BLAS operations. Our experiments on an SGI Origin2000
show that with $t=2$ or $4$ the new code offers better estimates
than \texttt{PDLACON} with a similar execution time. Moreover,
with $t>4$, estimates exact over $90\%$ of the time are achieved
with execution time growing much slower than $t$.",
keywords = ""
}

@techreport{Higham:2001:DDH,
author = "Nicholas J. Higham and Fran\c{c}oise Tisseur and
Paul M. Van Dooren",
title = "Detecting a Definite {Hermitian} Pair
and a Hyperbolic or Elliptic Quadratic Eigenvalue Problem,
and Associated Nearness Problems",
institution =  inst-MCCM,
type = t-NAREP,
number = "373",
month = feb,
pages = "18",
year = "2001",
URL =          "ftp://ftp.ma.man.ac.uk/pub/narep/narep373.ps.gz",
abstract = "An important class of generalized eigenvalue problems $Ax = \lambda Bx$ is those in which $A$ and $B$ are Hermitian and some
real linear combination of them is definite. For the quadratic
eigenvalue problem (QEP) $(\l^2 A + \l B + C)x = 0$ with
Hermitian $A$, $B$ and $C$ and positive definite $A$, particular
interest focuses on problems in which $(x^*Bx)^2 - 4 (x^*Ax)(x^*Cx)$ is one-signed for all nonzero $x$---for the
positive sign these problems are called hyperbolic and for the
negative sign elliptic. The important class of overdamped
problems arising in mechanics is a sub-class of the hyperbolic
problems. For each of these classes of generalized and quadratic
eigenvalue problems we show how to check that a putative member
has the required properties and we derive the distance to the
nearest problem outside the class. For definite pairs $(A,B)$ the
distance is the Crawford number, and we derive bisection and
level set algorithms both for testing its positivity and for
computing it. Testing hyperbolicity of a QEP is shown to reduce
to testing a related pair for definiteness. The distance to the
nearest non-hyperbolic or non-elliptic $n\times n$ QEP is shown
to be the solution of a global minimization problem with $n-1$
degrees of freedom. Numerical results are given to illustrate the
theory and algorithms.",
keywords = "Generalized eigenvalue problem, definite pair, Crawford number,
quadratic eigenvalue problem, hyperbolic system, elliptic system,
overdamped system, gyroscopic system, bisection algorithm, level
set algorithm, nearness problems"
}

@techreport{Higham:2001:MPP,
author = "Nicholas J. Higham and Fran\c{c}oise Tisseur",
title = "More on Pseudospectra for Polynomial Eigenvalue Problems
and Applications in Control Theory",
institution =  inst-MCCM,
type = t-NAREP,
number = "372",
month = jan,
pages = "21",
year = "2001",
URL =          "ftp://ftp.ma.man.ac.uk/pub/narep/narep372.ps.gz",
abstract = "Definitions and characterizations of pseudospectra are given for
rectangular matrix polynomials expressed in homogeneous form:
$P(\a,\b) = (\a^dA_d + \a^{d-1}\b A_{d-1} + \cdots + \b^d A_0)$.
It is shown that problems with infinite (pseudo)eigenvalues are
elegantly treated in this framework. For such problems
stereographic projection onto the Riemann sphere is shown to
provide a convenient way to visualize pseudospectra. Lower bounds
for the distance to the nearest non-regular polynomial and the
nearest uncontrollable $d$th order system (with equality for
standard state-space systems) are obtained in terms of
pseudospectra, showing that pseudospectra are a fundamental tool
for reasoning about matrix polynomials in areas such as control
theory. How and why to incorporate linear structure into
pseudospectra is also discussed by example.",
keywords = "polynomial eigenvalue problem, $\lambda$-matrix, matrix
polynomial, homogeneous form, pseudospectrum, stereographic
projection, Riemann sphere, nearest non-regular polynomial,
nearest uncontrollable system, structured perturbations"
}

@techreport{Tisseur:2001:BEM,
author = "Nicholas J. Higham and Fran\c{c}oise Tisseur",
title = "Bounds for Eigenvalues of Matrix Polynomials",
institution =  inst-MCCM,
type = t-NAREP,
number = "371",
month = jan,
pages = "16",
year = "2001",
URL =          "ftp://ftp.ma.man.ac.uk/pub/narep/narep371.pdf",
abstract = "Upper and lower bounds are derived for the absolute values of the
eigenvalues of a matrix polynomial (or $\lambda$-matrix). The
bounds are based on norms of the coefficient matrices and involve
the inverses of the leading and trailing coefficient matrices.
They generalize various existing bounds for scalar polynomials
and single matrices. A variety of tools are used in the
derivations, including block companion matrices, Gershgorin's
theorem, the numerical radius, and associated scalar polynomials.
Numerical experiments show that the bounds can be surprisingly
sharp on practical problems.",
keywords = "polynomial eigenvalue problem, $\lambda$-matrix, matrix
polynomial, block companion matrix, Gershgorin's theorem,
}

@techreport{Tisseur:2000:SQEP,
author = "Fran\c{c}oise Tisseur and Karl Meerbergen",
title = "The Quadratic Eigenvalue Problem",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 370",
pages = "49",
month = nov,
year = "2000",
URL =          "ftp://ftp.ma.man.ac.uk/pub/narep/narep370.ps.gz",
abstract = "We survey the quadratic eigenvalue problem,
treating its many applications, its
mathematical properties, and a variety of numerical
solution techniques.  Emphasis is given to exploiting
both the structure of the matrices in the problem
(dense, sparse, real, complex, Hermitian,
skew-Hermitian) and the spectral properties of the
problem. We classify the available choices of methods
and catalogue available software.",
keywords = "quadratic eigenvalue problem, eigenvalue, eigenvector,
$\lambda$-matrix, matrix polynomial, second-order
differential equation, overdamped system, gyroscopic
system, linearization, backward error, pseudospectrum,
condition number, Krylov methods, Arnoldi method, Lanczos
method, Jacobi-Davidson method"
}

@techreport{Higham:2000:CNC,
author = "Nicholas J. Higham",
title = "Computing the Nearest Correlation Matrix---{A} Problem from Finance",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 369",
month = oct,
pages = 14,
year = 2000,
note = "Revised July 2001",
URL =          "http://www.maths.man.ac.uk/~nareports/narep369.pdf",
abstract = "Given a symmetric matrix what is the nearest correlation matrix,
that is, the nearest symmetric positive semidefinite matrix with
unit diagonal? This problem arises in the finance industry, where
the correlations are between stocks. For distance measured in the
Frobenius norm, we characterize the solution using convex
analysis. We show how the modified alternating projections method
can be used to compute the solution. In the finance application
the original matrix has many zero or negative eigenvalues; we
show that the nearest correlation matrix has correspondingly many
zero eigenvalues and that this fact can be exploited in the
computation.",
keywords = "correlation matrix, positive semidefinite matrix, nearness
problem, convex analysis, alternating projections method"
}

@techreport{Paul:2000:DES,
author = "Christopher Paul",
title = "Designing Efficient Software for Solving Delay Differential
Equations",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 368",
month = oct,
pages = 13,
year = 2000,
note = "To appear in Journal of Computational Applied Mathematics",
URL = "http://www.maths.man.ac.uk/~nareports/narep368.ps.gz",
abstract = "In this paper the efficient implementation of numerical software
for solving delay differential equations is addressed. Several
strategies that have been developed over the past 25 years for
improving the efficiency of delay differential equation solvers
are described. Of particular interest is a new method of
automatically constructing the network dependency graph used in
tracking derivative discontinuities."
}

@techreport{Baker:2000:DAE,
author = "Christopher Baker and Christopher Paul and Hongjiong Tian",
title = "Differential Algebraic Equations with After-Effect",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 367",
month = oct,
pages = 20,
year = 2000,
note = "Submitted to Journal of Computational Applied Mathematics",
URL = "http://www.maths.man.ac.uk/~nareports/narep367.ps.gz",
abstract = "We consider the numerical solution of delay differential
algebraic equations - they are differential algebraic equations
with after-effect, or constrained delay differential equations.
The general semi-explicit form of the problem consists of a set
of delay differential equations along with a set of constraints
that may involve retarded arguments. Even simply-stated problems
of this type can give rise to difficult analytical and numerical
problems. The more tractable examples can be shown to be
equivalent to systems of delay or neutral delay differential
equations. We shall explore, in particular through examples, some
of the complexities and obstacles that can arise when solving
these problems."
}

@techreport{MCCM:2001:AR,
author =       "Manchester Centre for Computational Mathematics",
title =        "Annual Report: {January--December} 1999",
type =         "Numerical Analysis Report",
number =       "No. 364",
institution =  inst-MCCM,
pages =        "24",
month =        mar,
year =         "2001",
URL =          "http://www.maths.man.ac.uk/~nareports/narep364.ps.gz",
}

@techreport{Dedieu:2001:PTH,
author = "Jean-Pierre Dedieu and Fran\c{c}oise Tisseur",
title = "Perturbation Theory for Homogeneous Polynomial Eigenvalue Problems",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 363",
pages = "21",
month = mar,
year = "2001",
URL =          "ftp://ftp.ma.man.ac.uk/pub/narep/narep363.ps.gz",
abstract = "We consider polynomial eigenvalue problems
$P(A,\alpha,\beta)x=0$ in which the matrix polynomial is
homogeneous in the eigenvalue $(\alpha,\beta)\in\mathbb{C}^2$.
In this framework infinite eigenvalues are on the same
footing as finite eigenvalues.  We view the problem in
projective spaces to avoid normalization of the
eigenpairs. We show that a polynomial eigenvalue problem
is well-posed when its eigenvalues are simple. We define
the condition numbers of a simple eigenvalue
$(\a,\beta)$ and a corresponding eigenvector $x$ and
show that the distance to the nearest ill-posed problem
is equal to the reciprocal of the condition number of
the eigenvector $x$.  We describe a bihomogeneous Newton
method for the solution of the homogeneous polynomial
eigenvalue problem.",
keywords = "polynomial eigenvalue problem, matrix polynomial,
}

@techreport{Silvester:00:EPB,
author = "David Silvester and Milan Mihajlovic",
title = "A black-box multigrid preconditioner for the biharmonic equation",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 362",
month = aug,
year = 2002,
note = "Submitted to BIT",
URL =          "http://www.maths.man.ac.uk/~nareports/narep362.ps.gz",
abstract = "We examine the convergence characteristics of a preconditioned
Krylov subspace solver applied to the linear systems arising from low-order
mixed finite element approximation of the biharmonic problem. The key
feature of our approach is that the preconditioning can be realized using
any black-box'' multigrid solver designed for the discrete Dirichlet
Laplacian operator. This leads to preconditioned systems having an
eigenvalue  distribution consisting of a tightly clustered set together
with a small number of outliers. Numerical results show that the
performance of the methodology is competitive with that of specialized
fast iteration  methods that have been developed in the context of
biharmonic problems."
}

@techreport{Baker:2000:CTM,
author = "Christopher T. H. Baker and Evelyn Buckwar",
title = "Continuous $\Theta$-methods for the Stochastic Pantograph
Equation",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 361",
month = may,
year = 2000,
URL =          "http://www.maths.man.ac.uk/~nareports/narep361.ps.gz",
abstract = "We consider a stochastic version of the pantograph equation:
\begin{eqnarray}\nonumber dX(t)  &= & \Big\{ a X(t)+ b X (qt)
\Big\} dt + \Big\{\sigma_1 + \sigma_2\ X(t) + \sigma_3 X
(qt)\Big\}dW(t), \$-1mm]\nonumber && \label{spanto} \\[-1mm]\nonumber X(0) &= & X_0, \nonumber \end{eqnarray} for  t \in [0,T], a given Wiener noise W and 0 < q < 1. This is an example of an It\^{o} stochastic delay differential equation with unbounded memory. We give the necessary analytical theory for existence and uniqueness of a strong solution of (\ddagger), and of strong approximations to the solution obtained by a continuous extension of the \Theta-Euler scheme (\Theta \in [0,1]). We establish {\mathcal O}(h^\half) mean-square convergence of approximations obtained using a bounded mesh of uniform step h, rising in the case of additive noise only to {\mathcal O}(h). Illustrative numerical examples are provided.", keywords = "Stochastic delay differential equation, continuous \Theta-method, mean-square convergence." } @techreport{Davies:2000:ACM, author = "Philip I. Davies and Nicholas J. Higham and Fran\c{c}oise Tisseur", title = "Analysis of the {Cholesky} Method with Iterative Refinement for Solving the Symmetric Definite Generalized Eigenproblem", institution = inst-MCCM, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 360", month = jun, pages = 21, year = "2000", URL = "http://www.maths.man.ac.uk/~nareports/narep360.pdf", abstract = "A standard method for solving the symmetric definite generalized eigenvalue problem Ax = \lambda Bx, where A is symmetric and B is symmetric positive definite, is to compute a Cholesky factorization B = LL^T (optionally with complete pivoting) and solve the equivalent standard symmetric eigenvalue problem C y = \l y where C = L^{-1} A L^{-T}. Provided that a stable eigensolver is used, standard error analysis says that the computed eigenvalues are exact for A+\dA and B+\dB with \max( \normt{\dA}/\normt{A}, \normt{\dB}/\normt{B} ) bounded by a multiple of \kappa_2(B)u, where u is the unit roundoff. We take the Jacobi method as the eigensolver and give a detailed error analysis that yields backward error bounds potentially much smaller than \kappa_2(B)u. To show the practical utility of our bounds we describe a vibration problem from structural engineering in which B is ill conditioned yet the error bounds are small. We show how, in cases of instability, iterative refinement based on Newton's method can be used to produce eigenpairs with small backward errors. Our analysis and experiments also give insight into the popular Cholesky--QR method, in which the QR method is used as the eigensolver. We argue that it is desirable to augment current implementations of this method with pivoting in the Cholesky factorization.", keywords = "symmetric definite generalized eigenvalue problem, Cholesky method, Cholesky factorization with complete pivoting, Jacobi method, backward error analysis, rounding error analysis, iterative refinement, Newton's method, LAPACK, MATLAB." } @techreport{Tisseur:2000:SPP, author = "Fran\c{c}oise Tisseur and Nicholas J. Higham", title = "Structured Pseudospectra for Polynomial Eigenvalue Problems, with Applications", institution = inst-MCCM, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 359", month = apr, pages = 22, year = "2000", URL = "http://www.maths.man.ac.uk/~nareports/narep359.pdf", abstract = "Pseudospectra associated with the standard and generalized eigenvalue problems have been widely investigated in recent years. We extend the usual definitions in two respects, by treating the polynomial eigenvalue problem and by allowing structured perturbations of a type arising in control theory. We explore connections between structured pseudospectra, structured backward errors, and structured stability radii. Two main approaches for computing pseudospectra are described. One is based on a transfer function and employs a generalized Schur decomposition of the companion form pencil. The other, specific to quadratic polynomials, finds a solvent of the associated quadratic matrix equation and thereby factorizes the quadratic \lambda-matrix. Possible approaches for large, sparse problems are also outlined. A collection of examples from vibrating systems, control theory, acoustics and fluid mechanics is given to illustrate the techniques.", keywords = "polynomial eigenvalue problem, \lambda-matrix, matrix polynomial, pseudospectrum, stability radius, backward error, transfer function, quadratic matrix equation, solvent, structured perturbations, Orr--Sommerfeld equation" } @techreport{Higham:2000:EPA, author = "Nicholas J. Higham", title = "Evaluating {Pad\'e} Approximants of the Matrix Logarithm", institution = inst-MCCM, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 358", month = mar, pages = 10, year = 2000, URL = "http://www.maths.man.ac.uk/~nareports/narep358.pdf", abstract = "The inverse scaling and squaring method for evaluating the logarithm of a matrix takes repeated square roots to bring the matrix close to the identity, computes a Pad\'e approximant, and then scales back. We analyze several methods for evaluating the Pad\'e approximant, including Horner's method (used in some existing codes), suitably customized versions of the Paterson--Stockmeyer method and Van Loan's variant, and methods based on continued fraction and partial fraction expansions. The computational cost, storage, and numerical accuracy of the methods are compared. We find the partial fraction method to be the best method overall and illustrate the benefits it brings to a transformation-free form of the inverse scaling and squaring method recently proposed by Cheng, Higham, Kenney and Laub. We comment briefly on how the analysis carries over to the matrix exponential.", keywords = "matrix logarithm, Pad\'e approximation, inverse scaling and squaring method, Horner's method, Paterson--Stockmeyer method, continued fraction, partial fraction expansion" } @techreport{Tisseur:2000:SSHE, author = "Fran\c{c}oise Tisseur", title = "Stability of Structured {Hamiltonian} Eigensolvers", institution = inst-MCCM, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 357", pages = "23", month = feb, year = "2000", URL = "ftp://ftp.ma.man.ac.uk/pub/narep/narep357.ps.gz", abstract = "Various applications give rise to eigenvalue problems for which the matrices are Hamiltonian or skew-Hamiltonian and also symmetric or skew-symmetric. We define structured backward errors that are useful for testing the stability of numerical methods for the solution of these four classes of structured eigenproblems. We introduce the symplectic quasi-QR factorisation and show that for three of the classes it enables the structured backward error to be efficiently computed. We also give a detailed rounding error analysis of some recently developed Jacobi-like algorithms of Fa{\ss}bender, Mackey \& Mackey for these eigenproblems. Based on the direct solution of 4\times 4, and in one case 8\times 8, structured subproblems these algorithms produce a complete basis of symplectic orthogonal eigenvectors for the two symmetric cases and a symplectic orthogonal basis for all the real invariant subspaces for the two skew-symmetric cases. We prove that, when the rotations are implemented using suitable formulae, the algorithms are strongly backward stable and we show that the QR algorithm does not have this desirable property", keywords = "Backward error, Hamiltonian, skew-Hamiltonian, symmetric-Hamiltonian, skew-symmetric skew-Hamiltonian, structure-preserving, rounding error, symplectic, quaternion rotation" } @techreport{Sturgeon:2000:SSD, author = "J. A. Sturgeon and R. M. Thomas and I. Gladwell", title = "Solving a Singular DAE Model of Unconfined Detonation", type = "Numerical Analysis Report", number = "No. 356", institution = inst-mccm, address = inst-mccm:adr, month = apr, year = "2000", pages = "27", URL = "http://www.maths.man.ac.uk/~nareports/narep356.ps.gz", abstract = "We consider a simplified model of an unconfined one-dimensional detonation problem, giving a brief survey of the history of the problem and of its numerical solution. This problem with its mathematical features is typical of those solved commercially by ICI plc, and the specific values used for the chemical constants in the example are typical of those of interest. Unfortunately, not all obvious methods work well, because of the singular nature of the problem at the Chapman-Jouguet shock front. We concentrate on shooting methods for the detonation problem based on backward differentiation formula integrators, and present a new analysis which explains how these methods work. Finally, we outline some possibilities for further work, including discussing a more general detonation problem, previous solutions and potential future solution methods", keywords = "detonation problem, differential-algebraic equations, backward differentiation formulae, shooting methods, numerical analysis, Chapman-Jouguet shock front" } @techreport{Langlois:1999:ALC, author = "Philippe Langlois", title = "Automatic Linear Correction of Rounding Errors", institution = inst-MCCM, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 355", month = nov, pages = "24", year = "1999", URL = "http://www.maths.man.ac.uk/~nareports/narep355.ps.gz", abstract = "A new automatic method to correct the first-order effect of floating point rounding errors on the result of a numerical algorithm is presented. A correcting term and a confidence threshold are computed using automatic differentiation, computation of elementary rounding error and running error analysis. Algorithms for which the accuracy of the result is not affected by higher order terms are identified. The correction is applied to the final result or to sensitive intermediate results. The properties and the efficiency of the method are illustrated with a sample numerical example.", keywords = "automatic error analysis, rounding error, floating point arithmetic" } @techreport{Davies:1999:NSG, author = "Philip I. Davies and Nicholas J. Higham", title = "Numerically Stable Generation of Correlation Matrices and their Factors", institution = inst-MCCM, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 354", month = nov, pages = "15", year = "1999", URL = "http://www.maths.man.ac.uk/~nareports/narep354.pdf", abstract = "Correlation matrices---symmetric positive semidefinite matrices with unit diagonal---are important in statistics and in numerical linear algebra. For simulation and testing it is desirable to be able to generate random correlation matrices with specified eigenvalues (which must be nonnegative and sum to the dimension of the matrix). A popular algorithm of Bendel and Mickey takes a matrix having the specified eigenvalues and uses a finite sequence of Given rotations to introduce 1s on the diagonal. We give improved formulae for computing the rotations and prove that the resulting algorithm is numerically stable. We show by example that the formulae originally proposed, which are used in certain existing Fortran implementations, can lead to serious instability. We also show how to modify the algorithm to generate a rectangular matrix with columns of unit 2-norm. Such a matrix represents a correlation matrix in factored form, which can be preferable to representing the matrix itself, for example when the correlation matrix is nearly singular to working precision.", keywords = "random correlation matrix, Bendel--Mickey algorithm, eigenvalues, singular value decomposition, test matrices, forward error bounds, relative error bounds, IMSL, NAG Library, Jacobi method" } @techreport{Cheng:1999:ALM, author = "Sheung Hun Cheng and Nicholas J. Higham and Charles S. Kenney and Alan J. Laub", title = "Approximating the Logarithm of a Matrix to Specified Accuracy", institution = inst-MCCM, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 353", pages = "14", month = nov, year = "1999", URL = "http://www.maths.man.ac.uk/~nareports/narep353.ps.gz", abstract = "The standard inverse scaling and squaring algorithm for computing the matrix logarithm begins by transforming the matrix to Schur triangular form in order to facilitate subsequent matrix square root and Pad\'e approximation computations. A transformation-free form of this method is presented that exploits incomplete Denman--Beavers square root iterations and aims for a specified accuracy. The error introduced by using approximate square roots is accounted for by a novel splitting lemma for logarithms of matrix products. The number of square root stages and the degree of the final Pad\'e approximation are chosen to minimize the computational work. This new method is attractive for high-performance computation since it uses only the basic building blocks of matrix multiplication, LU factorizati on and matrix inversion.", keywords = "matrix logarithm, Pad\'e approximation, inverse scaling and squaring method, matrix square root, Denman--Beavers iteration" } @techreport{Silvester:99:EPL, author = "David Silvester and Howard Elman and David Kay and Andrew Wathen", title = "Efficient Preconditioning of the Linearized {Navier-Stokes} Equations", institution = inst-MCCM, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 352", month = oct, year = "1999", note = "Solicited by the Journal of Computational and Applied Mathematics", URL = "http://www.maths.man.ac.uk/~nareports/narep352.ps.gz", abstract = "We outline a new class of robust and efficient methods for solving subproblems that arise in the linearization and operator splitting of Navier-Stokes equations. We describe a very general strategy for preconditioning that has two basic building blocks; a multigrid V-cycle for the scalar convection-diffusion operator, and a multigrid V-cycle for a pressure Poisson operator. We present numerical experiments illustrating that a simple implementation of our approach leads to an effective and robust solver strategy in that the convergence rate is independent of the grid, robust with respect to the time-step, and only deteriorates very slowly as the Reynolds number is increased." } @techreport{Ford:1999:HDN, author = "Neville J. Ford and Volker Wulf", title = "How Do Numerical Methods Perform for Delay Differential Equations Undergoing a {Hopf} Bifurcation?", institution = inst-MCCM, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 351", pages = "10", month = sep, year = "1999", URL = "http://www.maths.man.ac.uk/~nareports/narep351.ps.gz", abstract = "In this report we consider the numerical solution of delay differential equations that undergo a Hopf bifurcation. We show that, under suitable conditions on the choice of method, the Hopf bifurcation is maintained in the approximation, and that the class of bifurcation is also preserved.", keywords = "Delay differential equations, numerical methods, Hopf bifurcation" } @techreport{Ford:1999:Charroots, author = "Neville J. Ford", title = "Numerical Approximation of the Characteristic Values for a Delay Differential Equation", institution = inst-MCCM, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 350", pages = "11", month = sep, year = "1999", URL = "http://www.maths.man.ac.uk/~nareports/narep350.ps.gz", abstract = "In this report we consider the characteristic roots of numerical approximations of delay differential equations and their relationship to the characteristic values of the original problem. We are able to show that, as h -> 0, the characteristic roots are approximated to the order of the method.", keywords = "Delay differential equations, characteristic roots, numerical methods, stability coefficient" } @techreport{Rihan:99:SAP, author = "Christopher T. H. Baker and Fathalla A. Rihan", title = "Sensitivity Analysis of Parameters in Modelling With Delay-Differential Equations", institution = inst-MCCM, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 349", month = aug, pages = 22, year = 1999, URL = "http://www.maths.man.ac.uk/~nareports/narep349.ps.gz", abstract = "Many problems in bioscience for which observations are reported in the literature can be modelled by suitable functional differential equations incorporating a delay, parameterized by parameters  p_1,p_2,\dots,p_L. Given such observations (which usually contain error or 'noise'), we may determine the parameters \{p_i\} by fitting the model equations to the data and optimizing a measure of best fit. Our aim in this paper is to produce a new method to estimate (i) the sensitivity of the state variables to the parameter estimates \{ p_i\}, (ii) the sensitivity of the parameter estimates to the observations and (iii) the nonlinearity effects for delay differential models. The sensitivity of the parameter estimate to the observation is {\it low} if the sensitivity of the state variable to the parameter estimate is {\it high}. Sensitivity coefficients are used to determine the covariance matrix of parameter estimates and hence to determine the standard deviations. Numerical results are used to illustrate the results.", keywords = "Sensitivity analysis, parameter estimates, neutral delay differential equation, time-lag, nonlinearity effect." } @techreport{Norburn:99:FAS, author = "Sean Norburn and David Silvester ", title = "Fourier Analysis of Stabilised {Q1-Q1} Mixed Finite Element Approximation", institution = inst-MCCM, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 348", month = sep, year = 1999, note = "Submitted to the SIAM Journal on Numerical Analysis", URL = "http://www.maths.man.ac.uk/~nareports/narep348.ps.gz", abstract = "We use Fourier analysis to investigate the instability of an equal-order mixed finite element approximation method for elliptic incompressible flow equations. The lack of stability can be attributed to the fact that the associated discrete LBB (Ladyzhenskaya-Babuska-Brezzi) constant tends to zero as the mesh size is reduced. We develop a stabilisation approach that is appropriate to the periodic setting, and deduce optimal choices of the associated stabilisation parameter." } @techreport{Higham:99:NAQ, author = "Nicholas J. Higham and Hyun-Min Kim", title = "Numerical Analysis of a Quadratic Matrix Equation", institution = inst-MCCM, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 347", month = aug, pages = 26, year = 1999, URL = "http://www.maths.man.ac.uk/~nareports/narep347.pdf", abstract = "The quadratic matrix equation AX^2 + BX + C = 0 in n\times n matrices arises in applications and is of intrinsic interest as one of the simplest nonlinear matrix equations. We give a complete characterization of solutions in terms of the generalized Schur decomposition and describe and compare various numerical solution techniques. In particular, we give a thorough treatment of functional iteration methods based on Bernoulli's method. Other methods considered include Newton's method with exact line searches, symbolic solution and continued fractions. We show that functional iteration applied to the quadratic matrix equation can provide an efficient way to solve the associated quadratic eigenvalue problem (\lambda^2 A + \lambda B + C)x=0.", keywords = "quadratic matrix equation, solvent, generalized Schur decomposition, scaling, functional iteration, Bernoulli's method, Newton's method, exact line searches, continued fractions, quadratic eigenvalue problem", } @techreport{Tisseur:99:NMFPA, author = "Fran\c{c}oise Tisseur", title = "Newton's Method in Floating Point Arithmetic and Iterative Refinement of Generalized Eigenvalue Problems", institution = inst-MCCM, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 346", pages = "26", month = aug, year = "1999", URL = "http://www.maths.man.ac.uk/~nareports/narep346.ps.gz", abstract = "We examine the behaviour of Newton's method in floating point arithmetic, allowing for extended precision in computation of the residual, inaccurate evaluation of the Jacobian and unstable solution of the linear systems. We bound the limiting accuracy and the smallest norm of the residual. The application that motivates this work is iterative refinement for the generalized eigenvalue problem. We show that iterative refinement by Newton's method can be used to improve the forward and backward errors of computed eigenpairs.", keywords = " Newton's method, generalized eigenvalue problem, iterative refinement, backward error, forward error, rounding error analysis, limiting accuracy, limiting residual" } @techreport{Baker:99:INA, author = "Christopher T. H. Baker and Evelyn Buckwar", title = "Introduction to the Numerical Analysis of Stochastic Delay Differential Equations", institution = inst-MCCM, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 345", month = aug, year = 1999, URL = "http://www.maths.man.ac.uk/~nareports/narep345.ps.gz", abstract = "We consider the problem of strong approximations of the solution of stochastic delay differential equations of It\^{o} form \[ \begin{array}{rcl} dX(t) &= & f(t, X(t),X (t- \tau) ) \; dt + g(t, X(t),X (t- \tau) )\; dW(t), \ t \in [0,T] \\ X(t) &= & \Psi(t), \qquad t \in [-\tau,0], \end{array}$
for given $f, g$, Wiener noise $W$ and given $\tau > 0$,
with a prescribed initial function $\Psi$. We indicate the
nature of the equations of interest and the problems
that arise in a discussion of their numerical treatment,
provide necessary background material, and give a
convergence proof for an Euler-Maruyama scheme. Some
illustrative numerical examples are provided."
}

@techreport{MCCM:2000:AR,
author =       "Manchester Centre for Computational Mathematics",
title =        "Annual Report: {January--December} 1998",
type =         "Numerical Analysis Report",
number =       "No. 344",
institution =  inst-MCCM,
pages =        "24",
month =        mar,
year =         "2000",
URL =          "http://www.maths.man.ac.uk/~nareports/narep344.ps.gz",
}

@techreport{Rihan:99:RUD,
author = "Christopher T. H. Baker and Gennadii A. Bocharov and
Fathalla A. Rihan",
title = "A Report on the Use of Delay Differential Equations
in Numerical Modelling in the Biosciences",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 343",
month = aug,
pages = 45,
year = 1999,
URL =          "http://www.maths.man.ac.uk/~nareports/narep343.ps.gz",
abstract = "We review the application of numerical techniques to investigate
mathematical models of phenomena in the biosciences using delay
differential equations. We show that there are {\em prima facie}
reasons for using such models: $(i)$ they have a richer mathematical
framework (compared with ordinary differential equations) for the
analysis of biosystems dynamics, $(ii)$ they display better
consistency with the nature of the underlying processes and predictive
results. We now have suitable computational techniques to treat
numerically the emerging models for the biosciences.",
keywords = "Delay differential equations, dynamical systems,
biological systems, deterministic and stochastic models,
numerical modelling, parameter estimation, optimization."
}

@techreport{Edwards:1999:EBS,
author = "John T. Edwards and Jason A. Roberts",
title = "On the Existence of Bounded Solutions of Discrete Analogues for a
Nonlinear Integrodifferential Equation",
institution = inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 342",
pages = "14",
month = sep,
year = "1999",
URL =          "http://www.maths.man.ac.uk/~nareports/narep342.ps.gz",
abstract = "This work is concerned with the boundedness and  ultimately the
asymptotic stability of solutions of discrete analogues of a
nonlinear integrodifferential equation with convolution kernel.
In previous work, the stability analysis of the nonlinear
equation required the assumption that bounded solutions exist. We
prove in this work that, for a certain choice of discretization
of the equation, all solutions are                   bounded (and
hence asymptotically stable) provided the equation's parameters
fall within certain ranges.",
keywords = "Volterra integrodifferential equations, numerical methods,
difference equations, stability analysis"
}

@techreport{Higham:99:BAM,
author = "Nicholas J. Higham and Fran\c{c}oise Tisseur",
title = "A Block Algorithm for Matrix 1-Norm Estimation,
with an Application to 1-Norm Pseudospectra",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 341",
month = apr,
pages = 32,
year = 1999,
URL =          "http://www.maths.man.ac.uk/~nareports/narep341.ps.gz",
abstract = "The matrix 1-norm estimation algorithm used in LAPACK and various
other software libraries and packages has proved to be a valuable
tool. However, it has the limitations that it offers the user no
control over the accuracy and reliability of the estimate and
that it is based on level~2 BLAS operations. A block
generalization of the 1-norm power method underlying the
estimator is derived here and developed into a practical
algorithm applicable to both real and complex matrices. The
algorithm works with $n\times t$ matrices, where $t$ is a
parameter. For $t=1$ the original algorithm is recovered, but
with two improvements (one for real matrices and one for complex
matrices). The accuracy and reliability of the estimates
generally increases with $t$ and the computational kernels are
level~3 BLAS operations for $t>1$. The last $t-1$ columns of the
starting matrix are randomly chosen, giving the algorithm a
statistical flavour. As a by-product of our investigations we
identify a matrix for which the 1-norm power method takes the
maximum number of iterations. As an application of the new
estimator we show how it can be used to efficiently approximate
1-norm pseudospectra.",
keywords = "matrix 1-norm, matrix norm estimation, matrix condition number,
condition number estimation, $p$-norm power method, 1-norm
pseudospectrum, LAPACK, level 3 BLAS"
}

@techreport{Ford:1999:NHBClass,
author = "Neville J. Ford and Volker Wulf",
title = "Numerical {Hopf} Bifurcation for a Class of
Delay Differential Equations",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 340",
pages = "17",
month = feb,
year = "1999",
URL =          "http://www.maths.man.ac.uk/~nareports/narep340.ps.gz",
abstract = "In this report we consider discretisation of parameter-dependent
delay differential equations. We show that, if the delay
differential equation undergoes a Hopf bifurcation, then the
discrete scheme undergoes a Hopf bifurcation of the same type.
The results of this report extend our previous results for the
delay logistic equation to a wider class of problems.",
keywords = "Delay differential equations, Hopf bifurcation, numerical methods,
stability coefficient"
}

@techreport{Higham:99:SQM,
author = "Nicholas J. Higham and Hyun-Min Kim",
title = "Solving a Quadratic Matrix Equation by
{Newton's} Method with Exact Line Searches",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 339",
month = jan,
pages = 18,
year = 1999,
URL =          "http://www.maths.man.ac.uk/~nareports/narep339.pdf",
abstract = "We show how to incorporate exact line searches into Newton's
method for solving the quadratic matrix equation $AX^2 + BX + C = 0$, where $A$, $B$ and $C$ are square matrices. The line searches
are relatively inexpensive and improve the global convergence
properties of Newton's method in theory and in practice. We also
derive a condition number for the problem and show how to compute
the backward error of an approximate solution.",
keywords = "quadratic matrix equation, solvent, Newton's method, generalized
Sylvester equation, exact line searches, quadratic eigenvalue
problem, condition number, backward error"
}

@techreport{Davies:99:GTM,
author = "Philip I. Davies and Nicholas J. Higham",
title = "Generating Test Matrices for the One- and Two-Sided {Jacobi} Methods",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 338",
month = jan,
pages = "13",
year = 1999,
URL =          "http://www.maths.man.ac.uk/~nareports/narep338.ps.gz",
abstract = "For a symmetric positive definite matrix the relative error in
the eigenvalues computed by the two-sided Jacobi method is
bounded in terms of the condition number of the matrix scaled to
have unit diagonal. Similarly, for a general matrix the relative
error in the singular values computed by the one-sided Jacobi
method is bounded in terms of the condition number of the matrix
scaled to have rows or columns of unit-2 norm. We show how to
generate random matrices having these scalings and given
eigenvalues and singular values, respectively. For the two-sided
Jacobi method we apply an algorithm of Bendel and Mickey for
generating random correlation matrices, with an improved formula
for the rotations. We show how to modify the algorithm to
generate matrices for the one-sided Jacobi method. Using these
test matrices we show empirically that the forward error bounds
for the one- and two-sided Jacobi methods are sharp.",
keywords = "Jacobi method, correlation matrix, eigenvalues, singular value
decomposition, test matrices, forward error bounds, relative
error bounds"
}

@techreport{Paul:99:TDDDE,
author = "Christopher A.H. Paul",
title = "The Treatment of Derivative Discontinuities in Differential
Equations",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 337",
month = jan,
pages = 11,
year = "1999",
URL =          "http://www.maths.man.ac.uk/~nareports/narep337.ps.gz",
abstract = "The assumption of sufficiently smooth derivatives underlies much
of the analysis of numerical methods for both ordinary and delay
differential equations. However, derivative discontinuities can
arise in ordinary differential equations and usually do arise in
delay differential equations. In this paper, I review two of the
approaches for treating derivative discontinuities, tracking''
using the discontinuity tracking equations and detection and
location'' using defect error control. I conclude that neither
approach is ideal when it comes to the treatment of derivative
discontinuities in delay differential algebraic equations.",
keywords = "differential equations, discontinuities, numerical solution"
}

@techreport{Higham:99:NSM,
author = "Nicholas J. Higham",
title = "A New \texttt{sqrtm} for \textsc{Matlab}",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 336",
month = jan,
pages = 11,
year = "1999",
URL =          "http://www.maths.man.ac.uk/~nareports/narep336.ps.gz",
abstract = "\Matlab's function \t{sqrtm} computes a square root of a matrix.
We propose a replacement for the \t{sqrtm} in \Matlab~5.2 that is
more accurate and returns useful information about the stability
and conditioning of the problem.",
keywords = "matrix square root, \Matlab, Schur decomposition, condition
number, stability"
}

@techreport{Tshelametse:98:ANT,
author = "Ronald Tshelametse and Jack Williams",
title = "A New Termination Criterion for Nonlinear Iterations in {ODE} Codes",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 334",
month = oct,
pages = 23,
year = "1998",
URL =          "http://www.maths.man.ac.uk/~nareports/narep334.ps.gz",
abstract = "Many codes that solve systems of stiff initial value problems for
ordinary differential equations implement implicit linear
multistep methods. This yields nonlinear equations which are
commonly solved using some form of Newton iteration. For linear
multistep methods we apply the theory of Dorsselaer and Spijker
\cite{dorss94} to justify a new stopping criterion for the Newton
iteration. The ideas are tested on the \textsc{Matlab} ODE suite
of Shampine and Reichelt \cite{shamp97} and are shown to lead to
a significant improvement in overall efficiency for some
problems.",
keywords = "stiff problems, IVPs, ODEs, BDFs, Newton method"
}

@techreport{Higham:98:NAS,
author = "Nicholas J. Higham",
title = "Notes on Accuracy and Stability of Algorithms in
Numerical Linear Algebra",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 333",
month = aug,
pages = 42,
year = "1998",
URL =          "http://www.maths.man.ac.uk/~nareports/narep333.ps.gz",
note = "To appear in proceedings of EPSRC Numerical Analysis Summer School,
Leicester University, July 1998."
}

@techreport{Tisseur:1998:BEC,
author = "Fran\c{c}oise Tisseur",
title = "Backward Error and Condition of Polynomial Eigenvalue Problems",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 332",
pages = "20",
month = aug,
year = "1998",
URL =          "http://www.maths.man.ac.uk/~nareports/narep332.ps.gz",
abstract = "We develop normwise backward errors and condition
numbers for the polynomial eigenvalue problem (PEP).
The standard way of dealing
with the PEP is to reformulate it as a generalized
eigenvalue problem (GEP). For
the special case of the quadratic eigenvalue problem (QEP), we
show that solving the QEP by applying the QZ algorithm to a
corresponding GEP can be backward unstable.  The QEP can
be reformulated as a GEP in many ways. We investigate
the sensitivity of a given eigenvalue to perturbations
in each ofthe GEP formulations and identify which
formulations are to be preferred for large and small
eigenvalues respectively.",
keywords = "Polynomial eigenvalue problem,
problem, backward error, condition number"
}

@techreport{Barlow:1998:MAB,
author = "Jesse L. Barlow",
title = "More Accurate Bidiagonal Reduction for Computing the Singular Value
Decomposition",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 331",
month = aug,
year = "1998",
URL =          "http://www.maths.man.ac.uk/~nareports/narep331.ps.gz",
abstract = "Bidiagonal reduction is the preliminary stage for the fastest
stable algorithms for computing the singular value decomposition.
However, the best error bounds on bidiagonal reduction methods
are of the form $A + \delta A = UBV^T ,$ $\normtwo{\delta A } \leq \macheps f(n) \normtwo{A}$ where $B$ is bidiagonal, $U$
and $V$ are orthogonal, $\macheps$ is machine precision, and
$f(n)$ is a modestly growing function of the dimensions of $A$.
A Givens-based bidiagonal reduction procedure is proposed that satisfies
$A + \delta A = U (B + \delta B ) V^T,$
where $\delta B$ is bounded {\em componentwise} and $\delta A$
satisfies a tighter columnwise bound. Thus the routine obtains
more accurate singular values for matrices that have poor column
scaling or those arising from rank revealing decompositions.",
keywords = ""
}

@techreport{Cheng:1998:NDP,
author = "Sheung Hun Cheng and Nicholas J. Higham",
title = "The Nearest Definite Pair for the
{Hermitian} Generalized Eigenvalue Problem",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 325",
pages = "15",
month = may,
year = "1998",
URL =          "http://www.maths.man.ac.uk/~nareports/narep325.ps.gz",
abstract = "The generalized eigenvalue problem $Ax = \lambda Bx$ has special
properties when $(A,B)$ is a Hermitian and definite pair. Given a
general Hermitian pair $(A,B)$ it is of interest to find the
nearest definite pair having a specified Crawford number $\delta > 0$. We solve the problem in terms of the inner numerical radius
associated with the field of values of $A+iB$. We show that once
the problem has been solved it is trivial to rotate the perturbed
pair $(A+\dA,B+\dB)$ to a pair $(\widetilde{A},\widetilde{B})$
for which $\lambda_{\min}(\widetilde{B})$ achieves its maximum
value $\delta$, which is a numerically desirable property when
solving the eigenvalue problem by methods that convert to a
standard eigenvalue problem by inverting $B$''. Numerical
examples are given to illustrate the analysis.",
keywords = "Nearest definite pair, Crawford number, Hermitian pair,
generalized eigenvalue problem, field of values, inner numerical
}

@techreport{Higham:1998:QRF,
author = "Nicholas J. Higham",
title = "{QR} Factorization with Complete Pivoting
and Accurate Computation of the {SVD}",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 324",
month = sep,
pages = 26,
year = "1998",
URL =          "http://www.maths.man.ac.uk/~nareports/narep324.ps.gz",
abstract = "A new algorithm of Demmel et al.\ for computing the singular
value decomposition (SVD) to high relative accuracy begins by
computing a rank-revealing decomposition. Demmel et al.\ analyse
the use of Gaussian elimination with complete pivoting (GECP) for
computing the rank-revealing decomposition. We investigate the
use of QR factorization with complete pivoting (that is, column
pivoting together with row sorting or row pivoting) as an
alternative to GECP, since this leads to a faster SVD algorithm.
We derive a new componentwise backward error result for
Householder QR factorization and combine it with the theory of
Demmel et al.\ to show when high relative accuracy in the
computed SVD can be expected. An a posteriori error bound is
derived that gives useful estimates of the relative accuracy of
the computed singular values. Numerical experiments confirm the
theoretical predictions.",
keywords = "QR factorization, Householder matrix, row pivoting, row sorting,
column pivoting, complete pivoting, backward error analysis,
singular value decomposition, relative accuracy"
}

@techreport{Ford:1998:NHB,
author = "Neville J. Ford and Volker Wulf",
title = "Numerical {Hopf} Bifurcation for the Delay Logistic Equation",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 323",
pages = "20",
month = apr,
year = "1998",
URL =          "http://www.maths.man.ac.uk/~nareports/narep323.ps.gz",
abstract = "This report is concerned with the approximate solution of the
delay logistic equation close to the Hopf bifurcation point. We
show that, for sufficiently small step sizes the bifurcation
point of the original problem is approximated by a bifurcation
point in the discrete problem and that furthermore, the nature of
the bifurcation is preserved.",
keywords = "Delay differential equations, Hopf bifurcation, numerical
methods, stability coefficient"
}

@techreport{Ford:1998:BLP,
author = "Neville J. Ford and Volker Wulf",
title = "The Use of Boundary Locus Plots in the Identification of Bifurcation
Points in Numerical Approximation of Delay Differential Equations",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 322",
pages = "15",
month = mar,
year = "1998",
URL =          "http://www.maths.man.ac.uk/~nareports/narep322.ps.gz",
abstract = "We use the boundary locus method as a tool to determine precise
parameter values at which discrete analogues of delay
differential equations may have a Hopf bifurcation. We prove that,
for strictly stable multistep methods, the boundary of the
stability region for the discrete equation approximates the true
stability region to the order of the method and that Hopf
bifurcation points in the discrete scheme approximate the
bifurcation points for the original equation to the same order.",
keywords = "Delay differential equations, Hopf bifurcations,
Boundary locus method"
}

@techreport{Cox:1998:BEB,
author = "Anthony J. Cox and Nicholas J. Higham",
title = "Backward Error Bounds for Constrained Least Squares Problems",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 321",
pages = "20",
month = apr,
year = "1998",
URL =          "http://www.maths.man.ac.uk/~nareports/narep321.ps.gz",
abstract = "We derive an upper bound on the normwise backward error of an
approximate solution to the equality constrained least squares
problem $\min_{Bx=d}\|b-Ax\|_2$. Instead of minimizing over the
four perturbations to $A$, $b$, $B$ and $d$, we fix those to $B$
and $d$ and minimize over the remaining two; we obtain an
explicit solution of this simplified minimization problem. Our
experiments show that backward error bounds of practical use are
obtained when $B$ and $d$ are chosen as the optimal normwise
relative backward perturbations to the constraint system, and we
find that when the bounds are weak they can be improved by direct
search optimization. We also derive upper and lower backward
error bounds for the problem of least squares minimization over a
sphere: $\min_{\|x\|_2\le\alpha}\|b-Ax\|_2$.",
keywords = "Equality constrained least squares problem, least squares
minimization over a sphere, null space method, elimination
method, method of weighting, backward error, backward stability"
}

@techreport{MCCM:1998:AR,
author =       "Manchester Centre for Computational Mathematics",
title =        "Annual Report: {January--December} 1997",
type =         "Numerical Analysis Report",
number =       "No. 320",
institution =  inst-MCCM,
pages =        "20",
month =        mar,
year =         "1998",
URL =          "http://www.maths.man.ac.uk/~nareports/narep320.ps.gz",
}

@techreport{Cox:1998:RWB,
author = "Anthony J. Cox and Nicholas J. Higham",
title = "Row-Wise Backward Stable Elimination Methods
for the Equality Constrained Least Squares Problem",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 319",
pages = "18",
month = mar,
year = "1998",
URL =          "http://www.maths.man.ac.uk/~nareports/narep319.ps.gz",
abstract = " It is well known that the solution of the equality constrained
least squares (LSE) problem $\min_{Bx=d}\|b-Ax\|_2$ is the limit
of the solution of the unconstrained weighted least squares
problem $$\min_x\left\| \bmatrix{ \mu d \cr b } - \bmatrix{\mu B \cr A } x \right\|_2$$ as the weight $\mu$ tends to infinity,
assuming that $\bmatrix{B^T & A^T \cr}^T$ has full rank. We
derive a new method for the LSE problem by applying Householder
QR factorization with column pivoting to this weighted problem
and taking the limit analytically, with an appropriate rescaling
of rows. The method obtained is a type of direct elimination
method. We adapt existing error analysis for the unconstrained
problem to obtain a row-wise backward error bound for the
method. The bound shows that, provided row pivoting or row
sorting is used, the method is well-suited to problems in which
the rows of $A$ and $B$ vary widely in norm. As a by-product of
our analysis, we show that the standard elimination method for
solving the LSE problem satisfies a row-wise backward error
bound of precisely the same form as the new method. We
illustrate our results with numerical tests.",
keywords = "Constrained least squares problem, weighted least squares
problem, Householder QR factorization, Gaussian elimination,
elimination method, rounding error analysis, backward stability,
row pivoting, row sorting, column pivoting"
}

@techreport{Tisseur:1998:PDC,
author = "Fran\c{c}oise Tisseur and Jack J. Dongarra",
title = "Parallelizing the Divide and Conquer Algorithm for
the Symmetric Tridiagonal Eigenvalue Problem on
Distributed Memory Architectures",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 317",
pages = "28",
month = mar,
year = "1998",
URL =          "http://www.maths.man.ac.uk/~nareports/narep317.ps.gz",
abstract = "We present a new parallel implementation of a divide and
conquer algorithm for computing the spectral
decomposition of a symmetric tridiagonal matrix on
distributed memory architectures.  The implementation we
develop differs from other implementations in that we use a two
dimensional block cyclic distribution of the data, we
use the L{\"o}wner theorem approach to compute orthogonal
eigenvectors, and we introduce permutations before the
back transformation of each rank-one update in order to
make good use of deflation.  This algorithm yields the
first scalable, portable and numerically stable parallel
divide and conquer eigensolver.  Numerical results
confirm the effectiveness of our algorithm.  We compare
performance of the algorithm with that of the QR
algorithm and of bisection followed by inverse iteration
on an IBM SP2 and a cluster of Pentium PII's.",
keywords = "Divide and conquer, symmetric eigenvalue problem,
tridiagonal matrix, rank-one modification, parallel algorithm,
ScaLAPACK, LAPACK, distributed memory architecture"
}

@techreport{kay:97:APE,
author = "David Kay and David Silvester",
title = "A-Posteriori Error Estimation for Stabilised Mixed Approximations of
the {Stokes} Equations.",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 316",
month = dec,
year = 1997,
note = "Submitted to the SIAM Journal on Scientifc Computing",
URL =          "http://www.maths.man.ac.uk/~nareports/narep316.ps.gz",
abstract = "The practical implementation  of the lowest-order \plpo\ (linear
velocity, constant pressure) finite element method for the
addressed in this work. Three different types of a-posteriori
error indicator are introduced and each is shown to be equivalent
to the discretisation  error. Our numerical results show that
these indicators can be used to drive an adaptive refinement
process which is  specially tailored to create grids which
conform to the requirements of the local stabilisation. It is
also shown that the indicators provide an effective method for
detecting local singularities in the flow solution."
}

@techreport{Baker:1997:MMI2,
author = "Christopher T. H. Baker and Gennadii A. Bocharov and
Christopher A. H. Paul",
title = "Mathematical Modelling Of The Interleukin-2 T-Cell System: A
Comparative Study Of Approaches Based On Ordinary And Delay
Differential Equations",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 314",
pages = "17",
month = nov,
year = "1997",
URL =          "http://www.maths.man.ac.uk/~nareports/narep314.ps.gz",
abstract = "Cell proliferation and differentiation phenomena are key issues
in immunology, tumour growth and cell biology. We study the
kinetics of cell growth in the immune system using mathematical
models formulated in terms of ordinary and delay differential
equations. We study how the suitability of the mathematical
models depends on the nature of the cell growth data and the
types of differential equations by minimizing an objective
function to give a best-fit parameterized solution. We show that
mathematical models that incorporate a time-lag in the cell
division phase are more consistent with certain reported data.
They also allow various cell proliferation characteristics to be
estimated directly, such as the average cell-doubling time and
the rate of commitment of cells to cell division. Specifically,
we study the Interleukin-2-dependent cell division of
phytohemagglutinin stimulated T-cells -- the model of which can
be considered to be a general model of cell growth. We also
review the numerical techniques available for solving delay
differential equations and calculating the least-squares
best-fit parameterized solution.",
keywords = "Cell proliferation, Interleukin-2, mathematical modelling,
parameter estimation, time-lag"
}

@techreport{Baker:1997:MATL,
author = "Christopher T. H. Baker and Gennadii A. Bocharova and
Christopher A. H. Paul and Fathalla A. Rihan",
title = "Modelling and Analysis of Time-Lags in Cell Proliferation",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 313",
pages = "16",
month = nov,
year = "1997",
URL =          "http://www.maths.man.ac.uk/~nareports/narep313.ps.gz",
abstract = "In this paper, we present a systematic approach for obtaining
qualitatively and quantitatively correct mathematical models of
some biological phenomena with time-lags. Features of our
approach are the development of a hierarchy of related models
and the estimation of parameter values, along with their
non-linear biases and standard deviations, for sets of
experimental data.
We demonstrate our method of solving parameter estimation
problems for neutral delay differential equations by analyzing
some models of cell growth that incorporate a time-lag in the
cell division phase. We show that these models are more
consistent with certain reported data than the classic
exponential growth model. Although the exponential growth model
provides estimates of some of the growth characteristics, such
as the population-doubling time, the time-lag growth models can
additionally provide estimates of: (i) the fraction of cells
that are dividing, (ii) the rate of commitment of cells to cell
division, (iii) the initial distribution of cells in the cell
cycle, and (iv) the degree of synchronization of cells in the
(initial) cell population.",
keywords = "Cell proliferation, mathematical modelling, time-lag, neutral
delay differential equation"
}

@techreport{Ford:1997:SDA,
author = "Neville J. Ford and John T. Edwards and Jason A. Roberts and
Leonid E. Shaikhet",
title = "Stability of a difference analogue for a nonlinear
integro differential equation of convolution type",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 312",
pages = "15",
month = oct,
year = "1997",
URL =          "http://www.maths.man.ac.uk/~nareports/narep312.ps.gz",
abstract = "This paper is concerned with the stability of discrete
analogues of nonlinear integro-differential equations.
We apply the general method of Lyapunov functional
construction to give sufficient conditions for
the asymptotic stability of the zero solution",
keywords = "Integro-differential equations, numerical stability,
Lyapunov stability"
}

@techreport{Ford:1997:NVI,
author = "Neville J. Ford and Christopher T. H. Baker and Jason A. Roberts",
title = "Nonlinear {Volterra} integro-differential equations -
stability and numerical stability of theta methods",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 311",
pages = "15",
month = oct,
year = "1997",
URL =          "http://www.maths.man.ac.uk/~nareports/narep311.ps.gz",
abstract = "This paper considers whether existing stability results
for a non-linear integro-differential equation are
preserved when a solution is found using numerical
approximations. We choose to analyse theta quadrature methods
and we apply the direct method of Lyapunov for the stability analysis",
keywords = "Integro-differential equations, numerical stability, Lyapunov stability"
}

@techreport{Ford:1997:QBA,
author = "Neville J. Ford and John T. Edwards and Kurt Frischmuth",
title = "Qualitative behaviour of approximate solutions of some {Volterra}
integral equations with non-{Lipschitz} nonlinearity",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 310",
pages = "17",
month = oct,
year = "1997",
URL =          "http://www.maths.man.ac.uk/~nareports/narep310.ps.gz",
abstract = "In some recent work the  long-term behaviour of numerical
solutions of the nonlinear convolution equation has been
considered. In the present paper, we consider examples
which do not satisfy classical  conditions guaranteeing
existence and uniqueness  of  the  exact solutions and
suitable behaviour  of  approximate solutions. We are able
to give a strengthened  version of a theorem given
originally by Corduneanu for certain kernels and
nonlinearities. On that basis  we consider how
certain numerical codes may fail. We explore and test methods
known to preserve properties of the solutions in the linear
case",
keywords = "Volterra integral equations, non-Lipschitz kernel, numerical
methods"
}

author = "John T. Edwards and Jason A. Roberts and Neville J. Ford",
title = "A comparison of {Adomian's} decomposition method and {Runge Kutta}
methods for approximate solution of some predator prey
model equations",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 309",
pages = "17",
month = oct,
year = "1997",
URL =          "http://www.maths.man.ac.uk/~nareports/narep309.ps.gz",
for example, by Bellman in 1986 as a significant powerful method
for solving systems of nonlinear equations. A recent paper by
Olek indicates that ADM is a powerful tool in the solution
of some predator prey model equations and that it compares
favourably with alternative approximate methods (both
analytical and numerical). In this report we investigate
whether this property holds more generally",
keywords = "Adomian's Decomposition Method, Predator Prey equations"
}

@techreport{Higham:1997:SBF,
author = "Nicholas J. Higham",
title = "Stability of Block $\mathrm{LDL^T}$ Factorization of a
Symmetric Tridiagonal Matrix",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 308",
month = sep,
pages = "8",
year = "1997",
URL = "http://www.maths.man.ac.uk/~nareports/narep308.ps.gz",
abstract = "For symmetric indefinite tridiagonal matrices, block
$\mathrm{LDL^T}$ factorization without interchanges is shown to
have excellent numerical stability when a pivoting strategy of
Bunch is used to choose the dimension (1 or 2) of the pivots.",
keywords = "tridiagonal matrix, symmetric indefinite matrix, diagonal
pivoting method, $\mathrm{LDL^T}$ factorization, growth factor,
numerical stability, rounding error analysis, LAPACK,
LINPACK."
}

@techreport{Norburn:1997:SSM,
author = "Sean Norburn and David Silvester",
title = "Stabilised Vs Stable Mixed Methods for Incompressible Flow",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 307",
month = sep,
year = "1997",
note = "Submitted to Computer Methods in Applied Mechanics and Engineering",
URL = "http://www.maths.man.ac.uk/~nareports/narep307.ps.gz",
abstract = "The accuracy of low-order mixed finite element methods for the
incompressible (Navier-)Stokes equations is investigated in this
work. Some numerical experiments suggest that the lowest-order
stabilised $P_1$--$P_0$ method (linear velocity, constant
pressure) is more efficient than the alternative a-priori stable
non-conforming Crouzeix-Raviart ($P_{-1}$--$P_0$) approach. The
relative accuracy of stabilised $P_1$--$P_0$ and $P_1$--$P_1$
(linear velocity, continuous linear pressure) is also assessed
herein."
}

@techreport{Cox:1997:ASN,
author = "Anthony J. Cox and Nicholas J. Higham",
title = "Accuracy and Stability of the Null Space Method for Solving the
Equality Constrained Least Squares Problem",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 306",
pages = "20",
month = aug,
year = "1997",
URL = "http://www.maths.man.ac.uk/~nareports/narep306.ps.gz",
abstract = "The null space method is a standard method for solving the linear
least squares problem subject to equality constraints (the LSE
problem). We show that three variants of the method, including
one used in LAPACK that is based on the generalized QR
factorization, are numerically stable. We derive two perturbation
bounds for the LSE problem: one of standard form that is not
attainable, and a bound that yields the condition number of the
LSE problem to within a small constant factor. By combining the
backward error analysis and perturbation bounds we derive an
approximate forward error bound suitable for practical
computation. Numerical experiments are given to illustrate the
sharpness of this bound.",
keywords = "Constrained least squares problem, null space method,
rounding error analysis, condition number,
generalized QR factorization, LAPACK"
}

@techreport{Higham:1997:SIM,
author = "Nicholas J. Higham",
title = "Stable Iterations for the Matrix Square Root",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 305",
month = apr,
pages = "20",
year = "1997",
URL = "http://www.maths.man.ac.uk/~nareports/narep305.pdf",
abstract = "Any matrix with no nonpositive real eigenvalues has a unique
square root for which every eigenvalue lies in the open right
half-plane. A link between the matrix sign function and this
square root is exploited to derive both old and new iterations
for the square root from iterations for the sign function. One
new iteration is a quadratically convergent Schulz iteration
based entirely on matrix multiplication; it converges only
locally, but can be used to compute the square root of any
nonsingular $M$-matrix. A new Pad\'e iteration well suited to
parallel implementation is also derived and its properties
explained. Iterative methods for the matrix square root are
notorious for suffering from numerical instability. It is shown
that apparently innocuous algorithmic modifications to the Pad\'e
iteration can lead to instability, and a perturbation analysis is
given to provide some explanation. Numerical experiments are
included and some advice is offered on the choice of iterative
method for computing the matrix square root.",
keywords = "Matrix square root, matrix logarithm, matrix sign function,
$M$-matrix, symmetric positive definite matrix, Pad\'e
approximation, numerical stability, Newton's method, Schulz
method"
}

@techreport{MCCM:1997:AR,
author =       "Manchester Centre for Computational Mathematics",
title =        "Annual Report: {January--December} 1996",
type =         "Numerical Analysis Report",
number =       "No. 304",
institution =  inst-MCCM,
pages = "16",
year =         "1997",
URL =          "http://www.maths.man.ac.uk/~nareports/narep304.ps.gz",
}

@techreport{Braconnier:1997:CBE,
author = "Thierry Braconnier and Fran\c{c}oise Chaitin-Chatelin",
title = "Chaotic Behavior for Eigensolvers Applied on
Highly Nonnormal Matrices in Finite Precision",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 303",
pages = "15",
month = may,
year = "1997",
URL = "http://www.maths.man.ac.uk/~nareports/narep303.ps.gz",
abstract = "In this paper, we provide examples of nonnormal matrices where
the power method is shown to be backward unstable. Investigating
the behavior of the computed sequence of Rayleigh quotients, we
show via two measures of chaos, i.e, the computation of the
Hausdorff dimension (more precisely of the correlation dimension)
and the computation of the Lyapunov exponent, that the computed
dominant eigenvalue lies on an attractor which has a dimension
larger than the dimension predicted by the theoretical analysis
and which increases with the nonnormality of the matrix.",
keywords = "Backward stability, Hausdorff dimension, correlation dimension,
Lyapunov exponents, eigenvalues, departure from normality."
}

@techreport{Williams:1997:LHC,
author = "Jack Williams and Z. Kalogiratou",
title = "The Local {Haar} Condition in Parameter
Estimation for Second Order Ordinary Differential Equations",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 302",
month = feb,
pages = "17",
year = "1997",
URL = "http://www.maths.man.ac.uk/~nareports/narep302.ps.gz"
}

@techreport{Cox:1997:SHQ,
author = "Anthony J. Cox and Nicholas J. Higham",
title = "Stability of {Householder QR} Factorization
for Weighted Least Squares Problems",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 301",
pages = "17",
month = feb,
year = "1997",
URL = "http://www.maths.man.ac.uk/~nareports/narep301.ps.gz",
abstract = "For least squares problems in which the rows of the coefficient
matrix vary widely in norm, Householder QR factorization (without
pivoting) has unsatisfactory backward stability properties.
Powell and Reid showed in 1969 that the use of both row and
column pivoting leads to a desirable row-wise backward error
result. We give a reworked backward error analysis in modern
notation and prove two new results. First, sorting the rows by
decreasing $\infty$-norm at the start of the factorization
obviates the need for row pivoting. Second, row-wise backward
stability is obtained for only one of the two possible choices of
keywords = "Weighted least squares problem, Householder matrix, QR
factorization, row pivoting, column pivoting, backward error
analysis"
}

@techreport{Usman:1996:ESP,
author = "Anila Usman and George Hall",
title = "Equilibrium States for Predictor-Corrector Methods",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 300",
month = dec,
pages = "35",
year = "1996",
URL = "http://www.maths.man.ac.uk/~nareports/narep300.ps.gz"
}

@techreport{Baker:1996:VCT,
author = "Christopher T. H. Baker and Arsalang Tang",
title = "Volterra Centennial Meetings---Invited talks at {Arlington and
Tempe}",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 299",
month = nov,
pages = "28",
year = "1996",
URL = "http://www.maths.man.ac.uk/~nareports/narep299.ps.gz"
}

@techreport{Higham:1996:FCS,
author = "Nicholas J. Higham",
title = "Factorizing Complex Symmetric Matrices with
Positive Definite Real and Imaginary Parts",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 298",
month = nov,
pages = "12",
year = "1996",
URL =          "http://www.maths.man.ac.uk/~nareports/narep298.ps.gz",
abstract = "Complex symmetric matrices whose real and imaginary parts are
positive definite are shown to have a growth factor bounded by 2
for LU factorization. This result adds to the classes of matrix
for which it is known to be safe not to pivot in LU
factorization. Block $\mathrm{LDL^T}$ \fact\ with the pivoting
strategy of Bunch and Kaufman is also considered, and it is shown
that for such matrices only $1\times 1$ pivots are used and the
same growth factor bound of 2 holds, but that interchanges that
destroy band structure may be made. The latter results hold
whether the pivoting strategy uses the usual absolute value or
the modification employed in LINPACK and LAPACK.",
keywords = "Complex symmetric matrices, LU factorization, diagonal pivoting
factorization, block $\mathrm{LDL^T}$ factorization,
Bunch--Kaufman pivoting strategy, growth factor, band matrix,
LINPACK, LAPACK"
}

@techreport{Higham:1996:SBE,
author = "Desmond J. Higham and Nicholas J. Higham",
title = "Structured Backward Error and Condition of
Generalized Eigenvalue Problems",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 297",
month = nov,
pages = "24",
year = "1996",
URL =          "http://www.maths.man.ac.uk/~nareports/narep297.pdf",
abstract = "Backward errors and condition numbers are defined and evaluated
for eigenvalues and eigenvectors of generalized eigenvalue
problems. Both normwise and componentwise measures are used.
Unstructured problems are considered first, and then the basic
definitions are extended so that linear structure in the
coefficient matrices (for example, Hermitian, Toeplitz or band
structure) is preserved by the perturbations.",
keywords = "Generalized eigenvalue problem, backward error, condition number,
structured matrices"
}

@techreport{Barlow:1996:GGE,
author = "Jesse L. Barlow and Hongyuan Zha",
title = "Growth in {Gaussian} Elimination, Orthogonal Matrices, and the
{Euclidean} Norm",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 296",
month = sep,
year = "1996",
URL =          "http://www.maths.man.ac.uk/~nareports/narep296.ps.gz",
abstract = "It is shown that maximal growth for Gaussian elimination as
measured in the Euclidean norm is achieved by orthogonal
matrices. A new, sharp bound on that growth is given.",
keywords = "LU factorization, orthogonal invariance, triangular matrix"
}

@techreport{Higham:1996:MIM,
author = "Nicholas J. Higham and Sheung Hun Cheng",
title = "Modifying the Inertia of Matrices Arising in Optimization",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 295",
month = sep,
pages = "17",
year = "1996",
URL =          "http://www.maths.man.ac.uk/~nareports/narep295.ps.gz",
abstract = "Applications in constrained optimization (and other areas)
produce symmetric matrices with a natural block $2\times 2$
structure. An optimality condition leads to the problem of
perturbing the $(1,1)$ block of the matrix to achieve a specific
inertia. We derive a perturbation of minimal norm, for any
unitarily invariant norm, that increases the number of
nonnegative eigenvalues by a given amount, and we show how it can
be computed efficiently given a factorization of the original
matrix. We also consider an alternative way to satisfy the
optimality condition based on a projection approach. Theoretical
tools developed here include an extension of Ostrowski's theorem
on congruences and some lemmas on inertias of block $2\times 2$
symmetric matrices.",
keywords = "Inertia, optimization, nonlinear programming, unitarily
invariant norm"
}

@techreport{Higham:1996:TLA,
author = "Nicholas J. Higham",
title = "Testing Linear Algebra Software",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 294",
month = aug,
pages = "14",
year = "1996",
URL =          "http://www.maths.man.ac.uk/~nareports/narep294.ps.gz",
abstract = "How can we test the correctness of a computer implementation of
an algorithm such as Gaussian elimination, or the QR algorithm
for the eigenproblem? This is an important question for program
libraries such as LAPACK, that are designed to run on a wide
range of systems. We discuss testing based on verifying known
backward or forward error properties of the algorithms, with
particular reference to the test software in LAPACK. Issues
considered include the choice of bound to verify, computation of
the backward error, and choice of test matrices. Some examples of
bugs in widely used linear algebra software are described.",
keywords = "Mathematical software, testing, linear algebra, backward error,
LAPACK.",
note = "To appear in
{\em The Quality Of Numerical Software: Assessment and Enhancement},
Ronald F. Boisvert, editor,
Chapman and Hall, 1997, ISBN 0-412-80530-8"
}

@techreport{Braconnier:1996:FAF,
author = "Thierry Braconnier",
title = "Fvpspack: {A Fortran} and {PVM} Package to Compute the Field of
Values and Pseudospectra of Large Matrices",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 293",
month = aug,
year = "1996",
pages = "34",
URL =          "http://www.maths.man.ac.uk/~nareports/narep293.ps.gz",
abstract = "The field of values and pseudospectra are tools which yield
insight into the spectral behavior of a matrix. For large sparse
matrices, both sets can be efficiently computed using a Lanczos
type method. Since both computations can be done in a natural
parallel way, we have developed a package including Fortran and
PVM routines. Experiments show that the PVM codes achieve
excellent speed-ups and efficiency.",
keywords = "Field of values, pseudospectra, PVM, parallel algorithm,
speed-up."
}

@techreport{Baker:1996:NAV,
author =       "Christopher T. H. Baker",
title =        "Numerical Analysis of {Volterra} Functional and Integral
Equations---{State} of the Art",
type =         "Numerical Analysis Report",
number =       "No. 292",
institution =  inst-MCCM,
month =        sep,
year =         "1996",
URL =          "http://www.maths.man.ac.uk/~nareports/narep292.ps.gz",
abstract =     "Volterra equations are, broadly speaking, evolutionary
equations incorporating memory. We indicate where they arise
in practice, discuss their significant features (notably
smoothness properties of the solutions), describe numerical
methods for their approximate solution, and indicate the
mathematical foundations of numerical procedures.",
}

@techreport{Smith:1996:IAL,
author = "Antony Smith and David Silvester",
title = "Implicit Algorithms and their Linearisation for the
Transient Incompressible {Navier-Stokes} Equations",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 290",
month = oct,
year = "1996",
note = "Submitted to the IMA Journal of Numerical Analysis",
URL =          "http://www.maths.man.ac.uk/~nareports/narep290.ps.gz",
abstract = "The issue of appropriate time discretisation methods for the
unsteady incompressible Navier-Stokes equations is considered
from a practical perspective here. Conventional implicit
time-stepping algorithms are not feasible for long time
simulations since they inherit the quadratic nonlinearity of the
steady-state equations. As a result linearised versions of the
pure'' algorithms are analysed herein. These have similar
stability properties and comparable accuracy to the underlying
nonlinear methods. "
}

@techreport{Cheng:1996:MCA,
author = "Sheung Hun Cheng and Nicholas J. Higham",
title = "A Modified {Cholesky} Algorithm Based on a Symmetric Indefinite
Factorization",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 289",
pages = "18",
month = apr,
year = "1996",
URL =          "http://www.maths.man.ac.uk/~nareports/narep289.pdf",
abstract = "Given a symmetric and not necessarily positive definite matrix
$A$, a modified Cholesky algorithm computes a Cholesky
factorization $P(A+E)P^T = R^T R$, where $P$ is a permutation
matrix and $E$ is a perturbation chosen to make $A+E$ positive
definite. The aims include producing a suitable small-normed $E$
and making $A+E$ reasonably well conditioned. Modified Cholesky
factorizations are widely used in optimization. We propose a new
modified Cholesky algorithm based on a symmetric indefinite
factorization computed using a new pivoting strategy of Ashcraft,
Grimes and Lewis. We analyze the effectiveness of the algorithm,
both in theory and practice, showing that the algorithm is
competitive with the existing algorithms of Gill, Murray and
Wright, and Schnabel and Eskow. Attractive features of the new
algorithm include easy-to-interpret inequalities that explain the
extent to which it satisfies its design goals, and the fact that
it can be implemented in terms of existing software.",
keywords = "Modified Cholesky factorization, optimization, Newton's method,
symmetric indefinite factorization"
}

@techreport{Higham:1996:RDN,
author = "Nicholas J. Higham",
title = "Recent Developments in Dense Numerical Linear Algebra",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 288",
pages = "26",
month = apr,
year = "1996",
URL =          "http://www.maths.man.ac.uk/~nareports/narep288.ps.gz",
abstract = "We survey recent developments in dense numerical linear algebra,
covering linear systems, least squares problems and
eigenproblems. Topics considered include the design and analysis
of block, partitioned and parallel algorithms, condition number
estimation, componentwise error analysis, and the computation of
practical error bounds. Frequent reference is made to LAPACK, the
state of the art package of Fortran software designed to solve
linear algebra problems efficiently and accurately on
high-performance computers.",
keywords = "Numerical linear algebra, linear system,
least squares problem, eigenvalue problem, LAPACK",
note = "Revised August 1996.
To appear in
{\em The State of the Art in Numerical Analysis},
I. S. Duff and G. A. Watson, editors,
Oxford University Press",
}

@techreport{Blank:1996:NTD,
author = "Luise Blank",
title = "Numerical Treatment of Differential Equations of Fractional Order",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 287",
pages = "18",
month = mar,
year = "1996",
URL = "http://www.maths.man.ac.uk/~nareports/narep287.ps.gz",
abstract = "A brief introduction to derivatives of fractional order is given
and results on the qualitative behaviour of solutions of basic
linear differential equations of fractional order are cited.
The collocation approximation with polynomial splines is applied
to differential equations of fractional order and the systems of
equations characterizing the numerical solution are determined.
In particular, the weight matrices resulting from the fractional
derivative of the spline are deduced and decomposed for
numerical implementation. The main result of this paper is the
simplification of the numerical method under a specific
smoothness condition on the chosen splines. Numerical results
conclude the work and show that the expected qualitative
behaviour is given by collocation approximation.",
keywords = "Differential equations, fractional order, numerical treatment"
}

@techreport{MCCM:1996:AR,
author =       "Manchester Centre for Computational Mathematics",
title =        "Annual Report: {January--December} 1995",
type =         "Numerical Analysis Report",
number =       "No. 286",
institution =  inst-MCCM,
month = mar,
pages = "18",
year =         "1996",
URL =          "http://www.maths.man.ac.uk/~nareports/narep286.ps.gz",
}

@techreport{Ford:1995:PTV,
author = "Neville J. Ford and Christopher T. H. Baker",
title = "Preserving transience in numerical solutions of {Volterra}
integral equations of convolution type",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 284",
pages = "13",
month = nov,
year = "1995",
URL = "http://www.maths.man.ac.uk/~nareports/narep284.ps.gz",
abstract = "We explore the property of transience in solutions of a nonlinear
integral equation of convolution type under discretisation. We
show that, under suitable conditions on the kernel function and
with the choice of a numerical method based on a positive
quadrature, transient solutions are approximated by transient
solutions. We show how the analysis can be extended to a wider
class of kernels",
keywords = "Volterra integral equations, numerical methods, transient solution, stability
}

@techreport{Paul:1995:UGA,
author = "Christopher A.H. Paul",
title = "A User-Guide to {Archi}---An Explicit {Runge-Kutta} Code for Solving
Delay and Neutral Differential Equations",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 283",
pages = "20",
month = dec,
year = "1995",
URL =          "http://www.maths.man.ac.uk/~nareports/narep283.ps.gz",
abstract = "Archi is based on the fifth-order Dormand and Prince explicit
Runge-Kutta method with a fifth-order Hermite interpolant due to
Shampine. The choice of a Hermite interpolant was influenced by:
1. The conclusion of Gladwell, Shampine, Baca and Brankin that
for non-monotonic solutions, quintic Hermite interpolants are
required in order to preserve the shape of a solution.
2. In order to maintain the asymptotic correctness of a $p$-th
order local error estimator, the local error in the approximation
scheme used to evaluate delay terms must be a least order $p+1$.
Archi has also been designed to track the propagation of
discontinuities in a solution using the discontinuity tracking
theory developed by Wille and Baker, and to solve vanishing-lag
DDEs using either extrapolation or the predictor-corrector
approach developed by Baker and Paul. In addition, Archi can
solve a limited class of delay-integro-differential equations,
specifically those equations that have a sufficiently smooth
by Paul.",
keywords = "Delay differential equations, numerical solution, user-guide"
}

@techreport{Blank:1995:STR,
author = "Luise Blank",
title = "Stability-Type Results for Collocation Methods for {Volterra}
Integral Equations",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 282",
pages = "19",
month = dec,
year = "1995",
URL = "http://www.maths.man.ac.uk/~nareports/narep282.ps.gz",
abstract = "Stability behaviour is investigated for polynomial spline
collocation applied to Volterra integral equations with special
emphasis on a weakly singular test equation. The characterization
of the corresponding stability domain is related to stability
results for solutions of finite recursion relations, as they
occur for methods applied to ordinary differential equations.
Estimates of the stability domains are established and, as a
consequence, conditions for $A(\pi/2)$-stability are obtained.
Exploiting these estimates, collocation approximations with
constant, linear, and continuous splines are considered in more
depth and several numerical illustrations are presented.",
keywords = "Stability, Volterra integral equation, weakly singular kernel,
collocation method"
}

@techreport{Braconnier:1995:IOB,
author = "Thierry Braconnier",
title = "Influence of Orthogonality on the Backward Error
and the Stopping Criterion for Krylov Methods",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 281",
month = dec,
year = "1995",
pages = "30",
URL =          "http://www.maths.man.ac.uk/~nareports/narep281.ps.gz",
abstract = "Many algorithms for solving linear systems, least squares
problems or eigenproblems need to compute an orthonormal basis.
The computation is commonly performed using a QR factorization
computed using the classical or the modified Gram-Schmidt
algorithm, the Householder algorithm, the Givens algorithm or the
Gram-Schmidt algorithm with iterative reorthogonalization. For
linear systems, it is well-known that the backward stability of
the process depends crucially on the algorithm used for the QR
factorization. For the eigenproblem, although textbooks warn
users about the possible instability of eigensolvers due to loss
of orthogonality, few theoretical results exist. In this paper we
show that the loss of orthogonality of the computed basis can
affect the reliability of the computed eigenpair. We also show
that the classical residual norm $\norm{Ax-\lambda x}$ and the
specific computed one using a Krylov method can differ because of
the loss of orthogonality of the computed basis of the Krylov
subspace. We give a bound which quantifies the difference between
both residuals in terms of the loss of orthogonality.",
keywords = "QR factorization, backward error, eigenvalues, Krylov methods"
}

@techreport{Braconnier:1995:CFV,
author = "Thierry Braconnier and Nicholas J. Higham",
title = "Computing the Field of Values and Pseudospectra
Using the Lanczos Method with Continuation",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 279",
month = nov,
year = "1995",
pages = "20",
URL =          "http://www.maths.man.ac.uk/~nareports/narep279.ps.gz",
abstract = "The field of values and pseudospectra are useful tools for
understanding the behaviour of various matrix processes. To
compute these subsets of the complex plane it is necessary to
estimate one or two eigenvalues of a large number of parametrized
Hermitian matrices; these computations are prohibitively
expensive for large, possibly sparse, matrices, if done by use of
the QR algorithm. We describe an approach based on the Lanczos
method with selective reorthogonalization and Chebyshev
acceleration that, when combined with continuation and a shift
and invert technique, enables efficient and reliable computation
of the field of values and pseudospectra for large matrices. The
idea of using the Lanczos method with continuation to compute
pseudospectra is not new, but in experiments reported here our
algorithm is faster and more accurate than existing algorithms of
this type.",
keywords = "Field of values, pseudospectra, Lanczos method, Chebyshev
acceleration, continuation."
}

@techreport{Higham:1995:IRL,
author = "Nicholas J. Higham",
title = "Iterative Refinement and {LAPACK}",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 277",
month = sep,
year = "1995",
URL =          "http://www.maths.man.ac.uk/~nareports/narep277.ps.gz",
abstract = "The technique of iterative refinement for improving the computed
solution to a linear system was used on desk calculators and
computers in the 1940s and has remained popular. In the 1990s
iterative refinement is well supported in software libraries,
notably in LAPACK. Although the behaviour of iterative refinement
in floating point arithmetic is reasonably well understood, the
existing theory is not sufficient to justify the use of fixed
precision iterative refinement in all the LAPACK routines in
which it is implemented. We present analysis that provides the
theoretical support needed for LAPACK. The analysis covers both
mixed and fixed precision iterative refinement with an arbitrary
number of iterations, makes only a general assumption on the
underlying solver, and is relatively short. We identify some
remaining open problems.",
keywords = "iterative refinement, rounding error analysis, backward error,
condition number, LAPACK"
}

@techreport{Higham:1995:TMT,
author = "Nicholas J. Higham",
title = "The {Test Matrix Toolbox} for \textsc{Matlab} (Version 3.0)",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 276",
month = sep,
year = "1995",
URL =          "http://www.maths.man.ac.uk/~nareports/narep276.ps.gz",
abstract = "We describe version 3.0 of the Test Matrix Toolbox for
Matlab~4.2. The toolbox contains a collection of test matrices,
routines for visualizing matrices, routines for direct search
optimization, and miscellaneous routines that provide useful
additions to Matlab's existing set of functions. There are 58
parametrized test matrices, which are mostly square, dense,
nonrandom, and of arbitrary dimension. The test matrices include
ones with known inverses or known eigenvalues; ill-conditioned or
rank deficient matrices; and symmetric, positive definite,
orthogonal, defective, involutary, and totally positive matrices.
The visualization routines display surface plots of a matrix and
its (pseudo-) inverse, the field of values, Gershgorin disks, and
two- and three-dimensional views of pseudospectra. The direct
search optimization routines implement the alternating directions
method, the multidirectional search method and the Nelder--Mead
simplex method. We explain the need for collections of test
matrices and summarize the features of the collection in the
toolbox. We give examples of the use of the toolbox and explain
some of the interesting properties of the Frank matrix and magic
square matrices. The leading comment lines from all the toolbox
routines are listed.",
keywords = "test matrix, Matlab, pseudospectrum, visualization, Frank matrix,
magic square matrix, random matrix, direct search optimization"
}

@techreport{MCCM:1995:AR,
author =       "Manchester Centre for Computational Mathematics",
title =        "Annual Report: {January--December} 1994",
type =         "Numerical Analysis Report",
number =       "No. 275",
institution =  inst-MCCM,
year =         "1995",
URL =          "http://www.maths.man.ac.uk/~nareports/narep275.ps.gz",
}

@techreport{Kalogiratou:1995:BCA,
author = "Z. Kalogiratou and Jack Williams",
title = "Best {Chebyshev} Approximation by Families of {ODEs} with a
Fixed Initial Condition",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 274",
month = jul,
year = "1995",
URL =          "http://www.maths.man.ac.uk/~nareports/narep274.ps.gz",
abstract = "Best Chebyshev approximation of real-valued data is treated by
approximations obtained from solutions of parameter dependent
initial value problems in ordinary differential equations, in
which the initial conditions are specified. The problem is a
modified form of the problem treated by the authors in
\cite{JH1}, \cite{JH2}, and \cite{JH3}. The fixing of the initial
condition at $t=a$ requires that approximation be carried out on
a set which excludes this point. Numerical examples illustrating
the theory are presented."
}

@techreport{Silvester:1995:FRS,
author = "David Silvester and Andrew Wathen ",
title = "Fast and Robust Solvers for Time-Discretised
Incompressible {Navier-Stokes} Equations",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 273",
month = jul,
year = "1995",
note = "To appear in proceedings of the 16th Dundee Bienniel
Numerical Analysis Conference",
URL =          "http://www.maths.man.ac.uk/~nareports/narep273.ps.gz",
abstract = "We consider the design of robust and efficient methods for
solving the Navier-Stokes equations governing laminar flow of a
viscous incompressible fluid. Three fundamental issues will be
assessed: the discretisation of the convective transport terms,
the weak enforcement of the incompressibility constraint in a
mixed finite element setting, and the solution of the indefinite
(Stokes-like) systems arising at each time-level. Our aim is to
prescribe a framework for adaptive error control. The essential
ingredients are: an unconditionally stable time-discretisation;
natural'' spatial discretisations which are (inf-sup) stable;
and a fast iterative solution strategy generating iterates
converging monotonically in an appropriate norm (which mimics the
dissipation inherent in the continuous system). A distinguishing
feature of our methodology is the use of multigrid
preconditioning to accelerate the convergence of our Krylov
subspace iterative solver. We motivate this with some analysis
showing that the  contraction rate is bounded away from unity
independently of the choice of the mixed finite element method
and the subdivision parameter. Analysis and implementation of
pure'' multigrid methods seems to be relatively complicated and
(discretisation-) method dependent by comparison. "
}

@techreport{Williams:1995:EUS,
author = "Jack Williams",
title = "Existence and Uniqueness of Solutions
of the Algebraic Equations in the {BDF} Methods",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 272",
month = jul,
year = "1995",
URL =          "http://www.maths.man.ac.uk/~nareports/narep272.ps.gz",
abstract = "The paper presents a new result for the existence and uniqueness
of solutions of the system of nonlinear algebraic equations
obtained from the application of the BDF methods to stiff initial
value problems in ordinary differential equations. Expressed in
terms of the local truncation error, starting errors in the BDF
and contractivity properties of the differential equation, the
result is used to explain how with a certain form of local error
control existence and uniqueness is guaranteed. Although the form
of local error control is a theoretical procedure it relates well
to the error control procedures used in practical codes and
provides insight into how at a typical step the BDF yields
solution values. A potential difficulty in trying to compute low
accuracy numerical solutions is also identified."
}

@techreport{Baker:1995:BNS,
author =       "C. T. H. Baker and C. A. H. Paul and D. R. Wille",
title =        "A Bibliography on the Numerical Solution of Delay
Differential Equations",
type =         "Numerical Analysis Report",
number =       "No. 269",
institution =  inst-MCCM,
month =        jun,
year =         "1995",
URL =          "http://www.maths.man.ac.uk/~nareports/narep269.ps.gz",
abstract =     "The aim of this bibliography is to provide an introduction
to papers and technical reports in the field of delay
differential equations and related differential equations.
In addition to the title, authors and reference of an
article, we provide the abstract which, if the article has
previously appeared as a technical report, comes from the
published paper unless indicated by a +'. The main interest
in this bibliography derives from the references to early
papers and technical reports in the field, as nowadays
on-line search facilities (such as BIDS in the U.K.) provide
to keep the bibliography up-to-date, most immediate effort
will be invested in extending the collection of earlier
references -- as it is far from being exhaustive. The
up-to-date bibliography is only available by anonymous ftp,
since it is hoped to update it regularly.",
keywords =     "Delay differential equations, numerical solution,
stability, applications"
}

@techreport{Hall:1995:ESM,
author =       "George Hall",
title =        "Equilibrium States for Multistep Methods",
type =         "Numerical Analysis Report",
number =       "No. 268",
institution =  inst-MCCM,
month =        jun,
year =         "1995",
URL =          "http://www.maths.man.ac.uk/~nareports/narep268.ps.gz",
abstract =     "When the stepsize in non-stiff ODE codes is restricted by
stability, an uneven pattern of stepsizes with many step
rejections is frequently observed. Results analysing this
behaviour have been obtained for Runge-Kutta methods, leading
to several papers attempting to improve stepsize control. It
is shown here that a similar analysis can be carried out for
multistep methods. The explicit Adams 2-step method is used
to illustrate the techniques required.",
keywords =     "Multistep methods, stepsize selection, stability"
}

@techreport{Baker:1995:PPE,
author =       "Christopher T. H. Baker and Christopher A. H. Paul",
title =        "Pitfalls in Parameter Estimation for Delay Differential
Equations",
type =         "Numerical Analysis Report",
number =       "No. 267",
institution =  inst-MCCM,
month =        may,
year =         "1995",
URL =          "http://www.maths.man.ac.uk/~nareports/narep267.ps.gz",
abstract =     "Given a set of data
$\{U(\gamma_{i}) \approx u(\gamma_{i};p^{*})\}$
corresponding to the delay differential equation (DDE)
$u'(t;p) = f(t,u(t;p),u(\alpha(t;p);p);p)$ for
$t \ge t_{0}(p)$, the basic problem addressed here is that
of calculating the vector $p^{*}$. (We also consider neutral
differential equations.) Most approaches to parameter
estimation calculate $p^{*}$ by minimizing a suitable
objective function that is assumed to be sufficiently smooth.
In this paper, we use the derivative discontinuity tracking
theory for DDEs to analyze how jumps can arise in the
derivatives of a natural objective function. These jumps can
occur when estimating parameters in lag functions and
estimating the position of the initial point, and as such are
not expected to occur in parameter estimation problems for
ordinary differential equations.",
keywords =     "Parameter estimation, delay differential equations,
derivative discontinuities"
}

@techreport{Baker:1995:SND,
author = "C. T. H. Baker and  G. A. Bocharov and C. A. H. Paul and
A. A. Romanyukha",
title = "A Short Note on Delay Effects in Cell Proliferation",
type = "Numerical Analysis Report",
number = "266",
institution =  inst-MCCM,
month = may,
year = "1995",
URL = "http://www.maths.man.ac.uk/~nareports/narep266.ps.gz",
abstract = "We show that a growth model for cell proliferation that
incorporates a time-lag in the cell division phase is more
consistent with certain reported data than the classical
exponential growth model. Although both models provide
estimates of some of the growth characteristics, such as
the population doubling-time, the time-lag growth model
additionally provides estimates of: $(i)$ the cell-doubling
time, $(ii)$ the fraction of the cells that are dividing,
$(iii)$ the rate of commitment to cell division and, $(iv)$
the initial distribution of cells in the cell cycle.",
keywords = "Parameter estimation, delay differential equations,
mathematical modelling"
}

@techreport{Higham:1995:SDP,
author = "Nicholas J. Higham",
title = "Stability of the Diagonal Pivoting Method with Partial Pivoting",
institution = inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 265",
month = jul,
year = "1995",
URL = "http://www.maths.man.ac.uk/~nareports/narep265.ps.gz",
abstract = "LAPACK and LINPACK both solve symmetric indefinite linear systems
using the diagonal pivoting method with the partial pivoting
strategy of Bunch and Kaufman (1977). No proof of the stability
of this method has appeared in the literature. It is tempting to
argue that the diagonal pivoting method is stable for a given
pivoting strategy if the growth factor is small. We show that
this argument is false in general, and give a sufficient
condition for stability. This condition is not satisfied by the
partial pivoting strategy, because the multipliers are unbounded.
Nevertheless, using a more specific approach we are able to prove
the stability of partial pivoting, thereby filling a gap in the
body of theory supporting LAPACK and LINPACK.",
keywords = "symmetric indefinite matrix, diagonal pivoting method,
$\mathrm{LDL^T}$ factorization, partial pivoting, growth factor,
numerical stability, rounding error analysis, LAPACK, LINPACK."
}

@techreport{Elman:1995:FNI,
author = "Howard Elman and David Silvester",
title = "Fast Nonsymmetric Iterations and Preconditioning for
{Navier-Stokes} Equations",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 263",
month = jun,
year = "1995",
note = "To appear in SIAM Journal of Scientific Computing",
keywords = "Preconditioning, Krylov subspace iterations,
Navier-Stokes equations, convection-diffusion problems",
URL =      "http://www.maths.man.ac.uk/~nareports/narep263.ps.gz",
abstract = "Discretization and linearization of the steady-state
Navier-Stokes equations gives rise to a nonsymmetric indefinite
linear system of equations. In this paper, we introduce
preconditioning techniques for such systems with the property
that the eigenvalues of the preconditioned matrices are bounded
independently of the mesh size used in the discretization. We
confirm and supplement these analytic results with a series of
numerical experiments indicating that Krylov subspace iterative
methods for nonsymmetric systems display rates of convergence
that are independent of the mesh parameter. In addition, we show
that preconditioning costs can be kept small by using iterative
methods for some intermediate steps performed by the
preconditioner."
}

@techreport{Silvester:95:SMF,
author = "David Silvester",
title = "Stabilised Mixed Finite Element Methods",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 262",
month = aug,
year = "1995",
note = "To appear in Incompressible Flow and the Finite Element
Method, by P. M. Gresho",
URL = "http://www.maths.man.ac.uk/~nareports/narep262.ps.gz",
abstract = "In this article we review a variety of techniques for stabilising
mixed finite element methods for solving incompressible flow
problems. Our emphasis is on low order finite element
approximations, typically based on piecewise linear velocity with
equal order or piecewise constant pressure. The optimal choice of
the stabilisation/regularisation parameter is discussed in
detail. "
}

@techreport{Ford:1994:QBS,
author = "Neville J. Ford and Christopher T. H. Baker",
title = "Qualitative behaviour and stability of solutions of discretised
non-linear {Volterra} integral equations of convolution type",
institution =  inst-MCCM,
type = "Numerical Analysis Report",
number = "No. 261",
pages = "11",
month = aug,
year = "1994",
URL =          "http://www.maths.man.ac.uk/~nareports/narep261.ps.gz",
abstract = "We provide a discrete analogue of a theorem given by Corduneanu
for a nonlinear Volterra equation of convolution type. We are
able to demonstrate the stability behaviour of some simple
quadrature rules applied to an example equation and show that the
qualitative behaviour of solutions is preserved under the
discretisation for certain families of quadrature rule",
keywords = "Volterra integral equations, numerical methods, stability analysis"
}

@techreport{Paul:1994:PPC,
author =       "Christopher A. H. Paul",
title =        "Performance and Properties of a Class of Parallel Continuous
Explicit Runge-Kutta Methods for Ordinary and Delay
Differential Equations equations",
type =         "Numerical Analysis Report",
number =       "No. 260",
institution =  inst-MCCM,
month =        oct,
year =         "1994",
URL =          "http://www.maths.man.ac.uk/~nareports/narep260.ps.gz",
abstract =     "In this paper we examine the parallel implementation of the
iterated continuous extensions (ICEs) of Baker \& Paul. We
indicate how to construct arbitrarily high-order ICEs, and
discuss some of the strategies for choosing the free'
parameters of the methods. The numerical results that we
present suggest that ICEs implemented in parallel provide an
easy and efficient way of writing dense-output ordinary
differential equation codes and delay differential equation
codes, even in the case that the lag vanishes and/or is
state-dependent.",
keywords =     "Ordinary differential equations, delay differential
equations, parallel algorithms, numerical performance,
stability"
}

@techreport{Griffiths:1994:UME,
author =       "D. F. Griffiths and D. J. Silvester",
title =        "Unstable modes of the {$Q_1$--$P_0$} element",
type =         "Numerical Analysis Report",
number =       "No. 257",
institution =  inst-MCCM,
month =        sep,
year =         "1994",
note =         "Submitted to SIAM J. Numer. Anal.",
URL =          "http://www.maths.man.ac.uk/~nareports/narep257.ps.gz",
abstract =     "In this paper the unstable eigenmodes of $Q_1$--$P_0$
velocity/pressure finite element approximation for
incompressible flow problems are characterised. It is
shown that the inf-sup stability constant is $O(h)$ in
two dimensions and $O(h^2)$ in three dimensions. The
basic tool for our analysis is the method of modified
equations which is applied to finite difference
representations of the underlying finite element
equations. The asymptotic estimates are confirmed and
supplemented by numerical experiments.",
keywords =     "mixed finite elements, inf-sup stability,
incompressible flow"
}

@techreport{Baxter:1994:NEI,
author =       "B. J. C. Baxter and C. A. Micchelli",
title =        "Norm Estimates for the {$l^2$}-Inverses of
Multivariate {Toeplitz} Matrices",
type =         "Numerical Analysis Report",
number =       "No. 255",
institution =  inst-MCCM,
month =        jul,
year =         "1994",
URL =          "http://www.maths.man.ac.uk/~nareports/narep255.ps.gz",
abstract =     "We bound the $l^2$-norms of inverses of certain
Toeplitz matrices arising from P\'olya frequency
functions.",
keywords =     "Radial basis functions, Toeplitz forms, Polya
frequency functions",
}

@techreport{Baxter:1994:RNR,
author =       "B. J. C. Baxter and N. Sivakumar and J. D. Ward",
title =        "Regarding the {$p$}-Norms of Radial Basis
Interpolation Matrices",
type =         "Numerical Analysis Report",
number =       "No. 254",
institution =  inst-MCCM,
month =        jul,
year =         "1994",
URL =          "http://www.maths.man.ac.uk/~nareports/narep254.ps.gz",
abstract =     "Radial basis functions are linear combinations of
translates of a given radially symmetric function. In
this paper we provide some upper bounds on
$\|A^{-1}\|_p$ for $p \ge 1$, where $A = (\phi(x_j - x_k))_{j,k=1}^n$ and $\phi$ is the radial basis
function. We also study some new applications of total
positivity in this context.",
keywords =     "Radial basis functions, p-norms, total positivity",
}

@techreport{Wille:1994:ESC,
author =       "David Richard Will\'e",
title =        "Experiments in Stepsize Control for {Adams} Linear
Multistep Methods",
type =         "Numerical Analysis Report",
number =       "No. 253",
institution =  inst-MCCM,
month =        oct,
year =         "1994",
URL =          "http://www.maths.man.ac.uk/~nareports/narep253.ps.gz",
abstract =     "In this paper we consider stepsize selection in one class
of Adams linear multistep methods for ordinary differential
equations. In particular, the exact form of the local error
for a variable step method is considered and a class of new
direct approximations proposed. The implications of this
approach are then discussed and illustrations provided with
numerical results.  Intended applications include differential
equations with derivative discontinuities and initial
stepsize strategies.",
keywords =     "Adams methods, linear multistep methods, error estimators,
stepsize control",
}

author =       "I. Gladwell and J. Wang and R. M. Thomas",
title =        "Iterations and Predictors for Second Order Problems",
type =         "Numerical Analysis Report",
number =       "No. 252",
institution =  inst-MCCM,
month =        jul,
year =         "1994",
URL =          "http://www.maths.man.ac.uk/~nareports/narep252.ps.gz",
abstract =     "This paper is concerned with the application of
implicit Runge-Kutta methods to initial value problems
consisting of second order ordinary differential
equations without first derivative terms. In
particular, we concentrate on implementing modified
Newton iterations and providing predictors for these
iterations.",
keywords =     "Ordinary differential equations, initial value
problems, implicit Runge-Kutta methods, predictors,
iteration schemes",
}

@techreport{MCCM:1994:AR,
author =       "Manchester Centre for Computational Mathematics",
title =        "Annual Report 1993",
type =         "Numerical Analysis Report",
number =       "No. 251",
institution =  inst-MCCM,
month =        may,
year =         "1994",
URL =          "http://www.maths.man.ac.uk/~nareports/narep251.ps.gz",
}

@techreport{Thomas:1994:FNT,
author =       "R. M. Thomas and T. E. Simos and G. V. Mitsou",
title =        "A Family of {Numerov}-Type Exponentially Fitted
Predictor-Corrector Methods for the Numerical
Integration of the Radial {Schr{\"o}dinger} Equation",
type =         "Numerical Analysis Report",
number =       "No. 249",
institution =  inst-MCCM,
month =        jul,
year =         "1994",
URL =          "http://www.maths.man.ac.uk/~nareports/narep249.ps.gz",
abstract =     "A family of predictor-corrector exponential
Numerov-type methods is developed for the numerical
integration of the one-dimensional Schr{\"o}dinger
equation. The formula considered contains certain free
parameters which allow it to be fitted automatically to
exponential functions. The new methods are very simple
and integrate more exponential functions than both the
well known fourth order Numerov type exponentially
fitted methods and the sixth algebraic order
Runge-Kutta type methods. Numerical results also
indicate that the new methods are much more accurate
than the other exponentially fitted methods mentioned
above.",
keywords =     "Schr{\"o}dinger equation, predictor-corrector methods,
Numerov-type methods, exponentially-fitted methods,
resonance problem",
}

@techreport{Baker:1994:INS,
author =       "Christopher T. H. Baker and Christopher A. H. Paul and
David R. Wille",
title =        "Issues in the Numerical Solution of Evolutionary Delay
Differential Equations",
type =         "Numerical Analysis Report",
number =       "No. 248",
institution =  inst-MCCM,
month =        apr,
year =         "1994",
URL =          "http://www.maths.man.ac.uk/~nareports/narep248.ps.gz",
abstract =     "Delay differential equations are of sufficient
importance in modelling real-life phenomena to merit
the attention of numerical analysts. In this paper, we
discuss key features of delay differential equations
and consider the main issues to be addressed when
constructing robust numerical codes for their solution.
We provide an introduction to the existing literature,
and in particular we indicate the approaches adopted by
the authors at Manchester.",
keywords =     "Numerical solution, delay differential equations",
}

@techreport{Wille:1994:NSE,
author =       "David R. Wille",
title =        "New Stepsize Estimators for Linear Multistep Methods",
type =         "Numerical Analysis Report",
number =       "No. 247",
institution =  inst-MCCM,
month =        mar,
year =         "1994",
URL =          "http://www.maths.man.ac.uk/~nareports/narep247.ps.gz",
abstract =     "The selection of appropriate stepsizes for linear
multistep methods requires in general the solution, or
approximation, of a non-trivial polynomial root-finding
problem. Here we show how this can be reduced to the
integration of a simple ODE. Results are presented for
formulae using both error-per-step' and
error-per-unit-step' controls.",
keywords =     "Linear multistep methods, error estimators, stepsize
control",
}

@techreport{Hanby:1994:CCS,
author =       "R. F. Hanby and D. J. Silvester and J. W. Chew",
title =        "A Comparison of Coupled and Segregated Iterative
Solution Techniques for Incompressible Swirling Flow",
type =         "Numerical Analysis Report",
number =       "No. 246",
institution =  inst-MCCM,
month =        mar,
year =         "1994",
URL =          "http://www.maths.man.ac.uk/~nareports/narep246.ps.gz",
abstract =     "In many popular solution algorithms for the
incompressible Navier-Stokes equations the coupling
between the momentum equations is neglected when the
linearized momentum equations are solved to update the
velocities. This is known to lead to poor convergence
in highly swirling flows where coupling between the
radial and tangential momentum equations is strong.
Here we propose a coupled solution algorithm in which
the linearized momentum and continuity equations are
solved simultaneously. Comparisons between the new
method and the well-known SIMPLEC method are
presented.",
keywords =     "Navier-Stokes equations; finite differences;
unsymmetric linear systems; Krylov subspace methods",
}

@techreport{Hall:1994:NSS,
author =       "George Hall",
title =        "A New Stepsize Strategy for {Runge-Kutta} codes",
type =         "Numerical Analysis Report",
number =       "No. 245",
institution =  inst-MCCM,
month =        mar,
year =         "1994",
URL =          "http://www.maths.man.ac.uk/~nareports/narep245.ps.gz",
abstract =     "When the stepsize in Runge-Kutta codes is restricted
by stability, an uneven pattern of stepsizes with many
step rejections is frequently observed. A modified
strategy is proposed to smooth out this type of
behaviour. Several new estimates for the dominant
eigenvalue of the Jacobian are derived. It is shown
that such estimates can be used, in a strictly
controlled way, to improve the stepsize strategy. Some
numerical evidence is presented to show that the
modified strategy is effective on a set of widely used
test problems.",
keywords =     "Runge-Kutta, stepsize selection, stability",
}

@techreport{Paul:1994:TSF,
author =       "Christopher A. H. Paul",
title =        "A Test Set of Functional Differential Equations",
type =         "Numerical Analysis Report",
number =       "No. 243",
institution =  inst-MCCM,
month =        feb,
year =         "1994",
URL =          "http://www.maths.man.ac.uk/~nareports/narep243.ps.gz",
abstract =     "This document is a collection of delay, neutral and
functional differential equations (mostly with
analytical solutions) taken from the mathematical
literature. Also included is the original source of the
equation, and other useful information. Further test
equations are welcomed, and will appear in the
anonymous ftp version of this document. The reference
solutions quoted were obtained using an explicit
fifth-order continuous Runge-Kutta method based on the
fifth-order Dormand and Prince method with a Hermite
interpolant.",
keywords =     "Test set, functional differential equations",
}

title =        "The {KSR1}---{A} Numerical Analyst's Perspective",
type =         "Numerical Analysis Report",
number =       "No. 242",
institution =  inst-MCCM,
month =        dec,
year =         "1993",
URL =          "http://www.maths.man.ac.uk/~nareports/narep242.ps.gz",
abstract =     "The Kendall Square Research KSR1 is a virtual shared
memory MIMD computer. We give a description of the KSR1
aimed at a numerical analyst who wishes to use the KSR1
for research. The basic architecture of the KSR1 is
described. Parallel constructs and language extensions
provided in KSR Fortran are discussed, and numerical
software issues are also considered. We describe our
practical experiences in using the KSR1 and draw some
conclusions.",
keywords =     "Kendall Square Research KSR1, Virtual shared memory,
ALLCACHE memory system, KSR Fortran, parallel
constructs, LAPACK, BLAS",
}

@techreport{Higham:1994:SCP,
author =       "Nicholas J. Higham",
title =        "A Survey of Componentwise Perturbation Theory in
Numerical Linear Algebra",
type =         "Numerical Analysis Report",
number =       "No. 241",
institution =  inst-MCCM,
month =        feb,
year =         "1994",
note =         "In W.~Gautschi, editor, {\em {Mathematics of Computation}
1943--1993: {A} Half Century of Computational Mathematics},
volume~48 of {\em Proceedings of Symposia in Applied
Mathematics}, pages 49--77. American Mathematical Society,
Providence, RI, USA, 1994.",
URL =          "http://www.maths.man.ac.uk/~nareports/narep241.ps.gz",
abstract =     "Perturbation bounds in numerical linear algebra are
traditionally derived and expressed using norms. Norm
bounds cannot reflect the scaling or sparsity of a
problem and its perturbation, and so can be unduly
weak. If the problem data and its perturbation are
measured componentwise, much smaller and more revealing
bounds can be obtained. A survey is given of
componentwise perturbation theory in numerical linear
algebra, covering linear systems, the matrix inverse,
matrix factorizations, the least squares problem, and
the eigenvalue and singular value problems. Most of the
results described have been published in the last five
years.",
keywords =     "Componentwise perturbation bounds, componentwise
backward error, linear systems, matrix inverse, matrix
factorizations, least squares problem, underdetermined
system, eigenvalue problem, singular value problem",
}

@techreport{Higham:1993:PSV,
author =       "Nicholas J. Higham and Pythagoras Papadimitriou",
title =        "Parallel Singular Value Decomposition via the Polar
Decomposition",
type =         "Numerical Analysis Report",
number =       "No. 239",
institution =  inst-MCCM,
month =        oct,
year =         "1993",
URL =          "http://www.maths.man.ac.uk/~nareports/narep239.ps.gz",
abstract =     "A new method is described for computing the singular
value decomposition (SVD). It begins by computing the
polar decomposition and then computes the spectral
decomposition of the Hermitian polar factor. The method
is particularly attractive for shared memory parallel
computers with a relatively small number of processors,
because the polar decomposition can be computed
efficiently on such machines using an iterative method
developed recently by the authors. This iterative polar
decomposition method requires only matrix
multiplication and matrix inversion kernels for its
implementation and is designed for full rank matrices;
thus the proposed SVD method is intended for matrices
that are not too close to being rank-deficient. On the
Kendall Square KSR1 virtual shared memory computer the
new method is up to six times faster than a
parallelized version of the LAPACK SVD routine,
depending on the condition number of the matrix.",
keywords =     "singular value decomposition, polar decomposition,
numerical stability, LAPACK, level 3 BLAS, Kendall
Square Research KSR1 computer",
}

@techreport{Wille:1994:SID,
author =       "David R. Wille and Christopher T. H. Baker",
title =        "Some Issues in the Detection and Location of
Derivative Discontinuities in Delay Differential
Equations",
type =         "Numerical Analysis Report",
number =       "No. 238",
institution =  inst-MCCM,
month =        mar,
year =         "1994",
URL =          "http://www.maths.man.ac.uk/~nareports/narep238.ps.gz",
abstract =     "The detection and location of derivative
discontinuities is a central issue in the design of
robust solvers for delay differential problems. In this
paper we discuss a number of particular features
associated with one common discontinuity tracking
strategy before addressing a particular problem arising
for a class of state-dependent problems.",
keywords =     "Delay differential equations, derivative
discontinuities",
}

@techreport{Higham:1993:TMT,
author =       "Nicholas J. Higham",
title =        "The {Test Matrix Toolbox} for {MATLAB}",
type =         "Numerical Analysis Report",
number =       "No. 237",
institution =  inst-MCCM,
month =        dec,
year =         "1993",
URL =          "http://www.maths.man.ac.uk/~nareports/narep237.ps.gz",
abstract-1 =   "We describe version 2.0 of the Test Matrix Toolbox for
Matlab 4. The toolbox contains a collection of test
matrices, routines for visualizing matrices, and
miscellaneous routines that provide useful additions to
Matlab's existing set of functions. There are 58
parametrized test matrices, which are mostly square,
dense, nonrandom, and of arbitrary dimension. The test
matrices include ones with known inverses or known
eigenvalues; ill-conditioned or rank deficient
matrices; and symmetric, positive definite, orthogonal,
defective, involutary, and totally positive matrices.",
abstract-2 =   "We describe version 2.0 of the Test Matrix Toolbox for
Matlab The visualization routines display surface plots
of a matrix and its (pseudo-) inverse, the field of
values, Gershgorin disks, and two- and
three-dimensional views of pseudospectra. We explain
the need for collections of test matrices and summarize
the features of the collection in the toolbox. We give
examples of the use of the toolbox and explain some of
the interesting properties of the Frank and Pascal
matrices, and of random, magic square and companion
matrices. The leading comment lines from all the
toolbox routines are listed.",
keywords =     "test matrix, Matlab, pseudospectrum, visualization,
Frank matrix, Pascal matrix, companion matrix, magic
square matrix, random matrix",
note = "*** This report has been withdrawn and is replaced by
Report No. 276. ***"
}

@techreport{Higham:1993:SPT,
author =       "Nicholas J. Higham",
title =        "Stability of Parallel Triangular System Solvers",
type =         "Numerical Analysis Report",
number =       "No. 236",
institution =  inst-MCCM,
month =        jul,
year =         "1993",
URL =          "http://www.maths.man.ac.uk/~nareports/narep236.ps.gz",
abstract =     "Several parallel algorithms have been proposed for
solution of triangular systems. The stability of four
of them is analysed here: a fan-in algorithm, a block
elimination method, a method based on a factorized
power series expansion of the matrix inverse, and a
method based on a divide and conquer matrix inversion
technique. We derive new forward error and residual
bounds, including an improvement on the bounds of Sameh
and Brent for the fan-in algorithm. We identify a
forward error bound that holds not only for all the
methods described here, but for any triangular equation
solver that does not rely on algebraic cancellation;
among the implications of the bound is that any such
method is extremely accurate for certain special types
of triangular system.",
keywords =     "triangular system, matrix inversion, parallel
algorithms, fan-in operation, numerical stability,
rounding error analysis",
}

@techreport{Higham:1993:MPFb,
author =       "Nicholas J. Higham and Philip A. Knight",
title =        "Matrix Powers in Finite Precision Arithmetic",
type =         "Numerical Analysis Report",
number =       "No. 235",
institution =  inst-MCCM,
month =        aug,
year =         "1993",
URL =          "http://www.maths.man.ac.uk/~nareports/narep235.ps.gz",
abstract =     "If $A$ is a square matrix with spectral radius less
than 1 then $A^k \to 0$ as $k \to \infty$, but the
powers computed in finite precision arithmetic may or
may not converge. We derive a sufficient condition for
$fl(A^k) \to 0$ as $k\to\infty$ and a bound on
$\norm{fl(A^k)}$, both expressed in terms of the Jordan
canonical form of $A$. Examples show that the results
can be sharp. We show that the sufficient condition can
be rephrased in terms of a pseudospectrum of $A$ when
$A$ is diagonalizable, under certain assumptions. Our
analysis leads to the rule of thumb that convergence or
divergence of the computed powers of $A$ can be
expected according as the spectral radius computed by
any backward stable algorithm is less than or greater
than 1.",
keywords =     "matrix powers, rounding errors, Jordan canonical form,
nonnormal matrices, pseudospectrum",
}

@techreport{MCCM:1993:AR,
author =       "Manchester Centre for Computational Mathematics",
title =        "Annual Report 1992",
type =         "Numerical Analysis Report",
number =       "No. 233",
institution =  inst-MCCM,
month =        jun,
year =         "1993",
URL =          "http://www.maths.man.ac.uk/~nareports/narep233.pdf",
}

@techreport{Paul:1993:FCP,
author =       "Christopher A. H. Paul",
title =        "A {FORTRAN} 77 Code for Plotting Stability Regions",
type =         "Numerical Analysis Report",
number =       "No. 230",
institution =  inst-MCCM,
month =        may,
year =         "1993",
URL =          "http://www.maths.man.ac.uk/~nareports/narep230.ps.gz",
abstract =     "The standard technique for obtaining the stability
regions of numerical methods for ordinary differential
equations (ODEs) and delay differential equations
(DDEs) is the boundary-locus algorithm (BLA). However,
in the case of the DDE $y'(t) = \lambda y(t) + \mu y(t-\tau)$ for $\lambda$ and $\mu$ real, the BLA often
fails to map out the stability region correctly. In
this paper we give a FORTRAN 77 listing of an
alternative stability boundary algorithm, which can be
used for both implicit and explicit numerical methods
applied to both ODEs and DDEs.",
keywords =     "stability regions, boundary-locus algorithm,
Runge-Kutta methods",
}

@techreport{Baker:1993:GCT,
author = "Christopher T. H. Baker and Christopher A. H. Paul",
title = "A Global Convergence Theorem for a Class of Parallel Continuous
Explicit Runge-Kutta Methods and Vanishing Lag Delay Differential
Equations",
institution = "University of Manchester",
type = "Numerical Analysis Report",
number = "No. 229",
month = may,
year = "1993",
keywords = "parallel continuous explicit Runge-Kutta methods, iterated
continuous extensions, delay differential equations, vanishing
lag",
URL = "http://www.maths.man.ac.uk/~nareports/narep229.ps.gz",
abstract = "Iterated continuous extensions (ICEs) are continuous explicit
Runge-Kutta methods developed for the numerical solution of
evolutionary problems in ordinary and delay differential
equations (DDEs). ICEs have a particular r\^{o}le in the
explicit solution of DDEs with vanishing lags. They may be
regarded as parallel continuous explicit Runge-Kutta (PCERK)
methods, as they allow advantage to be taken of parallel
architectures. ICEs can also be related to a collocation
method.
The purpose of this paper is to provide a theorem giving the
global order of convergence for variable-step implementations
of ICEs applied to state-dependent DDEs with and without
vanishing lags. Implications of the theoryfor the
implementation of this class of methods are discussed and
demonstrated. The results establish that our approach allows
the construction of PCERK methodswhose order of convergence is
restricted only by the continuity of the solution."
}

@techreport{Higham:1993:PAC,
author =       "Nicholas J. Higham and Pythagoras Papadimitriou",
title =        "A Parallel Algorithm for Computing the Polar
Decomposition",
type =         "Numerical Analysis Report",
number =       "No. 226",
institution =  inst-MCCM,
month =        jun,
year =         "1993",
URL =          "http://www.maths.man.ac.uk/~nareports/narep226.ps.gz",
abstract =     "The polar decomposition $A=UH$ of a rectangular matrix
$A$, where $U$ is unitary and $H$ is Hermitian positive
semidefinite, is an important tool in various
applications, including aerospace computations, factor
analysis and signal processing. We consider a $p$th
order iteration for computing $U$ that involves $p$
independent matrix inversions per step and which is
hence very amenable to parallel computation. We show
that scaling the iterates speeds convergence of the
iteration but makes the iteration only conditionally
stable, with the backward error typically $\kappa_2(A)$
times bigger than the unit roundoff. In our
implementation of the iteration on the Kendall Square
Research KSR1 virtual shared memory MIMD computer we
take $p$ to be the number of processors ($p \le 16$ in
our experiments). Our code is found to be significantly
faster than two existing techniques for computing the
polar decomposition: one a Newton iteration, the other
based on the singular value decomposition.",
keywords =     "polar decomposition, singular value decomposition,
numerical stability, LAPACK, level 3 BLAS, Kendall
Square Research KSR1 computer",
}

@techreport{Higham:1993:MSD,
author =       "Nicholas J. Higham",
title =        "The Matrix Sign Decomposition and Its Relation to the
Polar Decomposition",
type =         "Numerical Analysis Report",
number =       "No. 225",
institution =  inst-MCCM,
month =        apr,
year =         "1993",
URL =          "http://www.maths.man.ac.uk/~nareports/narep225.ps.gz",
abstract-1 =   "The sign function of a square matrix was introduced by
Roberts in 1971. We show that it is useful to regard
$S=\sign(A)$ as being part of a matrix sign
decomposition $A=SN$, where $N = (A^2)^{1/2}$. This
decomposition leads to the new representation $\sign(A) = A(A^2)^{-1/2}$. Most results for the matrix sign
decomposition have a counterpart for the polar
decomposition $A=UH$, and vice versa.",
abstract-2 =   "To illustrate this, we derive best approximation
properties of the factors $U$, $H$ and~$S$, determine
bounds for $\norm{A-S}$ and $\norm{A-U}$, and describe
integral formulas for $S$ and $U$. We also derive
explicit expressions for the condition numbers of the
factors $S$ and $N$. An important equation expresses
the sign of a block $2\times2$ matrix involving $A$ in
terms of the polar factor $U$ of $A$. We apply this
equation to a family of iterations for computing $S$ by
Pandey, Kenney and Laub, to obtain a new family of
iterations for computing $U$. The iterations have some
attractive properties, including suitability for
parallel computation.",
keywords =     "matrix sign function, polar decomposition,
conditioning, best approximation, properties, parallel
iteration",
}

@techreport{Higham:1992:SPI,
author =       "Nicholas J. Higham and Alex Pothen",
title =        "Stability of the Partitioned Inverse Method for
Parallel Solution of Sparse Triangular Systems",
type =         "Numerical Analysis Report",
number =       "No. 222",
institution =  inst-U-MANCHESTER,
month =        oct,
year =         "1992",
note =         "To appear in {\em SIAM J. Sci. Comput.}.",
URL =          "http://www.maths.man.ac.uk/~nareports/narep222.ps.gz",
abstract =     "Several authors have recently considered a parallel
method for solving sparse triangular systems with many
right-hand sides. The method employs a partition into
sparse factors of the product form of the inverse of
the coefficient matrix. It is shown here that while the
method can be unstable, stability is guaranteed if a
certain scalar that depends on the matrix and the
partition is small, and that this scalar is small when
the matrix is well-conditioned. Moreover, when the
partition is chosen so that the factors have the same
sparsity structure as the coefficient matrix, the
backward error matrix can be taken to be sparse.",
keywords =     "sparse matrix, triangular system, substitution
algorithm, parallel algorithm, rounding error analysis,
matrix inverse",
}

@Article{Higham:1993:FPB,
author =       "Nicholas J. Higham and Philip A. Knight",
title =        "Finite Precision Behavior of Stationary Iteration for
Solving Singular Systems",
journal =      "Linear Algebra and Appl.",
volume =       "192",
pages =        "165--186",
year =         "1993",
URL =          "http://www.maths.man.ac.uk/~nareports/narep215.ps.gz",
abstract-1 =   "A stationary iterative method for solving a singular
system $Ax=b$ converges for any starting vector if
$\lim_{i\to\infty}G^i$ exists, where $G$ is the
iteration matrix, and the solution to which it
converges depends on the starting vector. We examine
the behaviour of stationary iteration in finite
precision arithmetic. A perturbation bound for singular
systems is derived and used to define forward stability
of a numerical method. A rounding error analysis
enables us to deduce conditions under which a
stationary iterative method is forward stable or
backward stable.",
abstract-2 =   "A stationary iterative method for solving a singular
system The component of the forward error in the null
space of $A$ can grow linearly with the number of
iterations, but it is innocuous as long as the
iteration converges reasonably quickly. As special
cases, we show that when $A$ is symmetric positive
semi-definite the Richardson iteration with optimal
parameter is forward stable and, if $A$ also has unit
diagonal and property A, then the Gauss-Seidel method
is both forward and backward stable. Two numerical
examples are given to illustrate the analysis.",
keywords =     "stationary iteration, singular system, Drazin inverse,
Jacobi method, Gauss-Seidel method, successive
over-relaxation, Richardson iteration, error analysis,
numerical stability, Neumann problem",
mynote =       "Special issue for Proceedings of Workshop on
Computational Linear Algebra in Algebraic and Related
Problems",
}

@techreport{Paul:1992:FEV,
author =       "Christopher A. H. Paul",
title =        "A fast, efficient, very low storage, adaptive
type =         "Numerical Analysis Report",
number =       "No. 213",
institution =  inst-U-MANCHESTER,
month =        may,
year =         "1992",
URL =          "http://www.maths.man.ac.uk/~nareports/narep213.ps.gz",
abstract =     "In this report we review the types of quadrature
scheme that are available for the evaluation of
which avoids the storage requirements of normal
adaptive schemes. Finally we provide a comparison of
all the types of quadrature scheme considered for a set
of test integrals.",
}

@techreport{Paul:1992:ERK,
author =       "Christopher A. H. Paul and Christopher T. H. Baker",
title =        "Explicit {Runge-Kutta} Methods for the Numerical
Solution of Singular Delay Differential Equations",
type =         "Numerical Analysis Report",
number =       "No. 212",
institution =  inst-U-MANCHESTER,
month =        apr,
year =         "1992",
URL =          "http://www.maths.man.ac.uk/~nareports/narep212.ps.gz",
abstract =     "In this paper we are concerned with the development of
an explicit Runge-Kutta scheme for the numerical
solution of delay differential equations (DDEs) where
one or more delay lies in the current Runge-Kutta
interval. The scheme presented is also applicable to
the numerical solution of {\em Volterra functional
equations} (VFEs), although the theory is not covered
in this paper. We also derive the stability equations
of the scheme for the ODE $y'(t) = \lambda y(t)$, and
the DDE $y'(t) = \lambda y(t) + \mu y(t-\tau)$, where
the delay $\tau$ and the Runge-Kutta stepsize $H$ are
both constant. In the case of the DDE, we consider the
two distinct cases: $(i)$ $\tau \ge H$, and $(ii)$
$\tau < H$. We evaluate the performance of the scheme
by solving several types of {\em singular} DDE and a
VFE.",
keywords =     "explicit Runge-Kutta, singular delay differential
equations",
}

@Article{Higham:1993:PTB,
author = "Nicholas J. Higham",
title = "Perturbation Theory and Backward Error for {$AX-XB=C$}",
journal = "BIT",
volume = 33,
pages = "124-136",
year = "1993",
URL = "http://www.maths.man.ac.uk/~nareports/narep211.ps.gz",
abstract = "Because of the special structure of the equations $AX-XB=C$ the
usual relation for linear equations $\mbox{backward error} = \mbox{relative residual}$'' does not hold, and application of the
standard perturbation result for $Ax=b$ yields a perturbation
bound involving ${\rm sep}(A,B)^{-1}$ that is not always
attainable.  An expression is derived for the backward error of
an approximate solution $Y$; it shows that the backward error can
exceed the relative residual by an arbitrary factor.  A sharp
perturbation bound is derived and it is shown that the condition
number it defines can be arbitrarily smaller than the ${\rm sep}(A,B)^{-1}$-based quantity that is usually used to measure
sensitivity.  For practical error estimation using the residual
of a computed solution an LAPACK-style'' bound is shown to be
efficiently computable and potentially much smaller than a
sep-based bound. A Fortran~77 code has been written that solves
the Sylvester equation and computes this bound, making use of
LAPACK routines.",
keywords = "Sylvester equation, Lyapunov equation, backward error,
perturbation bound, condition number, error estimate, LAPACK."
}

@techreport{Baker:1992:ESA,
author =       "Christopher T. H. Baker and John C. Butcher and
Christopher A. H. Paul",
title =        "Experience of {STRIDE} applied to delay differential
equations",
type =         "Numerical Analysis Report",
number =       "No. 208",
institution =  inst-U-MANCHESTER,
month =        feb,
year =         "1992",
URL =          "http://www.maths.man.ac.uk/~nareports/narep208.ps.gz",
abstract =     "STRIDE is intended as a robust adaptive code for
solving {\em initial value problems} for {\em ordinary
differential equations} (ODEs). The acronym STRIDE
stands for {\bf ST}able {\bf R}unge-Kutta {\bf
I}ntegrator for {\bf D}ifferential {\bf E}quations. Our
purpose here is to report on its adaptation for the
numerical solution of a set of delay and neutral
differential equations.",
keywords =     "STRIDE, delay and neutral differential equations",
}

@InProceedings{Higham:1993:MPFa,
author =       "Nicholas J. Higham and Philip A. Knight",
editor =       "Carl D. Meyer and Robert J. Plemmons",
booktitle =    "Linear Algebra, Markov Chains, and Queueing Models",
title =        "Componentwise Error Analysis for Stationary Iterative
Methods",
volume =       "48",
publisher =    pub-SV,
pages =        "29--46",
year =         "1993",
series =       "IMA Volumes in Mathematics and its Applications",
URL =          "http://www.maths.man.ac.uk/~nareports/narep206.ps.gz",
asbtract =     "How small can a stationary iterative method for
solving a linear system $Ax=b$ make the error and the
residual in the presence of rounding errors? We give a
componentwise error analysis that provides an answer to
this question and we examine the implications for
numerical stability. The Jacobi, Gauss-Seidel and
successive over-relaxation methods are all found to be
forward stable in a componentwise sense and backward
stable in a normwise sense, provided certain conditions
are satisfied that involve the matrix, its splitting,
and the computed iterates. We show that the stronger
property of componentwise backward stability can be
achieved using one step of iterative refinement in
fixed precision, under suitable assumptions.",
keywords =     "stationary iteration, Jacobi method, Gauss-Seidel
method, successive over-relaxation, error analysis,
numerical stability",
}

@techReport{Paul:1991:SBR,
author =       "Christopher A. H. Paul and Christopher T. H. Baker",
title =        "Stability boundaries revisited --- {Runge-Kutta}
methods for delay differential equations",
type =         "Numerical Analysis Report",
number =       "No. 205",
institution =  inst-U-MANCHESTER,
month =        dec,
year =         "1991",
URL =          "http://www.maths.man.ac.uk/~nareports/narep205.ps.gz",
abstract =     "We discuss the practical determination of stability
regions when various Runge-Kutta methods, combined with
continuous extensions, are applied with a fixed
stepsize $h$ to the linear delay differential test
equation $y'(t) = \lambda y(t) +\mu y(t-\tau) \quad (t \ge 0)$, with fixed delay $\tau$. The delay $\tau$ is
not limited to an integer multiple of the stepsize
$h$.",
keywords =     "delay differential equation, continuous extension,
stability",
}

@techReport{Paul:1991:DDD,
author =       "Christopher A. H. Paul",
title =        "Developing a Delay Differential Equation Solver",
type =         "Numerical Analysis Report",
number =       "No. 204",
institution =  inst-U-MANCHESTER,
month =        dec,
year =         "1991",
URL =          "http://www.maths.man.ac.uk/~nareports/narep204.ps.gz",
abstract =     "We discuss briefly various phenomena to be tested when
designing a robust code for delay differential
equations: choice of interpolant; tracking of
discontinuities; vanishing delays; and problems arising
from floating-point arithmetic.",
keywords =     "delay differential equations, interpolation,
discontinuities, vanishing delay, real arithmetic",
}

@InProceedings{Higham:1990:HAG,
author =       "Nicholas J. Higham",
editor =       "D. F. Griffiths and G. A. Watson",
booktitle =    "Numerical Analysis 1989, Proceedings of the 13th
Dundee Conference",
title =        "How Accurate is {Gaussian} Elimination?",
volume =       "228",
publisher =    pub-LONGMAN,
pages =        "137--154",
year =         "1990",
series =       "Pitman Research Notes in Mathematics",
URL =          "http://www.maths.man.ac.uk/~nareports/narep173.ps.gz",
abstract =     "J. H. Wilkinson put Gaussian elimination (GE) on a
sound numerical footing in the 1960s when he showed
that with partial pivoting the method is stable in the
sense of yielding a small backward error. He also
derived bounds proportional to the condition number
$\kappa(A)$ for the forward error $\|x-\widehat x\|$,
where $\widehat x$ is the computed solution to $Ax=b$.
More recent work has furthered our understanding of GE,
largely through the use of componentwise rather than
normwise analysis. We survey what is known about the
accuracy of GE in both the forward and backward error
senses. Particular topics include: classes of matrix
for which it is advantageous not to pivot; how to
estimate or compute the backward error; iterative
refinement in single precision; and how to compute
efficiently a bound on the forward error.",
keywords =     "Gaussian elimination, partial pivoting, rounding error
analysis, backward error, forward error, condition
number, iterative refinement in single precision,
growth factor, componentwise bounds, condition
estimator",
}

@InProceedings{Higham:1990:ACD,
author =       "Nicholas J. Higham",
editor =       "M. G. Cox and S. J. Hammarling",
booktitle =    "Reliable Numerical Computation",
title =        "Analysis of the {Cholesky} Decomposition of a
Semi-definite Matrix",
publisher =    pub-OUP,
pages =        "161--185",
year =         "1990",
URL =          "http://www.maths.man.ac.uk/~nareports/narep128.ps.gz",
abstract =     "Perturbation theory is developed for the \Chol\ of an
$n \times n$ symmetric \psd\ matrix $A$ of rank~$r$.
The matrix $W=\All^{-1}\A{12}$ is found to play a key
role in the \pert\ bounds, where $\All$ and $\A{12}$
are $r \times r$ and $r \times (n-r)$ submatrices of
$A$ respectively. A backward error analysis is given;
it shows that the computed Cholesky factors are the
exact ones of a matrix whose distance from $A$ is
bounded by
$4r(r+1)\bigl(\norm{W}+1\bigr)^2u\norm{A}+O(u^2)$,
where $u$ is the unit roundoff. For the \cpiv\ strategy
it is shown that $\norm{W}^2 \le {1 \over 3}(n-r)(4^r- 1)$, and empirical evidence that $\norm{W}$ is usually
small is presented. The overall conclusion is that the
\Cholalg\ with \cpiv\ is stable for semi-definite
matrices. Similar \pert\ results are derived for the
\QR\ with \colpiv\ and for the LU decomposition with
\cpiv. The results give new insight into the
reliability of these decompositions in rank
estimation.",
keywords =     "Cholesky decomposition, positive semi-definite matrix,
perturbation theory, backward error analysis, QR
decomposition, rank estimation, LINPACK",
}