Assignment 6, Math 575A

Part I Matlab Section:
MATLAB has special functions to deal with polynomials. Using these commands is usually recommended, since they make the code easier to write and understand and are usually more efficient. In this HW assignment you should try to use MATLAB polynomial commands (and avoid for loops) as much as possible.

The polynomial $P(x) = 2x^2+2x-4$ and $Q(x) = x^2-6$ are represented in MATLAB by:

       P = [2 2 -4];
       Q = [1 0 -6];

$P(4.7)$ is evaluated, for example, using:

  polyval(P,4.7)

You can plot $Q(x)$ in the interval $[-6,   6]$ using:

  x = [-6:0.1:6];
  plot(x,polyval(Q,x))

The polynomial $S(x) = P(x)+Q(x)$ is calculated using

   S = P+Q;

(but addition or subtraction of polynomials with different degrees takes somewhat more effort). Multiplication is easy and the degrees do not have to be equal. The multiplication $T(x) = P(x)*Q(x)$ is represented by

   T = conv(P,Q);

Additional useful commands are prod, roots (finds the roots of a polynomial) and poly (constructs a polynomial with specified roots). For vectors, roots and poly are inverse functions of each other, up to ordering, scaling, and roundoff error. For more info look at the table of the MATLAB Primer and use the online help or ask for technical assistance.

In this assignment you are asked to hand in many plots. You can save paper by combining several plots on one page using the subplot command. Use the help feature to find out how to use subplot

Part II Theory Part:

1. Find the best approximation in the $L_2$-norm of $e^x$ on $x\in[0,1]$. using polynomials of at most degree $5$. Compare this norm to the norm of the polynomial obtained by computing the Taylor series of degree $5$ about $0$.
2. Construct the Lagrange interpolating polynomials for the following functions and find a bound for the absolute error on the interval $[x_0,x_n]$:
   1. $f(x) = e^{2 x} \cos 3 x$ where $x_0=0, x_1=0.3, x_2 = 0.6$, where $n=2$.
   2. $f(x) = \cos x + \sin x$, where $x_0=0, x_1=0.25, x_2 =0.5, x_3=1.0$, where $n=3$.
3. Find the polynomial of least degree that interpolates these sets of data:
   1. $(x_i,y_i) = (3,5), (7,-1)$
   2. $(x_i,y_i) = (3,12), (7,146), (1,2), (2,1)$
4. In class we discussed piece-wise linear interpolation and introduced the ``hat functions'' $B_i(x)$ ( $i=1,2,3,\ldots,n$). These were claimed to form a basis for ${\mathcal S}_1^0(\Delta)$, where $\Delta \equiv a=x_1 < x_2 < x_3 \ldots < x_{n-1} < x_n = b$. Show that these indeed form a basis.
5. Show that the Chebyshev polynomials, which are defined over $x \in [-1,1]$, have the following properties:
   1. $T_n(x) = cos(n \theta)$, where $\theta = arc cos( x)$.
   2. $\int_{-1}^{1} T_n(x) T_m(x) \frac{1}{\sqrt{1-x^2}} dx = \frac{\pi}{2} \delta_{n,m}$, where $\delta_{n,m}$ is the Kroneker delta function.
   3. They satisfy the ordinary differential equation
      
      $$(1-x^2) y'' - x y' + n^2 y = 0$$
      
      

   4. Show that $(1-x^2) T_n'(x) = n [ T_{n-1}(x) - x T_n(x)]$ and $2 T_n(x) T_m(x) = T_{n+m}(x) + T_{n-m}(x)$, for $n \ge m$.
   5. Plot a unit circle, centered at $0$. Project the location of equally spaced points on this circle to the horizontal axis $t$, and convince yourself that these lie at locations
      
      $$t_k = -cos( \frac{(2 k-1) \pi}{2 n}),$$
      
      
      with $k=1, 2, \ldots, n.$. These are called the Chebyshev interpolation points.

Part III Experimental Part: For the following exercises you can adapt the following code, which implements ``Neville's Algorithm'' for the generation of the Lagrange polynomials.

%Lagrange Polynomial Algorithm by 
%Neville Interpolation at a vector of %points.

% coded by Rachel Labes, Oct 2002


clc;
clear;
format long
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Input
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

n = input('Enter the number of points to interpolate at: ');
x = input('Enter the vector of points to interpolate at: ');
d = input('Enter the vector of function values of interpolation points: ');
t = input('Enter the value to approximate the function at: ');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Neville Method
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Q = zeros(n,n);	%Initializes an n x n matrix
Q(:,1) = d';	%Enters the vector of 
                % function values as the first column of Q.
Q;					

