Unformatted Attachment Preview
Math 128A, Spring 2019
Problem Set 07
1 (a) Find an exact formula for the cubic polynomial P3 (x) = x3 + · · · such that
Z 1
P3 (x)q(x)dx = 0
−1
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
1
Z
q(x)dx =
−1
3
X
wj q(xj )
j=1
exactly whenever q is a polynomial of degree 5.
(d) Given any real numbers a < b, find exact formulas for points yj ∈ [a, b] and weights
uj > 0 such that
Z b
3
X
q(x)dx =
uj q(yj )
a
j=1
whenever q is a polynomial of degree 5.
(e) Explain why each of the three factors in the error estimate
Z
b
f (x)dx −
a
3
X
uj f (yj ) = C6 f
(6)
Z
b
(ξ)
(y − y1 )2 (y − y2 )2 (y − y3 )2 dy
a
j=1
is inevitable and determine the exact value of the constant C6 .
(f) Use your code ectr.m to evaluate
Z 1
E6 =
(x − x1 )2 (x − x2 )2 (x − x3 )2 dx
−1
to 3-digit accuracy. Use parameters r = [x1 , x2 , x3 ].
2 Implement, debug and test a MATLAB function pleg.m of the form
function p = pleg(t, n)
% t: evaluation point
% n: degree of polynomial
This function evaluates a single value Pn (t) of the monic Legendre polynomial Pn of degree
n, at evaluation point t with |t| ≤ 1. Here P0 = 1, P1 (t) = t and Pn is determined by the
recurrence
Pn (t) = tPn−1 (t) − cn Pn−2 (t)
for n ≥ 2, where cn = (n − 1)2 /(4(n − 1)2 − 1). Be sure to iterate forward from n = 0 rather
than recurse backward from n, and do not generate any new function handles. Test that
your function gives the right values for small n where you know Pn .
1
Math 128A, Spring 2019
Problem Set 07
3 Implement a MATLAB function gaussint.m of the form
function [w, t] = gaussint( n )
% n: Number of Gauss weights and points
which computes weights w and points t for the n-point Gaussian integration rule
Z
1
f (t)dt ≈
−1
n
X
wj f (tj ).
j=1
(a) Find the points tj to as high precision as possible, by applying your code bisection.m
to pleg.m. Bracket each tj initially by the observation that the zeroes of Pn−1 separate the
zeroes of Pn for every n. Thus the single zero of P1 = t 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 find all the zeroes of P1 , P2 , . . . ,
Pn−1 in the process of finding all the zeroes of Pn .
(b) Find the weights wj to as high precision as possible by applying your code ectr.m to
Z
1
Lj (t)2 dt
wj =
−1
where Lj is the jth Lagrange basis polynomial for interpolating at t1 , t2 , . . . , tn .
(c) For 1 ≤ n ≤ 20, test that your weights and points integrate monomials f (t) = tj
exactly for 0 ≤ j ≤ 2n − 1.
4 (a) Show that
Z
1
−x
x dx =
0
∞
X
n−n
n=1
(b) Use the sum in (a) to evaluate the integral in (a) to 12-digit accuracy.
(c) Evaluate the integral in (a) by ectr.m to 1, 2, and 3-digit accuracy. Estimate how
many function evaluations will be required to achieve p-digit accuracy for 1 ≤ p ≤ 12.
Explain the agreement or disagreement of your results with theory.
5 (a) Write, test and debug an adaptive 3-point Gaussian integration code gadap.m of the
form
function [int, abt] = gadap(a, b, f, r, tol)
% a,b: interval endpoints with a < b
% f: function handle f(x, r) to integrate
% r: parameters for f
% tol: User-provided tolerance for integral accuracy
% int: Approximation to the integral
% abt: Endpoints and approximations
2
Math 128A, Spring 2019
Problem Set 07
Build a list abt = {[a1 , b1 , t1 ], . . . , [an , bn , tn ]} of n intervals [aj , bj ] and approximate integrals
Rb
tj ≈ ajj f (x, r)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 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 infinite loops and floating-point issues as you see fit and briefly
justify your design decisions in comments.
R1
(b) Approximate the integral 0 x−x dx using your code from (a). Tabulate the total
number of function evaluations required to obtain p-digit accuracy for 1 ≤ p ≤ 12. Compare
your results with the results and estimates for endpoint-corrected trapezoidal integration
obtained in problem 4.
3