Numerical Mathematics 1 PDF Lecture Notes

Numerical Mathematics 1 Prof. Dr. Sabine Le Borne Institute of Mathematics Hamburg University of Technology October 11, 2024 Contents 1 Introduction 5 2 Finite precision arithmetic...

Numerical Mathematics 1 Prof. Dr. Sabine Le Borne Institute of Mathematics Hamburg University of Technology October 11, 2024 Contents 1 Introduction 5 2 Finite precision arithmetic 6 2.1 Machine numbers............................. 7 2.1.1 IEEE standard.......................... 7 2.1.2 Rounding to machine numbers.................. 8 2.1.3 Computer arithmetic....................... 9 2.1.4 Cancellation............................ 9 2.2 Condition of a problem and stability of an algorithm......... 10 2.2.1 Norms............................... 11 2.2.2 Landau symbols.......................... 12 2.2.3 Condition number of a problem................. 13 2.2.4 Stability of an algorithm..................... 15 2.2.5 Stopping criteria......................... 17 3 Linear systems of equations 19 3.1 Gaussian elimination........................... 19 3.2 Condition of a system of linear equations................ 23 3.3 Cholesky decomposition......................... 25 3.4 Elimination with Givens rotations.................... 27 4 Interpolation 29 4.1 Polynomial interpolation......................... 29 4.1.1 Direct interpolation........................ 30 4.1.2 Lagrange polynomials...................... 30 4.1.3 Interpolation error........................ 31 4.1.4 Barycentric formula........................ 33 2  CONTENTS 4.1.5 Newton’s interpolation formula................. 34 4.1.6 Interpolation with orthogonal polynomials........... 37 4.1.7 Connection between different polynomial interpolation schemes 38 4.1.8 Extrapolation........................... 39 4.2 Piecewise interpolation with polynomials................ 40 4.2.1 Classical cubic splines...................... 40 4.2.2 Piecewise linear splines...................... 45 4.2.3 Spline curves........................... 46 4.3 Trigonometric interpolation....................... 46 4.3.1 Trigonometric polynomials.................... 47 4.3.2 Fast Fourier transform (FFT).................. 50 4.3.3 Convolutions using FFT..................... 51 5 Nonlinear equations 55 5.1 Scalar nonlinear equations........................ 55 5.1.1 Bisection.............................. 55 5.1.2 Fixed point iteration....................... 56 5.1.3 Convergence rates......................... 58 5.1.4 Construction of one-step-methods................ 59 5.1.5 Multi-step iteration methods................... 62 5.2 Nonlinear systems of equations..................... 64 5.2.1 Multivariate Taylor expansion.................. 64 5.2.2 Fixed point iteration....................... 66 5.2.3 Newton’s method......................... 68 5.2.4 Quasi-Newton Methods...................... 69 6 Least squares problems 72 6.1 Linear least squares problems and the normal equations........ 72 6.2 Singular value decomposition (SVD).................. 74 6.3 Condition of the linear least squares problem.............. 79 6.4 Solving the linear least squares problem using QR factorization.... 80 6.4.1 Method of Householder...................... 81 6.4.2 Method of Givens......................... 83 6.4.3 Gram-Schmidt orthogonalization................ 84 3 CONTENTS 6.5 Nonlinear least squares problems.................... 85 6.5.1 Newton’s method......................... 85 6.5.2 Gauss-Newton method...................... 87 6.5.3 Levenberg-Marquardt algorithm................. 88 7 Eigenvalue problems 93 7.1 Review of the theory of eigenvalue problems.............. 93 7.2 Power methods.............................. 97 7.3 QR algorithm............................... 99 7.4 Computation of the Singular Value Decomposition........... 103 8 Differentiation 106 8.1 Finite difference approximation..................... 106 8.2 Algorithmic differentiation........................ 113 9 Quadrature 115 9.1 Newton-Cotes rules............................ 116 9.2 Error of Newton-Cotes rules....................... 118 9.3 Composite rules.............................. 119 9.4 Romberg integration........................... 121 9.5 Gauss quadrature............................. 122 9.6 Adaptive quadrature........................... 126 4 Chapter 1 Introduction Numerical Mathematics, often also called Scientific Computing, deals with the devel- opment and analysis of techniques (algorithms) to exploit computers in the solution of technical and scientific problems. These lecture notes are based on the book “Scientific Computing: An introduction using Maple and MATLAB” by W. Gander, M. J. Gander and F. Kwok (Springer Texts in Computational Acience and Engineering 11, 2014), available as an E-book in the TUHH library. 5 Chapter 2 Finite precision arithmetic We are used to compute within the set of real numbers, R. Computers, however, have only finite storage and hence in most cases use truncated approximations to real numbers. We will begin with an elementary example to illustrate an effect of such truncations before formally introducing floating point numbers and finite precision arithmetic. Example 2.1 (Numerical differentiation). Given a differentiable function f : R → R, its derivative at x0 is defined by f (x0 + h) − f (x0 ) f 0 (x0 ) = lim. h→0 h We may compute an approximation of the derivative as f 0 (x0 ) ≈ f (x0 +h)−f h (x0 ) for some small h ∈ R. For f (x) = ex , x0 = 1, hi := 10−i , i = 0,..., 15 we obtain the following approxima- tion errors f (1 + hi ) − f (1) r(h) := f 0 (1) − hi i 1+hi = e − 10 (e − e) = e − 10i · e(ehi − 1) We observe that f (x0 +h)−f h (x0 ) 9 f 0 (x0 ). On the computer, the smallest error is obtained for h ≈ 10−8. The reason is that a computer uses finite precision arithmetic which causes rounding errors that affect the computations. In some cases, rounding errors are negligible, in other cases, they can have catastrophic effects. 6  2.1. MACHINE NUMBERS 2.1 Machine numbers A computer can only store a finite set of numbers, hence most real numbers in R = (−∞, ∞) must be approximated by the closest machine number. Most computers use the so-called floating point representation: For B: base; m: mantissa; e: exponent; K∈ N: bias of exponent; m=D.DD... D e=DD... D with digits D ∈ {0,... , B − 1}, let e −1 `X a = ±m · B eN −K e where eN = (e`e −1... e1 e0 )N := ek B k ∈ N k=0 with given lengths (number of digits) `m and `e for the mantissa and exponent. Convention: If e a 6= 0, then the first digit in the mantissa m is nonzero. We denote the set of floating point numbers with base B, mantissa length `m , exponent length `e and bias K by M = M(B, `m , `e , K). (Some textbooks use a leading 0 in the mantissa, i. e. m=0.DD... D and decrease the bias K by 1, yielding the same set of machine numbers in slightly different representation.) 2.1.1 IEEE standard We now assume the base B = 2 and bias K := 2`e −1 − 1. Single precision: 32 bits, S,E,M ∈ {0, 1}, M(2, 24, 8, 127) (K = 27 − 1 = 127) S E E E E E E E E MM M =b (−1)S · 2e−127 · 1.m 8 bits 32-8-1=23 bits sign exponent e mantissa m Since the first digit in the mantissa is always 1 (for any nonzero number), it is not stored explicitly, hence 23 bits are sufficient for mantissa length `m = 24. Example 2.2. 0 1 0 0 0 0 0 0 1 0 1 0 0... b + 1.01 · 2129−127 0 = + 27 26 25 24 23 22 21 20 = 1.01 · 22 128 64 32 16 8 4 2 1 = 101 b · 22 + 0 · 21 + 1 · 20 =1 =5 normal numbers: 0  e  28 − 1 = 255 (all except e=0...0 or e=1...1) subnormal numbers: e = 0, in particular e = 0, m = 0 ⇒ e a := 0 exceptions: e = 255: m 6= 0: ea :=NaN (“Not a Number”) m = 0: ea :=Inf (“Infinity”) 7 CHAPTER 2. FINITE PRECISION ARITHMETIC Double precision: 64 bits, S, E, M ∈ {0, 1}, M(2, 53, 11, K), K = 210 − 1 = 1023 SEE EM M b (−1)S · 2e−1023 · 1.m = sign 11 bits 64-11-1=52 bits normal numbers: e 6= 0, e 6= 211 − 1 = 2047 2.1.2 Rounding to machine numbers We will introduce rounding for single precision numbers - a generalization to double precision or other formats is straightforward. Let a ∈ R \ {0} be a real number between the smallest positive and largest machine number: a = 1.a1 a2...a23 a24... · 2e−127 , ai ∈ {0, 1}. We round to (the largest of) the nearest machine number by defining  1.a1...a23 · 2e−127 : a24 = 0, a := f l(a) = (1.a1... a23 + 0. 0... 0 1) · 2e−127 : a24 = 1 e  | {z } 22 zeros where “rounding up” might result in an increase of the exponent (if m = 1.a1 a2... a23 = 1.1... 1). The rounding error is bounded as follows: Absolute rounding error:.. 0} 1 · 2e−127 |a − f l(a)| ≤ 0. 0|.{z 23 zeros Relative rounding error: |a−f l(a)| |a| ≤ 0. 0|.{z b 21 · 2−23 =: 21 epssingle.. 0} 1= 23 zeros Relative rounding error in double precision: ≤ 1 2 · 2−52 =: 21 epsdouble ≈ 1 2 · 2.2 · 10−16 Machine precision eps: historically: smallest positive machine number e a ∈ M such that 1 + e a 6= 1 spacing of machine numbers between 1 and B (if B = 2, then eps= 2−`m +1 , which is twice the smallest relative rounding error) colloquially: log10 eps is the “number of accurate digits” upper bound for (twice) the relative rounding error 8  2.1. MACHINE NUMBERS 2.1.3 Computer arithmetic + a, eb ∈ M, the exact result c := e Given two machine numbers e a −· eb may not be a / machine number and will have to be rounded to e c = f l(c). We define ⊕ + e e a − eb) =: e a b := f l(e c. ·  / c − c be the absolute rounding error and Let ra := e ra r := c be the relative rounding error. Then there holds: (i) |r| ≤ eps. (ii) c · r = ra = e c−c ⇔ cr + c = e c ⇔ c(1 + r) = e c + ⊕ ⇔ (e a −· eb)(1 + r) = e a eb /  |r| exact computer ≤ eps arithmetic arithmetic Standard model of computer arithmetic: (⊗: computer operation, ×: exact operation) a, b ∈ M : a ⊗ b = (a × b)(1 + r) for some |r| ≤ eps. 2.1.4 Cancellation Cancellation is loosely translated as “losing precision in subtracting similar num- bers”. Example 2.3. a = 1.23456789 round to 5 digits: e a = 1.2346 relative rounding error: | a−e a a |≤ 0.0001 1.23... ≤ 10−4 (4 correct digits) b = 1.23468 round to 5 digits: eb = 1.2347 relative rounding error ≤ 10−4 Compute the difference: b 1.23468000 eb 1.2347 - a 1.23456789 - ea 1.2346 0.00011211 0.0001 only 1 correct digit in the computed difference, 3 digits of accuracy are “lost”. 9 CHAPTER 2. FINITE PRECISION ARITHMETIC Recall Example 2.1 to compute an approximation of a derivative by a difference formula: Subtraction of similar numbers, f (x0 + h) − f (x0 ), leads to cancellation, i. e., rounding errors that limit the accuracy that may be obtained. 2.2 Condition of a problem and stability of an algorithm The condition number measures how sensitive a problem P is to small changes in the input data x (e. g. through rounding or measurement errors), assuming exact arithmetic. exact input x −−−−−−−−−−−−−−−−−→ exact output y = P(x) problem solution using exact arithmetic e −−−−−−−−−−−−−−−−−→ inexact output ye = P(e inexact input x x) If x ≈ xe but y 6≈ ye (w. r. t. their relative errors), then the problem is called ill- conditioned. The condition number describes the factor by which a relative error in the input might be multiplied to result in an upper bound for the relative error in the output. The stability of an algorithm describes the propagation of rounding errors through- out the individual steps of the algorithm using finite precision arithmetic. inexact input x e solution algorithm ye = P(e x) −−−−−−−−−−−−−−−−−−−−−−−→ using exact arithmetic exact input x algorithm with several steps ye = P(x) e b −−−−−−−−−−−−−−−−−−−−−−−−−→ and intermediate rounding An algorithm is called stable if it computes an approximate solution of relative error comparable to the exact solution to rounded input, i. e., ye ≈ e ye. Before we introduce formal definitions of conditioning and stability in §2.2.3 and §2.2.4, resp., we introduce norms in §2.2.1 in order to measure errors of quantities other than scalars. To describe the asymptotic behaviour of errors and complexities, we introduce Landau symbols in §2.2.2. 10  2.2. CONDITION OF A PROBLEM AND STABILITY OF AN ALGORITHM 2.2.1 Norms Definition 2.4: Vector norm A vector norm is a function k · k : Rn → R that satisfies 1) kxk > 0 ∀x 6= 0. 2) kαxk = |α| · kxk ∀α ∈ R, x ∈ Rn. 3) kx + yk ≤ kxk + kyk ∀x, y ∈ Rn (triangle inequality). Example 2.5. 1) 2-norm, spectral norm, Euclidian norm: v u n uX kxk2 := t |xi |2 i=1 2) Maximum norm, infinity norm: kxk∞ := max |xi | 1≤i≤n 3) 1-norm: n X kxk1 := |xi | i=1 Definition 2.6: Matrix norm A matrix norm is a function k · k : Rm×n → R that satisfies 1) kAk > 0 ∀A 6= 0. 2) kαAk = |α| · kAk ∀α ∈ R, A ∈ Rm×n. 3) kA + Bk ≤ kAk + kBk ∀A, B ∈ Rm×n (triangle inequality). Example 2.7. Let A ∈ Am×n. 1) If k · kRn , k · kRm are given vector norms, then kAk := sup kAxkRm kxkRn =1 is a matrix norm, called induced matrix norm or operator norm. For m = n and k · kRm = k · kRn = k · k2 , one obtains the spectral norm that satisfies p kAk2 = ρ(AT A) 11 CHAPTER 2. FINITE PRECISION ARITHMETIC where ρ(AT A) denotes the largest eigenvalue of AT A. For m = n and k · kRn = k · kRm = k · k∞ , one obtains the infinity or maximum row sum norm n X kAk∞ = max |aij |. 1≤i≤n j=1 For m = n and k · kRn = k · kRm = k · k1 , one obtains the 1-norm or maximum column sum norm n X kAk1 = max |aij |. 1≤j≤n i=1 2) The Frobenius norm is defined by v u m X n uX kAkF = t |aij |2. i=1 j=1 Remark 2.8. Given a finite dimensional vector space V (e. g. Rn or Rn×m ), all norms are equivalent, i. e., there exist constants m, M > 0 such that m · kxka ≤ kxkb ≤ M kxka ∀x ∈ V for two different norms k · ka , k · kb. 2.2.2 Landau symbols Landau symbols are introduced to describe the asymptotic behavior of functions as x → L for L ∈ R ∪ {±∞}. 12  2.2. CONDITION OF A PROBLEM AND STABILITY OF AN ALGORITHM Definition 2.9: Big-O, little-o Let f, g be two functions, let L ∈ R ∪ {±∞}. 1) We write f (x) = Ox→L (g(x)) or f (x) ∈ Ox→L (g(x)) if there exists a constant C ∈ R such that |f (x)| ≤ C · |g(x)| for all x in a neighborhood of L. 2) We write f (x) = ox→L (g(x)) or f (x) ∈ ox→L (g(x)) if |f (x)| lim = 0. x→L |g(x)| If the limit point L is clear from the context, we simply write O(g(x)) or o(g(x)). Example 2.10. 1) If we analyze the computational complexity of an algorithm depending on the problem size n (e. g. matrix vector multiplication of A ∈ Rn×n and x ∈ Rn ), we have L = ∞ (and obtain complexity O(n2 )). 2) If we analyze the (approximation) accuracy with respect to a parameter h (e. g. in the approximation of the derivative, see Example 2.1), we have L = 0 (and obtain O(h) in exact arithmetic). 2.2.3 Condition number of a problem Definition 2.11: Condition number The condition number κ of a problem P : Rn → Rm is the smallest number κ ∈ R such that |b xi − xi | kP(bx) − P(x)k ≤ ∀1≤i≤n =⇒ ≤ κ + o(). |xi | kP(x)k | {z } | {z } (componentwise) relative relative output error input error 13 CHAPTER 2. FINITE PRECISION ARITHMETIC Example 2.12 (condition of +, −, ·, /). Let x1 , x2 ∈ R, 1 , 2 ∈ R with |i | ≤  and bi := xi (1 + i ), i. e., xbix−x x i i = |i | ≤ . 1) Addition: P(x1 , x2 ) = x1 + x2 with x1 , x2 > 0, x = (x1 , x2 ), x b = (b x1 , x b2 ). |P(b x) − P(x)| x1 (1 + 1 ) + x2 (1 + 2 ) − (x1 + x2 ) = |P(x)| x1 + x2 x1 1 + x2 2 |x1 | + |x2 | xi >0 = ≤· = , x1 + x2 |x1 + x2 | i. e., addition is well-conditioned with κ = 1. 2) Subtraction: P(x1 , x2 ) = x1 − x2 with x1 , x2 > 0, x1 6= x2. |P(b x) − P(x)| x1 (1 + 1 ) − x2 (1 + 2 ) − (x1 − x2 ) = |P(x)| x1 − x2 x1 1 − x2 2 |x1 | + |x2 | = ≤·. x1 − x2 |x1 − x2 | Hence, κ = |x|x11|+|x 2| −x2 | may become very large if x1 ≈ x2 , i. e., subtraction of numbers of similar sizes is ill-conditioned (Cancellation! See Example 2.3). 3) Multiplication: P(x1 , x2 ) = x1 · x2 P(b x) − P(x) x1 (1 + 1 ) · x2 (1 + 2 ) − x1 x2 = P(x) x1 x2 = |1 + 2 + 1 2 | ≤ 2 + o(), i. e., multiplication is well-conditioned with κ = 2. 1 4) Division: P(x) = x 1 1 P(b x) − P(x) x(1+1 ) − x 1 − (1 + 1 ) |1 | (∗) = 1 = = ≤  + o(), P(x) x 1 + 1 |1 + 1 | i. e., division is well-conditioned with κ = 1. x (*) This inequality follows from a Taylor expansion of f (x) = 1+x about x0 = 0. Example 2.13 (divided differences). We assume that the function values are rounded to machine numbers, i. e., we compute divided differences for slightly perturbed input values fb(x0 + h) = f (x0 + h)(1 + 1 ), fb(x0 ) = f (x0 )(1 + 2 ). 14  2.2. CONDITION OF A PROBLEM AND STABILITY OF AN ALGORITHM Then there holds fb(x0 +h)−fb(x0 ) h − f (x0 +h)−f h (x0 ) f (x0 +h)−f (x0 ) h f (x0 + h)(1 + 1 ) − f (x0 )(1 + 2 ) − (f (x0 + h) − f (x0 )) = f (x0 + h) − f (x0 ) 1 f (x0 + h) − 2 f (x0 )  |f (x0 + h)| + |f (x0 )| 2|f (x0 )| = ≤ · ≈·. h(f (x0 +h)−f (x0 )) h f (x0 +h)−f (x0 ) h|f 0 (x0 )| h h Hence, the problem of computing divided differences has a condition number κ = 2|f (x0 )| h|f 0 (x0 )| , i. e., divided differences become highly ill-conditioned as h → 0. Theorem 2.14: Condition of a differentiable function evaluation The condition number of a function evaluation P : R → R, x 7→ f (x) for a 0 (x) differentiable function f at x is given by κ = xff (x). Proof. Excercise. Remark 2.15. For some problems, it is also suitable to replace the componentwise relative input error in Definition 2.11 of the condition number by the relative error with respect to some norm, i. e., kx−bxk kxk ≤ ε. 2.2.4 Stability of an algorithm An algorithm for solving a problem P : Rn → R is a sequence of (elementary) operations (which may all be considered as (sub-) problems themselves), P(x) = fn (fn−1 (...f2 (f1 (x))...)). In general, there exist different algorithms for a given problem. 1 Example 2.16. Problem P: Compute P(x) = x(1+x). Algorithm 1: f1 (x) = 1 + x f2 (x, y) = x · y f3 (x) = x−1 ⇒ P(x) = f3 (f2 (f1 (x), x)) = (x · (1 + x))−1 Algorithm 2: f1 (x) = x1 f2 (x) = 1 + x 1 1 f3 (x, y) = x − y ⇒ P(x) = f3 (f1 (x), f1 (f2 (x))) = x − x+1 Since only subtractions may be problematic in view of cancellation, we focus on these. Algorithm 1: Possible cancellation if x ≈ −1, otherwise unproblematic. Algorithm 2: Possible cancellation if x ≈ −1 OR if x1 ≈ 1+x 1 , i. e., if |x| is large. Hence we should always use Algorithm 1. 15 CHAPTER 2. FINITE PRECISION ARITHMETIC There exist two types of stability definitions, forward and backward stability. We are typically interested in forward stability, but often this is very hard to prove. Backward stability might be easier to show, and this is helpful since it implies forward stability. Definition 2.17: Forward stability An algorithm for a problem P is called forward stable if there exists a (not too large) constant C such that κ(f1 ) ·... · κ(fn ) ≤ C · κ(P). Interpretation: Relative (rounding) errors are enlarged by the steps of the algorithm not much worse (up to a constant C) than a relative input error solving P in exact arithmetic. Definition 2.18: Backward stability An algorithm for a problem P is called backward stable if the result yb obtained with input x is the exact result for slightly perturbed input xb, i. e., yb = P(x) b |b xi −xi | and yb = P(bx) with componentwise relative error bounds |xi | ≤ C · eps for a (not too large) constant C and machine precision eps. Here, P(x) b denotes the algorithm with rounding errors at substeps. Interpretation: The result of a computation on the computer is the same as the result of an exact computation using slightly perturbed input data. Lemma 2.19. If an algorithm is backward stable, then it is forward stable. Proof. Assume P b is backward stable, i. e., for yb := P(x) b b with yb = P(b there exists x x) xbi −xi and xi ≤ C · eps. Let y := P(x) be the exact result. Then it follows that kb y − yk kP(b x) − P(x)k |b xi − xi | = ≤ κ(P) · max kyk kP(x)k i |xi | | {z } ↑ relative error Definition of a produced by condition number the algorithm of a problem ≤ κ(P) · C · eps = C · κ(P) · eps. | {z } “unavoidable” output error, even if exact arithmetic had been used. 16  2.2. CONDITION OF A PROBLEM AND STABILITY OF AN ALGORITHM Example 2.20. Computation of an inner product 2 X T x y= xi y i f or x, y ∈ R2. i=1   ( x1 (1 + 1 ) · y1 (1 + 2 ) ) (1 + η1 ) + (x2 (1 + 3 ) · y2 (1 + 4 ))(1 + η2 ) (1 + η3 ) | {z } | {z } | {z } | {z } rounding rounding rounding rounding input x1 input y1 result of result of multipl. addition b1 · yb1 + x =x b2 · yb2 with x b1 = x1 (1 + 1 )(1 + η1 ), yb1 = y1 (1 + 2 )(1 + η3 ), x b2 = x2 (1 + 3 )(1 + η2 ), yb2 = y2 (1 + 4 )(1 + η3 ). Hence, since xbix−x i i = |(1 + ∗ )(1 + η∗ ) − 1| = |∗ + η∗ + ∗ η∗ | ≤ 2eps + o(eps), the backward stability condition is satisfied with C = 2. 2.2.5 Stopping criteria An important problem when computing approximate solutions using a computer, and hence finite precision arithmetic, is to decide when the approximation is accurate enough. This might be especially difficult when the exact solution is unknown. (For testing purposes, problems might be designed such that the exact solution is indeed known.) Test successive approximations If we produce a sequence xk of approximations, a commonly used stopping criterion is to check the absolute or relative difference of two successive approximations: |xk+1 − xk | < tol or |xk+1 − xk | < tol · |xk+1 |. Intuitively, we stop the iteration if not much progress is made anymore. However, there is no guarantee that we are close to the exact solution - the algorithm might just converge very slowly (or not at all). Check the residual If we try to find a solution x to the general problem f (x) = 0, we may compute the residual rk := f (xk ) of the current iterate xk and measure its size in some norm. 17 CHAPTER 2. FINITE PRECISION ARITHMETIC Example 2.21. For a linear system Ax = b, we compute the residual of f (x) = b − Ax as rk = b − Axk. While the residual of the exact solution satisfies r = 0, it is hoped, but not guaran- teed, that a small residual rk = b − Axk = Ax − Axk = A(x − xk ) indicates a small error xk − x. 18 Chapter 3 Linear systems of equations In this chapter, we will review Gaussian elimination for the solution of a linear system of equations. We will extract the underlying LU factorization and discuss pivoting strategies. We will show how the condition number of the matrix influences the expected accuracy of the solution of the linear system. For the special case of symmetric positive definite systems, we introduce the Cholesky factorization. An alternative elimination which is especially useful for sparse matrices (those ma- trices with many zero entries) is based on Givens rotations. In this chapter, as a basic prerequisite we assume that A ∈ Rn×n is a square, non- singular matrix. 3.1 Gaussian elimination Gaussian elimination reduces a system to upper triangular form which is then easy to solve. The necessary steps can be written in compact form as P A = LU where P denotes a permutation matrix (row permutations for pivoting), L is a unit lower triangular matrix responsible for the elimination of the lower triangular part of A resulting in an upper triangular matrix U. If a decomposition (factorization) P A = LU is available, a system Ax = b (with A = P T LU ) is solved in two steps: 1) Ly = P b (forward substitution), 2) U x = y (backward substitution). 19 CHAPTER 3. LINEAR SYSTEMS OF EQUATIONS Gaussian elimination without pivoting “Ax = b” is a compact notation for the system of equations a11 x1 + a12 x2 + · · · + a1n xn = b1 ,... an1 x1 + an2 x2 + · · · + ann xn = bn. In the i’th step of Gaussian elimination, entries below the i’th diagonal are elimi- nated by subtracting a multiple of the i’th row/equation: 1: for k=(i+1):n do 2: new row/eqn k = old row/eqn k - aaki ii · row/eqn i; 3: (in Python notation: A[k, :] = A[k, :] − A[k, i]/A[i, i] ∗ A[i, :]) 4: end for Algorithm 1: i’th step of Gauss elimination Gaussian elimination with pivoting If aii = 0, then the process breaks down due to a division by zero. We try to avoid such breakdowns (and more generally, division by small numbers) by introducing partial pivoting, i. e., we reorder the remaining equations (rows) before each elimi- nation step so that the entry with the largest absolute value in the current column i determines the row to be used to eliminate the next unknown. Unless the matrix is singular, the theory guarantees that there exists a non-zero en- try in the current column (see Theorem 3.3). Since rounding errors might produce small, non-zero entries in the singular case, we consider a pivot element to be zero if it is smaller than 10−4 · kAk1. Example 3.1.   x2 +3x3 = 5 ) 0 1 3 5 2x1 −x2 +x3 = 5 ⇒  2 −1 1 5  −3x1 +5x2 −7x3 = −6 −3 5 −7 −6   −3 5 −7 −6 1) Swap rows 1&3:  2 −1 1 5  (pivoting) 0 1 3 5   −3 5 −7 −6 2) Eliminate entries:  0 7 − 11 1  below a11 = −3 3 3 0 1 3 5   −3 5 −7 −6 3) Eliminate entries:  0 7 − 11 1  below a22 = 73 3 3 32 32 0 0 7 7 20  3.1. GAUSSIAN ELIMINATION Subsequent backward substitution yields the solution x1 = 3, x2 = 2, x3 = 1. In order to derive the LU -factorization, we first assume that no pivoting is necessary, i. e., the equations (rows) are already in the preferable order. The steps to eliminate entries below the i’th diagonal (see Algorithm 1) can be written in compact matrix notation as Li Ai−1 x = Li bi−1 (i−1) where Ai−1 := Li−1 ·... · L1 A with entries ajk and bi−1 := Li−1 ·... · L1 b are the updated matrix and right hand side after the previous (i − 1) elimination steps, and   1 ..  0.  ...... 0  ..  .    (i−1)  0 1  aji Li := ..  with `j,i := (i−1) , j = i + 1,... , n.  . −`i+1,i 1   aii ....   . 0.   ............  ......  0 ··· 0 −`n,i 0 ··· 0 1 Remark 3.2. i) Li is a unit lower triangular matrix.   1 ...       1  ii) Li is easy to invert: L−1 i = .... .   ..  ...    +`ji  ... 1 iii) L−1 i is again a unit lower triangular matrix. iv) The product of two unit lower triangular matrices is again unit lower triangular. In particular, there holds ··· 0   1 0.....   `21....  −1 −1 −1 L : = (Ln ·... · L1 ) = L1 ·... · Ln = ... . ...... 0 `n1 · · · `n,n−1 1 21 CHAPTER 3. LINEAR SYSTEMS OF EQUATIONS Since the lower triangular matrices Li transform A into upper triangular form, we obtain Ln... L2 L1 A = U | {z } =L−1 ⇐⇒ A = LU. With partial pivoting, i. e., permutations of the order of equations (rows in A), we obtain a triangular decomposition of P A where P is a permutation matrix. Theorem 3.3: LU decomposition Let A ∈ Rn×n be a non-singular matrix. Then there exists a permutation matrix P ∈ Rn×n such that P A = LU , where L and U are lower unit and upper triangular matrices obtained from Gaussian elimination. Proof. Thm 3.2. in. The advantage of storing the lower triangular factor L in addition to U lies in the fact that it can be reused for multiple right hand side vectors b: The decomposition P A = LU is computed only once (in O( 23 n3 ) operations), and multiple systems P A x = P bi , |{z} i = 1,... , k LU can be solved by Forward-Backward-Substitution, each one in O(n2 ). Theorem 3.4: Stability of Gauss elimination with partial pivoting Let A ∈ Rn×n be invertible and let L, b Ub be the numerically computed LU - factors using Gaussian elimination with partial pivoting. Then, for the ele- ments of Ab := P T · L b·U b , we have aij − aij | ≤ 2α · min(i − 1, j)eps + O(eps2 ) |b (k) where α := maxijk |b aij | (A b(k) is the computed matrix after k elimination steps). Gaussian elimination is stable in the sense that kA b − Ak∞ ≤ ρ p(n)epskAk∞ if the growth factor α ρ := max |aij | is not too large. Here, p denotes a cubic polynomial. Proof. Thm 3.3. in. 22  3.2. CONDITION OF A SYSTEM OF LINEAR EQUATIONS Remark 3.5. 1. For most matrices that occur in applications, the growth factor is acceptably small, i. e., Gauss elimination is backward stable in practice. 2. Backward and forward substitution can also be proven to be backward stable if the growth factor is not too large. 3. The computational complexity of Gaussian elimination amounts to 32 n3 +O(n2 ). 3.2 Condition of a system of linear equations Definition 3.6: Condition number of a matrix The condition number of a non-singular matrix A with respect to a matrix norm k · k is defined by K(A) := kAk · kA−1 k. Theorem 3.7: Condition of the solution of linear systems Consider two linear systems of equations Ax = b and Ab bx = bb satisfying aij = aij (1 + ij ), bbi = bi (1 + i ) b with |ij | ≤ e A , |i | ≤ e b. Then there exist A , b such that kA b − Ak ≤ A kAk, kbb − bk ≤ b kbk. (3.1) For the infinity or 1-norm, we have eA = A , eb = b , for other induced norms they differ only by a constant factor. If A is invertible and A · K(A) < 1 with condition number K(A), then there holds kb x − xk K(A) ≤ (A + b ). kxk 1 − A · K(A) kb x−xk The special case A b = A simplifies to kxk ≤ K(A) · b. Proof. From bb − b = Ab bx − Ax = (A b − A)b x − x) x + A(b we obtain   −1 b−x = A x −(A − A)b b x+b−b , b 23 CHAPTER 3. LINEAR SYSTEMS OF EQUATIONS and hence, using (3.1)   x − xk ≤ kA−1 k kA kb b − Akkb xk + kbb − bk ≤ kA−1 k(A kAkkb xk + b kbk). We estimate kb xk = kb x − x + xk ≤ kb x − xk + kxk, kbk = kAxk ≤ kAk · kxk, which leads to x − xk ≤ kA−1 k · kAk · (A (kb kb x − xk + kxk) + b kxk) = K(A) · (A + b )kxk + K(A) · A kb x − xk, and then K(A) kb x − xk ≤ · (A + b ) · kxk 1 − K(A)A which completes the proof. Remark 3.8. Theorem 3.7 implies that we have to expect an increase of the relative error by a factor up to K(A) in the solution. If A ≈ b ≈ eps (rounding error due to machine precision), the relative error of the solution may be of size eps · K(A). Example 3.9. a) Consider the linear system      1 1 x1 2 =. 1 0.999 x2 1.999 | {z } | {z } | {z } A x b
T
T Its exact solution is given by x1 x2 = 1 1. It can be checked that there holds      1 1 5 1.998 = , 1 0.999 −3.002 2.001002 | {z } | {z } x b b b i. e., a small change in the right hand side b leads to a much larger change in the solution x. Theorem 3.7 is applied as follows: There holds A = 0 (no perturbation in A) and kb − bbk∞ k(0.002, − 0.002002)T k∞ 0.002002 = T = = 0.001001 =: b. kbk∞ k(2, 1.999) k∞ 2 The condition number of A follows with   −1 −999 1000 A = =⇒ κ∞ (A) = kAk∞ kA−1 k∞ = 2·2000 = 4000. 1000 −1000 24  3.3. CHOLESKY DECOMPOSITION Theorem 3.7 provides an upper bound on the relative error of the type kx − x bk∞ ≤ κ∞ (A)b = 4000 · 0.001001 = 4.004. kxk∞ The actual relative error amounts to kx − x bk∞ k(−4, −4.002)T k∞ = = 4.002 kxk∞ k(1, 1)T k∞ which in this case is slightly smaller than the bound. b) Consider the matrix   21.6257 51.2930 1.5724 93.4650  5.2284 83.4314 37.6507 84.7163 A= 68.3400 3.6422 6.4801 52.5777 ,  67.7589 4.5447 42.3687 9.2995 the solution vector x := [1, 2, 3, 4]T and right hand side b := A · x. Solving x b := A\b (using MATLAB backslash) yields the numerical solution   0.999999999974... 1.999999999951... x b= 3.000000000039... ,  4.000000000032... i. e., a relative error of   2.6 kb x − xk 10−11  4.9 −11 ≈ ·  −3.9 ≈ 1.369 · 10.  kxk kxk −3.2 Theorem 3.7 guarantees (in exact arithmetic) a relative error smaller than K(A) · eps ≈ (6.01 · 105 ) · (2.2 · 10−16 ) ≈ 1.3 · 10−10 , i. e., about 10 correct digits in the numerical solution. Hence MATLAB backslash has computed the solution in a stable way. 3.3 Cholesky decomposition Definition 3.10: Symmetric positive definite matrices A matrix A ∈ Rn×n is called symmetric if AT = A, and positive definite if xT Ax > 0 ∀x 6= 0. Positive definite matrices occur frequently in applications and have special properties listed in the following lemma. 25 CHAPTER 3. LINEAR SYSTEMS OF EQUATIONS Lemma 3.11. Let A ∈ Rn×n be symmetric positive definite. i) If L ∈ Rm×n (m ≤ n) has rank m, then LALT ∈ Rm×m is symmetric positive definite. ii) Every principal submatrix of A is symmetric positive definite. iii) All diagonal entries of A are positive. iv) The largest entry of A in absolute value must be on the diagonal, and hence is positive. v) All eigenvalues of A are positive. Proof. Lemma 3.1. in. Theorem 3.12: Cholesky decomposition Let A ∈ Rn×n be symmetric positive definite. Then, in exact arithmetic, Gaussian elimination without pivoting does not break down, and only positive (i−1) pivots aii occur. Gaussian elimination results in A = LU with U = DLT for a diagonal matrix D = diag(d1 ,... , dn ) with only positive entries dj on the diagonal. 1 1 √ √ The upper triangular matrix R := D 2 LT with D 2 := diag( d1 ,... , dn ) leads to the so-called Cholesky factorization A = RT · R. Remark 3.13. The Cholesky factorization might be computed from an LU -factori- zation in 23 n3 + O(n2 ) operations. However, the factor R can also be computed directly in 13 n3 + O(n2 ) operations, i. e., in half the time (since only half of the n2 entries are needed):    r11 0 r11 · · · r1n A = RT R = ......... ...     r1n · · · rnn 0 rnn The entries r11 ,... , r1n of the first row of R are computed as follows: √ a1k r11 := a11 ; r1k =, k = 2,... , n. r11 Assuming entries of rows 1,... , j − 1 of R have already been computed, the next row rj,: follows from multiplication with the j’th unit vector ej from the left: eTj A = eTj RT R (aj1... ajn ) = (r1j... rjj 0... 0) · R j X = rkj · rk,:. |{z} k=1 k0 th row of R 26  3.4. ELIMINATION WITH GIVENS ROTATIONS Solving for the j’th row rj,: (all others are already known) yields j−1 X rjj · rj,: = aj,: − rkj · rk,: =: v T. k=1 The vector v can already be computed. Multiplication of both sides by ej gives 2 √ rjj = vj ⇒ rjj := vj. Now the remaining row rj,j+1:n of R can be computed: 1 T rj,j+1:n := √ vj+1:n. vj 3.4 Elimination with Givens rotations Gauss elimination transforms A into a triangular system by (successively) multi- plying with lower triangular systems, i. e., Ln... L1 A = U and then solving two | {z } L−1 triangular systems involving L, U. It may happen that L, U are worse conditioned than A. In this case, it would be preferable to transform A into triangular form using orthogonal transformations, i. e., QA = R with an orthogonal matrix Q (i. e., QT Q = I) and upper triangular matrix R and then solve A x=b |{z} QT Rx = b | · Q (from left) Rx = Qb ⇒ backward substitution For orthogonal matrices Q, there holds kQAk = kAk = kAQk when k · k is the Frobenius or Euclidian norm, hence K(R) = K(A). Givens rotation is one of several possibilities to compute such an orthogonal matrix Q. We illustrate the approach for a 2 × 2 matrix. Let   a11 a12 A=. a21 a22 We wish to find an orthogonal matrix     q11 q12 ∗ ∗ Q= such that QA = is upper triangular. q21 q22 0 ∗ 27 CHAPTER 3. LINEAR SYSTEMS OF EQUATIONS One way to obtain an orthogonal matrix Q results from a (plane) rotation setting   c −s Q= with s := sin α, c := cos α for some angle α. +s c ! One easily checks that QT Q = I. The requirement (QA)21 = 0 leads to s·a11 +c·a21 = 0 which is satisfied for a11 −a21 c := p , s := p 2. a211 + a221 a11 + a221 The 2 × 2 case generalizes to an n × n matrix where we use matrices i j   10 ··· ··· 0 ....  0.. .  .. 1    i  c ··· · · · −s   ....   . 1.   Qij,k =  .. . ....  . 1.     j s ··· ··· c    ..    1. ... .. . 0 0 ··· ··· 0 1 with −ajk aik s= q , c= q a2ik + a2jk a2ik + a2jk to eliminate an entry ajk , i. e. multiplication by Qij,k zeros out the entry ajk 6= 0 for a given k: a1k   ...    a   ik   ..  (Qij A)jk  = 0... 0 s 0... 0 c 0... ·  .  = s · aik + c · ajk = 0  a   jk  .  ..  ank for the above choice of s, c. The product Qij,k A will differ from A only in rows i and j while the product A(Qij,k )T will change only columns i and j of A. 28 Chapter 4 Interpolation Interpolation means to determine a function that satisfies certain given function values (of an otherwise unknown function) exactly. One has a choice of the type of functions to consider - the most common ones are polynomials. For periodic data, one may also choose trigonometric functions. Example 4.1. Interpolation of the function 1 f (x) = on [−5, 5] 1 + x2 with polynomials of degree 4 and 8 on an equidistant grid (left) and so-called Cheby- shev nodes (right). 1 1 0.8 0.8 0.6 0.6 0.4 0.4 0.2 0.2 0 0 -0.2 -0.2 -0.4 -0.4 -0.6 -0.6 -0.8 -0.8 -1 -1 -5 -4 -3 -2 -1 0 1 2 3 4 5 -5 -4 -3 -2 -1 0 1 2 3 4 5 4.1 Polynomial interpolation Polynomials are a common choice for interpolation since they are easy to evaluate, infinitely often differentiable and easy to integrate. Interpolation problem: Given n + 1 points (xi , yi ) ∈ R2 , i = 0,... , n, xi 6= xj ∀i 6= j, find a polynomial P ∈ Πn := {polynomials of degree ≤ n}, 29 CHAPTER 4. INTERPOLATION i. e., P (x) = a0 + a1 x +... + an−1 xn−1 + an xn , such that P (xi ) = yi ∀i = 0,... n. We will illustrate several approaches to solve the interpolation problem: i) direct approach using the Vandermonde matrix; ii) Lagrange polynomials; iii) Barycentric formula; iv) Newton’s formula; v) Orthogonal polynomials; vi) Aitken-Neville. 4.1.1 Direct interpolation Simply plugging in the n + 1 interpolation conditions leads to the following linear system of equations, ! P (xi ) = a0 + a1 xi +... + an xni = yi , i = 0,... , n, 1 x0 x20 · · · xn0      a0 y0 .....   a1   y1  . x1 x21      ..  = .. . (4.1)  .... ........  .  .  1 xn xn · · · xnn 2 an yn | {z } | {z } | {z } V := a:= y:= The matrix V is called Vandermonde matrix. If the nodes xi are pairwise distinct, then V is non-singular, i. e., the interpolation problem has a unique solution (see Theorem 4.2). However, if n is large or two nodes xi , xj are close to each other, V may become very ill-conditioned. Hence, this direct approach is not practical for a numerical solution of the interpolation problem. 4.1.2 Lagrange polynomials Let xi , i = 0,... , n, be n + 1 pairwise distinct nodes. We define the i’th Lagrange polynomial `i (x) ∈ Πn to satisfy the interpolation problem `i (xi ) = 1, `i (xj ) = 0 ∀j 6= i. Such a Lagrange polynomial is easily constructed, n Y x − xk `i (x) :=. k=0 xi − xk k6=i Now we can construct an interpolating polynomial that satisfies P (xi ) = yi , i = 0,... , n, as a weighted sum (linear combination) of Lagrange polynomials, n X P (x) := yi `i (x). i=0 30  4.1. POLYNOMIAL INTERPOLATION It is easy to check that P (x) satisfies the interpolation conditions. The following Theorem confirms that the interpolating polynomial using Lagrange polynomials is identical to the polynomial computed using the Vandermonde matrix in the direct approach, the two polynomials are just given in different representa- tions. Theorem 4.2: Existence and uniqueness of the interpolation poly- nomial Given n + 1 distinct nodes, the Vandermonde matrix in (4.1) is non-singular, i. e., there exists a unique interpolating polynomial P of degree ≤ n with P (xj ) = yj , j = 0,... , n. Proof. The construction of an interpolating polynomial P (x) using Lagrange poly- nomials proves the existence of an interpolating polynomial. To show unique- ness, assume that Q(x) is another interpolating polynomial of degree ≤ n. Then d(x) := P (x) − Q(x) is a polynomial of degree ≤ n satisfying d(xi ) = 0, i = 0,... , n. Since a non-zero polynomial of degree ≤ n can have at most n roots, we conclude that d(x) is the zero function, i. e., P (x) = Q(x). Example 4.3. Find a polynomial (of smallest possible degree) that interpolates the four points (−1, −4), (0, −1), (1, 2), (2, 11). Direct approach (using the Vandermonde matrix):      1 −1 1 −1 a0 −4 T 1  0 0 0 a1  =  −1  ⇒ (a0 , a1 , a2 , a3 ) = (−1,     2, 0, 1)T 1 1 1 1 a2   2  ⇒ P (x) = −1 + 2x + x3 1 2 4 8 a3 11 Lagrange polynomials: x · (x − 1)(x − 2) (x + 1)(x − 1)(x − 2) P (x) = − 4 · −1· (−1) · (−2) · (−3) 1 · (−1) · (−2) (x + 1) · x · (x − 2) (x + 1) · x · (x − 1) +2· + 11 · 2 · 1 · (−1) 3·2·1 Multiplying out all Lagrange polynomials and simplifying should lead to P (x) = −1 + 2x + x3. 4.1.3 Interpolation error In Theorem 4.4, we provide a representation of the interpolation error f (x) − Pn (x) for the polynomial Pn (x) that interpolates f (x) at distinct nodes xi , i = 0,... , n, which is later used to bound |f (x) − Pn (x)|. 31 CHAPTER 4. INTERPOLATION Theorem 4.4: Interpolation Error Let f ∈ C n+1 (R) with nodes x0 <... < xn. If Pn interpolates f at the nodes xj , then for every x ∈ R there exists a ξx ∈ R between the nodes x0 ,... , xn and x such that n 1 Y f (x) − Pn (x) = (x − xi ) f (n+1) (ξx ). (n + 1)! i=0 | {z } L(x):= Proof. For any node x = xk , the equality holds trivially (zero on both sides). Let x 6= xk , k = 0,... , n, be fixed and define the node polynomial n Y L(t) := (t − xi ) ∈ Πn+1 i=0 as well as the function F ∈ C n+1 (R), f (x) − Pn (x) F (t) = f (t) − Pn (t) − c · L(t) with c = ∈ R. L(x) Then there holds F (xk ) = 0 ∀k = 0,... , n as well as F (x) = 0. Therefore, F has at least n + 2 distinct zeros. By Rolle’s Theorem, the continuity of F 0 implies that there is at least one zero of F 0 between any two zeros of F. Thus, F 0 has at least n + 1 zeros. If we continue the argument, we conclude that F 00 has at least n zeros, F 000 has at least n − 1 zeros,... F (n+1) has at least 1 zero. Let ξ be such that F (n+1) (ξ) = 0, i. e., F (n+1) (ξ) = f (n+1) (ξ) − Pn(n+1) (ξ) − c · L(n+1) (ξ) = 0. (n+1) Pn (t) is polynomial of degree ≤ n, hence its (n+1)st derivative is zero: Pn (ξ) = 0. st L(t) is polynomial of degree n + 1, hence its (n + 1) derivative is (n + 1)!. We obtain 0 = F (n+1) (ξ) = f (n+1) (ξ) − c · (n + 1)!. Solving for “f (x) − Pn (x)” (contained in the constant c) yields the result. Remark 4.5. i) Theorem 4.4 yields a pointwise error bound. A more general re- sult is obtained using the maximum of f (n+1) (x) in the interval [a, b] of interest: 1 |f (x) − Pn (x)| ≤ · max |L(x)| · max |f (n+1) (x)| (n + 1)! x∈[a,b] x∈[a,b] 1 = · kL(x)k∞,[a,b] · kf (n+1) (x)k∞,[a,b] (n + 1)! for all x ∈ [a, b] (with xk ∈ [a, b] ∀k = 0,... , n). 32  4.1. POLYNOMIAL INTERPOLATION ii) The derivative f (n+1) may not be available which makes it difficult to obtain a bound on the interpolation error. iii) The node polynomial L(x) may become large if x lies near the boundary of the interval, especially when many nodes are used (→ see Example 4.1). An upper bound is given by kL(x)k∞,[a,b] ≤ (b − a)n+1. 4.1.4 Barycentric formula Lagrange polynomials are typically only used in theoretical considerations since they are rather expensive to evaluate (O(n2 ) operations). An alternative is the barycentric formula which requires only O(n) operations to evaluate the interpolation polyno- mial (after a one-time initial setup with costs O(n2 )). The formula is derived from the Lagrange interpolation polynomial as follows: Define coefficients 1 λi = (xi − x0 ) ·... · (xi − xi−1 )(xi − xi+1 ) ·... · (xi − xn ) 1 = Qn , i = 0,... , n. (xi − xk ) k=0 k6=i Let x 6= xi , i = 0,... , n. Then there holds n n X Y x − xj Pn (x) = yi i=0 x − xj j=0 i j6=i | {z } Lagrange representation Xn Y n = yi · λi · (x − xj ) i=0 j=0 j6=i n n Y X λi = (x − xj ) · · yi. j=0 i=0 x − xi If we apply the interpolation formula to y0 = y1 =... = yn = 1, we obtain the relation n n n Y X λi Y 1 1 = Pn (x) = (x − xj ) · =⇒ (x − xj ) = P n j=0 i=0 x − xi j=0 λi x−xi i=0 which leads to the formula n λi P x−xi · yi i=0 Pn (x) = n. λi P x−xi i=0 Once the coefficients λi are precomputed (in O(n2 )), we can evaluate Pn (x) for any x 6= xk in O(n). 33 CHAPTER 4. INTERPOLATION 4.1.5 Newton’s interpolation formula Newton’s formula is useful if we already have an interpolating (Newton) polynomial for n nodes and want to update it when more information/points become available. For the previous approaches, we cannot use the previously computed polynomial but would have to start from scratch. For pairwise distinct nodes xi , i = 0,... , n, we define Newton polynomials π0 (x) := 1, k−1 Y πk (x) := (x − xj ) ∈ Πk , k = 1,... , n, j=0 and represent a polynomial in the form Pn (x) = d0 π0 (x) + d1 π1 (x) +... + dn πn (x). The coefficients di could be successively computed by inserting the interpolation conditions, ! y0 = Pn (x0 ) = d0 + d1 · 0 +... + dn · 0 =⇒ d0 = y0 , ! y1 − d0 y1 = Pn (x1 ) = d0 + d1 · (x1 − x0 ) + d2 · 0 +... + dn · 0 =⇒ d1 = , x1 − x0... ! yn = Pn (xn ) = d0 + d1 · (xn − x0 ) + d2 (xn − x0 )(xn − x1 ) +... + dn (xn − x0 ) ·... · (xn − xn−1 ) n−1 P n−1 P yn − dj πj (xn ) yn − dj πj (xn ) j=0 j=0 =⇒ dn = =. (xn − x0 )... (xn − xn−1 ) πn (xn ) If a new point (xn+1 , yn+1 ) is added, we simply add dn+1 πn+1 (x) to the interpolating polynomial and compute the new coefficient dn+1 from the condition Pn+1 (xn+1 ) = yn+1. There exists a more efficient way to compute the coefficients di based on so-called divided differences which are defined recursively, f [xi ] := f (xi ) = yi , i = 0,... , n, f [xi+1 ,... , xi+k ] − f [xi ,... , xi+k−1 ] f [xi , xi+1 ,... , xi+k ] :=. xi+k − xi One typically writes all divided differences in a triangular array: 34  4.1. POLYNOMIAL INTERPOLATION x0 f [x0 ] x1 f [x1 ] — f [x0 , x1 ]...... — f [x1 , x2 ] — f [x0 , x1 , x2 ]............ xn f [xn ] — f [xn−1 , xn ]...... f [x0 ,... , xn ] It turns out that the coefficients di in the Newton polynomial are the top diagonal of the table of divided differences, dk = f [x0 ,... , xk ] i. e., n X k−1 Y Pn (x) = f [x0 ,... , xk ] · (x − xq ). | {z } k=0 q=0 dk | {z } Newton polyn. πk (x) The identity dk = f [x0 ,... , xk ] follows from the Aitken Lemma. Lemma 4.6 (Aitken). Given points (xi , yi ) ∈ R2 , i = 0,... , n, xi 6= xj for i 6= j, let P ∈ Πn be the interpolation polynomial satisfying P (xi ) = yi ∀i = 0,... , n. Let P[0,n−1] , P[1,n] ∈ Πn−1 be interpolation polynomials satisfying P[0,n−1] (xi ) = yi ∀i = 0,... , n − 1, P[1,n] (xi ) = yi ∀i = 1,... , n. Then there holds P[1,n] (x) · (x − x0 ) − P[0,n−1] (x) · (x − xn ) P (x) =. xn − x0 P[1,n] (x)·(x−x0 )−P[0,n−1] (x)·(x−xn ) Proof. Pe(x) := xn −x0 ∈ Πn satisfies the interpolation conditions since −y0 (x0 − xn ) Pe(x0 ) = = y0 , xn − x0 yn (xn − x0 ) Pe(xn ) = = yn , xn − x0 yj (xj − x0 ) − yj (xj − xn ) Pe(xj ) = = yj , j = 1,... , n − 1. xn − x0 Since the interpolating polynomial is unique, it follows that Pe(x) = P (x). We can now show by induction that the divided difference f [xi ,... , xi+k ] equals the (highest order) coefficient dk of the interpolating polynomial Pi,i+k (x) that satisfies Pi,i+k (xj ) = yj , j = i,... , i + k. 35 CHAPTER 4. INTERPOLATION For k = 0, this is true since a (constant) polynomial that interpolates a single point (xi , yi ) is given by Pi,i (x) = yi = f [xi ]. We now assume that the proposition is true for any interpolating polynomial through k consecutive points. By Aitkens Lemma, it follows that highest coeff. in Pi+1,i+k − highest coeff. in Pi,i+k−1 highest coeff. in Pi,i+k = xi+k − xi f [xi+1 ,... , xi+k ] − f [xi ,... , xi+k−1 ] = = f [xi ,... , xi+k ]. xi+k − xi Hence, the coefficient dn for the highest monomial xn in the Newton interpolation polynomial is the divided difference f [x1 ,... , xn ] − f [x0 ,... , xn−1 ] dn = f [x0 ,... , xn ] =. xn − x0 Remark 4.7. The table of divided differences contains all coefficients that are necessary to write down a Newton interpolation polynomial in Πk for any subset of consecutive interpolation points (xi , yi ), (xi+1 , yi+1 ),... , (xi+k , yi+k ). Evaluation of the Newton polynomial via the Horner scheme in O(n): Pn (x) = d0 + (x − x0 )(d1 + (x − x1 )(d2 +... + (x − xn−1 )dn )...). (4.2) Lemma 4.8. For some ξ between the nodes x, x0 ,... , xn , we have f (n+1) (ξ) f [x0 ,... , xn , x] =. (n + 1)! Proof. Let Pn be the interpolating Newton polynomial for (xi , f (xi )), i = 0,... , n. If we add another point, (x, f (x)), we obtain the Newton polynomial n Y Pn+1 (t) = Pn (t) + f [x0 ,... , xn , x] · (t − xj ) j=0 which satisfies n Y f (x) = Pn+1 (x) = Pn (x) + f [x0 ,... , xn , x] · (x − xj ) j=0 n Y =⇒ f (x) − Pn (x) = f [x0 ,... , xn , x] (x − xj ). | {z } j=0 interpolation error of Pn at x Comparison with the interpolation error in Theorem 4.4 implies that there exists a 1 ξ such that f [x0 ,... , xn , x] = (n+1)! f (n+1) (ξ). 36  4.1. POLYNOMIAL INTERPOLATION Lemma 4.9 (Symmetry of Divided Differences). The divided difference f [xi , xi+1 ,... , xi+k ] is a symmetric function of its arguments, i. e., f [xi+π(0) , xi+π(1) ,... , xi+π(k) ] = f [xi , xi+1 ,... , xi+k ] for any permutation π : {0,... , k} → {0,... , k}. Proof. Intuitiv: f [xi ,... , xi+k ] is the coefficient for xk in the polynomial that interpolates (xi , f (xi )),... , (xi+k , f (xi+k )). This coefficient is independent of the ordering of the nodes. The Aitken-Lemma can be used to evaluate the interpolant at z ∈ R without formally computing the interpolating polynomial explicitly. For i ≤ j, let Tij be the evaluation of the interpolating polynomial for points (xi , yi ), (xi+1 , yi+1 ),... , (xj , yj ) at z ∈ R. By the Aitken-Lemma, there holds T00 = y0 ,... , Tnn = yn , Ti+1,j (z − xi ) − Ti,j−1 (z − xj ) Tij = , i < j. (4.3) xj − xi The formula in (4.3) is called Aitken-Neville recurrence. Hence we compute all values in the table T00 T11 — T01... — T12 — T02......... Tnn — Tn−1,n —... T0n to obtain T0n = Pn (z). 4.1.6 Interpolation with orthogonal polynomials We define an inner product for functions, n X hpj , pk i = pj (xi ) · pk (xi ) i=0 for n + 1 pairwise distinct nodes xi , i = 0,... , n. 37 CHAPTER 4. INTERPOLATION Definition 4.10: Orthogonal polynomials A set of polynomials {pj }nj=0 is said to be orthogonal with respect to h·, ·i if pj has degree j for 0 ≤ j ≤ n, hpj , pk i = 0, j 6= k. We define the norm of a polynomial by kpk2 = hp, pi. Theorem 4.11 Let p−1 (x) ≡ 0, p0 (x) ≡ 1, β = 0 and pk+1 (x) = (x − αk+1 )pk (x) − βk pk−1 (x) hxpk ,pk i kpk k2 with αk+1 = kpk k2 and βk = kpk−1 k2. Then the polynomials pj (x) are orthog- onal. Proof. By induction, “Lanczos Algorithm”, see Theorem 4.4 in. Given an orthogonal set of polynomials {pj }nj=0 with kpj k = 1, the interpolation problem Find coefficients bj in P (x) = b0 p0 (x) +... + bn pn (x) such that P (xi ) = yi is easily solved,      p0 (x0 )... pn (x0 ) b0 y0 ....  ..  =! .. . .. .  .  p0 (xn )... pn (xn ) bn yn | {z } | {z } | {z } Q:= b y The matrix Q is orthogonal, i. e., QT Q = I and hence b = QT y, i. e., n X bi = pi (xj ) · yj. j=0 4.1.7 Connection between different polynomial interpolation schemes All discussed schemes compute the same interpolating polynomial - it is uniquely determined. The representations and algorithms to compute it differ since different bases are used for the vector space of polynomials up to degree n: Πn = {a0 + a1 x +... + an xn | ai ∈ R}. 38  4.1. POLYNOMIAL INTERPOLATION direct approach, Vandermonde: {1, x, x2 ,... , xn } Lagrange polynomials: {`0 (x), `1 (x),... , `n (x)} Newton interpolation: {π0 (x), π1 (x),... , πn (x)} orthogonal polynomials: {p0 (x), p1 (x),... , pn (x)} Computing the interpolation polynomial for data (xi , yi ), i = 0,... , n, we have n X n X n X n X i Pn (x) = ai x = yi `i (x) = di πi (x) = bi pi (x). i=0 i=0 i=0 i=0 4.1.8 Extrapolation Extrapolation is the same as interpolation, only the interpolant is evaluated at z outside the range of the nodes xj , i. e., z < xj ∀j or z > xj ∀j. W. l. o. g. we assume z = 0 (since we may shift the independent variable, x0 = x − z) and compute approximations using nodes hi > 0. The sequence {Pn (0)}, given by the values on the diagonal of the Aitken-Neville scheme, consists of evaluations of the interpolating polynomial through (hi , f (hi )), i = 0,... , n, at z = 0. The Aitken-Neville recurrence (4.3) simplifies to −hi Ti+1,j + hj Ti,j−1 Tij = (i ≤ j). hj − hi If we choose hi = h0 · 2−i , the recurrence becomes −2−i Ti+1,j + 2−j Ti,j−1 2j−i Ti+1,j − Ti,j−1 Tij = =. 2−j − 2−i 2j−i − 1 Such an extrapolation converges to f (0) if f indeed has an asymptotic expansion f (h) = a0 + a1 h +... + ak hk + Rk (h) with |Rk (h)| < Ck · hk+1. If f is known to have an asymptotic expansion of the form (only even exponents) f (h) = a0 + a2 h2 + a4 h4 +... one would interpolate ((hi )2 , f (hi )), i = 0,... , n, to obtain the recurrence formula −(hi )2 Ti+1,j + (hj )2 Ti,j−1 Tij =. (hj )2 − (hi )2 Again using hi = h0 · 2−i , we obtain 4j−i Ti+1,j − Ti,j−1 Tij =. 4j−i − 1 39 CHAPTER 4. INTERPOLATION 4.2 Piecewise interpolation with polynomials The interpolating polynomial may not always produce the result one would like. One (counter-) example was Runge’s function in Example 4.1, another one follows here. Example 4.12. For the given data x 1 2.5 3 5 13 18 20 , y 2 3 4 5 7 6 3 we compute an interpolating polynomial of degree 6 (left) and simply connect the points by straight lines (right). Interpolating polynomial of degree 6 Piecewise linear interpolant 14 14 12 12 10 10 8 8 6 6 4 4 2 2 0 0 -2 -2 -4 -4 0 2 4 6 8 10 12 14 16 18 20 0 2 4 6 8 10 12 14 16 18 20 The polynomial in the interval [5, 20] does not seem to be helpful while the piecewise linear function has only positive values (as anticipated) but is no longer differentiable at the interpolation nodes. In the following, we will introduce spline interpolation - piecewise polynomials that are differentiable (as many times as possible) at the nodes where the intervals connect. 4.2.1 Classical cubic splines Given nodes x1 < x2 <... < xn , the goal is to determine n − 1 polynomials Pi (x) = ai + bi x + ci x2 + di x3 for every interval [xi , xi+1 ], i = 1,... , n − 1. Since we have n − 1 intervals with 4 unknowns ai , bi , ci , di per interval, we have 4(n − 1) degrees of freedom/unknowns, i. e., we may prescribe as many conditions for the spline. The obvious conditions are the interpolating conditions, Pi (xi ) = yi and Pi (xi+1 ) = yi+1 , which are 2(n − 1) conditions in total. We will discuss two possibilities to prescribe the remaining 2(n − 1) conditions: 40  4.2. PIECEWISE INTERPOLATION WITH POLYNOMIALS 1. Defective spline functions: Prescribe conditions for approximate derivatives, Pi0 (xi ) = yi0 and Pi0 (xi+1 ) = yi+1 0 , where yi0 , yi+1 0 are approximations of derivatives at xi , xi+1 computed as difference quotients from the given data. 2. Classical cubic splines: Prescribe 2(n − 2) conditions 0 Pi−1 (xi ) = Pi0 (xi ), 00 Pi−1 (xi ) = Pi00 (xi ) for i = 2,... , n − 1 (all inner nodes), and 2 additional boundary conditions. Defective spline functions To begin, suppose that xi , yi , yi0 , i = 1,... , n, are given. Let hi = xi+1 − xi , i = 1,... , n − 1. We need to determine n − 1 cubic polynomials Pi that satisfy Pi (xi ) = yi , Pi (xi+1 ) = yi+1 , Pi0 (xi ) = yi0 , Pi0 (xi+1 ) = yi+1 0. One possibility to solve the problem is to set up Pi (x) = ai + bi x + ci x2 + di x3 , Pi0 (x) = bi + 2ci x + 3di x2 , and insert the conditions which leads to the system      1 xi x2i x3i ai yi 1 xi+1 x2i+1 x3i+1   bi  yi+1      2   =  0  (4.4)  0 1 2xi 3xi ci yi 2 0 0 1 2xi+1 3xi+1 di yi+1 for the coefficients ai , bi , ci , di. The unknown derivatives yi0 may be estimated through difference quotients using neighboring points: