Computational methods and applications (AMS 147)
Homework 4 – Due Sunday February 24 Please submit to CANVAS a .zip file that includes the following Matlab functions:
poly least squares.m test least squares.m
For the Extra Credit 1 and 2, scan your notes into one PDF file scan.pdf, and attach it to your submission (as a separate file, not as part of the .zip file)
Exercise 1 (50 points) Write a Matlab function poly least squares.m that implements the least squares method to approximate a data set in terms of a polynomial model of degree M. The function should be of the form
function [a,err] = poly least squares(xi,yi,M)
Input:
xi: vector of nodes xi=[xi(1) … xi(N)]
yi: vector of data points yi=[yi(1) … yi(N)] corresponding to [xi(1) … xi(N)] M: degree of the polynomial model
ψ(x) =a(1)+a(2)x+a(3)x2 +···+a(M+1)xM (1) a: vector of coefficients representing the polynomial (1), i.e., a=[a(1) … a(M+1)]
err: error between the model and the data in the 2-norm
N
err=[yi −ψ(xi)]2. (2)
i=1
Hint: You can compare the results of your least squares implementation with the output of the Matlab functions polyfit() and polyval() (see the Matlab/Octave documentation for details).
Exercise 2 (50 points) Use the function you coded in Exercise 1 to determine the least squares polynomial ap- proximants of the attached data set DJI 2014 2019.dat representing the daily opening value of the Dow Jones Industrial Average from 1/1/14 to 2/17/19. To this end, write a Matlab/Octave function test least squares.m of the form
function [x,p1,p2,p4,p8] = test least squares()
Output:
x: vector of 1000 evenly-spaced evaluation nodes in [0, 1] (including endpoints)
p1, p2, p4, p8: least squares polynomial models (1) of the DJI value with degrees M = 1, 2, 4, 8, respectively,
Output:
1
evaluated at x.
The function should also return one figure with the plots of the DJI opening prices {x(i),y(i)}i=1,2,… (from the data file DJI 2014 2019.dat) (in blue) and the least-squares polynomial models p1, p2, p4 and p8 you computed above (in red).
Hint: To load the DJI data in Matlab/Octave use the command load (see the Matlab/Octave documentation for further details).
Extra Credit 1 (10 points) Show that you can integrate exactly any polynomial p3(x) of degree 3 in [0,1] by integrating the parabola Π2p3(x) that interpolates p3(x) at{0, 1/2, 1}, i.e.,
1 1
p3(x)dx = Π2p3(x)dx. (3)
00
Extra Credit 2 (10 points) Let xj = j/n (j = 0, …, n), be a set of n+1 evenly spaced points in [0, 1]. We have seen in class that we can approximate any smooth function f : [0, 1] → R with an interpolatory (not-a-knot) cubic spline s3(x). This means that we can approximate the integral of f(x) by integrating the corresponding spline s3(x) (which can be done analytically). The error associated with such approximation can bounded as
1
f(x)dx−
0
1 0
5 ′′′′ s3(x)dx≤ max |f
(x)|, (4)
384n4 x∈[0,1]
(see the lecture noteLECTURE8piecewiseinterpolationandsplines.pdfposted in CANVAS). Derive a similar error estimate for quadratic cubic splines and compare the upper bound you obtain with the one of the Simpson rule1. To this end, note that if we interpolate the function f(x) at the n+2 nodes x0 = 0, xn = 1 and xj = (xj+1 + xj )/2 (j = 0, .., n − 1), then the spline s2(x) is uniquely defined, and we have the error estimate
max |f(x) − s2(x)| ≤ 1 max |f′′′(x)| . (5) x∈[0,1] 7n3 x∈[0,1]
Which method between Simpson and quadratic spline interpolation converges faster to the integral of f(x) in [0, 1] as we increase the number of points n?
1Note that in both cases you approximate the function f(x) locally with a polynomial of degree 2. The difference between spline interpolation and simple piecewise polynomial interpolation (used in the Simpson rule) is that s2(x) is C1([0,1]), while the classical piecewise polynomial interpolant is only C0([0, 1]).
2