Math 400 Assignment #3 -- Numerical Integration and Differentiation 1. For the function g x e x 2 2 Due: 6am April 09, 2016 use the centered difference formula with h 101 , 102 , , 1020 to create a table of estimates of g 1.4 . Include a column for the relative error of each estimate. 2. Using the data below estimate h 1.4 as accurately as you can using the centered difference formula with Richardson extrapolation. Also find the Lagrange interpolating polynomial for the data and evaluate L 1.4 . Compare your results. x h x 1.0 1.1 1.00000000 0.90483742 1.2 1.3 0.81873075 0.74081822 1.4 1.5 0.67032005 0.60653066 1.6 1.7 0.54881164 0.49658530 1.8 0.44932896 For more work with the data table estimate 1.8 1 h x dx using Romberg integration. 3. Estimate each of the following integrals to seven decimal place accuracy (absolute error less than 107 ) using Gaussian quadrature, the trapezoidal method, and Romberg integration. For each method state how many points were needed to achieve the desired accuracy. Also provide a graph of each function over an interval that includes the interval of integration. a) c) 12 0 1 0 sin x dx Sinc 2 x dx , b) 2 0.5 e x 2 dx 2 sin 2 x where Sinc 2 x 2 x 1 if x 0 if x 0 4. Find the exact values of the Gauss points 0 x1 x2 x3 4 and weights for n 3. Use the derivation given in class, but work over the interval [0, 4]. The Gauss points for a given value of n will be the zeros of the polynomial of degree n that is obtained by using Gram-Schmidt on the standard power function basis for the polynomials, and requiring each of the orthogonal basis polynomials to take on the value 1 when evaluated at x 1. An Introduction to Numerical Analysis Based on Lectures by David Sklar and Lecture Notes taken by Lauren Rath March 2, 2016 Parts I-VI Draft 3.04 1 Copyright These notes are Lauren Rath's undergraduate Applied Mathematics Project (Math 696 and 697). They began with Lauren's class notes as a student in David Sklar's Numerical Analysis course (Math 400) during the Spring 2014 semester. Over the Fall 2014 and Spring 2015 semesters Lauren created a rough draft using the LYX Document Processsor as a front end A to the L TEX typesetting system. The current draft is the result of multiple cycles of gure creation, editing, rewriting, and re-rewriting. It will provide the bulk of the reading material for the Numerical Analysis students at San Francisco State in Spring 2016. It is our hope that these notes will prove helpful to students interested in learning the basics of numerical analysis. Lauren Rath David Sklar 2 Contents 1 Systems of Linear Equations 1.1 5 Direct Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.1.1 10 1.1.2 Pivoting Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 1.1.2.1 Partial Pivoting . . . . . . . . . . . . . . . . . . . . . . . . . 13 1.1.2.2 Scaled Partial Pivoting . . . . . . . . . . . . . . . . . . . . . 13 1.1.3 Summary of Gaussian Elimination Algorithms . . . . . . . . . . . . . 14 1.1.4 Operation Counts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 1.1.5 A Brief Comment on Determinants . . . . . . . . . . . . . . . . . . . 18 1.1.6 Solving Multiple Systems with the Same Coecient Matrix . . . . . . 18 Iterative Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 1.2.1 Norms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 1.2.1.1 Vector Norms . . . . . . . . . . . . . . . . . . . . . . . . . . 25 1.2.1.2 1.2 Naive Gaussian Elimination . . . . . . . . . . . . . . . . . . . . . . . Matrix Norms . . . . . . . . . . . . . . . . . . . . . . . . . . 27 1.2.2 The Method of Iterative Renement . . . . . . . . . . . . . . . . . . 29 1.2.3 The Gauss-Seidel Method . . . . . . . . . . . . . . . . . . . . . . . . 30 1.2.4 Error Bounds and Condition Numbers . . . . . . . . . . . . . . . . . 32 1.2.5 Relaxation and Some Theory of Iterative Methods . . . . . . . . . . . 34 2 Eigenvalues and Eigenvectors 2.1 37 Dominant Eigenvalues and The Power Method . . . . . . . . . . . . . . . . . 2.1.0.1 The Power Method . . . . . . . . . . . . . . . . . . . . . . . 3 Interpolation 3.1 37 39 40 The Lagrange Interpolating Polynomial . . . . . . . . . . . . . . . . . . . . . 41 3.1.1 Errors in Lagrange Interpolation . . . . . . . . . . . . . . . . . . . . . 44 3.1.1.1 A Proof of Theorem 1.5 (Optional) . . . . . . . . . . . . . . 45 Hermite Interpolation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 3.2.1 Errors in Hermite Interpolation . . . . . . . . . . . . . . . . . . . . . 51 3.3 Linear Splines . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 3.4 Cubic Splines . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 3.5 A Brief Introduction to 2-D Spline Interpolation . . . . . . . . . . . . . . . 61 3.5.1 Bilinear Splines . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 3.5.2 Bicubic Splines . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 3.2 3 4 Numerical Dierentiation 66 4.1 The Central Dierence Formula . . . . . . . . . . . . . . . . . . . . . . . . . 66 4.2 Richardson's Extrapolation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 5 Numerical Integration, Estimating 5.1 b a f (x) dx 70 The Trapezoidal Rule . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 5.1.1 Big O notation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 5.2 Romberg Integration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 5.3 Hermite Quadrature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 5.3.1 Errors in Hermite Quadrature . . . . . . . . . . . . . . . . . . . . . . 77 5.4 Gaussian Quadrature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 5.5 Linear Algebra Review: Real Inner Product Spaces . . . . . . . . . . . . . . 80 5.6 Derivation of Gaussian Quadrature . . . . . . . . . . . . . . . . . . . . . . . 1 1 5.6.1 A Proof that 1 [Li (x)]2 dx = 1 Li (x) dx . . . . . . . . . . . . . . . 84 5.6.2 A Theorem on the Roots of Orthogonal Polynomials . . . . . . . . . 90 Summary of Numerical Integration Methods . . . . . . . . . . . . . . . . . . 91 5.7 6 Numerical Solutions to Ordinary Dierential Equations 89 93 6.1 Euler's Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 6.2 Modied Euler Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 6.3 RK-4 - The Basic Steps . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 6.4 Higher order Taylor Methods . . . . . . . . . . . . . . . . . . . . . . . . . . 101 6.4.1 6.5 Higher order Taylor Methods . . . . . . . . . . . . . . . . . . . . . . 101 Systems of First Order and Higher Order Ordinary Dierential Equations . . 103 6.5.1 6.5.2 RK-4 For n = 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 6.5.3 6.6 The Euler Method For n = 2 . . . . . . . . . . . . . . . . . . . . . . . 103 Second and Higher Order Equations . . . . . . . . . . . . . . . . . . . 104 Runge - Kutta Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 6.6.1 6.6.2 A General Class of Runge-Kutta Methods and a Subclass Containing RK-4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 6.6.3 Deriving Runge-Kutta Methods . . . . . . . . . . . . . . . . . . . . . 109 6.6.4 6.7 Background Materials for Deriving Runge-Kutta Methods . . . . . . 105 The (Optional) Messy Details for RK-4 . . . . . . . . . . . . . . . . . 111 Concluding Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 4 1 Systems of Linear Equations Many biological, economic, chemical, and physical phenomena are modeled using sophisticated mathematics. Linear and nonlinear systems of algebraic, ordinary and partial dierential equations arise frequently. While we can often prove that exact solutions exist, we can very rarely nd one. We generally nd approximate solutions by numerical methods. A key part of many of the methods for nding approximate solutions involves solving systems of linear equations. Our rst mathematical topic will be Methods for Solving Linear Systems. By Linear System we mean a system of n linear algebraic equations in n unknowns. The methods fall into two classes: Direct Methods, which include versions Gaussian elimination and matrix factorizations; and Iterative Methods, which generate sequences of approximate solutions that (hopefully) approach a good enough approximate solution. 1.1 Direct Methods We begin with an important preliminary example. Example 1.1. Consider the 2 2 system x1 + x2 = 2 x1 + x2 = 1 (1.1) where is a small nonnegative constant. Setting up the augmented matrix and row-reducing we get 1 1 2 1 1 2 (1.2) 1 1 0 1 12 using back substitution, we get x2= 12 1 and x 1 = 2 x 2 We give this solution a name 2 x2 x1 SOL1 = = 1 2 x2 1 Next consider the same system, but with the order of the equations reversed. x1 + x2 = 1 x1 + x2 = 2 (1.3) (1.4) where is a small positive constant. This system is mathematically equivalent (has the same solution as) the rst system. Setting up the augmented matrix and row-reducing we get 1 1 1 1 (1.5) 1 1 0 1 1 1 2 2 5 Solving, using back substitution, we get 2 x2 = 1 1 1 1 and x1 = (1 x2 ) We also give this solution a name SOL2 = x1 x2 1 (1 x 2 ) 1 = 2 1 1 (1.6) It can be shown, by algebraically simplifying both solutions, that SOL1 and SOL2 are alge1 braically equivalent. Clearly the linear system has a solution near for small values of 1 . More precisely, the true solution vector x= x1 x2 1 1 as 0 Using a computer to evaluate the solutions for smaller and smaller values of , the numerical 1 0 data below indicate that as 0, SOL1 , but SOL2 . 1 1 When = 101 SOL1 = 1.1111111111111111 0.8888888888888889 SOL2 = 1.1111111111111112 0.8888888888888889 SOL1 = 1.000010000100001 0.999989999899999 SOL2 = 1.000010000096196 0.999989999899999 SOL1 1.000000000100000 0.999999999999999 SOL2 = 1.000000082740371 0.999999999999999 = 105 = 1010 6 Notice that the the second solution for x 1 seems to contain random digits in the last eight decimal places. Continuing with even smaller values of we have = 1015 SOL1 = 1.000000000000009 0.999999999999999 SOL2 = 0.999200722162641 0.999999999999999 SOL1 = 1.000000000000000 0.999999999999999 SOL2 = 2.220446049250313 0.999999999999998 SOL1 = 1.000000000000000 1.000000000000000 SOL2 = 0.000000000000000 1.000000000000000 SOL1 = 1.000000000000000 1.000000000000000 SOL2 = 0.000000000000000 1.000000000000000 = 1016 = 1017 = 1020 The fact that SOL1 and SOL2 are mathematically equivalent does not mean that they are numerically equivalent. For a given value of we do dierent arithmetic in the evaluation of the dierent solutions. SOL1 is numerically stable SOL2 is NOT numerically stable The errors in the computation become particularly drastic in SOL2 due to the way in which the oating point hardware of modern computers stores numbers. The arithmetic of oating point numbers used in calculating SOL1 and SOL2 , as well as their storage, is determined by the standard titled the \"IEEE Standard for Floating-Point Arithmetic\". The standard species arithmetic storage formats, rounding rules, and arithmetic operations (and much more). The default format for oating point arithmetic in most computer languages is binary64 which uses sixty-four binary bits to represent a oating point number. Before looking at aspects of the IEEE standard, let's look at a computation using four digit decimal arithmetic with rounding. 7 Example 1.2. (Source: [BF]) Consider the 2 2 system 0.003x 1 + 59.14x 2 = 59.17 5.291x 1 6.130x 2 = 46.78 (1.7) which has exact solution x 1 = 10 and x 2 = 1. Solving the system using Gaussian elimination with four signicant digit arithmetic the augmented matrix is 0.003000 59.14 59.17 5.291 6.130 46.78 (1.8) (Note: To write a number to four signicant digits it is often convenient to rst write it in scientic notation. For example 0.003 = 3.000 103 = 0.003000.) Our rst pivot a11 is quite small in comparison with the other entries in the matrix. This leads to the large multiplier m21 = 5.291 = 1763.66 . . . 0.003000 which, to four signicant digits, rounds to 1764. Row reducing and rounding to 4 digits of precision at each step, we have row2 m21row1 = = 5.291 1764 3 104 6.130 1764 (59.14) 46.78 1764 (59.17) 5.291 5292 103 6.130 104300 46.78 104400 To four signicant digits the row reduced matrix is .003000 59.14 59.17 .0010 104300 104400 (1.9) (Note, had we used higher precision to do the arithmetic and then rounded to 4 signicant .003000 59.14 59.17 digits, we would have gotten .) .0000 104400 104400 and we see that error has been introduced by the small pivot .003. Solving using backsubstitution, we get x2 = 59.17 59.14 (1.001) 104400 = 1.001 and x 1 = = 9.713 104300 .003000 The exact solution of this system is x 1 = 10 and x 2 = 1. The small error of .001 in x 2 was 59.14 20, 000 during the back-substitution introducing a large error into multiplied by .003000 the approximation for x 1 . If you begin the row reduction by rst interchanging the rows, and still work to 4 signicant digits, you should get the solution x 1 = 10.00 and x 2 = 1.000. 8 It is clear from Examples 1.1 and 1.2 that how we organize a computation, how much precision we are working with, and how we round all aect the accuracy of our approximations. To get consistent results on dierent computer systems it is important that precision and rounding are handled according to widely accepted standards. The IEEE Standard for Floating-Point Arithmetic (see the Wikipedia pages on this standard) was developed for this purpose. The current version is IEEE 754 2008. It provides ve basic formats for numerical representation. The default format for oating point arithmetic in most computer languages today is binary64 which uses sixty-four binary bits to represent a oating point number, giving about 16 signicant decimal digits precision. Each format includes a nite number of digits, which can be base 2 (binary) or base 10 (decimal). A oating point number is described by three integers 1. s is 0 or 1, and it determines the sign of the number 2. c is called the signicand or coecient 3. q is an exponent. The numerical value of a normal number is given by (1) s cbq For example, in decimal format the number 46.875 is described by (1) 0 46875 (10) 3 with s = 0, c = 46875, b = 10, and q = 3. For a more interesting binary example consider the number with b = 2, s = 1, c = 1.01110111002 , and q = 001012 . Stored in a 16 bit computer memory word in binary16 oating point format, it would look like 1 1 0 1 0 0 0 1 1 1 0 1 1 1 0 0 15 14 13 12 11 10 9 8 7 6 5 4 3 2 1 0 where the second row numbers the 16 bits from 0 to 15 (right to left). Bit 15 is the sign bit s = 1, bits 14 through 10 contain the biased exponent in base 2 with bit 14 being the most signicant. The 10 bits 9 through 0 represent the binary number in a normalized form 10 1.b9 b8 b7 b6 b5 b4 b3 b2 b1 b0 = 1 + i=1 b10i 2i . (Note: In binary64 bit 63 is the sign bit, 62 to 52 the exponent bits, and 51 to 0 the signicand bits.) (1) s = (1) 1 = 1 q = 101002 bias = 101002 24 1 = 20 15 = 5 b9 b8 b7 b6 b5 b4 b3 b2 b1 b0 = 0111011100 9 which corresponds to c = 1.01110111002 = 1 + 22 + 23 + 24 + 26 + 27 + 28 And the represented number (in base ten) is (1) s cbq = (1) 1 + 22 + 23 + 24 + 26 + 27 + 28 25 = 25 + 23 + 22 + 21 + 21 + 22 + 23 = (32 + 8 + 4 + 2 + .5 + .25 + .125) = 46.875 The binary and decimal representations in our examples are exact. Not all numbers, not even all rational numbers, can be represented exactly in nite precision oating format. 1.1.1 Naive Gaussian Elimination Naive Gaussian elimination is an algorithm for solving systems of linear equations. We review it through an example which we also use to introduce a very explicit notation for describing the steps of the algorithm. Example 1.3. Consider the 3 3 system x1 4x 1 2x 2 8x 3 = 8 2x 2 + x3 = 0 + 5x 2 + 9x 3 = 9 We can think of the system as the matrix equation Ax = b such that 0 2 8 x 1 1 2 1 , x = x 2 , b = A= 4 5 9 x3 8 0 9 We form the augmented matrix 0 2 8 8 a 11 a12 a13 a14 A : b = A = 1 2 1 0 = a21 a22 a23 a24 4 5 9 9 a31 a32 a33 a34 and nd the rst row with a non-zero entry in the rst column, in our case R2 . We then interchange R1 and R2 . Now let A(1) 1 2 1 0 a (1) a (1) a (1) a (1) 11 12 13 14 (1 (1 (1 (1 = 0 2 8 8 = a21) a22) a23) a24) 4 5 9 9 a (1) a (1) a (1) a (1) 31 32 33 34 10 (1 Next we row-reduce the matrix until it is in row echelon form. As a21) = 0, we only need to (1 eliminate a31) in the rst column. We perform the operation R3 (1 a31) (1 a11) R1 R3 and A(1) is reduced to A(2) Next we compute R3 1 2 1 0 a (1) a (1) a (1) a (1) 11 12 13 14 (2 (2 (2 = 0 2 8 8 = 0 a22) a23) a24) 0 3 13 9 0 a (2) a (2) a (2) 32 33 34 (2 a32) (2 a22) R2 , and nd A(3) 1 2 1 = 0 2 8 0 0 1 0 8 3 Once the matrix is in row echelon form, we can use back substitution to solve for x. 3 =3 1 1 1 x 2 = (8 + 8x 3 ) = (8 + 8 (3)) = 16 2 2 1 x 1 = (2x 2 x 3 ) = 1 (2 (16) 3) = 29 1 x3 = The solution vector is 29 x = 16 3 It can be checked by computing r = b Ax which should be equal (or very close) to 0. We now formalize our notation. We begin by noting A(1) corresponds to the matrix A (1) (2) (possibly after a row interchange) with a non-zero pivot a11 , A corresponds to the partially (2 reduced matrix with a non-zero pivot a22) , etc. A(1) (1) (1) (1) (1) a 11 a12 a13 a14 = a (1) a (1) a (1) a (1) 21 22 23 24 (1) (1) (1) (1) a 31 a32 a33 a34 (1) (1) (1) (1) a 11 a12 a13 a14 A(2) = 0 a (2) a (2) a (2) 22 23 24 0 a (2) a (2) a (2) 32 33 34 11 (1) (1) (1) (1) a 11 a12 a13 a14 A(3) = 0 a (2) a (2) a (2) 22 23 24 (3) (3 0 0 a33 a34) The matrix superscript is incremented by one after any rows interchanges, but before any (2 other row operations are performed. For instance, in the previous example a33) = 13. In general, (1) (1) (1) a 11 a12 a1 k 0 a (2) a (2) 22 2k . . .. .. . . . . . . = . (k) . 0 ak k . . . . . . . . . . (k) 0 0 a kn ( ( a11) a11) 1 n n+ (2) (2) a2n a2n+1 . . . . . . (k) (k) a kn a kn+1 . . . . . . (k) (k) ann ann+1 If A is non-singular, the row-reduced matrix is (1) (1) (1) a 11 a12 a1 k 0 a (2) a (2) 22 2k . . .. .. . . . . . . (n) = A . .. (k) . . ak k . . .. . . . 0 ( ( ) a11) a11n+1 n (2) (2) a2n a2 n+1 . . . . . . (k) (k) a kn a k n+1 . . .. . . . . . (n) (n) 0 ann an n+1 A(k) (k) with a k k (k) with a k k 0. 0 for k = 1, . . . , n. Solving by Back Substitution we have 1 (n) xn = ann+1 (n) ann 1 (n ) (n x n1 = an11n+1 an11n) x n (n1) an1n1 . . . xk = . . . x1 1 (k) ak k n (k) a k n+1 (k) ak j x j j=k+1 1 = a (1) a11 1 n+1 . . . 12 n ( a11j) x j i=2 1.1.2 Pivoting Methods As we saw in Examples 1.1 and 1.2, naive Gaussian elimination may not produce accurate results when one entry is small relative to other entries in the system. Pivoting methods address this issue. We will cover two such methods here, Partial Pivoting, and Scaled Partial Pivoting. 1.1.2.1 Partial Pivoting Partial pivoting is used to avoid small pivots. The method adds the following step to the (k) naive Gaussian elimination algorithm: To nd the pivot a k k nd the rst row p such that k p n such that a pk = max {|aik |}. Exchange row k with row p before performing any k in eliminations. With this method the row-reduction of system (1.1) would proceed as follows x1 + x2 = 1 x1 + x2 = 2 1 1 1 1 2 1 1 2 1 1 Since |a21 | = max {|a11 | , |a21 |} = 1 we exchanged rows 1 and 2. The solution then follows the steps of the rst solution in Example 1.1. 1.1.2.2 Scaled Partial Pivoting Scaled partial pivoting is used when partial pivoting would not chose the largest pivot relative to the entries in its row. The method builds on partial pivoting by placing the largest entry, scaled relative to the entries in its row, in the pivot position. The specic steps are as follows, rst nd the largest element in each row without regard to sign si = max j[1,n] ai j Next to nd the rst pivot march down the rst column and nd the least integer p for which a p1 |ai1 | = max sp si i[1,n] Exchange row 1 with row p (the scale factors follow their rows on the interchanges). Continue with columns 2 through n 1. This method would prevent loss of precision problems with the system 10x 1 + 1018 x 2 = 1018 x1 + x2 = 2 |a11 | |a21 | = 1017 , and = 1, s1 s2 the pivot is a21 , and we exchange rows 1 and 2. Note the scalings are calculated only once, and are interchanged with their rows. which is modied from Example 1.1. Since |s1 | = 1018 , |s2 | = 1, 13 1.1.3 Summary of Gaussian Elimination Algorithms Naive Gaussian Elimination is summarized in the following algorithm: Input: An n (n + 1) augmented matrix, A b . Output: An approximate solution vector x = [x 1 x 2 . . . x n ]T to the system Ax = b or a message: \"No unique solution\". Step 1: (Elimination Step) For k = 1, . . . . , n 1, do steps 2 through 6. Step 2: Find the \"naive\" pivot row p, such that p is the least integer such that k p n and a pk 0. If no such p exists STOP Return message: \"No unique solution\" Step 3: (Row exchange step) If p k interchange rows p and k. Steps 46: Use pivot row to \"zero\" out the k th column. (n Step 7: If ann 1) = 0, STOP Return message: \"No unique solution\" (n) (n Step 8: Set ann = ann 1) . Step 9: (Back substitution step) Set x n = i = n 1, n 2, . . . , 1. (n) ann+1 (n) ann and x i = (i) ain+1 n 1 ai(i) x j j=i+ j (i) aii , for Step 10: Output x = [x 1, x 2, . . . . , x n ]T . Done Partial Pivoting replaces Step 2 in the naive Gaussian algorithm with the following: Step 2: Find the \"partial pivot row\