Math 128A, Spring 2019 Problem Set 07

Math 128A, Spring 2019 Problem Set 07
1 (a) Find an exact formula for the cubic polynomial P3(x) = x
3 + · · · such that
Z 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
Z 1
−1
q(x)dx =
X
3
j=1
wjq(xj )
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
a
q(x)dx =
X
3
j=1
ujq(yj )
whenever q is a polynomial of degree 5.
(e) Explain why each of the three factors in the error estimate
Z b
a
f(x)dx −
X
3
j=1
ujf(yj ) = C6f
(6)(ξ)
Z 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 ectr.m to evaluate
E6 =
Z 1
−1
(x − x1)
2
(x − x2)
2
(x − x3)
2
dx
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) − cnPn−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
−1
f(t)dt ≈
Xn
j=1
wjf(tj ).
(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
wj =
Z 1
−1
Lj (t)
2
dt
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) = t
j
exactly for 0 ≤ j ≤ 2n − 1.
4 (a) Show that
Z 1
0
x
−x
dx =
X∞
n=1
n
−n
(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
tj ≈
R bj
aj
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.
(b) Approximate the integral R 1
0
x
−xdx 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