for i = 2:n
   for j = 2:i
      Q(i,j) =   ...
     [(t - x(i - 1))*Q(i,j-1) - (t - x(i))*Q(i-1,j-1)]/(x(i) - x(i-1));
   end
end
%Q(i,j) calculates the polynomial approximation for each matrix point.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Output%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Q		%Displays the complete matrix Q.

P = input('Enter the point in the matrix to be outputed: ');
disp('The value of the polynomial is approximately: ')
disp(P)		%Displays the function approximation for that Lagrange
		%polynomial.

1. For $P$ and $Q$ defined in the introduction evaluate $T = P*Q$ by hand and compare with the result you get from MATLAB. Plot $P$, $Q$, $T$ and the x-axis (e.g. using grid on) on the same graph. What do you observe about the roots of the three polynomials?

2. Use appropriate Lagrange interpolation polynomials of degree one, two three and four to approximate each of the following:

   1. $f(8.4)$ if $f(8) = 16.63553$, $f(8.1) = 17.61549$, $f(8.3) = 17.56492$, $f(8.6) = 18.50515$, $f(8.7) = 18.82091$.
   2. $f(-1/3)$ if $f(-1) = 0.10000000$, $f(-0.75) = -0.07181250$, $f(-0.5) = -0.02475000$, $f(-0.25) = 0.33493750$, $f(0) = 1.10100000$

3. To approximate the function $f(x) = \sqrt{x}$, we will use the points (1,1) (4,2) and (9,3).
   1. Write the formula for the Lagrange Polynomial $P_2(x)$ that interpolates these three points.
   2. Write a .m function file that implements $P_2(x)$ using the MATLAB polynomial commands.
   3. To see how well the approximation is:
      1. Plot $f$ and $P_2$ in the interval $[0.01,   12]$.
      2. Plot $f-P_2$ in the interval $[0.01,   12]$.
      3. Plot $abs((f-P_2)./f)$ in the interval $[0.01,   12]$.
   4. Using the plots determine where the approximation is better/worse.
   5. Why do you get these results? Do they agree with the formula derived in class for the error bound?
   6. Add the point (0,0) and repeat (a)-(c). Do you see any improvement? Deterioration? Where? Why?

4. In this exercise we will pretend that MATLAB stores the values of the function $e^x$ only for $x = 0, 0.1, 0.2, \dots$. Use these values and the appropriate Lagrange polynomials of degrees zero, one, two three and four to approximate $e^{0.743}$. How does the absolute error change?

5. We would like to approximate the function $h(x) = 1/(1+x^2)$ in the interval $[-5  5]$ using equally spaced x-values.
   1. Write a program that evaluates the corresponding Lagrange polynomials using $n$ points.
   2. Plot on one page (using subplot) four plots of the function versus the interpolation polynomial using 3,7,11 and 15 points.
   3. What do you observe?
   4. For each case, plot also the absolute error $\vert h-P_n\vert$ versus the theoretical error bound $\vert\prod_{j=0}^n (x-x_j)\vert$. You may want to plot them initially on separate graphs. Then, come up with a meaningful way to compare the plots using just one graph.
   5. What do you observe?

6. In this exercise you will analyze the problems that arise when the interpolation polynomial is evaluated using the Vandermonde matrix. We will find the interpolation polynomial $P_n$ for the nice smooth function sin (x) and see that as $n$ increases problems arise.

   To do so, run the following program (available as Vandermun.m on the class home page)

   %  Comparison of the interpolation polynomial  Pn
   %  for sin(x) in the interval [1 2] with sin(x)
   %  Pn(x) = c_1 x^n + ..+ c_n x + c_{n+1}
   %
     for i = 0:3
        Nx = 10*3^i;
        dx = 1/Nx;
        x = [1:dx:2]';                 
   % Find the interpolation polynomial using
   % the Vandermonde matrix
        V = vander(x);
        C = V\sin(x);  % Warning: use \, not /
   % Plot the results at grid values other than the 
   % ones used in the interpolation.
        y = x+0.1*dx*rand(size(x));
        subplot(2,2,i+1)
        plot(y,polyval(C,y)-sin(y))
        xlabel('x')
        ylabel('Pn(x)-f(x)')
        title(['Nx = ',num2str(Nx)])
     end
    figure(gcf);
   

   Do not submit the plots or the program, but answer the following questions:

   1. Using the error formula derived in class, show that $\vert\sin(x)-P_n(x)\vert\rightarrow0$ as $n\rightarrow\infty$ for $x\in[0,2]$.
   2. Does this agree with your plots?
   3. In practice, is larger $n$ good or bad?
   4. What is the source of the problem with large $n$? (Hint: Run the program again but use the grid points $x$ rather than $y$ in the plot command. Also note the MATLAB error messages.)