Homework 10 Solution

Homework 10 Math 471, Fall 2007 Assigned: Friday, November 16, 2007 Due: Monday, December 3, 2007 Include a cover page Clearly label all plots using title, xlabel, ylabel, legend Use the subplot command to compare multiple plots Include printouts of all Matlab code, labeled with your name, date, section, etc. (1) (Piecewise Linear Interpolation) Use the error theorem on p. 390 to compute an upper bound on the pointwise error from interpolating f(x) = e x2 on the interval [ 4;4] with piecewise linear functions at the points 4, 2, 0, 2, and 4. maxx2[ 4;4]je x2 s(x)j 18h2 maxx2[ 4;4]jf00(x)j= 1: (2) (Cubic Spline Interpolation) P. 403 #10, P. 404 #15 P. 403 #10a j aj bj cj dj 0 0.500 1.49 0.286 -0.153 1 1.43 2.18 0.579 -0.153 2 2.64 2.64 0.350 -0.330 3 4.01 2.75 -0.144 -0.330 P. 403 #10b j aj bj cj dj 0 0.500 1.50 0.756 -0.106 1 1.43 2.18 0.600 -0.175 2 2.64 2.64 0.334 -0.287 3 4.01 2.76 -0.0961 -0.478 P. 404 #15 j aj bj cj dj 0 0.500 1.72 0 0.523 1 1.43 2.11 0.784 -0.297 2 2.64 2.67 0.338 -0.426 3 4.01 2.69 -0.301 0.200 (3) (Richardson Extrapolation) P. 454 #11 (a) Adding the Taylor series for f(x0 +h) and f(x0 h) gives f(x0 +h) +f(x0 h) = 2f(x0) +h2f00(x0) + h 4 12f (4)(x0) + h8 2880f (8)(x0) + h8 2880[f (8)( ) f(8)(x0)]: Thus 1 2 f00(x0) = f(x0 +h) 2f(x0) +f(x0 h)h2 h 2 12f (4)(x0) h4 360f (6)(x0) h6 2880f (8)(x0) +o(h6): (b) h O(h2) O(h4) O(h6) 0.5 2.255252 0.25 2.062826 2.049998 0.125 2.015645 2.012500 2.010000 (4) (Newton{Cotes). Suppose that f is a function with four continuous derivatives on the interval [a;b]. Recall that the error bound for the composite trapezoidal rule T(h) with panel width h is T(h) Rba f(x) dx = (b a)h212 f00( ) for some 2[a;b]. The error in the composite Simpson?s rule S(h) with panel width h is S(h) Rba f(x) dx = (b a)h4180 f(4)( ) for some (di erent) 2[a;b]. (a) Let f(x) = e x sinx. For the composite trapezoidal rule and the composite Simpson?s rule, nd the number of panels n required to integrate f on the interval [0;2 ] with error at most 10 4. Recall that h = 2 =n. How many function evaluations are required in each case? The rst ve derivatives of f are f0(x) = e x(cosx sinx) f00(x) = 2e x cosx f000(x) = 2e x(sinx+ cosx) f(4)(x) = 4e x sinx f(5)(x) = 4e x(sinx cosx) To develop the bound for the composite trapezoidal rule, we need to calculate the maximum value of jf00(x)j on the interval [0;2 ]. The extremum must occur at x = 0, x = 2 , or at a point where f000(x) = 0. A short calculation establishes that the maximum occurs when x = 0, so maxx2[0;2 ]jf00(x)j= 2: To ensure that the error is no more than 10 4, the error bound requires 2 h2 12 maxx2[0;2 ]jf 00(x)j 10 4: Solving this inequality for the panel width h yields h q 3 10 4 : Since h = 2 =n, we conclude that n 2 p 3 10 4; or n 643. To develop the bound for Simpson?s rule, we need to calculate the maximum value of jf(4)(x)j on the interval [0;2 ]. The extremum must occur at x = 0, x = 2 , or at a 3 point where f(5)(x) = 0. A short calculation establishes that the last case is in force, and the extremal value of f(4)(x) occurs at x = =4. Thus max x2[0;2 ] jf(4)(x)j= 2p2 e =4: To ensure that the error does not exceed 10 4, the error bound requires 2 h4 180 maxx2[0;2 ]jf (4)(x)j 10 4: Solving this inequality for the panel width h yields h h 45p2 10 4 2 e =4 i1=4 : Since h = 2 =n, we determine that n 2 " 2 e =4 45p2 10 4 #1=4 ; or n 29. Simpson?s rule requires an even number of panels, so we must increase n to 30. Nevertheless, by using Simpson?s rule, we obtain a drastic reduction in the number of panels required. (b) Using the number of panels determined in the last part of the problem, use each rule to approximate the integral numerically. Compare your results with the actual value of the integral. For your information, an antiderivative of f isR e x sinxdx = 0:5 e x(sinx+ cosx): Using the antiderivative, which normally would not be available, we determine that the actual value of the integral is Z 2 0 e x sinxdx = 0:5 e x(sinx+ cosx) 2 x=0 = 0:5 0:5 e 2 = 0:49906628: First, we apply the composite trapezoidal rule with n = 643 panels, which corresponds with a panel width h = 2 =643. So T(h) = h2 " f(0) + 2 642X i=1 f(ih) +f(2 ) # = 0:49905834: The actual error is T(h) R2 0 f(x) dx = 7:94 10 6; which is about one order of magnitude better than we required. The number of function evaluations involved is 644. The composite Simpson?s rule requires an even number of panels. The number n = 30 is acceptable, which corresponds with a panel width of h = 2 =30. So S(h) = h3 " f(0) + 4 15X i=1 f((2i 1)h) + 2 14X i=1 f(2ih) +f(2 ) # = 0:49904472: 4 The actual error is S(h) R2 0 f(x) dx = 2:16 10 5; which is about an order of magnitude better than we required. The number of function evaluations involved is 31, which is far superior to the composite trapezoidal rule. (5) (Gaussian Quadrature). The Gaussian quadrature rule for n = 6 has the following nodes and weights: n xn wn 1 0:93246951 0:17132449 2 0:66120939 0:36076157 3 0:23861918 0:46791393 4 0:23861918 0:46791393 5 0:66120939 0:36076157 6 0:93246951 0:17132449 (a) Verify that the Gaussian quadrature rule with six nodes correctly integrates the poly- nomials 1;x;x2;:::;x11 on the interval [ 1;1] with error less than 10 7 (due to nite precision issues). Verify that the rule fails for x12. Note that the Gaussian quadrature rule succeeds for any odd-degree monomial because the nodes are symmetric about the origin. Therefore, we need only check the even- degree monomials. Function Estimated Actual Error 1 1.99999998 2.00000000 2:0 10 8 x2 0.66666666 0.66666667 9:4 10 9 x4 0.39999999 0.40000000 6:8 10 9 x6 0.28571428 0.28571429 7:8 10 9 x8 0.22222221 0.22222222 8:9 10 9 x10 0.18181817 0.18181818 9:5 10 9 x12 0.15310807 0.15384615 7:4 10 4 Although the failure is not dramatic, it is clear that the quadrature rule is no longer yielding precise estimates when the degree of the monomial increases to 12. (b) Apply the Gaussian quadrature rule with six nodes to approximateR 2 0 e x sinxdx: Make sure that you rst change the domain of integration to [ 1;1]. What is the actual error? How many function evaluations are performed? First, we change variables in the integral:R 2 0 e x sinxdx = R1 1 e (x+1) sin( (x+ 1)) dx: We apply the Gaussian quadrature rule to the integrand g(x) = e (x+1) sin( (x+1)) to obtain G6(g) = 0:49907053956816: The actual error is G6(g) R1 1g(x) dx = 4:26 10 6: The number of function evaluations is just six. Compare this with the composite trapezoidal rule and the composite Simpson?s rule. 5 (6) (Romberg Integration) (a) Following the directions in Section 6.7 of the book, write a Matlab program to perform Romberg integration of an arbitrary input function f on the interval [a;b]. It should halt when the estimated error declines below a speci ed threshold epsilon. It should report the estimated value of the integral and the number of function evaluations performed. The rst line will read function [value, fnevals] = romberg( f, a, b, epsilon ), You can de ne a function f to pass to your code using the inline command. As an example, f(x) = e x sinx is de ned via f = inline( ?exp(-x) * sin(x)? ); % [value, fnevals] = romberg( f, a, b, epsilon ) % % Uses Romberg integration to evaluate the definite % integral of the function f on the interval [a, b] % It halts when the estimated error declines below % epsilon % % Returns value, the estimated integral % and fnevals, the number of function evaluations function [value, fnevals] = romberg( f, a, b, epsilon ), if (a >= b), error( ?The endpoints do not satisfy a < b.? ); end % Initializations errorEst = Inf; % Trivial initial error estimate row = 1; % Row number n = 1; % Number of panels h = abs(b - a); % The panel width % Calculate first entry of the first row table( row, 1 ) = h/2 * ( f(a) + f(b) ); % Calculate successive rows % Stop when the error estimate drops below epsilon % We also require that the panels be somewhat narrow % before we stop. This is because the extrapolation formulae % and the error estimates only work for "small" panel width 6 while ( or( h >= 0.5, errorEst >= epsilon) ), row = row + 1; % Increment the row number n = 2*n; % Double the number of panels h = h/2; % Halve the panel width % Calculate the first entry in the row. % Note that we reuse the prior estimate % to save function evaluations update = 0; for j = 1 : 2 : (n - 1), update = update + f( a + j*h ); end table( row, 1 ) = table( row-1, 1 ) / 2 + h * update; % Extrapolate for k = 2 : row, table( row, k ) = ... (4^(k-1) * table( row, k-1 ) - table( row-1, k-1 ) ) ... / (4^(k-1) - 1); end % Estimate error errorEst = abs(table(row, row) - table(row-1, row-1)) / 2^(row-1); end % Uncomment the next line to have a look at the table of estimates % table value = table( row, row ); % Estimated integral fnevals = n + 1; % Number of evaluations (b) Use your code to estimate the value ofR 2 0 e x sinxdx with an error tolerance of 10 4. How many function evaluations are performed? What is the actual error? How does this compare with the composite trapezoidal rule? >> f = inline( ?exp(-x)*sin(x)? ); >> [val, fnevals] = romberg( f, 0, 2*pi, 1e-4 ) 7 table = -0.00000000 0.00000000 0.00000000 0.31242555 0.41656740 0.44433856 0.44884307 0.49431557 0.49949878 0.50037435 0.48630564 0.49879317 0.49909167 0.49908522 0.49908016 val = 0.49908016 fnevals = 17 The error using Romberg integration is 1:38 10 5, which again is about one order of magnitude better than we required. The number of function evaluations is 17, which easily beats the composite trapezoidal rule and the composite Simpson?s rule. The powerful performance is not surprising since these two rules generate the rst two columns of the table. Since the integrand is in nitely di erentiable, the Euler{ Maclaurin summation formula can be extended inde nitely, and repeated extrapolation will improve the estimates of the integral until oating-point errors begin to interfere. Note that the rst entries in the table are zero because the integrand is inconveniently zero at the endpoints and midpoint of the domain of integration. This observation illustrates that numerical quadrature algorithms can easily be fooled, so it is important to make sure that your answers are reasonable. (c) Use your code to solve Problems 15 and 17 on page 504. >> f = inline( ?1/sqrt(1+x^4)? ); >> [val, fnevals] = romberg( f, 0, 1, 5e-7 ) val = 0.92703736 fnevals = 17 >> f = inline( ?sin(x^2)? ); >> [val, fnevals] = romberg( f, 0, 1, 5e-7 ) val = 0.31026830 fnevals = 17