M439 Laboratory on Lagrange Polynomials in MATLAB

Download the following two files to your N:\m439 folder:

lagran.m

1.  MATLAB stores polynomials by storing their coefficients in a row vector ( 3 x 1 matrix)
Hence the polynomial  x2 + 2x +3  would be stored as the vector [1 2 3]

» P = [1 2 3]
P =
     1     2     3

2.  MATLAB has a several routines that allows one to manipulate polynomials.
We can evaluate polynomials at a particular point with the command polyval
The following evaluates the polynomial with coefficients stored in P above at the value x = 2:

» polyval(P,2)
ans =
    11

3. The conv command can be used to multiply two polynomial.  For example, to multiply the above polynomial by (x + 1) -- i.e. to compute (x2 + 2x + 3) (x + 1):

» Q = [1 1]
Q =
     1     1

» conv(P,Q)
ans =
     1     3     5     3
( The result is thus x3+ 3x2 + 5x + 3. )

4.  The MATLAB Symbolic Toolbox also has some commands for manipulating polynomials symbolically.  The poly2sym command converts a vector of coefficients to a symbolic polynomial.

» syms x
» poly2sym(P,x)
ans =
x^2+2*x+3

5.  To find the  polynomial of degree n that passes through n + 1 data points we can use polyfit.
» help polyfit

To find the polynomial passing through the points in Example 4.4, page 200:
» X = [1 2 3 5]
» Y = [1.06 1.12 1.34 1.78]
» P = polyfit(X,Y,3)
P =
   -0.0200    0.2000   -0.4000    1.2800
 
 

Thus the polynomial fitting these points is

-0.02x3 +0.2x2-0.4x +1.28.
You can see that MATLAB produces the correct coefficients agreeing with your text.

The  function polyfit uses the linear algebra method of solving a system of linear equations to find the coefficients.  We are using the  Lagrange method which will give the same result, but is much less subject to computer round-off error in the calculations.

Constructing the Lagrange Interpolating Polynomial


6.  In order to construct the Lagrange Coefficients for the Lagrange Polynomial in MATLAB, we can use the built-in function poly, which constructs a polynomial with given roots.
Enter the following to construct a polynomial with roots 1 and 2 for example:

» poly([1 2])
ans =
     1    -3     2

Thus this is the polynomial (x-1)(x-2) = x2 -3x +2 which has roots 1 and 2.

Example: (Exercise 2a, Section 4.3, page 218).   Lagrange interpolating polynomial of degree two (quadratic) approximating f(x) = 1 + 2/x  at points x0 = 1, x1 = 2, and x3 = 2.5.

Looking at our formula for the Lagrange coefficient polynomial L2,0 we can see that we need the numerator to be the polynomial with roots x1 = 2, and x2  = 2.5:  (x -2)(x-2.5)  The denominator is the constant (x0 -x1)*(x0-x2) = (1 -2)*(1-2.5).

7.  Assuming that we want to store the Lagrange coefficient polynomials in the 3x3 array L  (with the 1st row being the coefficients for L2,0,  the 2nd row being the coefficients for L2,1, and the third row being the coefficients for L2,2 ), we proceed as follows:

» L(1,:)= poly([2 2.5])/((1 - 2)*(1 - 2.5))
L =
    0.6667   -3.0000    3.3333

» L(2,:)= poly([ 1 2.5])/((2 - 1)*( 2 - 2.5))
L =
    0.6667   -3.0000    3.3333
   -2.0000    7.0000   -5.0000

» L(3,:)= poly([1 2])/((2.5 - 1)*(2.5 - 2))
L =
    0.6667   -3.0000    3.3333
   -2.0000    7.0000   -5.0000
    1.3333   -4.0000    2.6667

8.  The final Lagrange polynomial is y0*L2,0 (x)+ y1*L2,1(x)+y2*L2,2(x).
Since y0 = f(x0) = 1+2/1 = 3, and similarly y1 = 3, y2 = 3.3, we compute the polynomial P as

» P = 3*L(1,:) + 3*L(2,:) + 3.3*L(3,:)

P =
    0.4000   -1.2000    3.8000

9.  To see this formatted as a polynomial in x enter:

» pretty(poly2sym(P))

10.  To evaluate the polynomial at 1.5  we use:

» polyval(P,1.5)

ans =
    2.9000

11. Comparing the true value of f(x) = x + 2/x at 1.5
» 1.5 +2/1.5

ans =
    2.8333

Similarly for x = 1.2:
» polyval(P,1.2)

» 1.2 + 2/1.5

12.  We can plot both of these on a graph to compare the two:
First create the symbolic polynomial from the coefficients in P:
» SP = poly2sym(P)

SP =
 2/5*x^2-6/5*x+19/5

