M439 Laboratory on Interpolation with  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 (1,1.06), (2, 1.12), (3, 1.34), (5, 1.78) use:

» 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.

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: Lagrange interpolating polynomial of degree two 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)   Should get 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  Compare 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.


 

Laboratory Exercise to Turn In:  

Consider the ongoing example from the e-text on the upward velocity of a rocket, data given below.

 

Time t in s.

Velocity V(t)  m/s

10

227.04

15

362.78

20

527.35

22.5

602.97

 

 

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.
 

a)   Use the lagran function to compute the third degree interpolating polynomial for the above data points.    Set up your X, Y vectors   ( X is time values,  Y is velocity values ).

b)   Write an m-file function called integral.m.  This will start with the heading
             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 antiderivative (indefinite integral) of the polynomial stored in C. 

c)  Call the function with the command D = integral(C).

 

d)  Use polyval on D to evaluate indefinite integral at the endpoints in order to compute the distance covered from the record from time 11 s. to time 16 s.

e)   (Review the Mean Value Theorem for Integrals):   The average of a continuous function over an interval is computed as the integral of the function (your approximating polynomial for temperature)  over the limits of the interval  divided by the length of the interval.     Compute the average velocity over the interval from 11s. to time 16 s. using your answer to d.


f  )  Print the data points  above and the polynomial you found in part a) on the same plot.  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 original data points in blue stars
» hold on;                Retain the graph for next plot also.
» ezplot(SP,10,25)