VIEWS: 104 PAGES: 10 CATEGORY: Legal POSTED ON: 1/14/2010
CS3911 Intro. to Numerical Methods with Fortran Final Solutions 1 CS3911 Introduction to Numerical Methods with Fortran Final Exam Solutions 1. Accuracy and Reliability (a) [20 points] Let be a very small positive number (i.e., ≈ 0). Do the following two problems: (1) solve the following system of linear equations with and without partial pivoting, assuming you use ﬁnite precision arithmetics; and (2) ﬁnd the major reason or reasons that can explain the diﬀerence(s) between the two solutions. You have to clearly state your ﬁndings with a convincing argument. Just stating a “reason” such as “it is because of cancelation” or “overﬂow” will receive zero point. 1 x 1 · = 1 1 y 2 Answer: Without pivoting, one multiplies −1/ to the ﬁrst equation and adds the result to the second. This yields the following: 1 x 1 · = 0 1 − 1/ y 2 − 1/ Since is small, 1/ is large. As a result, -1/ is the dominating term in both 1 − 1/ and 2 − 1/ . In other words, 1 − 1/ ≈ −1/ and 2 − 1/ ≈ −1/ . Hence, the equation (1 − 1/ )y = 2 − 1/ numerically becomes 1 1 − y≈− Therefore, we have y = 1. Plugging y = 1 into the ﬁrst equation x + y = 1 yields x = 0. This is certainly incorrect because x = 0 and y = 1 do not satisfy the second equation x + y = 2. With pivoting, swapping the equations transforms the system to the following: 1 1 x 2 · = 1 y 1 Multiplying the ﬁrst equation by − and adding the result to the second yields: 1 1 x 2 · = 0 1− y 1−2 Since is small, it does not contribute much to 1 − and 1 − 2 . As a result, 1 − ≈ 1 and 1 − 2 ≈ 1, and the second equation (1 − )y = 1 − 2 numerically becomes y = 1. Plugging y = 1 into the ﬁrst equation gives x = 1. Obviously, x = y = 1 is the correct numerical solution if is very small. The reason for the non-pivoting solution to go wrong is rounding error in computing 1 − 1/ and 2 − 1/ if 1/ is very large. In this case, rounding error makes both terms nearly equal to −1/ . CS3911 Intro. to Numerical Methods with Fortran Final Solutions 2 Consequently, we have y = 1 which is still correct. But, in solving x from x + y = 1 with y = 1, we have x = 0/ = 0. It does not require an extremely small to cause this problem. For example, on a 7-signiﬁcant-digit computer, if = 0.00000001 we have 1/ = 100, 000, 000. Then, 1 − 1/ = −99, 999, 999 and 2 − 1/ = −99, 999, 998. Both would be rounded to 7-digit and the result is 0.1 × 109. This example clearly illustrates that pivoting is necessary. 2. Systems of Linear Equations (a) [10 points] Answer the following two questions: (1) deﬁne the meaning of a diagonal dominant matrix A = [ai,j ]n×n ; and (2) what is the advantage(s) for a matrix to be diagonal dominant? Answer: A n × n matrix A = [ai,j ] is diagonal dominant if for every row i the following condition holds: n |ai,j | < |ai,i | j=1,j=i In other words, for each row, the sum of the absolute values of all oﬀ-diagonal entries is less than the absolute value of the diagonal element. If a matrix is diagonal dominant, Gauss elimination, LU-decomposition included, does not need pivoting, and Jacobi’s and Gauss-Seidel’s methods always converge. 3. Interpolation and Approximation (a) [10 points] Given a polynomial of degree n as follows, Pn (x) = a0 + a1 x + a2x2 + a3x3 + · · · + an xn suggest an eﬃcient way to evaluate Pn (x) with n multiplications. You should provide an algorithm and discuss its correctness and complexity. A method that does not achieve n multiplications receives zero point. Answer: The given polynomial can be rewritten in a nested form as follows: Pn (x) = a0 + (a1 + (a2 + (a3 + (· · · (an−1 + an x)x) · · ·)x)x) For example, P4 (x) = a0 + a1 x + a2 x2 + a3 x3 + a4 x4 = a0 + (a1 + (a2 + (a3 + a4x)x)x)x. In this way, only n multiplications are needed (i.e., O(n)). On the other hand, since ai xi requires i multiplications (i.e., i − 1 for xi and 1 for ai (xi)), the total number of multiplications with a direct evaluation of the power form is 0+1+2+· · ·+n = n(n+1)/2 (i.e., O(n2)). The following is a possible algorithm, which assumes that coeﬃcients ai ’s are in array a( ), x is in variable x, and Px has the result: Px = a(n) DO i = n-1, 0, -1 Px = a(i) + Px*x END DO Since this DO loop iterates n times, each of which uses exactly one multiplication, the order of complexity is O(n). One may suggest the following O(n) implementation: CS3911 Intro. to Numerical Methods with Fortran Final Solutions 3 Px = a(0) Power = x DO i = 1, n Px = Px + a(i)*Power Power = Power*x END DO Although it is an O(n) method, it uses two multiplications per iteration and the total number of multiplications is 2n. This is certainly slower than the nested form. (b) [20 points] Suppose ﬁve data points (x0 , y0) = (−2, 1), (x1, y1 ) = (−1, 1), (x2 , y2) = (0, 1), (x3, y3) = (1, 1) and (x4, y4) = (2, 25) are given. Find the Newton interpolating polynomial for the given data points with the divided diﬀerence method. You should show all computation steps. Only providing an answer and/or using a wrong method receives zero point. Answer: The Newton interpolating polynomial for ﬁve data points is P4 (x) = f [x0] + f [x0 , x1](x − x0) + f [x0 , x1, x2](x − x0)(x − x1 ) + f [x0, x1, x2, x3](x − x0)(x − x1 )(x − x2 ) + f [x0, x1, x2, x3, x4](x − x0)(x − x1 )(x − x2)(x − x3 ) The divided diﬀerence table is shown below: xi f[xi ] −2 1 −1 1 → 0 0 1 → 0 → 0 1 1 → 0 → 0 → 0 2 25 → 24 → 12 → 4 → 1 Therefore, we have f [x0 ] = 1, f [x0, x1] = 0, f [x0 , x1, x2] = 0, f [x0, x1, x2, x3] = 0 and f [x0, x1, x2, x3, x4] = 1, and the interpolating polynomial is P4 (x) = 1 + (x + 2)(x + 1)x(x − 1). 4. Numerical Diﬀerentiation: (a) [10 points] Use the 3-point backward diﬀerence method to compute the derivative of f (x) = ex sin(x) at x = 1 with ∆ = 0.01. You should show all computation steps. Only providing an answer and/or using a wrong method receives zero point. Answer: The 3-point backward diﬀerence method uses the following formula: 1 f (x) ≈ (3fi − 4fi−1 + fi−2 ) 2∆ Since ∆ = 0.01, fi = f (1) = 2.287355, fi−1 = f (1 − 0.01) = 2.249942 and fi−2 = f (1 − 0.02) = 2.212824, we have 1 f (1) ≈ (3 × 2.287355 − 4 × 2.249942 + 2.212824) = 3.756102 2 × 0.01 CS3911 Intro. to Numerical Methods with Fortran Final Solutions 4 (b) [20 points] Use Richardson’s extrapolation method to compute the derivative of f (x) = ex sin(x) at x = 1 with initial ∆ = 0.1. You should carry out the extrapolation steps until the computation for ∆ = 0.025 completes, and clearly indicate the desired answer. You should show all computation steps. Only providing an answer and/or using a wrong method receives zero point. Answer: Column 0 (i.e., d0,0, d1,0, d2,0, etc) uses 2-point central diﬀerence. Since the starting ∆ is 0.1 (i.e., row 0), the ∆’s for row 1 and row 2 are 0.05 = 0.1/2 and 0.025 = 0.5/2, respectively. The following shows d0,0, d1,0 and d2,0 using central diﬀerence with ∆ being 0.1, 0.05 and 0.025: f (1 + 0.1) − f (1 − 0.1) d0,0 = = 3.753307857 2 × 0.1 f (1 + 0.05) − f (1 − 0.05) d1,0 = = 3.755366227 2 × 0.05 f (1 + 0.025) − f (1 − 0.025) d2,0 = = 3.75587862 2 × 0.025 Richardson’s extrapolation uses the following formula: di,j − di−1,j di,j+1 = di,j + 4j+1 − 1 The following computes the desired result in d2,2: Row ∆ Column 0 Column 1 Column 2 0 0.1 d0,0 = 3.753307857 1 0.05 d1,0 = 3.755366227 → d1,1 = 3.7560523567 2 0.025 d2,0 = 3.75587862 → d2.1 = 3.7560494227 → d2,2 = 3.756049227 Therefore, the extrapolated derivative approximation is d2,2 = 3.756049227. 5. Numerical Integration: The following integral has an exact solution: 1 ex dx = e1 − e0 = e − 1 = 1.71828 . . . 0 (a) [20 points] Use the iterative 3-point Simpson method to compute the above integration. You should start with ∆ = 0.5, update the result in each iteration, and stop after the result of ∆ = 0.125 is obtained. You should show all computation steps. Only providing an answer and/or using a wrong method receives zero point. Answer: The 3-point Simpson method uses the following formula: m m ∆ (f0 + f2m ) + 4 f2i−1 + 2 f2i 3 i=1 i=1 The table below shows the computed results. In the table, “Previous Sum” is the sum of terms computed in the previous iteration, which becomes the sum of “even” terms in this iteration (i.e., m f2i ), and “This Odd Sum” is the sum of newly introduced terms i=1 in this iteration (i.e., m f2i−1 ). The sum of “Previous Sum” and “This Odd Sum” i=1 CS3911 Intro. to Numerical Methods with Fortran Final Solutions 5 becomes the “Previous Sum” for the next iteration. The desired result is 1.718284155. ∆ 0.5 0.25 0.125 f (0) + f (1) 3.718281828 3.718281828 3.718281828 0 1 0.125 1.133148453 0.25 1.284025417 0.375 1.454991415 0.5 1.648721271 0.625 1.868245957 0.75 2.117000017 0.825 2.398875294 1 2.718281828 Previous Sum 0 1.648721271 5.649746705 This Odd Sum 1.648721271 3.401025434 6.855261119 Integration 1.718861152 1.718318842 1.718284155 (b) [20 points] Use Romberg’s method to compute the above integration and use the trape- zoid method for column 0. You should carry out all steps until the computation for ∆ = 0.25 completes, and clearly indicate the desired answer. You should show all computation steps. Only providing an answer and/or using a wrong method receives zero point. Answer: Column 0 (i.e., r0,0, r1,0, r2,0, etc) uses the trapezoid method with 1, 2 and 4 subintervals. Since the starting interval is [0, 1], the length of subintervals for rows 0, 1 and 2 are 1, 0.5, and 0.25, respectively. The following shows r0,0, r1,0 and r2,0 using the trapezoid method with ∆ = 0.1, 0.05 and 0.025: f (0) + f (1) r0,0 = ∆ = 1.859149914 2 f (0) + f (1) r1,0 = ∆ + f (0.5) = 1.753931092 2 f (0) + f (1) r2,0 = ∆ + f (0.25) + f (0.5) + f (0.75) = 1.727221905 2 Romberg’s method uses the formula identical to that of Richardson’s extrapolation: ri,j − ri−1,j ri,j+1 = ri,j + 4j+1 − 1 The following computes the desired result in r2,2: Row ∆ Column 0 Column 1 Column 2 0 1 r0,0 = 1.859149914 1 0.5 r1,0 = 1.753931092 → r1,1 = 1.718861150 2 0.25 r2,0 = 1.727221905 → r2.1 = 1.718318843 → r2,2 = 1.718282689 Therefore, Romberg’s method yields r2,2 = 1.718282689. CS3911 Intro. to Numerical Methods with Fortran Final Solutions 6 6. ODE Initial Value Problems Suppose you are given an ODE y = xy with initial value y(0) = 1. (a) [10 points] Use Euler’s method to solve the given equation on [0, 4] with ∆ = 1, and ﬁll your results in the following table. Answer: Euler’s method uses y(x + ∆) = y(x) + ∆f (x, y) to compute the next y. The following table shows the computed results: Iteration x y Initial Value 0 1 1 1 1 2 2 2 3 3 6 4 4 24 Since f (x, y) = xy and ∆ = 1, Euler’s method uses y(x + 1) = y + xy, where y(0) = 1. The following shows the needed iterative steps. • Iteration 1: We have x = 0 and y = 1 from initial value. Plugging this into the above equation yields y(1) = 1 + 0 × 1 = 1. Thus, we have x = 1 and y = 1. • Iteration 2: Since we have x = 1 and y = 1, plugging them into the above equation yields y(2) = 1 + 1 × 2 = 2. Thus, we have x = 2 and y = 2. • Iteration 3: Plugging x = 2 and y = 2 into the above equation yields y(3) = 2 + 2 × 4 = 6. Hence, we have x = 3 and y = 6. • Iteration 4: Plugging x = 3 and y = 6 into the above equation yields y(4) = 6 + 3 × 6 = 24. Therefore, we have x = 4 and y = 24. (b) [15 points] Solve the given equation on [0, 3] using Ralston’s second-order Runge-Kutta method with ∆ = 1, and ﬁll your results in the following table. Answer: The second order Runge-Kutta method uses the following formula: y(x + ∆) = y(x) + ∆ [c1f (x, y) + c2f (x + p2∆, y + (a2,1f (x, y))∆)] Ralston’s method is a special case in which c1 = 1 , c2 = 2 and p2 = a2,1 = 3 . Note that 3 3 4 ∆ = 1. The following table shows the computed results: Iteration x y Initial Value 0 1 3 1 1 2 = 1.5 81 2 2 16 = 5.0625 2025 3 3 64 = 31.640625 CS3911 Intro. to Numerical Methods with Fortran Final Solutions 7 If you wish to do hand calculation, Ralston’s formula can be simpliﬁed to reduce com- putation eﬀort. Since f (x, y) = xy and ∆ = 1, Ralston’s formula can be rewritten as follows: 1 2 3 3 y(x + 1) = y + f (x, y) + f x + , y + f (x, y) 3 3 4 4 The f (x, y) in 1 f (x, y) and 3 f (x, y) should be replaced by xy. This yields the following: 3 4 1 2 3 3 y(x + 1) = y + xy + x + , y + xy 3 3 4 4 Since f (x, y) = xy (i.e., f (x, y) being the product of its two arguments), we have 3 3 3 3 f x + , y + xy = x + y + xy 4 4 4 4 As a result, Ralston’s formula becomes 1 2 3 3 y(x + 1) = y + xy + x+ y + xy 3 3 4 4 If you wish, the above can be simpliﬁed further, making calculations nearly eﬀortless. Expanding the product yields: 2 3 3 25 1 1 x+ y + xy = xy + x2 y + y 3 4 4 24 2 2 Putting this back to the simpliﬁed Ralston’s formula yields 1 2 3 3 y(x + 1) = y + xy + x + , y + xy 3 3 4 4 1 25 1 1 = y + xy + xy + x2 y + y 3 24 2 2 3 11 1 2 = y + xy + x y 2 8 2 1 11 = y 3 + x + x2 2 4 With this simpliﬁed form, calculations can be carried out easily. • Iteration 1: Initially we have x = 0 and y = 1. Plugging x and y into the above formula yields y(1) = 3 (i.e., x = 1 and y = 3 ). 2 2 • Iteration 2: When iteration 2 starts, we have x = 1 and y = 3 . Plugging x = 1 2 and y = 3 into Ralston’s formula yields y(2) = 81 (i.e., x = 2 and y = 81 ). 2 16 16 • Iteration 3: When iteration 3 starts, we have x = 2 and y = 81 . Plugging them 16 into Ralston’s formula gives y(3) = 2025 (i.e., x = 3 and y = 2025 ). This ﬁnishes 64 64 the third iteration. 7. Random Numbers (a) [10 points] Suppose we randomly draw a rectangle with side lengths a and b, both in the range of [0, 1). What is the probability of having a rectangle whose area is larger than of equal to 0.5? The following shows a sequence of pre-generated random numbers. Each number has ﬁve digits, and a decimal point is assumed to appear to the left of the ﬁrst digit. For example, random number 53479 should be interpreted as 0.53479. Therefore, all random numbers CS3911 Intro. to Numerical Methods with Fortran Final Solutions 8 are in the range of [0,1). Use Monte Carlo method with 12 runs to ﬁnd an approximation of this probability. You should show all computation steps. Only providing an answer and/or using a wrong method receives zero point. x1 to x8 53479 81115 98036 12217 59526 40238 40577 39351 x9 to x16 43211 69255 97344 70328 58116 91964 26240 44643 x17 to x24 83287 97391 92823 77578 66023 38277 74523 71118 Answer: Since the lengths are random numbers in [0, 1), the given table of random numbers can be used to estimate a probability. Rectangle Number Side a Side b Area = a × b Area ≥ 0.5? 1 0.53479 0.81115 0.43379 No 2 0.98036 0.12217 0.11977 No 3 0.59526 0.40238 0.23952 No 4 0.40577 0.39351 0.15967 No 5 0.43211 0.69255 0.29926 No 6 0.97344 0.70328 0.68460 Yes 7 0.58116 0.91964 0.53446 Yes 8 0.26240 0.44643 0.11714 No 9 0.83287 0.97391 0.81114 Yes 10 0.92823 0.77578 0.72010 Yes 11 0.66023 0.38277 0.25272 No 12 0.74523 0.71118 0.52999 Yes Of these 12 runs, ﬁve have an area larger than or equal to 0.5. Therefore, the estimated probability is 5/12 = 0.41666 . . . ≈ 0.42. In fact, this estimated probability is a bit too high. You may write a program to generate a large number of random rectangles and the result should be approximately 0.15. If you know more about probability theory, you should be able to determine this probability directly, because this is the probability of two uniform random numbers. We will not provide any discussion in this direction as this is not a topic in this course. 8. Eigenvalues and Eigenvectors (a) [15 points] Use the power method to ﬁnd the largest eigenvalue and its corresponding eigenvector of matrix A as shown below, and ﬁll the following table with your results. The initial value (i.e., iteration 0) is z = [1, 1]T , and you only do two iterations (i.e., iterations 2 and 3). 1 2 A= 4 3 Answer: This following table has the desired result: CS3911 Intro. to Numerical Methods with Fortran Final Solutions 9 Approx. Approx. Eigenvector Iteration Eigenvalue x y 0 1 1 1 3 1 7 7 = 0.428571428 1 33 17 2 7 = 4.714285714 33 = 0.515151515 1 Here is the details: • Iteration 1: Since the initial value is z = [1, 1]T , we have w = A · z = [3, 7]T . The maximum in w is 7, which is used to scale w. Thus, λ = 7 and z = w/7 = [ 3 , 1]T . 7 This completes iteration 1. • Iteration 2: Since z = [ 3 , 1]T , we have w = A · z = [ 17 , 33 ]T . The maximum in w 7 7 7 is 33 . We have λ = 33 and z = w/ 33 = [ 17 , 1]T . This completes iteration 2. 7 7 7 33 (b) [20 points] Use Jacobi method to ﬁnd all eigenvalues and their corresponding eigen- vectors of the following symmetric matrix A: √ √ 3 2 2 2 4 − 4 √ A= 4 2 9 4 3 4 √ 2 3 9 − 4 4 4 You should clearly provide all computation details, and match each eigen- value with its corresponding eigenvector. Otherwise, you will risk low grade. Additionally, you will receive zero point if you do not use Jacobi method. One more note: Unless you are very good at hand calculation, use your cal- culator to do this problem. Answer: Since the maximum of oﬀ-diagonal entries is a2,3 = 3/4, we have 2a2,3 2× 3 tan(2θ) = = 9 4 =∞ 9 a3,3 − a2,2 4 − 4 and θ = π/4. The rotation matrix of iteration 1 is 1 0 0 1 0 0 √ √ R1 = 0 cos(π/4) sin(π/4) = 0 2 2 2 2 √ √ 0 − sin(π/4) cos(π/4) 2 2 0 − 2 2 The rotated A is computed as follows: 3 1 2 2 0 A = RT 1 · A · R1 = 1 2 3 2 0 0 0 3 CS3911 Intro. to Numerical Methods with Fortran Final Solutions 10 Since A is not diagonal, it should be rotated again. In the second iteration, since the maximum oﬀ-diagonal entry is a1,2 = 1 , we have 2 2a1,2 2× 1 tan(2θ) = = 1 2 =∞ 1 a2,2 − a1,1 2 − 2 Again, θ = π/4 and the second rotation matrix is √ √ 2 2 cos(π/4) sin(π/4) 0 2 2 0 √ √ R2 = − sin(π/4) cos(π/4) 0 = − 2 2 0 2 2 0 0 1 0 0 1 The rotated A is computed as follows: 1 0 0 A = R T · A · R2 = 0 2 2 0 0 0 3 Since A is diagonal, eigenvalues of matrix A are 1, 2 and 3. The eigenvectors of A is computed as √ √ 2 2 1 0 0 cos(π/4) sin(π/4) 0 2 2 0 √ R1·R2 = 0 cos(π/4) sin(π/4) · − sin(π/4) cos(π/4) 0 = −1 1 2 2 2 2 √ 0 − sin(π/4) cos(π/4) 0 0 1 1 −1 2 2 2 2 The following table shows eigenvalues and their corresponding eigenvectors: Eigenvalues 1 2 3 √ √ 2 2 0 2 2 √ 1 Eigenvectors −2 1 2 2 2 √ 1 2 −1 2 2 2 More precisely, the corresponding eigenvector of the i-th eigenvalue is the i-th column of R1 · R2 .