13.  Plot the function f(x) = x + 2/x on the range x = .5 to 3
» ezplot('x + 2/x', [.5 3])

Plot the polynomial Lagrange interpolating polynomial SP on the same range on the same graph:
» hold on; ezplot(SP,[0 3])

The Code for Implementing an Algorithm for Computing the Lagrange Polynomial.

14.  Open the file lagran.m in the MATLAB editor
Note the interface for the m-file function lagran:

function [C,L]=lagran(X,Y)

%Input  - X is a vector that contains a list of abscissas
%       - Y is a vector that contains a list of ordinates
%Output - C is a matrix that contains the coefficents of
%         the Lagrange interpolatory polynomial
%       - L is a matrix that contains the Lagrange
%         coefficient polynomials

15.  Thus to call this function we set up the vectors X and Y with the x and y coefficients of the interpolating points.  Then call the function to return the interpolating polynomial in C and the Lagrange coefficients for that polynomial in L.  For our above example, it would be:

» X = [1 2 2.5]
» Y = [3 3 3.3]
» [C L] = lagran(X,Y)

We see that the same answer is returned as the one we computed step by step above.

Look at how the coefficients are computed in the body of the function

for k=1:n+1    %  Calculate each of n+1 Lagrange coefficient
   V=1;        %  Accumulate computations in V temporarily

   for j=1:n+1
                         % Multiply by (x - X(j))/(X(k) - X(j))
      if k~=j  %  Be sure to skip the k'th one
         V=conv(V,poly(X(j)))/(X(k)-X(j));
      end

   end
   L(k,:)=V;    %  Store Lagrange coefficient in kth row of L
end

The final step is to compute the interpolating polynomial:
C = Y*L
This uses matrix multiplication to compute C =  Y * L. It is not difficult to verify using the rules for matrix multiplication that this gives the correct polynomial.  For example, the first entry in C will be
Y(1)*L(1,1) + Y(2)*L(2,1) + Y(3)*L(3,1)  which is the correct coefficient of x^2 term.

16.  Exercise:  Use lagran to calculate the third degree Lagrange polynomial for cos(x) at the evenly spaced interpolating points for the x-values (abscissas) of 0.0, 0.4, 0.8, and 1.2  (see page 210-211, Example 4.7, part b), text for the computations and page 210 figure 4.12b)  comparing the graph of this polynomial to the graph for  function cos(x).

» X = [0.0 0.4 0.8 1.2]

» Y = cos(X)

» [C , L ] = lagran(X,Y)

Compare the graph of  y = cos(x) to your computed polynomial with coefficients stored in C with the following commands

» SP = poly2sym(C)
» clf
» ezplot('cos(x)',0,2)
» hold on
» ezplot(SP,0,2)

We see they look quite close at least through the interpolating range 0 to 1.2.

Use the error formula of Theorem 4.4 (see page 213 text) to compute an upper bound of the error for these nodes (interpolating points).
h = distance between evenly spaced nodes = .4 in this case.
M4 is the maximum of the absolute value of the  4th derivative of f(x) = cos(x)  on the interpolation range -- we can use 1 for this.   Thus our error should be <= (.4)^4/25:

» (.4)^4/24

Does this appear to be true? --Calculate some of the errors  with the following commands  (and compare to table 4.7, page 216)

» cos(.1)
» polyval(C,.1)
» cos(.1)-polyval(C,.1)

Laboratory Exercise to Turn In:   Page 220 Algorithms and Programs 2 a) b) c) 

Print your graph for c). Write on the print out what the approximating polynomial is returned by lagran for part a)
For part b)  also write on your printout, what the integral polynomial was, and your answer for the average.
Turn in by next lab period, Friday Feb 24th.

For part a)   Use the lagran function.  Set up your X, Y vectors   ( X is time values,  Y is temp values ).

For part c)  Use the following commands to plot your data points and the corresponding Lagrange polynomial stored in C
» SP = poly2sym(C)
» clf                         Clear any previous figures
» plot(X,Y,'r*')           Plot the original data points in blue starts
» hold on;                Retain the graph for next plot also.
» ezplot(SP,0,7)

For part b -- (See Theorem 1.10 -- Mean Value Theorem for Integrals, page 6).   Average temperature is the integral of the temperature function (your approximating polynomial for temperature)  over the limits of the time interval  divided by the length of the time interval.  

Thus you need to:
1.  Write an m-file function called integ.m.  This will start with the headng
             function D = integral(C)

The body of the function will create a new polynomial --coefficients stored in D, whose length will be one more than C (its last entry will be 0 for the constant term) and whose coefficients are those of the integral of the polynomial stored in C. 

2.  After calling the function with the command D = integral(C),  you can use polyval on D to evaluate indefinite integral at the endpoints.