Math 128A, Spring 2016 Problem Set 07 Question 1 (a) Find an exact formula for the cubic polynomial P3 (x) = x3 + such that 1 1 P3 (x)q(x)dx = 0 for any quadratic polynomial q. (b) Find exact formulas for the three roots x1 , x2 , x3 of the equation P3 (x) = 0. (c) Find exact formulas for the integration weights w1 , w2 , w3 such that 3 1 q(x)dx = 1 wj q(xj ) j=1 exactly whenever q is a polynomial of degree 5. (d) Given any real numbers a < b, nd exact formulas for points yj [a, b] and weights uj > 0 such that 3 b uj q(yj ) q(x)dx = a j=1 whenever q is a polynomial of degree 5. (e) Explain why each of the three factors in the error estimate 3 b uj f (yj ) = C6 f (6) () f (x)dx a j=1 b a (y y1 )2 (y y2 )2 (y y3 )2 dy is inevitable and determine the exact value of the constant C6 . (f) Use your code from Question 2 to evaluate 1 E6 = 0 (x x1 )2 (x x2 )2 (x x3 )2 dx to 3-digit accuracy. Question 2 (a) Write, test and debug an adaptive 3-point Gaussian integration code gadap.m of the form function [int, abt] = gadap(a, b, f, p, tol) % a,b: interval endpoints with a < b % f: function handle f(x, p) to integrate (p for user parameters) % tol: User-provided tolerance for integral accuracy % int: Approximation to the integral % abt: Endpoints and approximations 1 Math 128A, Spring 2016 Problem Set 07 Build a list abt = {[a1 , b1 , t1 ], . . . , [an , bn , tn ]} of n intervals [aj , bj ] and apb proximate integrals tj ajj f (x)dx, computed with 3-point Gaussian integration. Initialize with n = 1 and [a1 , b1 ] = [a, b]. At each step j = 1, 2, . . ., subdivide interval j into into left and right half-intervals l and r, and approximate the integrals tl and tr over each half-interval by 3-point Gaussian quadrature. If |tj (tl + tr )| > tol max(|tj |, |tl | + |tr |) add the half-intervals l and r and approximations tl and tr to the list. Otherwise, increment int by tj . Guard against innite loops and oating-point issues as you see t and briey justify your design decisions in comments. (b) Approximate the integral 01 xx dx using your code from (a). Measure the total number of function evaluations required to obtain 12-digit accuracy. Plot the accepted intervals. Compare your results with those obtained in the previous problem set by Romberg integration. 2 Math 128A, Spring 2016 Programming Project 02 For the following two problems, write and debug MATLAB codes and make sure they run with the test autograders from the course web page. Test them thoroughly on test cases of your own design. When you are convinced they work, submit your codes together with brief discussion of any design decisions brief comparison of test results with theory Part 1 Implement two MATLAB functions cleg.m and pleg.m of the form function c = cleg(n) % n: number of recurrence parameters function [ p, pp, ppp ] = pleg(x, c) % x: vector of evaluation points % c: n-vector of recurrence parameters These functions work together to evaluate all Legendre polynomials Pn of degree up to and including n, at an array of evaluation points xj with |xj | 1. Here P0 = 1, P1 (x) = x and Pn is determined by the recurrence Pn (x) = xPn1 (x) cn Pn2 (x) for n 2. (As a consequence, Pn (x) = Pn1 (x) + xPn1 (x) cn Pn2 (x) and similarly for Pn (x).) Each parameter cn is determined by cn = 1 1 xPn1 (x)Pn2 (x)dx 1 1 Pn2 (x)2 dx for n 2, so for example c2 = 1/3. Your code cleg should evaluate the integrals to 12digit accuracy with your code gadap.m from Problem Set 07. After cleg( n ) generates the constants c1 , c2 , . . . , cn . your code pleg should evaluate Pn (xj ), Pn (xj ), Pn (xj ) by recurrence. pleg determines n by checking the length of the parameter vector c. Part 2 Implement a MATLAB function gaussint.m of the form function [w, x] = gaussint( n ) % n: Number of Gauss weights and points which computes weights w and points x for the n-point Gaussian integration rule n 1 f (x)dx 1 wj f (xj ). j=1 1 Math 128A, Spring 2016 Programming Project 02 (a) Find the points xj by your code schroderbisection from Programming Project 1, modied as necessary. Bracket each xj initially by the observation that the zeroes of Pn1 separate the zeroes of Pn for every n. Thus the single zero of P1 = x separates the interval [1, 1] into two intervals, each containing exactly one zero of P2 . The two zeroes of P2 separate the interval [1, 1] into three intervals, and so forth. Thus you will nd all the zeroes of P1 , P2 , . . . , Pn1 in the process of nding all the zeroes of Pn . (b) Find the weights wj by applying gadap.m to 1 Lj (x)2 dx wj = 1 where Lj is the jth Lagrange basis polynomial for interpolating at x1 , x2 , . . . , xn . Code Submission: Upload the MATLAB les cleg.m , pleg.m , and gaussint.m and any supporting les to bCourses for your GSI to grade. 2 function [r, h] = schroderbisection(a, b, f, fp, fpp, t) % Function that combines the fast convergence of Schroder iteration for % multiple roots with the bracketing guarantee of bisection % % a: beginning of interval [a, b] % b: end of interval [a, b] % f: function handle f(x) to find a zero for % fp: function hanfle f'(x) % fpp: function handle f''(x) % t: user provided tolerance for interval width % the left endpoint of our interval left = a; % the right endpoint of our interval right = b; % the size of the interval range = right - left; % the matrix within which we will store the various values we acquire % as we progress through the iterations h = zeros(6, 1); % a counter to keep track of the number of iterations count = 1; % the function e as defined in the specification e = @(x) min(abs(f(right) - f(left))/8, abs(fpp(x)) * abs(right - left)^2); % we create two separate iteration functions - one for g+ and another for % g- which are simply the regular iteration function except with f(x) % replaced with (f(x) + e(p)) in the g+ case and (f(x) - e(p)) in the % g- case. gplus = @(x) x - (((f(x) + e(x)) * fp(x)) / (fp(x)^2 - ((f(x) + e(x)) * fpp(x)))); gminus = @(x) x - (((f(x) - e(x)) * fp(x)) / (fp(x)^2 - ((f(x) - e(x)) * fpp(x)))); while right - left >= t && right - left >= eps(left) h(1, count) = left; h(5, count) = right; % we define p at each step of the iteration through geometric % bisection if left <= 0 && right >= 0 if sign(f(realmin)) ~= sign(f(right)) p = realmin; else p = -realmin; end else p = sqrt(abs(left * right)) * sign(left); end h(2, h(3, h(4, h(6, count) count) count) count) = = = = p; gminus(p); gplus(p); f(p); % the five values from which we will select our new interval. We then % sort it and add it into a new list. values = [left, p, gminus(p), gplus(p), right]; sorted = sort(values); % if f(v) = 0 for any of the five values above, then we have found % the root and set both the left and right endpoint to be this value for i = 1:5 if f(sorted(i)) == 0 left = sorted(i); right = sorted(i); break end end % if the left and right endpoints are identical, then we have found % the root and we can exit the loop if left == right break end % we now iterate through adjacent pairs of values and select the first % pair where we observe a sign change - this interval will therefore % contain the root and since this list of values is sorted, we know % that this is the smallest interval containing the interval. We choose % to sort the list and iterate through it this way instead of iterating % through every pair and picking the smallest interval because if any % of the values are realmin or smaller than realmin, MATLAB will not be % able to distinguish between the two values and we could end up with a % bug where the code cannot distinguish between the two intervals. for i = 1:4 if sign(f(sorted(i))) ~= sign(f(sorted(i + 1))) left = sorted(i); right = sorted(i + 1); break; end end count = count + 1; end % now that we have found the required interval, we store the endpoints of % the interval and the other required values in our matrix if left <= 0 && right >= 0 p = 0; else p = sqrt(abs(left * right)) * sign(left); end h(1, h(2, h(3, h(4, h(5, h(6, count) count) count) count) count) count) = = = = = = left; p; gminus(p); gplus(p); right; f(p); % we compute the approximate to the root, r, as the arithmetic mean of % our left and right endpoint r = (left + right)/2; end