%%% ==================================================================== %%% BibTeX-file{ %%% author = "Nicholas J. Higham", %%% version = "1.00", %%% date = "March 10, 2014", %%% 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 = "nick.higham@manchester.ac.uk", %%% 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.manchester.ac.uk/~higham/narep/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"} @String{inst-MCCM:adr = "Manchester, England"} %%% 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-LONGMAN:adr = "Essex, UK"} @String{pub-OUP = "Oxford University Press"} @String{pub-OUP:adr = "Oxford, UK"} @String{pub-SV = "Spring{\-}er-Ver{\-}lag"} @String{pub-SV:adr = "Berlin, Germany~/ Heidelberg, 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, address = inst-MCCM:adr, type = t-NAREP, number = "No. 475", pages = "17", month = oct, year = "2005", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = t-NAREP, number = "No. 474", month = oct, year = 2005, note = "Submitted to ACM Transactions on Mathematical Software", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = t-NAREP, number = "No. 473", pages = "21", month = sep, year = "2005", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, pages = "25", month = sep, year = "2005", URL = "http://www.maths.manchester.ac.uk/~higham/narep/narep472.pdf" } @techreport{Bhatia:2005:IDM, author = "Rajendra Bhatia", title = "Infinitely Divisible Matrices", institution = inst-MCCM, address = inst-MCCM:adr, type = t-NAREP, number = "No. 471", pages = "36", month = aug, year = "2005", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = t-NAREP, number = "No. 470", pages = "6", month = jun, year = "2005", URL = "http://www.maths.manchester.ac.uk/~higham/narep/narep470.pdf" } @techreport{Berhanu:05:PBH, author = "Michael Berhanu", title = "Perturbation Bounds for Hyperbolic Matrix Factorizations", institution = inst-MCCM, address = inst-MCCM:adr, type = t-NAREP, number = "No. 469", pages = "24", month = jun, year = "2005", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = t-NAREP, number = "No. 468", pages = "10", month = jun, year = "2005", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = t-NAREP, number = "No. 467", pages = "16", month = apr, year = "2005", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = t-NAREP, number = "No. 466", pages = 26, month = apr, year = "2005", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = t-NAREP, number = "No. 465", pages = "25", month = apr, year = "2005", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = t-NAREP, number = "No. 464", pages = "32", month = apr, year = "2005", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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 innitely 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, address = inst-MCCM:adr, type = t-NAREP, number = "No. 463", pages = "77", month = mar, year = "2005", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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 matrices which are selfadjoint or skew-adjoint in a regular 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, address = inst-MCCM:adr, type = t-NAREP, number = "No. 462", pages = "24", month = mar, year = "2005", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = t-NAREP, number = "461", pages = "17", month = feb, year = "2005", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, Taylor series, Pad\'e approximation, Pad\'e approximant, 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, address = inst-MCCM:adr, type = t-NAREP, number = "No. 460", pages = "45", month = jan, year = "2005", URL = "http://www.maths.manchester.ac.uk/~higham/narep/narep460.pdf", abstract = "After a brief historical review and an account of the canonical forms attributed to Jordan and Kronecker, a systematic development is made of the simultaneous reduction of pairs of quadratic forms 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, address = inst-MCCM:adr, type = t-NAREP, number = "No. 459", pages = "23", month = jan, year = "2005", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = t-NAREP, number = "No. 458", pages = "17", month = jan, year = "2005", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = t-NAREP, number = "No. 457", pages = 51, month = oct, year = 2004, URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = t-NAREP, number = "No. 456", pages = 29, month = oct, year = 2004, URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = t-NAREP, number = "No. 455", pages = 12, month = oct, year = 2004, URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = t-NAREP, number = "No. 454", pages = 26, month = aug, year = 2004, URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = t-NAREP, number = "No. 453", month = aug, year = 2004, URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = t-NAREP, number = "No. 452", pages = "15", month = jul, year = "2004", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = t-NAREP, number = "No. 450", pages = 38, month = jun, year = 2004, url = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, pages = "25", month = may, year = "2004", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = t-NAREP, number = "No. 448", pages = 13, month = apr, year = 2004, url = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = t-NAREP, number = "No. 447", pages = 19, month = apr, year = 2004, url = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = t-NAREP, number = "446", pages = 29, month = mar, year = "2004", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = t-NAREP, number = "No. 445", pages = "17", month = feb, year = "2004", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = t-NAREP, number = "No. 444", pages = 27, month = mar, year = 2004, url = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = t-NAREP, number = "No. 443", pages = 29, month = mar, year = 2004, url = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = t-NAREP, number = "No. 442", pages = "31", month = feb, year = "2004", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = t-NAREP, number = "No. 441", pages = "33", month = jan, year = "2004", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = t-NAREP, number = "No. 440", pages = "10", month = dec, year = "2003", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = t-NAREP, number = "No. 439", month = dec, year = "2003", URL = "http://www.maths.manchester.ac.uk/~higham/narep/narep439.pdf" } @techreport{Shardlow:2003:NWE, author = "Shardlow, Tony", title = "Nucleation of Waves in Excitable Media by Noise", institution = inst-MCCM, address = inst-MCCM:adr, type = t-NAREP, number = "No. 438", month = dec, year = "2003", URL = "http://www.maths.manchester.ac.uk/~higham/narep/narep438.pdf" } @techreport{Shardlow:2003:NSS, author = "Shardlow, Tony", title = "Numerical Simulation of Stochastic {PDE}s for Excitable Media", institution = inst-MCCM, address = inst-MCCM:adr, type = t-NAREP, number = "No. 437", month = dec, year = "2003", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = t-NAREP, number = "No. 436", pages = "10", month = nov, year = "2003", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = t-NAREP, number = "No. 435", month = oct, year = "2003", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = t-NAREP, number = "No. 434", month = oct, pages = "11", year = "2003", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = t-NAREP, number = "No. 433", pages = "6", month = sep, year = "2003", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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:2004: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, address = inst-MCCM:adr, pages = "", month = nov, year = "2004", URL = "http://www.maths.manchester.ac.uk/~higham/narep/narep432.pdf", } @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, address = inst-MCCM:adr, type = t-NAREP, number = "No. 431", pages = 25, month = jan, year = 2004, url = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = t-NAREP, number = "No. 430", pages = 21, month = jan, year = 2004, url = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, pages = "19", month = jun, year = "2003", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, 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, address = inst-MCCM:adr, type = t-NAREP, number = "No. 427", pages = "30", month = may, year = "2003", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = t-NAREP, number = "No. 426", pages = "15", month = apr, year = "2003", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = t-NAREP, number = "No. 425", pages = "28", month = feb, year = "2003", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = t-NAREP, number = "No. 424", pages = "32", month = mar, year = "2003", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = t-NAREP, number = "No. 422", pages = "12", month = feb, year = "2003", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 421", pages = "30", month = Nov, year = "2004", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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 product, orthosymmetric, adjoint, factorization, symplectic, 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, address = inst-MCCM:adr, type = t-NAREP, number = "No. 420", pages = "24", month = feb, year = "2003", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = t-NAREP, number = "No. 419", pages = "36", month = feb, year = "2003", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = t-NAREP, number = "No. 418", month = jan, year = 2003, URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = t-NAREP, number = "No. 417", pages = 35, month = Dec, year = 2003, url = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = t-NAREP, number = "No. 416", pages = 49, month = dec, year = 2002, url = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, 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.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = t-NAREP, number = "No. 414", pages = 20, month = Nov, year = 2003, url = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = t-NAREP, number = "No. 413", pages = "11", month = oct, year = "2002", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = t-NAREP, number = "No. 412", month = sep, year = 2002, URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = t-NAREP, number = "No. 411", pages = 15, month = sep, year = 2002, URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = t-NAREP, number = "No. 410", pages = "19", month = aug, year = "2002", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = t-NAREP, number = "No. 409", pages = "20", month = sep, year = "2002", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = t-NAREP, number = "No. 408", pages = "15", month = sep, year = 2002, URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 407", pages = "18", month = aug, year = "2002", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, symmetric quadratic eigenvalue problem." } @techreport{Mihajlovic:2002:EPS, author = "Milan Mihajlovic and David Silvester", title = "Efficient Parallel Solvers for the Biharmonic Equation", institution = inst-MCCM, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 406", month = jul, year = 2002, note = "Submitted to Parallel Computing", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = t-NAREP, number = "No. 405", pages = "30", month = aug, year = "2002", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = t-NAREP, number = "No. 404", pages = "22", month = jul, year = "2002", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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 radius of convergence, but it can be adapted for certain other 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, address = inst-MCCM:adr, type = t-NAREP, number = "No. 403", pages = "18", month = jun, year = "2002", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, 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, address = inst-MCCM:adr, pages = "21", month = may, year = "2002", URL = "http://www.maths.manchester.ac.uk/~higham/narep/narep401.pdf", } @techreport{Baker:2002:PSD, author = "C. T. H. Baker and Y. Song", title = "Periodic Solutions of Discrete {Volterra} Equations", institution = inst-MCCM, address = inst-MCCM:adr, 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, address = inst-MCCM:adr, 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.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, 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, address = inst-MCCM:adr, 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, address = inst-MCCM:adr, 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, address = inst-MCCM:adr, 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, address = inst-MCCM:adr, 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, address = inst-MCCM:adr, 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, address = inst-MCCM:adr, 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, address = inst-MCCM:adr, pages = "23", month = mar, year = "2001", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, 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, address = inst-MCCM:adr, 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, address = inst-MCCM:adr, 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, address = inst-MCCM:adr, 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, address = inst-MCCM:adr, 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, address = inst-MCCM:adr, 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, address = inst-MCCM:adr, 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, address = inst-MCCM:adr, 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, address = inst-MCCM:adr, 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, address = inst-MCCM:adr, 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, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 379", month = feb, pages = 14, year = 2001, URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 378", month = feb, pages = 12, year = 2001, note = "Submitted to BIT", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, 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.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, 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, address = inst-MCCM:adr, month = oct, year = "2001", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, 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, address = inst-MCCM:adr, 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, address = inst-MCCM:adr, 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, address = inst-MCCM:adr, 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, numerical radius" } @techreport{Tisseur:2000:SQEP, author = "Fran\c{c}oise Tisseur and Karl Meerbergen", title = "The Quadratic Eigenvalue Problem", institution = inst-MCCM, address = inst-MCCM:adr, 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, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 369", month = oct, pages = 14, year = 2000, note = "Revised July 2001", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, 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.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, 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.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, pages = "24", month = mar, year = "2001", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, 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, quadratic eigenvalue problem,condition number" } @techreport{Silvester:00:EPB, author = "David Silvester and Milan Mihajlovic", title = "A black-box multigrid preconditioner for the biharmonic equation", institution = inst-MCCM, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 362", month = aug, year = 2002, note = "Submitted to BIT", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 361", month = may, year = 2000, URL = "http://www.maths.manchester.ac.uk/~higham/narep/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.manchester.ac.uk/~higham/narep/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.manchester.ac.uk/~higham/narep/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.manchester.ac.uk/~higham/narep/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.manchester.ac.uk/~higham/narep/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.manchester.ac.uk/~higham/narep/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.manchester.ac.uk/~higham/narep/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 $1$s 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.manchester.ac.uk/~higham/narep/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.manchester.ac.uk/~higham/narep/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.manchester.ac.uk/~higham/narep/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.manchester.ac.uk/~higham/narep/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.manchester.ac.uk/~higham/narep/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.manchester.ac.uk/~higham/narep/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.manchester.ac.uk/~higham/narep/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.manchester.ac.uk/~higham/narep/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.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, pages = "24", month = mar, year = "2000", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 343", month = aug, pages = 45, year = 1999, URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 342", pages = "14", month = sep, year = "1999", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 341", month = apr, pages = 32, year = 1999, URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 340", pages = "17", month = feb, year = "1999", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 339", month = jan, pages = 18, year = 1999, URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 338", month = jan, pages = "13", year = 1999, URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 337", month = jan, pages = 11, year = "1999", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 336", month = jan, pages = 11, year = "1999", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 334", month = oct, pages = 23, year = "1998", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 333", month = aug, pages = 42, year = "1998", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 332", pages = "20", month = aug, year = "1998", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, quadratic eigenvalue problem, generalized eigenvalue 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, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 331", month = aug, year = "1998", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 325", pages = "15", month = may, year = "1998", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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 radius, numerical radius" } @techreport{Higham:1998:QRF, author = "Nicholas J. Higham", title = "{QR} Factorization with Complete Pivoting and Accurate Computation of the {SVD}", institution = inst-MCCM, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 324", month = sep, pages = 26, year = "1998", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 323", pages = "20", month = apr, year = "1998", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 322", pages = "15", month = mar, year = "1998", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 321", pages = "20", month = apr, year = "1998", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, pages = "20", month = mar, year = "1998", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 319", pages = "18", month = mar, year = "1998", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 317", pages = "28", month = mar, year = "1998", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 316", month = dec, year = 1997, note = "Submitted to the SIAM Journal on Scientifc Computing", URL = "http://www.maths.manchester.ac.uk/~higham/narep/narep316.ps.gz", abstract = "The practical implementation of the lowest-order \plpo\ (linear velocity, constant pressure) finite element method for the steady-state incompressible (Navier-)Stokes equations is 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, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 314", pages = "17", month = nov, year = "1997", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 313", pages = "16", month = nov, year = "1997", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 312", pages = "15", month = oct, year = "1997", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 311", pages = "15", month = oct, year = "1997", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 310", pages = "17", month = oct, year = "1997", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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" } @techreport{Ford:1997:ADM, 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, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 309", pages = "17", month = oct, year = "1997", URL = "http://www.maths.manchester.ac.uk/~higham/narep/narep309.ps.gz", abstract = "Adomian's decomposition method (ADM) has been presented 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, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 308", month = sep, pages = "8", year = "1997", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, 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.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 306", pages = "20", month = aug, year = "1997", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 305", month = apr, pages = "20", year = "1997", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, pages = "16", year = "1997", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 303", pages = "15", month = may, year = "1997", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 302", month = feb, pages = "17", year = "1997", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 301", pages = "17", month = feb, year = "1997", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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 sign in the Householder vector.", 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, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 300", month = dec, pages = "35", year = "1996", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 299", month = nov, pages = "28", year = "1996", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 298", month = nov, pages = "12", year = "1996", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 297", month = nov, pages = "24", year = "1996", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 296", month = sep, year = "1996", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 295", month = sep, pages = "17", year = "1996", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 294", month = aug, pages = "14", year = "1996", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 293", month = aug, year = "1996", pages = "34", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, month = sep, year = "1996", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 290", month = oct, year = "1996", note = "Submitted to the IMA Journal of Numerical Analysis", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 289", pages = "18", month = apr, year = "1996", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 288", pages = "26", month = apr, year = "1996", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 287", pages = "18", month = mar, year = "1996", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, month = mar, pages = "18", year = "1996", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 284", pages = "13", month = nov, year = "1995", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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 analysis, positive quadrature" } @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, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 283", pages = "20", month = dec, year = "1995", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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 integrand, using the adaptive quadrature algorithm developed 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, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 282", pages = "19", month = dec, year = "1995", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 281", month = dec, year = "1995", pages = "30", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 279", month = nov, year = "1995", pages = "20", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 277", month = sep, year = "1995", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 276", month = sep, year = "1995", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, year = "1995", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 274", month = jul, year = "1995", ftpaddress = "ftp.ma.man.ac.uk", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 273", month = jul, year = "1995", note = "To appear in proceedings of the 16th Dundee Bienniel Numerical Analysis Conference", ftpaddress = "ftp.ma.man.ac.uk", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 272", month = jul, year = "1995", ftpaddress = "ftp.ma.man.ac.uk", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, month = jun, year = "1995", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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 access to the most recent publications. Although it is hoped 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, address = inst-MCCM:adr, month = jun, year = "1995", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, month = may, year = "1995", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, month = may, year = "1995", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 265", month = jul, year = "1995", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, 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.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, 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.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, type = "Numerical Analysis Report", number = "No. 261", pages = "11", month = aug, year = "1994", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, month = oct, year = "1994", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, month = sep, year = "1994", note = "Submitted to SIAM J. Numer. Anal.", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, month = jul, year = "1994", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, month = jul, year = "1994", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, month = oct, year = "1994", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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", } @techreport{Gladwell:1994:IPS, 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, address = inst-MCCM:adr, month = jul, year = "1994", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, month = may, year = "1994", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, month = jul, year = "1994", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, month = apr, year = "1994", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, month = mar, year = "1994", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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 Adams-Bashforth-Moulton and backward differentiation 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, address = inst-MCCM:adr, month = mar, year = "1994", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, month = mar, year = "1994", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, month = feb, year = "1994", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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", } @techreport{Papadimitriou:1993:KNA, author = "Pythagoras Papadimitriou", title = "The {KSR1}---{A} Numerical Analyst's Perspective", type = "Numerical Analysis Report", number = "No. 242", institution = inst-MCCM, address = inst-MCCM:adr, month = dec, year = "1993", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, 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.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, month = oct, year = "1993", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, month = mar, year = "1994", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, month = dec, year = "1993", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, month = jul, year = "1993", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, month = aug, year = "1993", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, month = jun, year = "1993", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, month = may, year = "1993", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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", address = "England", 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.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, month = jun, year = "1993", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-MCCM:adr, month = apr, year = "1993", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-U-MANCHESTER:adr, month = oct, year = "1992", note = "To appear in {\em SIAM J. Sci. Comput.}.", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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.manchester.ac.uk/~higham/narep/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 quadrature scheme", type = "Numerical Analysis Report", number = "No. 213", institution = inst-U-MANCHESTER, address = inst-U-MANCHESTER:adr, month = may, year = "1992", URL = "http://www.maths.manchester.ac.uk/~higham/narep/narep213.ps.gz", abstract = "In this report we review the types of quadrature scheme that are available for the evaluation of integrals. We present a new adaptive quadrature scheme 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.", keywords = "adaptive quadrature, minimal storage", } @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, address = inst-U-MANCHESTER:adr, month = apr, year = "1992", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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.manchester.ac.uk/~higham/narep/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, address = inst-U-MANCHESTER:adr, month = feb, year = "1992", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = pub-SV:adr, pages = "29--46", year = "1993", series = "IMA Volumes in Mathematics and its Applications", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-U-MANCHESTER:adr, month = dec, year = "1991", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = inst-U-MANCHESTER:adr, month = dec, year = "1991", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = pub-LONGMAN:adr, pages = "137--154", year = "1990", series = "Pitman Research Notes in Mathematics", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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, address = pub-OUP:adr, pages = "161--185", year = "1990", URL = "http://www.maths.manchester.ac.uk/~higham/narep/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", }