text:
Main Numerical Analysis: A Mathematical Introduction
Numerical Analysis: A Mathematical IntroductionMichelle Schatzman
This book provides professionals and students with a thorough understanding of the interface between mathematics and scientific computation. Ranging from classical questions to modern techniques, it explains why numerical computations succeed or fail. The book is divided into four sections, with an emphasis on the use of mathematics as a tool in determining the success rate of numerical methods. The text requires only a modest level of mathematical training, and is ideally suited for scientists and students in mathematics, physics and engineering.
You may be interested in
You can write a book review and share your experiences. Other readers will always be interested in your opinion of the books you've read. Whether you've loved the book or not, if you give your honest and detailed thoughts then people will find new books that are right for them.
OXFORD Numerical Analysis Numerical Analysis A Mathematical Introduction MICHELLE SCHATZMAN Directeur de Recherche au CNRS Universite Claude Bernard Lyon Translated by JOHN T A Y L O R Edinburgh Petroleum Services CLARENDON PRESS •OXFORD 2002 This book has been prin ted digitally and produced in a standard specification in order to ensure its continuing availability OX FOR D U N I V E R S I T Y PRESS Great Clarendon Street, Oxford 0X2 6DP Oxford University Press is a department o f the University o f Oxford. It furthers the University' s objective o f excellence in research, scholarship, and education by publishing worldwide in Oxford New York Auckland Cape Town Dar es Salaam Hong Kong Karachi Kuala Lumpur Madrid Melbourne Mexico City Nairobi New Delhi Shanghai Taipei Toronto With offices in Argentina Austria Brazil Chile Czech Republic France Greece Guatemala Hungary Italy Japan South Korea Poland Portugal Singapore Switzerland Thailand Turkey Ukraine Vietnam Oxford is a registered trade mark o f Oxford University Press in the UK and in certain other countries Published in the United States by Oxford University Press Inc., New York © Dunod & Masson Sciences, InterEdition, Paris 1991 English translation © Oxford University Press 2002 Translation o f Analyse Numerique by Michelle Schatzman originally published in French by Inter Editions, Paris 1991 Aide par le ministere fran<;ais charge de la culture The moral rights o f the author have been asserted Database right Oxford University Press (maker) Reprinted 2011 ISBN 978-0-19-850852-6 Preface Ne sois pas lascif et peureux Comme le lievre et l' amoureux. Mais que toujours ton cerveau soit la hase pleine qui congoit.1 Guillaume Apollinaire, Bestiaire (1920) Fairy tales Fair reader, the preface of a mathematical book is the right place for telling how it came about, and maybe spawn a few legends. Once upon a time, in a far away country2, there was a young research mathe matician who specialized in nonlinear partial differential equations. In 1984, this young mathematician was promoted to professor of numerical analysis, a subject that she mostly ignored. So, in order to teach, she learnt, a fact already known from the Talmud: ' Rav Hanina said, " Much have I learnt from my masters, more from my colleagues, but the most from my own students' "(Talmud of Babylon, Tractate Taanit, 6)3. You know how a fairy tale is composed: the heroine has to fall in love: I fell in love with numerical analysis4. And here is the result, or rather the second iteration of the result, since a first edition appeared in French in 1991. Sonke Adlung, who is an editor with Oxford University Press, and to whom I had been introduced by John Ball, thought that it would be nice to have this book translated and also revised. The heroine must go through hard times and meet a few fire-spitting dragons on her way, which would change her into stone just for the fun of it. I was moving xBe not lewd and fearful/ as the hare and the love-fool./ But let your brain ever be/ the hare-doe that conceives. 2We have to stick to the traditional format of fairy tales. 3*im> >*r>tonoi ,>rrcn»*im> n>nnai .' ircno nrnn :w>:n m -ion (.1 m u n ,'^22 TiTrtm .pIDO 4Some would consider that bad taste; theirs is the loss. v PREFACE VI so slowly at times that Sonke may well have thought that I had been turned into a piece of rock. Since I had agreed to update the 1991 edition, I had to rewrite parts of it, and add some material. This did not turn out to be easy. One reason is that my extra-scientific activities graduated from raising two children5 to chairing a research group of about 50 people, a task which includes dealing with five different administrations. Another reason is that I am still a research mathe matician, writing papers, advising students, and doing whatever the form of life called mathematics demands. My natural tendency is to consider that the most important mathematics are tomorrow' s, not yesterday' s or yesteryear' s. Though I should not admit it, I may have also fallen in love with my research group, MAPLY (Laboratoire de Mathematiques Appliquees de Lyon), and with its future; when I was born, there were trees; so now, I feel obligated to plant some. I did meet quite a few dragons in the forest6. This being a legend, fair reader, you have to remember that dragons are in the mind of the beholder. Acknowledgements The heroine needs to receive help from many quarters, and so did I. Sonke found a translator, John Taylor. John, you did a pretty good job, keeping as much as you could of the colloquialism of the original French style. Thank you, Sonke, for the idea which enabled me to conclude the task: that I should find myself someone who would not be so busy, who would be able to spot typos and mixed up indices7, and who would be knowledgeable enough to criticize or praise wherever applicable. And this is how Jean-Frangois Coulombel came in: in February 2000, he was sitting in one of my graduate courses, and one month later, he was bravely starting to push me forward, so that the slow motion would not be so slow. Thank you, Jean-Frangois, you did well, and that must not have been easy for you. And thank you again, Sonke, for never relenting before my procrastination. In the final runs, I shared an office with Stephane Descombes in a faraway 5Claude and Rene, you have built me. 6 Nel mezzo del cammin di nostra vita mi ritrovai per une selva oscura che la diritta via era smarrita. (In the middle of our life' s path/ I found myself in dark woods,/ the clear direct way being lost.) Dante, La Divina Commedia 71 guess that I became a mathematician because I have so much trouble with + and —, and with z, j, k, /, and m. And n too; so if there are still some errors, I am the one who put them in, and I apologize. PREFACE vii country8, and he was kind enough to provide a stern examination of the spline chapter. I also benefited from a rereading of the multistep chapter by Magali Ribot, who had a first-hand opportunity to learn about the fallibility of some thesis advisors. I have had many sources and I have been influenced in person or in print by the following authors: K. E. Atkinson [5], J. C. Butcher [13], P. G. Ciarlet [16], M. Crouzeix and A. L. Mignot [19], C. de Boor [20], G. H. Golub and C. F. van Loan [35], E. Hairer, S. P. Norsett, and G. Wanner [43,42], P. Henrici [45], E. Isaacson and H. B. Keller [51], A. Iserles [52], D. E. Knuth [38,53,54], H.-O. Kreiss [56], Y. Meyer [61], A. Ralston and P. Rabinowitz [68], R. D. Richtmyer and K. W. Morton [70], L. L. Schumaker [71], J. Stoer and R. Bulirsch [73], and H. S. Wilf [78]. To these, I am deeply indebted. I have now two institutions to thank: CNRS gave me the initial help in my research career; since 1995, it gave me the charge of planting mathematical trees; it told me, in very direct language, that one has to go forward and take responsibility, or get out of the way; it gave me a temporary research position which enabled me to remain a mathematician instead of turning into a full-time administrator. Finally, it gave me a permanent research position. The other institution is the Technion in Haifa, where I have been a frequent visitor since 1994. Some of the funding came from the binational Keshet/Arc-enciel binational agreement; I also got funding from the CNRS-MOSA. So, twice a year, I have some continuous time to myself, with only a few e-mails and faxes for French business. There is this wonderful library, from which I take out books on the account of my friend and colleague Koby Rubinstein, who from time to time gets a message telling him to return some. And the friends that I have there kindly let me be as bearish as I want; their hospitality is a blessing; I owe much to the kindness of Iris and Yehuda Pinchover. Finally, in a fairy tale, the heroine has the use of some life-saving tricks, spells or magic formulae. They are called T^X, SCILAB, GNUPLOT, and XFIG; they are all free software, and they are wonderful tools without which this book would simply not exist. Contents: a subjective approach Since this preface is so romantic, let me mention that love for a scientific subject is exactly the same as love for a human being: you always hate quite a few things in the loved one, but the balance looks good enough to keep the attachment going. So, let me talk subjectively about my subject. Looks are so important, and, alas, numerical analysis with all its burden of notations looks more like a heavy matron than a gracious ballerina. So what? Maybe a heavy matron can take care of a large brood of children. Better a heavy 8Israel, according to the age-old joke: this is far? far from where? viii PREFACE theory with many applications than an elegant one without offspring. Numerical analysis has to do with real life computation: an understanding of the floating number system and the machine arithmetic is essential. Chapter 1 gives an exercises guided tour to this world. It is highly recommended that you do the exercises. I observed that numerical analysis requires more maturity from mathemat ics students than other subjects at the same level of difficulty. A naive vision puts numerical analysis between two stools: the physicist or the engineer wants methods, and is satisfied with experimental numerics; the mathematician wants beautiful problems, and is not much interested in constructive solutions using a fallible and limited machine. However, I believe that numerical analysis sits on both stools, and has the best of both worlds: motivated problems where getting a solution and getting it fast can make a difference, and tools which can be, at the same time, elementary and powerful: not necessarily a contradiction according to the words of S. S. Abhyankar [1]. Chapter 2 gives a flavour of numerical analy sis, constructing the logarithm and the exponential from scratch, using methods which are the daily bread of numerical analysis, and which moreover generalize to situations where power series do not work. The daily life of numerical analysts includes much linear algebra; Chapter 3 summarizes some of the required knowledge, and adds to it the theory of block matrices. This makes up Part I. Part II describes polynomial approximation and piecewise polynomial approximation, in algebraic or trigonometric versions: inter polation and divided differences in Chapter 4, least-squares approximation in Chapter 5, and splines in Chapter 6. The recent surge of the use of splines in computer-aided geometrical design and image processing is one of the motiva tions for the spline chapter. It also turns out that splines are a nice generalization of Bernstein polynomials, and that they fit very well with two approaches: di vided differences and convex algorithms. I used to hate splines, a baseless and despicable prejudice. I hate them no longer. Chapter 7 is on Fourier series in one space dimension; it does not tackle any of the hard questions of which Fourier analysis is so replete. However, it treats easy and essential questions, including convolution and regularization, and it makes room for the Gibbs phenomenon, so important in applications. Chapter 8 is about quadrature: approximation by algebraic polynomials leads to the classical formulae for numerical integration; trigonometric approximation leads to the Euler-MacLaurin formula and to the fast Fourier transform (FFT), which may be the most important algorithm in scientific computation. The FFT is the prototype of recursive algorithms; it is the ancestor of multigrid and wavelet algorithms; it is the epitome of easy and powerful tools. Part III relates to numerical linear algebra. This part is important because operation counts are the limiting factor for any serious computation. Any sci entific computation program spends most of its time solving linear systems or approximating the solution of linear systems, even when trying to solve nonlinear PREFACE IX systems. Chapter 9 provides the direct methods for the resolution of linear systems of equations, with an emphasis on operation counts. Operation counts justify all of the acrobatics of iterative methods, treated in Chapter 11 after extra information on linear algebra is produced in Chapter 10. Chapter 12 relates orthogonality methods to the resolution of linear systems and introduces the QR decomposition. Part IV treats a selection of nonlinear or complex problems: resolution of linear equations and systems in Chapter 14, ordinary differential equations in Chapter 15, single-step schemes in Chapter 16 and multistep schemes in Chap ter 17, and introduction to partial differential equations in Section 18.4. It would have been my natural tendency to put in much more of these, but, since I decided that this book would avoid any functional spaces beyond spaces of continuous functions or Lebesgue spaces and, in particular, Sobolev spaces, there was little possibility to include more than a tiny introduction to partial differential equa tions. I tried to select a few important things which are accessible and attractive on an elementary level. I did the same for the ordinary differential equations part: I skipped the detailed analysis of the Runge-Kutta methods by trees, be cause it is long and difficult; however, I have given a full theory of the analysis and convergence of multistep schemes, because the use of appropriate norms makes it possible without tears. This book started as an elementary book; the revision put in some more advanced layers, but the layered structure remains; the less elementary parts are Sections 4.4 and 4.5, Chapter 6, Section 7.2, Section 8.6, and Chapters 17 and 18. A number of problems describe some classical algorithms together with some newer ones. Since mathematics is not a spectator sport, the more advanced parts are put into problems: the most exciting things can be found there. The prerequisites are linear algebra, calculus, and a tiny bit of Lebesgue the ory, which is used only in Chapter 5 on polynomial least-squares approximation, Chapter 7 on Fourier analysis, and Chapter 18 which introduces partial differ ential equations. I do not use the theory of distributions, though I disguise some of its ideas in the spline Chapter 6. Just a short word about notation: I decided not to use bold face for matrices or vectors, with very few exceptions. The reason is that I very often use block decomposition of matrices; if I decompose an n x n matrix into an (n —1) x (n —1), an ( n - l ) x l , a l x ( n - l ) , and a 1 x 1 block, what notation would make sense? I could not imagine an efficient answer, and so I dropped the bold faces. Computation and numerical analysis There is always a question about the role of computations and software in a numerical analysis book. I did not include any numerical software, or even quasi-programs which enable X PREFACE one to write one' s own software. The first reason is that good numerical software requires quite a bit of thought to be really efficient. Some really good scientific software, such as MATLAB can be bought; but there are also free products of quality, such as SC I LAB, distributed by INRIA, http//www-rocq.in r ia .fr / s c ila b , for use under several different operating systems. The second is that one must be sceptical of scientific computations: numerical analysis is concerned with the essence of the scientific method. But what is the foundation of such an attitude? It is conceivable if we know that the results are inaccurate, or plainly wrong. But how do we know that? Well, if we know that our mathematical equations are a good approximation of a natural phenomenon, and if the computations do not agree with the observations, then the scientific computation software must be guilty; or must it? Finding the guilty party may be a very difficult endeavour, because many factors may be involved: maybe the equations were wrong; maybe the numerical method was inappropriate; maybe the parameters of the method were badly chosen; maybe the software was incorrect. One of the purposes of numerical analysis is to find the specifically mathe matical factors which govern the success or failure of numerical computations. The third reason is that the successful numerical analyst destroys their own job by finding algorithms which are so clear and efficient that they can be safely implemented into software. Therefore, in order to stay in employment, he or she must keep finding new areas where existing software does not do the job. Mathematics provides the light which enables us to explore new territory. Mathematics is also a very cultural subject: for mathematicians, a hundred years old result can be as good as a new one, even if we do not use it for the purpose for which it had been initially crafted. In this book, I try to explain how to make your own light, and how to find your own way with it. This is obviously much harder—for you and for me—than inviting you to contemplate nice pictures without telling you how I might have come across them and how they fit together. Nevertheless, fair reader, do not believe that I have shyed away from com puting: I have only hidden it, since showing it would have led to a completely different book. So, in order to gain more understanding, try, dear reader, to program the most algorithms that you can think of. It may not be of the same high quality as the commercially or freely available software, but it will teach you much about the behaviour of computational methods, and about the difficulty of putting mathematics into code. Lyon September 2000 M. S. Contents I T h e entran ce fee 1 1 F loatin g n um bers 1.1 Counting in base /3 1.2 Expansion of the rational numbers in base /? 1.3 The machine representation of numbers 1.4 Summation of series in floating-point numbers 1.5 Even the obvious problems are rotten 1.6 Even the easy problems are hard 1.7 A floating conclusion 3 3 5 5 7 8 9 10 2 A flavour o f num erical analysis 2.1 Comparison of exponentials and powers 2.2 Convergence and divergence of classic series 2.3 Discrete approximation of the logarithm 2.4 Comparison of means 2.5 Elementary construction of theexponential 2.6 Exponentials of matrices 11 12 13 13 16 17 20 3 A lgebraic prelim in aries 3.1 Linear algebra refresher 3.1.1 The matrix of a linear mapping 3.1.2 The determinant 3.1.3 The fundamental theorem of linear algebra and its consequences 3.1.4 Eigenvalues and eigenvectors 3.1.5 Scalar products, adjoints, and company 3.1.6 Triangular matrices 3.2 Block matrices 3.2.1 Block decomposition of a linear mapping or matrix 3.2.2 Block multiplication 3.3 Exercises from Chapter 3 3.3.1 Elementary algebra 3.3.2 Block decomposition 25 25 25 27 xi 29 30 31 34 34 34 36 37 37 38 CONTENTS xii 3.3.3 3.3.4 3.3.5 3.3.6 3.3.7 II Graphs and matrices Functions of matrices Square roots, cosines, and sines of matrices Companion matrices and bounds of matrix powers The Kantorovich inequality P oly n om ia l and trig o n o m etric ap p rox im a tion 39 40 41 42 44 47 4 In terp ola tion and d iv id ed differences 4.1 Lagrange interpolation 4.1.1 The Lagrange interpolation problem 4.2 Newton' s interpolation 4.2.1 Newton' s basis is better than Lagrange' s basis 4.2.2 Integral representation of divided differences 4.3 Interpolation error 4.4 Hermite and osculating interpolation 4.5 Divided differences as operators 4.5.1 Finite differences on uniform grids 4.6 Exercises from Chapter 4 4.6.1 More on divided differences 4.6.2 Boundary value problem for an ODE 4.6.3 Extrapolation to the limit 51 51 51 54 54 57 59 64 67 70 72 72 74 75 5 Least-squares ap proxim ation for poly n om ia ls 5.1 Posing the problem 5.1.1 Least-squares is Pythagoras in many dimensions 5.1.2 Is it really calculable? 5.2 Orthogonal polynomials 5.2.1 Construction of orthogonal polynomials 5.2.2 Examples of orthogonal polynomials 5.2.3 Revival of special functions 5.2.4 Orthogonal polynomials and least-squares 5.3 Polynomial density: Bernstein polynomials 5.3.1 Modulus of continuity 5.3.2 Bernstein polynomials and Bernstein approximation 5.3.3 Application of Bernstein polynomials to graphics software: the Bezier curves 5.4 Least-squares convergence 5.5 Qualitative properties 5.6 Exercises from Chapter 5 5.6.1 Laguerre polynomials 5.6.2 Pade type and Pade approximations 77 77 77 80 81 82 83 87 88 88 89 90 94 96 97 99 99 101 CONTENTS xiii 6 Splines 6.1 Natural splines: the functional approach 6.1.1 Weak equality of functions 6.1.2 Weak integrals of functions 6.1.3 The space of natural splines 6.2 Numerics for cubic natural splines 6.3 Spaces of splines, 5-splines 6.3.1 Splines with distinct knots 6.3.2 The beautiful properties of 5-splines 6.3.3 Numerics with 5-splines 6.3.4 Using 5-splines to understand natural splines 6.3.5 5-splines in CAGD 6.4 Exercises from Chapter 6 6.4.1 Varied exercises on splines 6.4.2 Approximation by splines 6.4.3 Coincident knots 106 107 107 108 110 116 119 119 120 125 127 128 130 130 131 132 7 F ourier' s w orld 7.1 Trigonometric approximation andFourier series 7.1.1 Trigonometric polynomials 7.1.2 Integration of periodic functions 7.1.3 Least-squares approximation for trigonometric polynomials 7.1.4 Density of trigonometric polynomials in the space of continuous periodic functions 7.1.5 Convergence in the mean square of trigonometric approximation to continuous functions 7.1.6 Asymptotic behaviour of Fourier coefficients 7.1.7 Convergence of trigonometric approximation to Lj functions 7.1.8 Uniform convergence of Fourier series 7.2 Convolution and pointwiseconvergence 7.2.1 Convolution 7.2.2 Regularization 7.2.3 Constructive density results 7.2.4 Convolution and Fourier series 7.2.5 Convergence of Fourier series as a localphenomenon 7.2.6 Pointwise convergence of partial Fouriersums of absolutely continuous functions 7.2.7 Pointwise convergence of partial Fourier sums of piecewise absolutely continuous functions 7.2.8 Gibbs phenomenon 7.3 Exercises from Chapter 7 7.3.1 Elementary exercises on Fourier series 133 133 134 134 135 136 139 140 142 144 146 146 148 150 151 152 153 156 157 159 159 CONTENTS X IV 7.3.2 7.3.3 7.3.4 7.3.5 8 III 9 Fejer, La Vallee Poussin, and Poisson kernels There exists an integrable function whose Fourier coefficients decrease arbitrarily slowly to 0 The existence of sequences of numbers dk tending to 0 as |fc| tends to infinity which are not the Fourier coefficients of any integrable function Discrete least-squares approximation by trigonometric polynomials Q uadrature 8.1 Numerical integration 8.1.1 Numerical integration for dummies 8.2 The analysis of quadrature formulae 8.2.1 Order of a quadrature formula 8.2.2 On the practical interest of weighted formulae 8.2.3 Examples of simple formulae 8.2.4 Composite formulae 8.3 The Peano kernel and error estimates 8.3.1 Definition of the Peano kernel 8.3.2 Quadrature error in composite formulae 8.4 Gaussian quadrature 8.5 Periodic numerical integration 8.6 Bernoulli, Euler-MacLaurin 8.6.1 Detailed analysis of the trapezium formula 8.6.2 The Bernoulli polynomials 8.6.3 The Euler-MacLaurin formula 8.7 Discrete Fourier and fast Fourier transforms 8.7.1 Discrete Fourier transforms 8.7.2 Principle of the fast Fourier transform algorithm 8.7.3 FFT algorithm: decimation-in-frequency 8.8 Exercises from Chapter 8 8.8.1 Summation of series with Bernoulli numbers and polynomials 8.8.2 The Fredholm integralequation of the first kind 8.8.3 Towards Franklin' s periodic wavelets N u m erical linear algeb ra G a u ss' s w orld 9.1 Elimination without pivoting 9.1.1 Just elimination 9.1.2 Matrix interpretation of Gaussian elimination 9.2 Putting it into practice: operation counts 160 162 162 163 165 166 166 169 169 171 172 173 174 174 179 180 183 184 185 187 189 189 190 191 193 197 197 200 202 205 207 207 207 209 212 CONTENTS 9.2.1 9.2.2 9.2.3 9.2.4 The madness of Cram er' s rule Putting eliminationinto practice Operation counts for elimination Inverting a matrix: putting it into practice and the operation count 9.2.5 Do we need to invert matrices? 9.3 Elimination with pivoting 9.3.1 The effect of a small pivot 9.3.2 Partial pivoting and total pivoting: general description and cost 9.3.3 Aside: permutation matrices 9.3.4 Matrix interpretation of partial and total pivoting 9.3.5 The return of the determinant 9.3.6 Banded matrices 9.4 Other decompositions: LDU and Cholesky 9.4.1 The LDU decomposition 9.4.2 The Cholesky method 9.4.3 Putting the Cholesky method into practice and operation counts 9.5 Exercises from Chapter 9 9.5.1 Exercises on the rank of systems of vectors 9.5.2 Echelon matrices and least-squares 9.5.3 The conditioning of a linear system 9.5.4 Inverting persymmetric matrices xv 212 214 214 215 216 217 217 219 221 222 224 224 227 227 228 231 231 231 232 235 2371 0 10 T h eoretical interlude 240 10.1 The Rayleigh quotient 240 10.2 Spectral radius and norms 242 10.2.1 Spectral radius 242 10.2.2 Norms of vectors, operators, and matrices 243 10.3 Topology and norms 244 10.3.1 Topology refresher 244 10.3.2 Equivalence of norms 244 10.3.3 Linear mappings: continuity, norm 247 10.3.4 Subordinate norms 248 10.3.5 Examples of subordinate norms 250 10.3.6 The Frobenius norm is not subordinate 252 10.4 Exercises from Chapter 10 254 10.4.1 Continuity of the eigenvalues of a matrix with respect to itself 254 10.4.2 Various questions on norms 255 CONTENTS XV I 11 Iterations and recurren ce 11.1 Iterative solution of linear systems 11.2 Linear recurrence and powers of matrices 11.2.1 The spectrum of a finite difference matrix 11.3 Exercises from Chapter 11 11.3.1 Finite difference matrix of the Laplacian in a rectangle 11.3.2 Richardson' s and pre-conditioned Richardson' s methods 11.3.3 Convergence rate of the gradient method 11.3.4 The conjugate gradient 11.3.5 Introduction to multigrid methods 257 258 270 273 276 276 277 280 282 285 12 P y th a gora s'w orld 12.1 About orthogonalization 12.1.1 The Gram-Schmidt orthonormalization revisited 12.1.2 Paths of inertia 12.1.3 Topology for QR and Cholesky 12.1.4 Operation counts and numeric strategies 12.1.5 Hessenberg form 12.1.6 Householder transformations 12.1.7 QR decomposition by Householder transformations 12.1.8 Hessenberg form by Householder transformations 12.2 Exercises from Chapter 12 12.2.1 The square root of a Hermitian positive definite matrix 290 290 290 293 295 296 297 298 299 302 303 303 IV N onlinear p ro b lem s 13 S p ectra 13.1 Eigenvalues: the naive approach 13.1.1 Seekingeigenvalues and polynomialequations 13.1.2 The bisection method 13.2 Resonance and vibration 13.2.1 Galloping Gertie 13.2.2 Small vibrations 13.3 Power method 13.3.1 The straightforward case 13.3.2 Modification of the power method 13.3.3 Inverse power method 13.4 QR method 13.4.1 The algorithm and itsbasic properties 13.4.2 Convergence in a special case 13.4.3 Effectiveness of QR 13.5 Exercises from Chapter 13 13.5.1 Spectral pathology 305 307 311 311 312 313 313 314 316 316 319 321 322 322 324 327 328 328 CONTENTS 13.5.2 QR flow and Lax pairs xvii 328 14 N onlinear equ ation s and system s 14.1 From existence to construction 14.1.1 Existence and non-existence of solutions 14.1.2 Existence proofs translate into algorithms 14.1.3 A long and exciting history 14.1.4 An overview of existence proofs 14.2 Construction of several methods 14.2.1 The strictly contracting fixed point theorem 14.2.2 Newton' s method: geometric interpretation and examples 14.2.3 Convergence of Newton' s method 14.2.4 The secant method 14.2.5 The golden ratio and Fibonacci' s rabbits 14.2.6 Order of an iterative method 14.2.7 Ideas on the solution of vector problems 14.3 Exercises from Chapter 14 14.3.1 The Cardano formulae 14.3.2 Brouwer' s fixed point theorem in dimension 2 14.3.3 Comparison of two methods for calculating square roots 14.3.4 Newton' s method for finding the square roots of matrices 331 331 331 333 333 334 336 336 338 341 344 348 351 352 353 353 354 355 3571 5 15 S olvin g differential sy stem s 15.1 Cauchy-Lipschitz theory 15.1.1 Idea of the proof of existence for ODEs 15.1.2 Cauchy-Lipschitz existence theorem 15.1.3 Systems of order 1 and of order p 15.1.4 Autonomous and non-autonomous systems 15.2 Linear differential equations 15.2.1 Constant coefficient linear systems 15.2.2 Matrix exponentials 15.2.3 Duhamel' s formula 15.2.4 Linear equations and systems with variable coefficients 15.2.5 Gronwall' s lemma 15.2.6 Applications of Gronwall' s lemma 15.2.7 Smoother solutions 15.3 Exercises from Chapter 15 15.3.1 Lyapunov function for a 2 x 2 linear system 15.3.2 A delay differential equation 15.3.3 A second-order ordinary differential equation 362 362 362 363 366 367 367 367 369 371 372 377 379 380 382 382 382 384 xviii CONTENTS 16 Single-step sch em es 385 16.1 Single-step schemes: the basics 386 16.1.1 Convergence, stability, consistency 387 16.1.2 Necessary and sufficient condition of consistency 389 16.1.3 Sufficient condition for stability 390 16.2 Order of a one-step scheme 392 16.3 Explicit and implicit Euler schemes 395 16.3.1 The forward Euler scheme 395 16.3.2 Backwards Euler scheme 396 16.3.3 0-method 398 16.4 Runge-Kutta formulae 399 16.4.1 Examples of Runge-Kutta schemes 400 16.5 Exercises from Chapter 16 403 16.5.1 Detailed study of the 0-scheme 403 16.5.2 Euler scheme with variable step size and asymptotic error estimates 404 16.5.3 Numerical schemes for a delay differentialequation 407 16.5.4 Alternate directions 409 16.5.5 Numerical analysis of a second-orderdifferentialequation 412 17 Linear m u ltistep sch em es 17.1 Constructing multistep methods 17.1.1 Adams methods 17.1.2 The multistep methods of Adams 17.1.3 Backward differentiation 17.1.4 Other multistep methods 17.2 Order of multistep methods 17.2.1 The order is nice and easy for multistep methods 17.2.2 Order of some multistep methods 17.3 Stability of multistep methods 17.3.1 Multistep methods canbe very unstable 17.3.2 The stability theory for multistep methods 17.3.3 Stability of some multistep schemes 17.4 Convergence of multistep schemes 17.4.1 Initializing multistep methods 17.4.2 Solving in the implicit case 17.5 Exercises from Chapter 17 17.5.1 Short exercises 17.5.2 An alternative formulation of the ordercondition 17.5.3 Weak instability 17.5.4 Predictor-corrector methods 17.5.5 One-leg methods 414 414 415 416 419 420 421 421 423 424 424 426 432 432 432 433 433 433 434 434 435 437 CONTENTS xix 18 Tow ards partial differential equ ation s 439 18.1 The advection equation 439 18.1.1 The advection equation and its physical origin 439 18.1.2 Solving the advection equation 442 18.1.3 More general advection equations and systems 444 18.2 Numerics for the advection equation 446 18.2.1 Definition of some good and some bad schemes 446 18.2.2 Convergence of the scheme (18.2.6) 450 18.3 The wave equation in one dimension 455 18.3.1 Masses and springs 455 18.3.2 Elementary facts about the wave equation 458 18.3.3 A numerical scheme for the wave equation 459 18.4 The heat equation and separation of variables 461 18.4.1 Derivation of the heat equation 461 18.4.2 Seeking a particular solution by separation of variables 465 18.4.3 Solution by Fourier series 467 18.4.4 Relation between the heat equation and the discrete model 469 18.5 Exercises from Chapter 18 471 18.5.1 The eigenvectors of a strictly hyperbolic matrix 471 18.5.2 More on the upwind scheme 472 18.5.3 Fourier analysis of difference schemes for the advection equation 474 18.5.4 The Lax-Friedrichs scheme 474 18.5.5 The Lax-Wendroff scheme 475 18.5.6 Stability of the leap-frog scheme 475 18.5.7 Elementary questions on the wave equation 476 18.5.8 Generalized solutions for the advection equation 476 18.5.9 Advection-diffusion equation 478 R eferen ces 479 Index 485 Part I The entrance fee Most of numerical analysis uses recursive procedures: do the same thing again and again, until a reasonable degree of accuracy is reached. Even apparently algebraic problems admitting a solution in finite terms can be efficiently treated by iterative techniques, as we shall see in Chapter 11. Therefore, since numerical analysis involves a large number of machine computations, the error must be analysed. Numerical analysts are scared stiff of a phenomenon called instability: un controlled amplification of error. But where does error come from? It comes from the fact that we represent real numbers in a finite system, called floating numbers; therefore, any arithmetic operation leads to loss of precision. The idiosyncrasies of floating-point operations are examined in Chapter 1. Of course, instability is quite visible: usually one gets a computer message which says ' overflow' , and at this point the computation stops. What is impor tant is to find a cure for instability, and this requires a comprehensive understand ing of the mathematical methods which have been used for the computation. As presented here, numerical analysis is a part of mathematics, but it works on questions which are strongly related to the use of computers and to applica tions from other sciences. Therefore, numerical analysts create in their minds visions of mathematical objects which may be slightly different from the visions of other mathematicians. In particular, numerical analysis is about construct ing or approximating effective solutions. Of course, numerical analysts do use existence theorems, and they combine constructive and non-constructive infor mation. 1 2 I. THE ENTRANCE FEE I have tried to give a flavour of numerical analysis in Chapter 2 through concrete examples, which are probably well known to the reader. The fastest way of constructing the natural logarithm is to say that it is the integral of the function x i-» 1/x which vanishes at 1. This information can be made constructive: in order to compute the integral of x 1/x, we use ap proximate formulae, and these formulae enable us to construct from scratch the natural logarithm. In the same fashion, we can define the exponential as the reciprocal function of the logarithm, or using its famous entire series expansion. Unfortunately, neither of these methods leads to an efficient numerical construc tion. A good construction method for the exponential is the original method of Euler, based on a product formula; the exponential can be constructed from scratch with this method. It is of even more interest that such product for mulae are the core of numerical integration for ordinary and partial differential equations. Chapter 3 provides a review of standard results from linear algebra. Linear algebra is pervasive in this book. We cannot do anything without linear algebra, and we need a bit more than is usually taught at the elementary level. The only not completely standard feature is block matrices. 1 Floating numbers The beginning of a course on numerical analysis naturally includes some fairly ab stract considerations of how real numbers are represented in a computer, known as floating-point representation, and incidentally on systems of counting. To convince the reader that there is material here which is both surprising and thought-provoking, the following exercises are more effective than a long discourse. It is strongly recommended that the reader sharpens his or her mind by trying them. The level of mathematics and programming required is entirely elementary, though this does not spoil the fun of doing them. This chapter owes an enormous debt to the paper by G.E. Forsythe [30], which presents some striking examples of the calculation difficulties linked to floating-point numbers. 1.1. Counting in base (3 Exercise 1.1.1. Let /? be an integer greater than 1. Show that, for every integer n greater than or equal to 1, there exists a unique integer p and integers d*, 0 $ i ^ p, between 0 and /? — 1 inclusive, with dp ^ 0, such that (l.i.i) V The right-hand side of eqn (1.1.1) gives the representation of n in base /?, also denoted by n = dpdp—i • • • dido . Normally, we represent numbers in base f3 = 10 using the figures 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, and from an early age we use the result of Exercise 1.1.1 without question, at least in base (3 = 10. The choice of base 10 is linked to an anatomical peculiarity of the human species. We could also have counted in base 20, like the Mayans. 3 1. FLOATING NUMBERS 4 Nevertheless, in the history of humanity, counting with positional numerals is a relatively recent development which we owe to the Hindus and the Arabs. Arithmetic operations in roman numerals are very awkward. It is only thanks to positional numerals that efficient arithmetic algorithms could be developed, and these spread across Europe only from the twelfth century, finally triumphing in the eighteenth century. The reader who is interested in the history of numbers and systems of counting should consult, for example, [48] and [49], and their translation [50]. Other systems of counting have been devised by mathematicians, see, in particular, in [38, pp. 115-23]. Computers, which do not have ten fingers, count in base 2, with the figures 0 and 1 (binary or dyadic), in base 8 (octal), with the figures '0 to 7, and in base 16 (hexadecimal), with the figures "0 to "9, to which are added the letters "A to "F. The list of 128 ASCII characters is a list of standard characters corresponding to codes understood by all computers. Exercise 1.1.2. The ASCII characters are numbered from 0 to 127; the letter b has octal number 742. Give its decimal and hexadecimal numbers. What do we do with the fractional part of a number? By analogy with the representation of Exercise 1.1.1, consider expressions of the form x = j^ i——q where p and q are positive integers or zero and the di are integers from 0 to (3—1 inclusive. To fix p, we insist that dp ^ 0. We write x — dpdp— * d\d().d—j ... d—q , We can also take q = oo. What does this infinite sum mean? If we know the properties of real numbers, we claim that the sequence of rationals v rn = Pl(li i——n is increasing and bounded by p 1+ frdi. i=0 It is therefore a convergent sequence whose limit is the sum of the series with general term (/?*</*).<p. Warning: this course assumes that the reader is familiar with the properties of the real numbers. It is wise to revise them before continuing. A standard difficulty is the occurrence of a real number which can have two distinct representations in base /?, as we see in the following exercise: 1.3. THE MACHINE REPRESENTATION OF NUMBERS 5 Exercise 1.1.3. Let b = /3 — 1. Show that, in every base /?, 1 = O.bbbb---. How many reals are there which have two distinct representations in base /i? Few, as we will see in the exercise below: Exercise 1.1.4. Give the general form of the reals possessing two distinct repre sentations in base /?. Show that the set of all the reals for which there exists a base (3 in which the real number has two distinct representations is exactly the set of the rationals. 1.2. Expansion of the rational numbers in base In the decimal base, the quotient of two integers ' falls exactly'or ' does not fall exactly' , but after a certain point it is always periodic. The following exercises allow us to verify this periodicity result in any base /?. This section forms a short problem on elementary arithmetic. Before working on the general case, we will solve the following particular case: Exercise 1.2.1. Calculate the decimal expansion of 1/7. Exercise 1.2.2. Let m and n be two relatively prime integers such that m < n. Let ro = ra, and define d-j and r_j iteratively as being the quotient and the remainder, respectively, of the Euclidean division of /?r_j+i by n: /?r_j+1 = nd-j + r - j, 0^ < n. Show that, for every j ^ 1, 0 ^ d-j < (3. Exercise 1.2.3. Show that 0.d-\d-2d-$ • • •is the expansion of m/n in base (3. Exercise 1.2.4•Show that there exist two integers k and £ such that are equal. (Argue by contradiction.) and rt Exercise 1.2.5. Deduce from this that the expansion of m/n in base [3 is periodic from a certain point. Generalize this to the case m/n ^ 1. Exercise 1.2.6. When do the divisions m/n ' fall exactly'in base /?? Exercise 1.2.7. Give the expansion of 1/5 in bases 2, 8, and 16. (Observe care fully the result in base 16.) 1.3. The machine representation of numbers Several types of numbers can be represented in a computer. The machine integers are nothing special, the only point to note is that the set of integers which can be represented is finite, for example from -32768 to +32767. Recall that a bit of information is a binary digit, i.e., 0 or 1, and that a byte is a group of 8 bits. Since in this example there are exactly 216 different numbers, it is sufficient to have 16 bits or 2 bytes to represent them. 1. FLOATING NUMBERS 6 The floating-point numbers are more interesting. The set of these numbers is described by a base /?, a number of significant figures r, and two integers m and M. Every floating-point number is of the form s O.d-id-2 • • • d-r /3J with j between m and M, 1 ^ cLi < /?, and, for k > 1, 0 ^ d-k < /?. The letter s designates the sign of the number. We add the zero to the set of floating-point numbers. We thus obtain a subset F(/?,r, ra,M) of E which is formed from numbers whose expansion in base /? is finite. The expansion 0.d_id_2 • • • d-r is the mantissa of the number and j is its exponent. The normalization 1 ^ d-\ < /? is very important. It ensures that all the figures of all the nonzero floating-point numbers are significant. If the result of an operation is greater in absolute value than the largest floating-point number then the machine generally returns ' overflow' . If the result of a nonzero operation is rounded to zero then the machine returns ' underflow' . Exercise 1.3.1. Forsythe' s toy floating-point system: take /? = 2, r = 3, m = -1, and M — 2. Determine all the floating-point numbers and draw them as a scale on a straight line segment centred on zero. Does 0.1112~1 belong to this system? As we see from Exercise 1.3.1, there are many gaps between the floating point numbers. We therefore need a rounding function A which has the following properties: •A is defined for all E; •A leaves F(/3,r, m,M) invariant; •Let x E E. Let [/, /'] be the smallest interval containing x, and whose extremities are floating-point numbers. Then A(x) is equal to whichever of the numbers / and f is closest to x; •If x is equidistant from / and /', then A(x) is determined in a variety of ways, which can be dependent on the machine. Exercise 1.3.2. Calculate A( 1/3) in the toy system F(2,3, -1,2). Our machine can only recognize the numbers belonging to its floating-point system. We define arithmetic operations on the floating-point numbers by letting f Q f = A {f - f ) , s® r=A un, f ® f' = A ( f/ f) . Exercise 1.3.3. What are the results of the following operations in the toy floating-point system: 5 2 1.4. SUMMATION OF SERIES IN FLOATING-POINT NUMBERS 7 The program must be able to recognize the overflow and underflow situations. Exercise 1.3.4•Still in the toy system, calculate What can we say about the algebraic properties of the operations 0 and 0? Exercise 1.3.5. Simulate a small floating-point system in base fl = 10 by using a programming language such as PASCAL or FORTRAN and making a rounding procedure. Begin by fixing r = 3, m = -3, and M — 5. From this rounding procedure, program the floating-point arithmetic operations in F(10,3, —3,5). Return overflow and underflow when appropriate. A more elaborate version should allow the simulation of floating-point systems for various values of r, m, and M, with fl = 10. Passing to any base fl is more delicate, but represents an interesting exercise for those who can program. Remark 1.3.6. The unit of speed of a computer, used principally for scientific calculations, is the flops, or FLoating-point OPeration per Second. We would therefore talk about a machine calculating at 100 megaflops, that is 108 floating point operations per second. Traditionally, when we evaluate the efficiency of a scientific calculation algorithm, we only count the number of multiplications and divisions that it demands. It is more reasonable, with the current state of technology, to also count the additions and subtractions, since the relative time for multiplications and divisions has decreased. In the era of parallel machines, we cannot be content with evaluating speed in flops. We must also take account of the number of processors, and note that certain algorithms use the structure of the machine more efficiently than others. 1.4. Summation of series in floating-point numbers When we have a computer at our disposal we are tempted to calculate things that we did not know how to, or did not want to, calculate by hand. The preceding exercises have shown that operations on floating-point numbers suffer from significant arithmetic faults. We will see others which are more analytic. Let (un)n^o be a sequence of reals. We let E0 = j4(m0) and £n = £n_i 0 A (un) . This is therefore the sum of the numeric series whose general term is un. Exercise I.4.I. Show that if un tends to 0 as n tends to infinity, and if the sums £n stay below the level of overflow, the sequence of En is stationary from a certain point. 1. FLOATING NUMBERS 8 Exercise 1.4-2. We let un = l/n (harmonic series). Is this numeric series con vergent? Show that the partial sum of this series is equivalent to In n. Exercise 1.4-3. Show that, for n > /?r, the sequence of floating-point partial sums of the harmonic series is stationary. Find an upper bound for the partial sum thus obtained. Exercise 1.4-4•Calculate the sum of the harmonic series working in the sets F(10,3, -3,5) and in F(10,5, -5,5). This is the moment when the sequence of £n is stationary. Although all the series with terms tending to zero are convergent in a ma chine, their machine sum depends considerably on the floating-point system used. We would expect that a series with terms tending to zero is divergent if its sum depends on the machine. Exercise 1.4-5. The partial sums of a divergent series with positive terms exceeds, from a certain point, every given positive number. What must the minimum size of the mantissa be in base 10 so that the partial sums of the harmonic series exceed 100? On a computer doing 109 flops how much time, in years, is required to do this? (Take 1 year = 3 x 107 seconds and make suitable approximations.) Another peculiarity of the floating-point numbers is calculation instability. Recall that the series is convergent for all x £ E. Exercise 1.4-6. Program: calculate e5 and e~5, first using the exponential func tion of the chosen scientific programming language, and then in the floating-point numbers F(10,3, -3,5). Explain the difference between the relative error for e5 and the relative error for e~5. What happens if we reverse the summation order of the partial sums for e"5? What happens if we sum the positive and negative terms separately? 1.5. Even the obvious problems are rotten Consider the recurrence relation (1.5.1) ttn+i = (q+ -p. Exercise 1.5.1. Verify that the sequence, whose general term is the constant un — p/q, is a solution of eqn (1.5.1). Take q = 3 and p = 1 or p — 2, and program the above recurrence. What do you observe? Exercise 1.5.2. Take q — 4 and p — 1, 2, or 3, and do the same calculation. Are the phenomena the same? Exercise 1.5.3. Give an interpretation of the different behaviours observed. 1.6. EVEN THE EASY PROBLEMS ARE HARD 9 1.6. Even the easy problems are hard The case of solving the second degree equation (1.6.1) ax2 + bx + c = 0 allows us to see the difficulties linked to the orders of magnitude of the numbers with which we are working, and how the simplest nonlinear formulae can tie us up in knots. The roots of eqn (1.6.1) are given by the formulae ( 1.6.2) Xi -b + \Jb2 - 4ac 2a and —b — \]b2 —4ac 2a #2 Exercise 1.6.1. Write a program giving the roots of eqn (1.6.1), including the complex case and a test to ensure that a is not zero. Exercise 1.6.2. Let r be the maximum number of decimals representable in floating-point, and q the integer part of 1 + (r/2). We take a= 1, b = - 10\ c —1. Calculate the roots of eqn (1.6.1). Also calculate them by machine and compare the result. Exercise 1.6.3. To avoid the above difficulty, we can write the formulae (1.6.2) a little differently by noting that V y-'/z = 3/ - z Vv + Vz' Write the corresponding new formulae, which the program will apply if it detects the need after having made a test. Apply the new program to the following choices of coefficients, denoting the exponent of the largest power of 10 which can be represented by the machine by M, and denoting the integer part of 1 + (M/2) by n: a = 10"", (1.6.5) a = 1, Cr II 1 i—' O3 (1.6.4) c = -4 x 10", c = 10", cT 1 II -o a = 6 x 10", 6 = 5 x 10", 1 II CJ CO II <3 (1.6.3) 6 = 5, c = 1 - 10~r"1 Note that the case (1.6.3) can be solved by using a suitable scale. For the case (1.6.4) it is necessary to use more drastic means, for example, the change of variable y = l/x. Provide tests and modifications of the program to allow the solution of these problems. 1. FLOATING NUMBERS 10 Denote the smallest root of the equation (1.6.6) x2 - 2x+ 1 = e by X\(e) and the largest root by X2(e). Calculate the derivatives of x\ and X2 with respect to e, for e > 0. Explain the phenomena observed and describe the remedies that can be employed. 1.7. A floating conclusion If the floating-point numbers have so many faults, why not calculate with other representations of the numbers? In certain cases, it is wise to use exact rep resentations, such as rational numbers, or, more generally, a representation in some suitable field or ring. For applications arising directly from engineering science or from nature, such representations are often not convenient since they assume precise knowledge of the numerical data for the problem. In general, this knowledge is not accessible. There are, however, less obvious applications: if we seek the coefficients of a series that we know are rational a priori, it is natural to use an exact representation, that is, the rationals. More generally, formal calculation tools allow us to calculate the derivatives, and sometimes the integrals, of functions for which we possess an explicit expression, to find the explicit solutions of differential equations, and to give a large palette of tools which are later available for use in the heart of a scientific calculation program (it is possible for some formal calculation software to produce a procedure in FORTRAN, C, or some other language). It is particularly interesting if the pro gram relies on complex formulae which are difficult to check. In this case, the symbolic calculation tool can be remarkably effective, when it is reliable. But there is always the other side of the coin: the use of symbolic manipula tion software leads to extremely complex calculations more rapidly than the use of scientific calculation software, and is occasionally totally infeasible. Further more, this software is still new, and clearly less reliable than scientific calculation software. If, for example, we ask a formal calculation program to take the inte gral of a function which has distinct forms according to the interval considered, the program is frequently susceptible to giving the wrong answer. Similarly, for mal calculation programs are. often bad at simplifying complicated expressions. Finally, everyone who leaves the beaten track requires tailor-made code and there exist far fewer libraries of coherent programs than for scientific programming. For a numerical analyst, floating-point representation is the devil that we know. Formal calculation software is the devil that we know less well. Whatever happens, with these two devils, it is necessary to proceed with caution. 2 A flavour of numerical analysis Teaching numerical analysis, or, more generally, the analysis of calculation pro cedures, often has the reputation of being pointless. Chapter 1 will, I hope, have convinced the reader that it is useful to think about numerical methods. The present chapter introduces some mathematical techniques which are extremely common in numerical analysis. Since mathematics is certainly not a spectator sport, but more an activity where one only acquires skill and strength by doing it, studying proofs, and solving exercises, the format of Chapter 2 is similar to that of Chapter 1. It therefore consists of a self-guided visit in the garden of approximations of the continuous by the discrete. We place the discrete in the care of Don Knuth [38], and the continuous in that of the great Euler [25]. There are no surprising results in this chapter, just the comparison of powers and exponentials, and the construction of the logarithm and exponential func tions. It allows us to sum up results which are generally known to the reader, but by using an approach which is independent of the classical results of ele mentary real analysis. This approach rests on simple arithmetic identities and inequalities, and has a large place in familiar numerical analysis procedures, but is little used elsewhere. We must not be afraid of throwing ourselves ' in at the deep end'of indices and limits of sequences of functions. At the beginning, we will, perhaps, not feel completely at ease, but as we go along and use these tools, it will become beautiful and natural. 11 2. A FLAVOUR OF NUMERICAL ANALYSIS 12 2.1. Comparison of exponentials and powers There are numerous methods for proving the comparison theorems between ex ponentials and powers. The first exercises of this project will allow the reader to obtain these comparisons in a totally elementary manner, using only algebraic identities and simple inequalities. The results are eminently classical; the proofs guided by these exercises are perhaps not. We begin with some elementary identities. Exercise 2.1.1. Let m be a strictly positive integer. Let (using notation from [38]) (2.1.1) y(y+ 1 ) ■■■(y + m - 1),1. ym = Verify the following identities: (2.1.2) ym - -y ( l)m m (2.1.3) (y + z)m = J 2 CJm y ^ J z l 3=0 Hint: for identity (2.1.3), let F(y,z,m) = (y + z)m 3zJ. j =0 Calculate the difference F(y, 2 ,m )-F (y,2 ,m —1) with the aid of identity (2.1.2) and conclude by means of a recurrence on m. Exercise 2.1.2. Let x be a number (rational, real, or complex) which is different from 1. Calculate the sum n 5 (x, m, n) = ^ xkkm, k=1 using the fact that k ___ kw = m Y , j m~1 j =0 and changing the order of the summations. Deduce from this calculation that, if 0 ^ x < 1, then m\ 0 ^ S (x, m, n) ^ ----- — r v (1 — and (2.1.4) 0 nm xn+1 ^ m\ ( i^ y 2.3. DISCRETE APPROXIMATION OF THE LOGARITHM 13 Exercise 2.1.3. Show that for every y E ]x, 1[ there exists a number C such that, for every n ^ 0, we have the estimate nm xn ^ C yn. Show, by using eqn (2.1.3), that there exists a C' such that, for every n ^ 0 and for every N > n, N (2.1.5) x kkm ^ C ' ( n + l)w k=n+l Use the formula found in Exercise 2.1.2 and inequality (2.1.5) to calculate the limit of 5(x,m,n) as n tends to infinity. What can we say if x is a complex number of modulus strictly less that 1? The preceding exercises thus show, in an elementary way, that increasing exponentials dominate all power functions. 2.2. Convergence and divergence of classic series The divergence of the harmonic series is a well-known fact, but generally proved by comparison with the logarithm function. There follows a proof of this diver gence which is entirely independent of all knowledge of logarithms: Let (2.2.1) Hn = 1 + - + ... 4— • 2 n Exercise 2.2.1. Show that, for every n ^ 1, H' 2n —Hzn-l ^ —. Deduce the divergence of the harmonic series. The same technique can be used, suitably modified, to show convergence: Exercise 2.2.2. Show that, for every a > 1, the series oo E»-" n= 1 converges. Hint: do not delimit the packets exactly as in the preceding exercise. 2. 3. Discrete approximation of the logarithm This section gives a construction of the natural logarithm from scratch; we do not suppose that the reader has never seen a natural logarithm before, but we wish to lead him or her to explore known ground with new eyes. 2. A FLAVOUR OF NUMERICAL ANALYSIS 14 We call the largest integer less than or equal to the real x the ' floor'of x and denote it by We call the smallest integer greater than or equal to x the ' ceiling'of x and denote it by \x~\. These ideas differ from the widely known idea of the integer part, but are often rather easier to use. Recommendation: draw graphs of the floor, ceiling and integer part functions to see how they are different. We define the following functions for x ^ 1 : (2.3.1) L (x ,n) = J/fna.-| - Hn, (2.3.2) L (x, n) = Hynxj - ffn-i • Exercise 2.3.1. Show that L(x,n) is strictly positive if x > 1 and n ^ rio(x) = \{x — 1 )-1"|, and that L(x,n) is strictly positive if x ^ 1 and n ^ 2 . Exercise 2.3.2. Show that, for every real x > n, we have the inequalities 1 and every sufficiently large integer L (x, n) ^ L (x, 2n) ^ L (x, 2n) ^ L (x, n). Exercise 2.3.3. Show that the limits lim L(x, 2 n) n— >oc and lim L(x, 2 n) n— »oo exist and are equal to a function of x, which we will denote by L(x). Remark 2.3.4•Some simple numerical experiments with x = \ / 2 will convince the reader that the sequence of numbers (L(x, n))n>no^ is n°l monotonic and that, consequently, its convergence is not obvious. On the other hand, we have just shown that the sequence of numbers (L(x, 2n))2n>no^ is monotonic. We have thus extracted a subsequence whose convergence is easy to show. There are a lot of procedures for extracting subsequences. Here we have used an arith metic argument but we often call on a compactness argument. It would be a shame to be happy with only demonstrating the convergence of a subsequence: we will therefore show the convergence of the two sequences (L(x,n)) n>1 and (L(x,n))n>1, beginning with the convergence of (L(x,n) ) n > 1 when x is an in teger. Exercise 2.3.5. Show that, if p is an integer which is greater than or equal to then for every n ^ 1 L(p,n) ^ L (p, n -h 1). 2, Deduce from this that, for every integer p which is greater than or equal to the sequence (L(p, n)) converges to L(p). What is the value of L(l)? 2, The function L should be the Napierian logarithm. It is therefore necessary to verify that it has the additive property (2.3.3) We will first verify this on the integers: L(xy) = L(x) + L(y). 2.3. DISCRETE APPROXIMATION OF THE LOGARITHM Exercise 2.3.6. Show that, for all integers p, 15 and n greater than or equal to 2, L( pq,n) = (Lp, + ( , np). Deduce from this that L(pq) = L(p) + L(q), for all integers p and q greater than or equal to 1 . This property is immediately generalized as follows: Exercise 2.3.7. Let a: be a dyadic number with finite expansion: x — 2~lr > 1. Show that, for every integer p ^ 1 and for every integer n ^ /, we have L (% P ,2") = L (x, 2") + L (p,2n~lr), and deduce from Exercise 2.3.5 that L(px) = L(p) + L(x). Then show that L(x) — L(r) - IL(2) and, therefore, for all dyadic numbers x and y which are greater than or equal to 1 and have a finite expansion, L has the additive property (2.3.3). The function L has the additive property when it acts on dyadic numbers of finite expansion. How does it behave on reals greater than or equal to 1 ? A density and continuity argument gives the answer: Exercise 2.3.8. Show that, for all x and y such that 1 ^ x < y, and every integer n ^ n0(x), From this deduce the inequalities (2.3.4) -— - ^ L(y) V -^ x Show that eqn (2.3.3) holds for all reals x and y greater than or equal to 1 , by combining the continuity relation (2.3.4) and the additive property (2.3.3) already proved for dyadic numbers. Show that the function L is differentiable on [l,+oo). What is the value of its derivative? Now we are ready to show that every sequence L(x,n) converges to its limit: Exercise 2.3.9. Deduce from the inequalities (2.3.4) that and that, consequently, L ([xn\ + 1) - L (n+ 1) ^ - [znj) - (n). 2. A FLAVOUR OF NUMERICAL ANALYSIS 16 Use the properties (2.3.3) and (2.3.4), and the definition of the floor function to show that L (x)---- ^ L (x, n) ^ L (x ). n What can we say about L(x,n)? We have momentarily forgotten the logarithms of numbers less than they are treated in an analogous way: Exercise 2.3.10. Show that, for every x < we have the inequalities n- 1 1 1, but and for every sufficiently large n, - L (x, n) ^ —7— —- - L (x, 2n) 2n 1 1 _____ 1 - L (x, 2n) ^ £ |_xnj [2xn\ 2n - 1 1 —L (x, n). n —1 Exercise 2.3.11. Show that, if p and q are integers such that q < p, then L Conclude, by using the toolkit developed above, that L(x) — -L(l/x) for any strictly positive real x. Calculate the derivative of L at every positive real number. The function L thus obtained is the Napierian logarithm. It was introduced without calling on its definition as the integral of 1 /x which vanishes at x = 1 . Certainly, the integral definition has been hidden: without saying it, we approximated an integral with a finite sum by means of a rectangle method, allowing us to even drop the small pieces at the ends of the interval (draw a picture). It is recommended that the reader draw the areas defined by L(x,n) and L(x,n) for x > 1, and compare them with the area defined by f* dy/y. It will be useful to write, for example, 2.4. Comparison of means We are also going to construct the exponential function by elementary methods. To do this we will need classical inequalities between different means. The proof of these inequalities is very simple, and curiously dyadic—as the reader will appreciate. We recall some terminology: the arithmetic mean of a and b is equal to half of the sum of a and b; the geometric mean of a and b is the square root of the product of a and b; and the harmonic mean of a and b is the inverse of the arithmetic mean of 1/a and 1/b. 2.5. ELEMENTARY CONSTRUCTION OF THE EXPONENTIAL 17 Exercise 2-4-L Show that, for all positive real a and 5, we have the inequalities 2 (- + l) ^ V ^ b ^ - ( a + b). We are going to show that, for every integer n, and for every choice of strictly positive reals a\, a 2 ,..., an, we also have the following inequality: (2.4.1) / n^ y\ _1 \j=1 ai/ 1 n^ n \ Y ~ I^ y/ai «2 • • -a" ^ n Exercise 2.4-2. Show that inequality (2.4.1) is true when n is a positive integer power of 2 . Exercise 2-4-3. Let n < 2m. Show that we can deduce the second inequality in expression (2.4.1) for n from the corresponding inequality for 2m, on the condition that we choose b _ \ aj j ~\g if 1 ^ j ^ nifn + l and we take g to be the geometric mean of the aj. This elementary proof of the inequality between the geometric and the arith metic mean of n numbers can be found in [44, Chapter II, Section 5]. Exercise 2-4-4•Use a similar procedure to show the first inequality in expression (2.4.1) . 2.5. Elementary construction of the exponential The construction described below is entirely independent of the construction of the logarithm which was presented above. We will see in the differential equations chapter that this construction is a precursor to a lot of the ideas used for the numerical integration methods for ordinary differential equations. Just as the logarithm was obtained by an additive approximation of inte gration, the exponential will be obtained by a multiplicative approximation of integration. I owe the relation between the inequalities on the means and the approxi mation of the number e to the small book by P. P. Korovkin [55], and I have generalized this idea to the approximation of the exponential function. On the other hand, his proof of the inequalities on the means is more complicated than that given in the preceding section. Suppose that x is strictly positive, until indicated otherwise. We define the following functions: ( 7'\ ni 1 H--- ) m/ and _ / X \ —m E (x, m) = (1----- ) \ m/ 2. A FLAVOUR OF NUMERICAL ANALYSIS 18 It is proved in all calculus courses that E(*,m) and jE(-,m) converge to the exponential as m tends to infinity. Here, we shall see that these functions give a means to construct the exponential from scratch; moreover, we show that the exponential is the inverse function of the logarithm. Exercise 2.5.1. Let x > 0. Show that, for every integer m ^ E (x, m) 1 , we have E (x, m + 1 ). Hint: apply the second of the inequalities (2.4.1) with d\ = 1, 1 + x/ra, and n = m + 1 . = • • •= am+i = Exercise 2.5.2. In the same way, show that, if m > x, E (x, m) ^ E (x, m 4-1), using the first of the inequalities (2.4.1). Exercise 2.5.3. Show that, for every p > x which is sufficiently large and for every m greater than or equal to 1 , E(x,p) ^ E (x,m ). Exercise 2.5.4. Show that, for every m > x, we have _ _ 0 ^ E (x, m) —E_ (x, m) ^ — E (x, m ). m Exercise 2.5.5. Finally, show that, for every x > 0, the sequences (E(x,m)) and (i£(x,m)) m each have a limit as m tends to infinity and that these limits coincide. We will denote by E(x) the common limit of E(x, m) and E(x,m) as m tends to infinity. Exercise 2.5.6. Calculate E{0). Exercise 2.5.1. Show that, if p is greater than or equal to 1, then E{px) = E{x)p . Exercise 2.5.8. Deduce from the preceding result that, for every positive or zero rational p/q, we have Exercise 2.5.9. Show that, if x is positive or zero and if y is strictly greater than x, then (y ~ x) E (x, m) ( 1 4 - —) \ m/ ^ E(y,m) - E (x,m) ^ { y - x ) E ( y ,m ) ( l + —) \ m/ . 2.5. ELEMENTARY CONSTRUCTION OF THE EXPONENTIAL 19 From this, deduce that in the limit as m tends to infinity, (y —x) E(x) 4. E )- ) ^ ( (y (y). Show that E is differentiable for every x > 0 and right differentiable at 0. Cal culate E'(x). Exercise 2.5.10. Deduce from Exercises 2.5.8 and 2.5.9 that, for every x and y greater than or equal to 0, E satisfies the multiplicative property (2.5.1) E (x + y) = E (x) E (y) . From now on we are interested in the behaviour of E for x < 0. Exercise 2.5.11. If x is strictly negative, show that the two sequences (E(x,m))m>_ a xnd are well defined and converge monotonically to E(—x)~l . Then, show that the property (2.5.1) holds for all real x and y. The function E is therefore the exponential function. It remains to verify that it is the inverse of the logarithm. Exercise 2.5.12. Deduce from expression (2.3.4) that, for every x ^ 0, — — ^ L { l + x) ^ x . 1+ x Deduce from this inequality that — ^ «£(£(*,"•» and that consequently, for every x ^ 0, L(E(x)) = x. Show that this relation is still true for x < 0, and that we also have, for y > 0, E(L(y)) = y. Thus, we have constructed the exponential and logarithm functions by en tirely elementary methods and we have shown that one is the inverse of the other. Equally simple methods can be useful elsewhere— and this is what we are going to see in the next subsection. 2. A FLAVOUR OF NUMERICAL ANALYSIS 20 2.6. Exponentials of matrices Let AA.n be the set of real or complex square matrices with n rows and n columns. We will verify that M n can be equipped with a norm in the following way: let A— mid n M ill = max y^|Ajj| 1=1 Exercise 2.6.1. Show that the expression thus defined really is a norm, that is to say that, for all matrices A and B of M n and every scalar A, we have P||i ^ 0, P ||x = 0 positivity, if and only if A = 0, P + J9||i ^ P||i + ||£||i, the norm is positive definite, triangle inequality, ||A^4||i = |A| p||i, homogeneity of order 1 . The matrix norm that we have just defined is linked in a simple way to a certain n-component vector norm as follows: Exercise 2.6.2. If X\,X2,... ,xn are the coordinates of some vector x, we let Mi i=i ■ Show that we have the inequality (2.6.1) \Ax\i ^ ||A||i|x|i. Furthermore, show that, for every matrix A, we can always find a vector x / 0 such that |Ax|i = ||A||i|x|,. Then, show that |Ax|! P i l l = max x^O M l ' ( 2.6. 2) Exercise 2.6.3. Show that, if A and B are in M n, then ||A5||1 ^||A||1 ||5||1. Denote by 1 the identity matrix. Calculate ||l||i. Exercise 2.6.4- Let P(ra))m^ i be a sequence of matrices belonging to M n• Show that the sequence of norms P(m)||i tends to 0 if and only if each of the sequences {Aij(m))m^i tends to 0 for any i and j between 1 and n. 2.6. EXPONENTIALS OF MATRICES 21 Exercise 2.6.5. Show that the sequence of matrices A(m) converges if and only if it is a Cauchy sequence; in other words, if and only if, for every given number s > 0, we can find an integer M such that, for every m and p greater than or equal to M, ||A(m) —A(p)\\i is less than or equal to e. The result of the previous exercise states that M n is a complete space, i.e., all Cauchy sequences converge. The preceding ideas are sufficient to construct the exponential of a matrix A without using power series. We begin by showing an elementary criterion of invertibility: Exercise 2.6.6. Let A e M n be a matrix with corresponding norm strictly less than 1. Show, with the aid of inequality (2.6.1), that 1 — A is injective and, therefore, invertible. Also show, using eqn (2 .6 .2 ), that ||(1 - A)-1||i ^ (1 M U ) - 1. We will now show an elementary identity: Let C and D be matrices belonging to M n. Show that the following identity is valid for every integer p ^ 1 : (2.6.3) Cp - D p p-1 = ^ 2 °j (C - D) D p- l ~j . j =0 We now use the following approximations to the exponential: (2.6.4) E(A,m) and E (A,m) = H ") Exercise 2.6.7. Verify that the approximations (2.6.4) are well defined for every sufficiently large m. Compare ||£(,4,m)||i, ||E(A,m)||i, E{\\A\\i,m), and E{\\A\\i,m). Exercise 2.6.8. Show that, for every m greater than or equal to (2.6.5) \\E(A,2m) 1, ||i ^ Mill exP(Mlli) -E m Hint: apply eqn (2.6.3) with C = (1 -f A/2m)2 and D = 1 + A/m. This inequality allows us to show the convergence of a subsequence: Exercise 2.6.9. Verify that the sequence E(A, 2 m) has a limit, which will be denoted by E(A), as m tends to infinity. What can we say about the sequence E(A,2 m)? Exercise 2.6.10. Show that, even if A and B do not commute, we have (2-6.6) \\E(A) -E(B)\U ^ \\A - exp(max(P|| 1 ,||.B||1)). Hint: this is another application of identity (2 .6 .3 ). 2. A FLAVOUR OF NUMERICAL ANALYSIS 22 It remains to illustrate something of the multiplicative property of the expo nential: Exercise 2.6.11. Show that, if A and B commute, then (2.6.7) E {A + B) = E {A) E (B) = E { B ) E {A). Show that, if p is an integer greater than or equal to 1, then E(A)p = E(PA). Exercise 2.6.12. Prove that (2.6.8) P M - l - X I U < 2 M H ? e x p (HAH,). Hint: estimate \\E(A,2j ) - E(A, 2^"1)||i using expression (2.6.5) and sum with respect to j. The inequality (2.6.8) allows us to solve the remaining questions, and, in particular, the convergence of the sequence E(A,p) to E(A): Exercise 2.6.13. Show that, if p is an integer greater than or equal to 1 , \\E(A,p)-E(A)\\!^ 2exp(2 M i ) M | . P This allows us to verify the convergence of E(A1n) to E(A): Exercise 2.6.H. Show that there exists a constant C such that \\E(A,p)-E(A)\\1 ^ and estimate its value. This allows us to consider the real variable function which associates t with E(tA): Exercise 2.6.15. Show that the function t • -> E(tA) is infinitely differentiable and calculate all of its derivatives. We can also make some calculations of exponentials in practice: Exercise 2.6.16. Let A be a diagonal matrix. Calculate E(A). This was easy. Let us try a more complicated case: Exercise 2.6.17. Let A be a nilpotent matrix, that is, there exists p (less than or equal to n) such that A(p) is zero. Calculate exp(A). Show that, for any m, E(A, m) is a polynomial in A, whose degree is bounded independently of m. The case of nilpotent matrices is particularly curious, because we have an analogous behaviour for E(A,m): 2.6. EXPONENTIALS OF MATRICES 23 Exercise 2.6.18. Under the conditions of the preceding exercise, show that E(A,m) is also a polynomial in A whose degree is bounded independently of m, and find its limit as m tends to infinity. There are many ways of approximating a matrix exponential. Assume that the matrix A and the function B(t) are given such that ||J9(*) - (*^L)||i = te (t), where e(t) is a function of t which tends to 0 as t tends to 0 . Exercise 2.6.19. Show that there exists a constant C such that, for all m, Therefore, if we assume that e(t) is small compared with t, the approximation of E(A) by B(l/m)m is more precise than the approximation of E(A) by E(A,m) or by E(A,m). Exercise 2.6.20. For sufficiently small t, let (2.6.9) B(t) = ( l + (l - . Show that \\B(t)-E(tA)\\l = t3r](t), where 77 is a function of t which is bounded for sufficiently small t. What is the relationship between B(t/m), E(tA,m), and JE(tA,m)? Not only is this type of approximation more precise, but it allows us to obtain interesting information. Denote the Euclidean scalar product in W1 (or the Hermitian scalar product in C71) by (x,y). A matrix A is self-adjoint if, for any x and y, (Ax,y) = (x,Ay). It is skew-adjoint if, for any x and y, (Ax,y) = —(x,Ay). Finally, it is unitary if, for any x and y, (Ax, Ay) = (x,y). Exercise 2.6.21. Let A be a self-adjoint matrix. Show that E(A) is also selfadjoint. Let A be a skew-adjoint matrix. Show that the matrix B , defined by eqn (2.6.9), is unitary and, consequently, that E(A) is unitary. This is not all; we can estimate matrix exponentials when we have little information on the matrix itself: Exercise 2.6.22. Suppose that A is a self-adjoint matrix which is positive in the sense of quadratic forms, that is (2.6.10) (Ax,x) ^ 0, for any x in Rn (or C 1). Show that, for every x, we have 0 ^ (E (-A) x ,x) ^ (x ,x) . tffoL use E(A, n) rather than E(A, n) to obtain the answer. 2. A FLAVOUR OF NUMERICAL ANALYSIS 24 What use is such an example? It happens (and we will illustrate this in a particular case in the last chapter of this book) that a number of practical numerical analysis operators have the property (2.6.10), but their norm grows rapidly with the number of degrees of freedom, n, of the discretization; they arise from differential operators which have no chance of being bounded in a reasonable vector space. This does not prevent us from analysing them and doing the calculations which require them, even in practice. One property of the exponential is lost when we consider matrices: we still do not have E(A + B) = E{A)E(B). Below is an elementary example of this phenomenon. Let -(2 o) o)- Exercise 2.6.23. Calculate £(.4) and E(B). Is E{A)E{B) unitary? All of this section can be generalized without difficulty to the case of a unitary normed algebra, that is, a set A which has a vector space structure and, fur thermore, which possesses one multiplication which is distributive with respect to addition, and a norm for which this space is complete. 3 Algebraic preliminaries This chapter contains some elementary and some very elementary information. Normally, at level n one has not entirely assimilated the lessons of level n — 1. Therefore, the consultation of first year courses and books is strongly recom mended in case of difficulty. 3.1. Linear algebra refresher We assume the reader to be familiar with the ideas of a vector space over the field K, which could be the real field R or the complex field C, and linear mappings. We also assume the reader to be familiar with the ideas of linearly independent sets, spanning sets, dimension, and basis. 3.1.1. The matrix of a linear mapping Let V be a vector space of dimension n, and W a vector space of dimension m, over the field K. We choose a basis (v\,..., vn) of V and a basis (w\,..., wm) of W. Recall how we determine the matrix of the mapping / between these two bases: the images f(vj) of the vectors of the basis of V have the following decomposition on the basis Wji m f (vj) = i= 1 Every element x of V has the following unique decomposition: n x = J2 xi vjj =i By the linearity of /, [ J 2 xi vi ) = '52 ✓ n \j= l v n n m xi(fui) = 5Z ' 25 j= 1 ' j= l i= 3. ALGEBRAIC PRELIMINARIES 26 and, therefore, (3.1.1) The matrix of the mapping / between these two bases is therefore the following table of numbers belonging to K* ( An A21 Aw A22 • • • A in ^ *** A2n We frequently denote the element of the matrix A which is situated at the in tersection of the i-th row with the j-th column by Aij or (A)ij. The notational convention consists of first giving the row index and then the column index. The set of matrices with m rows, n columns and coefficients in IK, equipped with matrix addition and the multiplication by a scalar, forms a vector space denoted by A fm,n(K)- If m = n, this space is formed from square matrices and we denote it M n(lK). If the choice of the field IKis not important then we will denote it sim ply by M m,n or M n- The space .Mm>n(IK) therefore corresponds to the space of linear mappings from IKn to Km, once these two spaces are equipped with bases. Matrix multiplication corresponds to the composition of linear mappings. If B belongs to .Mm,n(lK) and A to -Mn,p(IK), we define the product AB by n In multiplication is an internal law. The identity matrix 1 is the identity in A1n(^); it will also be denoted by /, or by /n, if we wish to specify the dimension. A synonym of invertible is regular. A matrix which is not invertible is said to be singular. Remark 3.1.1 (Geometric Remark). If we identify V and Mn (or C71), and simi larly W and Rm (or C771), then we can say that the column vectors of the matrix A are the images of the vectors of the starting basis. We can perceive the matrix of the linear mapping as the assembly (in the right order) of the images of the basis vectors. Thus, the identity matrix is formed from the assembly of the canonical basis vectors of Kn . The above identification will be automatic in finite dimensions. Unless in dicated otherwise, we will identify linear form and row vector, linear mapping from the basis field IK to V and vector of V, and linear mapping from V to W \ and the matrix of this linear mapping. 3.1. LINEAR ALGEBRA REFRESHER 27 3.1.2. The determinant We recall the definition of the determinant as well as its principal properties. To do this, we identify families of n vectors of Kn = V and matrices formed from the assembly of these n vectors. The set An of alternate multilinear forms from Vn to K is formed of mappings /, from Vn to K, which have the following properties: / (V! + \v[,V2,...,Vn) = (f Vl, K", VA € K, / <rl),tV(2 ), • (v • • ,Va(n) )= where a is a permutation of n objects and e(cr) is its signature. It can be proved that An is a vector space of dimension 1. The determinant is the one of these alternate multilinear forms on Kn which has the value 1 on the canonical basis of V. The determinant of a matrix is the determinant of the family of its n column vectors. The determinant of a linear mapping from W to itself is the determinant of its matrix, provided that we choose the same basis before and after the mapping. It can be proved that this determinant does not depend on the chosen basis. The explicit expression for the determinant of a matrix is given by the for mula n (3.1.2) det tt i=l A = where the sum extends to all permutations a of n objects. We know that a square matrix A is invertible if and only if its determinant det A is not zero. Similarly, the determinant of a family of vectors is nonzero if and only if this family is a basis. We also know that the determinant of a product of matrices is the product of their determinants. However, this description is of little practical interest, since determinants are difficult to calculate. We happily calculate the determinant of a 2 x 2 numeric matrix, and we frequently solve a 2x 2 system by employing the Cramer formulae. For n ^ 4, it is already a bad idea to calculate a determinant with formula (3.1.2) . We will see in Chapter 9 how the tools used for solving linear systems give efficient methods for the numerical calculation of determinants. The determinant of a family of n vectors remains an interesting theoretical tool, if only because of its geometric significance. It is the volume of a paral lelepiped constructed on n vectors. Note that we allow this volume to be negative if the basis made out of these n vectors has the opposite orientation to the refer ence basis of the vector space. Consequently, the determinant of a matrix is the volume of the parallelepiped constructed from the column vectors of the matrix. e determinant therefore serves a crucial purpose when we make a change of 3. ALGEBRAIC PRELIMINARIES 28 variable in a multiple integration and, consequently, in the study of differential operators in geometry. For higher dimensions, the calculation of the determinant can be simplified by considering structure or symmetry. As an example of the exploitation of symmetry, we prove the following clas sical result, which concerns the determinant of Vandermonde: Lemma 3.1.2. The following identity holds: 1 1 Xi Xo rq 2 Z? X rj.n x x x0 l •• 1 z2 ' * xn -• rr2 x\ ■ ■ xn = 0^i<j^n \ ■ • Tn xn n ~ The determinant thus calculated is called the Vandermonde determinant. Proof. To verify this identity let P (x0,x i,...,«n) be the Vandermonde determinant. This is a polynomial in the n + 1 variables x o,xi,... ,xn. If two of the numbers Xi are identical, the determinant has two identical columns and, therefore, it vanishes. Consequently, P must be of the form Q (x0,X\, . . . ,Xfi) n ^ ' with Q being another polynomial. Moreover, examining the formula (3.1.2) indicates that each term of the determinant, given by Cram er' s formulae, is of global degree 1 + 2 + ... -f n = -n (n + 1). Consequently, P is a homogeneous polynomial of degree n(n + l)/2 with real coefficients and, therefore, Q is a real constant. It remains to calculate this con stant and, to do this, we note that, in the formula which gives the Vandermonde determinant, the monomial * 0 * l- "*n can only be obtained by taking a to be the identity permutation. The coefficient of this monomial in the determinant is 1. It remains to find the coefficient of this monomial in n ^ ■ o^i<j^n To obtain the term in x™ we must make xn appear n times and, therefore, we must take n —1 j[J[ (xn ~ %i) • i=0 3.1. LINEAR ALGEBRA REFRESHER 29 The coefficient of our monomial is therefore the coefficient of in JJ - . By induction, this is 1. □ We will come across this determinant again later when investigating the the ory of polynomial interpolation. To conclude this section, we recall the following two facts: •A nonzero matrix can have a zero determinant. For example, •Except in 1 dimension, the determinant of the sum of two matrices is not related to the sum of their determinants. The reader who is in doubt would do well to compare det(A + B) and det A + det J3, where B is given by 3.1.3. The fundamental theorem of linear algebra and its consequences The following theorem will be required frequently: Theorem 3.1.3. Let / be a linear mapping from a finite-dimensional space V to a finite-dimensional space W. Then, the dimension of the domain of / is equal to the sum of the dimension of the kernel of /, and the dimension of the image of/: dim Im / -I- dim ker / = dim V. o Here is an interesting exercise: Exercise 3.1.4. Prove, without using the determinant, that it is equivalent for a square matrix to be invertible, to have a left inverse, or to have a right inverse. Hint: if A has a left inverse B , then the kernel of A is reduced to 0 and the image of B is the whole space; use then the fundamental theorem of linear algebra, Theorem 3.1.3. We can immediately deduce from this the following corollary for the linear system Ax = 6, whose matrix A is square: 3. ALGEBRAIC PRELIMINARIES 30 Corollary 3.1.5. The following properties are equivalent: (i) For every vector 6, eqn (3.1.4) has at least one solution; (ii) There exists a vector b for which eqn (3.1.4) has at most one solution. The linearity of the problem implies that if we have uniqueness for a vector b then we have uniqueness for all vectors b. There is a simple criterion which guarantees that a square matrix is invertible: we say that an n by n matrix is strictly diagonally dominant if, for all j = 1,..., n, the following inequality holds: lajjl > Lemma 3.1.6. A strictly diagonally dominant matrix is invertible. Proof. We will show that the kernel of a strictly diagonally dominant matrix A is reduced to 0. We assume that there exists a vector x ^ 0 such that Ax vanishes, and let i be the index of the component of x with maximum absolute value. The i-th equation can be written AuXi = — ^ ^ AikXk. {k.k^i} Using the triangle inequality and the definition of i, this relation implies the inequality \Aa\ \x(\ ^ ^ ^ \x{\, {k.k^i} which, after division by |xj|, requires that \Au\ ^ ^ ^ {k.k^i} , contradicting the assumption on A. Therefore, due to the fundamental theorem of algebra, A is invertible. □ 3.1.4. Eigenvalues and eigenvectors Let / be a linear mapping from a vector space V of finite dimension n to itself. The spectrum of / is the complement of the set of complex numbers A for which / - XI is invertible. As we are in finite dimensions / — A/ is invertible if and only if its kernel is reduced to zero. Therefore, if A is in the spectrum of /, there exists a vector x (clearly nonzero) for which f (x) = Ax. 3.1. LINEAR ALGEBRA REFRESHER 31 In this case, A is an eigenvalue of /, and x an eigenvector of /. If we fix the basis of the space V , we can identify / with its matrix A in this basis, and thus we will refer to the eigenvalues of A and its eigenvectors. The characteristic polynomial of /, or of A, is P{X) = d et ( X I- A ) . The eigenvalues of A are the roots of its characteristic polynomial. It can be proved that every matrix is similar in the complex field to a triangular matrix of the form D + N, where D is a diagonal matrix and N is a nilpotent matrix (that is N n = 0) which commutes with D. 3.1.5. Scalar products, adjoints, and company Recall that a scalar product on a real (respectively, complex) vector space V is a positive definite bilinear (respectively, sesquilinear) form on V x V. On V = Mn the canonical bilinear scalar product is given by n (3.1.5) (x,x')v = j =i and on W = C™ the canonical sesquilinear scalar product is given by m (: t=l y,y')w= where yi is the complex number conjugate to yi. The Euclidean length of a vector x £ V is IMIv = \/(x,x)v , with an analogous definition in W . Given two finite-dimensional vector spaces V and IT, each equipped with a scalar product denoted by (*,*)v and respectively, the adjoint of a linear mapping / from V to W is the unique linear mapping f* from W to V such that (/ (*), y)w = (*> /* (y))v ' Vx g v, Vy g w. If V (respectively, IT) is identified with IKn (respectively, Km) equipped with its canonical basis, if the two spaces are equipped with their canonical scalar products, and if / has A £ M m,n as its matrix, then the matrix A* £ M n,m of /* is given by The transpose of a matrix A of m rows and n columns is the matrix AT of n rows and m columns whose coefficients are defined by (AT)ij = Aji. 3. ALGEBRAIC PRELIMINARIES 32 The transpose is the same as the adjoint in the real field; in the complex field, the adjoint is the conjugate transposed matrix. It is obvious that (a *) * = a , and that (3.1.6) (Ax, y)w = (x, A*y)v , Vx, Vy. Sometimes, the transpose of A is denoted by lA, which is somewhat incon sistent with the notation for the adjoint. Definition 3.1.7. A matrix A is said to be Hermitian if it is square and A* = A. We also say that it is self-adjoint. A matrix A is said to be skew-Hermitian if it is square and A* = —A. To be totally rigorous, it is not the same thing for a matrix to be Hermitian and to be self-adjoint: if we use a non-canonical scalar product, a matrix can be self-adjoint with respect to this scalar product without being Hermitian. As we do not treat mappings using these non-canonical scalar products in Kn in this book, we will use these two terms interchangeably without causing confusion. The spectral properties of Hermitian matrices are summarized by the follow ing assertion: Theorem 3.1.8. A Hermitian matrix is diagonalizable in an orthonormal basis and its eigenvalues are real. o If the matrix A has real coefficients and is Hermitian, it is said to be sym metric. In this case, the diagonalization basis can be taken to be real. A matrix A is said to be unitary if it is square and its inverse is equal to its adjoint: (3.1.7) A* A = AA*= I. The eigenvalues of A are complex with modulus 1 and A is diagonalizable in an orthonormal basis. If A is a unitary matrix with all real coefficients, we say that it is orthogonal. However, it should be noted that for an orthogonal matrix the eigenvectors are, in general, not real. To construct an orthonormal basis of eigenvectors of this matrix it is, therefore, necessary to use complex numbers. This idea will rarely be used in this course. We say that a linear mapping / from V to W conserves the Euclidean length if ll/(*)lk = IMIv, V* V. Lemma 3.1.9. Let / be a linear mapping with corresponding matrix A, which is assumed to be square. Then, A conserves Euclidean lengths if and only if it is unitary. 33 3.1. LINEAR ALGEBRA REFRESHER Proof. If A is unitary, we have, on using eqn (3.1.6), that (Ax, Ax) = (A*Ax,x) = (x,x). Therefore, A conserves lengths. Conversely, let A be a linear mapping matrix which conserves lengths. We then have the relation: (3.1.8) (Ax, Ax) = (x, x ) , Vx. To extract information from this relation we are going to pass to the polar form of the quadratic form x (Ax, Ax) - (x, x). We first of all remark that we have, on using the sesquilinearity of the scalar product, (Ax + Ay, Ax + Ay) - (Ax - Ay, Ax - Ay) = 4 3? (Ax, Ay), Vx, My. Similarly, (x + y, x -t- y) - (x - y, x - y) = 4 3? (x, y), Vx, My. It follows, by using eqn (3.1.8), that (3.1.9) 3i (Ax, Ay) = 3? (x, y), Vx, My. If we replace y by \y, the preceding relation becomes (3.1.10) 9 (Ax, Ay) = 3 (x, y), Vx, My. We deduce from eqns (3.1.9) and (3.1.10) that (3.1.11) (Ax, Ay) = (x, y), Vx, My. Using the definition of the adjoint we obtain ((A *A -I)x,y)= Q , Vx, My. which clearly implies that A*A = I. Since A is square, A* is also a right inverse of A. □ The Schur lemma states that, for every matrix A, there exists a unitary matrix U such that U~1AU is upper triangular. We call such an upper triangular matrix, which is unitarily equivalent to A, a Schur form. More generally, a square matrix is said to be normal if it commutes with its adjoint. Normal matrices are diagonalizable in an orthonormal basis. 34 3. ALGEBRAIC PRELIMINARIES 3.1.6. Triangular matrices Numerical analysts love triangular matrices because systems (3.1.4), where A is triangular, are very easy to solve. Suppose, for instance, that A is upper triangular. Then, the linear system in which we are interested is written a\\X i + ai2#2 + • • •+ Gl,n-l#n-l + n%n = ^1? 022^2 + • • •+ fl2,n-l^n-l + n%n = ^2? dn—l,n— l%n — 1 d n—l,n*^n ~ ^n—lj ® nn^n ~ ^n* Assuming that none of the are zero, we see that — X n —l — Xi = dnn 1 bn—1 Q'n—l,n% n &n—l,n—1 5 6l — 012^2 - ... - d in X n ------------------------ . d\\ Solving this is completely elementary and requires only a few operations, see Operation Counts 9.2.2 and 9.3.5. 3.2. Block matrices Block matrices generalize the concept of decomposition in coordinates. They come up almost systematically when we discretize differential equations in more than one variable. They are also very useful. 3.2.1. Block decomposition of a linear mapping or matrix Suppose that V, the domain space of the mapping / which we have already considered, decomposes into a direct sum of subspaces V ji v = m . | 3=1 Similarly, suppose that W, the image space of /, decomposes into a direct sum of subspaces Wf. M W = Q )W i. i—1 3.2. BLOCK MATRICES 35 Therefore, every element x of V decomposes on Vj in a unique way in the form N X— Xj . 3=1 The mapping x —► xj is the canonical projection of V to Vj and is denoted by Vj. The canonical injection J j is the mapping J j : Vj V, Xj -> Xj. In the same way, I* is the canonical injection from Wi to W, and Qi is the canonical projection from W to Wi. We can then write the linear mapping / by decomposing it on Vj and Wi as follows: N j=l and therefore the component of f(x) on Wi is given by N y i Q if {*?jxj) • j =i Let fij = Q' i °f °0j • Then, since Xj = VjX, M (3.2.1) N /(x) = 5 ^ 1 i ^ / y ( P j x). i=l j =1 Note the resemblance between eqns (3.2.1) and (3.1.1). Relation (3.2.1) is a generalization of eqn (3.1.1), and eqn (3.2.1) leads to eqn (3.1.1) provided that Vi = KVi, and similarly for the W j. By analogy with matrix notation, we can write the block decomposition of / in the form / / ll (3.2.2) /21 \/m 1 /12 flN ^ I f 2N ti f M2 * ■• I m n J 3. ALGEBRAIC PRELIMINARIES 36 We can obviously replace the fij by their matrices Aij, by fixing the bases of each of Vj and and obtaining the block decomposition of the matrix A of /: (3.2.3) A= ( An A21 A12 A22 • • • A in ^ *• * A2N \Ami Am2 AmnJ We say that the decompositions V = ® jVj and W = ® jWj equip / (or A) with a block structure. It so happens that certain problems impose a block structure. We are going to show that block decomposition makes matrix multiplication easier. 3.2.2. Block multiplication Suppose that / is a mapping from V to W which has a block structure given by eqn (3.2.2). In the same way, g is a mapping from W to X. If we suppose that L h=l with the canonical projections 7Zh and the canonical injections %/*, then g has the block structure (911 9\2 *•* 9im ^ 921 922 ••• 92M (3.2.4) \9l i 9L2 9lm J Clearly, 9hi — °9 ° > and, as in eqn (3.2.1), L 9(y) M ='^2'Hh'^2ghi(Q h=1 i= 1 Consequently, the composition g o / becomes L (:9 ° f) (y) = M E E 9hi h=1 L M =/i=l E n» E i=l L = h=l f ^ i=l M / M q<E V i'= l / N 53 9hi ( 53 A; i=l \j=l \ N Ej= lto (*#*)/ \ j 3') I ' / 3.3. 37 E X E R C IS E S F R O M C H A P T E R 3 since Qili> = 0 except when V = i. In this case is the identity of Wi. It follows that g o f has the block decomposition (g o f) hj with M (9° f) hj = ' E gh i°hi' h i= 1 This formula is obviously the generalization of the matrix product. If the matrix B of g is decomposed into matrix blocks Bhi, we obtain an analogous formula. To give an example, let A={Z Then, BA = ZZ)M d + B12A21 B\iA\2 -f B 12A 22 B n A is + B 12A 23) • It must be remembered from all of this that block multiplication is identical to normal matrix multiplication, provided that the dimensions of the blocks are compatible: block B^i must have as many columns as block Aij has rows. The Jordan decomposition of a square matrix A of order n makes use of the idea of block decomposition. Such a matrix has one Jordan form, that is, there exists one invertible matrix P such that P ~ 1AP = J, where /J(Xurn) 0 0 <7(^2>^2 ) • • • 0 0 0 \ J = 0 J(Ar_i,nr_i) 0 0 0 J (A^, J The Jordan blocks J(A,m) are m x m matrices of the form /A 0 J 1 A 0 1 • • • 0\ 0 (A, m) = \° A 1 0 \) We remark here that the Jordan form is numerically less stable than the Schur form, see [15]. 3.3. Exercises from Chapter 3 3.3.1. Elementary algebra Exercise 3.3.1. Let given by the formula Pb e a square matrix of order n with generic Pij = ^ij "h 3. ALGEBRAIC PRELIMINARIES 38 where a is a real number, Sij is the Kronecker delta, and m and r are integers between 1 and n. (i) Is the matrix P invertible? (ii) Calculate the inverse of P if it exists; (iii) Show that an endomorphism whose matrix is independent of basis is a scaling. To do this, use a change of basis whose transformation matrix is the form P. Exercise 3.3.2. Let A be a rectangular matrix of m rows and n columns. Show that, by applying a change of basis in the domain and image spaces, we can put A in the form n 0 • • • 0\ \ • •0 0 1 • (0) 0 • •V V> 0 • (0) (ob Reason in the following manner: •If the matrix A is identically zero there is nothing to do; •Suppose now that j\ is the index of the first column vector of A which is not identically zero. With the first change of basis we can move column ji to the position of column 1. Then, with a second change of basis, we replace the first column by a column containing only a 1 in its first row and zeros underneath; •The general case follows by induction. 3.3.2. Block decomposition Exercise 3.3.3. Let A be a square matrix of order n. Suppose that it has a block decomposition Aij of dimensions rij x n^, where 1 ^ i, j ^ m. We suppose, furthermore, that A is Hermitian and, therefore, diagonalizable in an orthonormal basis. Determine a sufficient condition such that the transformation matrix to this orthogonal basis has the same block structure. Exercise 3.3.4•Let A be a square matrix of order n. Suppose that it has a block decomposition Aij of dimensions n* x rij, where 1 ^ i, j ^ m. We say that A is block triangular if Aij = 0 for every index i and j such that i > j. Calculate the determinant of A as a function of the determinants of the blocks. Hint: begin with a matrix A decomposed into 2 x 2 blocks. There are at least two possible proofs of the result: 39 3.3. EXERCISES FROM CHAPTER 3 (i) By induction on the dimension of the diagonal blocks; (ii) Show that A is similar to a triangular matrix by a similarity which does not destroy the block structure. To complete the proof\ argue by induction on the number of blocks. Show by means of a counterexample that, in general, A\2 ) ± det (An) det (A22) - det (Ai2) det (A2i) . A22/ 3.3.3. Graphs and matrices A simple graph is defined as follows: we begin with a finite set X and a finite subset U of X x X. Suppose that U contains no element of the form (x,x). The graph is the pair G = (X, [/). If, for every (x,y) £ U, we also have (y,x) e 17, we say that the graph is symmetric. We will only consider symmetric graphs, which is the same as considering U as a set of subsets of two distinct elements of X. If u = (x,y) we say that x and y are the end-points of u. We say that (x, y) is an edge of the graph if (x, y) belongs to U. We say that two vertices x and y are adjacent if x and y are linked by an edge. Finally, we say that an edge u is incident to a vertex x towards the interior if u = (x, y) or if u = (y,x). The degree do{x) of a vertex x is the number of edges incident to x. We define the associated matrix of a graph to be the matrix A(G) such that Exercise 3.3.5. We say that a graph is properly coloured if we can colour each vertex in such a way that two adjacent vertices are always of different colours. It is fc-colourable if we can properly colour it with k colours. Show that if a graph is fc-colourable, its associated matrix has a block struc ture indexed by the colours. Exercise 3.3.6. The vertex-edge incidence matrix of a graph is the matrix R(G) defined by 1 if u has x as an end-point; («(c ))," ={ 0 otherwise. The degree matrix of G is given by Show that we have A(G) = R(G) R(G) t - D . 3. ALGEBRAIC PRELIMINARIES 40 Exercise 3.3.7. If G is a simple graph, we define a new graph S(G) by adding a new vertex £(u) for each edge u = (x ,y) and defining the new set of edges as the set of the (x,£(u)) and the (£(u),y). Show that the matrix associated to S(G) is given by Exercise 3.3.8. We define a matrix B in the following way: for every u G J7, consider a strictly positive number a(u) and let ^2{u:x is an end-point of u} a(u) if x = y \ if u = (x,y); otherwise. Exy — ^ a(u) 0 Show that the matrix thus defined is symmetric and either positive or zero. What is the dimension of its kernel? 3.3.4. Functions of matrices We begin by recalling some essential facts about square matrices. Every square matrix A of order n has a decomposition of the form A s T 4- N, where T and N commute together and with A and, furthermore, T is diagonalizable and N is nilpotent, that is, N n = 0. This fact is clear in the complex field; in the real field one has to be careful, and to proceed as follows: com plexify the problem and combine the projections on the generalized eigenspace corresponding to two complex conjugate eigenvalues. Exercise 3.3.9. Calculate the powers of A as functions of T and N. Exercise 3.3.10. Let Q be a polynomial of the variable x. Show that we have (3-3.1) q M) - o (T) + 'Q(T) N+ W 1 + " . + Q '"-1'm Exercise 3.3.11. Let / be a function of the variable x. Suppose that the domain of definition of / contains the spectrum of the operator A. If this spectrum is purely real, we will suppose that / is C n~l on an open subset of E containing this spectrum. If the spectrum of A contains points which are not real, we will suppose that / is n —1 times continuously differentiable with respect to the com plex variable x in a neighbourhood of the spectrum of A. The theory of analytic functions allows us to confirm that this condition is fulfilled if / is once contin uously differentiable, and at every point in the neighbourhood of the spectrum of A it has a convergent series expansion. Given a finite sequence we denote by diag(a«) the diagonal matrix whose elements are a i ,..., a^, in this 41 3.3. EXERCISES FROM CHAPTER 3 order. We extend the definition (3.3.1) by defining, for a diagonal matrix B, a diagonal matrix (3.3.2) f(B) = diag (/(&*)), for C = P~1B P a diagonalizable matrix, (3.3.3) f ( C ) = P - l f( B ) P , and, for any matrix A, (3.3.4) f" (T )N 2 f (A) = f (T) + /' (T) TV+ 7 ^ /• ("-!) m ATn"1 + ... + *...(n^ j ---- . Show that the definition (3.3.3) does not depend on the matrix P which is employed for the diagonalization. Exercise 3.3.12. Show that if g i

subset: commoncrawl