E5 Matlab 2

In this lab you will again be using MatLab.


To Turn in:

By Tuesday 11/22/2005 turn in the answers to the questions in pink, one assignment per group.  Make sure everybody's name is on the report.  You may choose to work alone, but if you choose to work in a group, make sure that everybody in the group gets approximately the same amount of time at the keyboard.


Writing Programs in MatLab (m-files)

Last week you typed in commands to get results from MatLab.  However, you will often want to repeat a series of commands.  The simplest way to do this is through the use of m-files, or MatLab scripts.

To open a script start MatLab and then go to File->New->M-file.  In the window that opens, enter (or cut and paste) the following commands.

t=0:5; %Define a time vector.
y=[0.59 2.98 6.16 9.09 11.91 15.36]; %Define some data values.
plot(t,y,'r+');
xlabel('Time');
ylabel('Data');
title('A plot of some data');

Save your file on your network user folder (usually the "H:" drive).  Note that the filename can have no blank spaces in it.  In the editor window menu, go to Debug->Run.   You will get a dialog box asking if you want to change the working directory to the one on you user folder.  Answer "OK."

A plot should occur that is similar to the one shown below.

Let's add a little bit more to the plot.  Edit the script file so that is is as shown, and then run it.

t=0:5;      %Define a time vector.
y=[0.59  2.98  6.16  9.09  11.91 15.36];    %Define some data values.
plot(t,y,'r+',[-1 6],[-3 18],'b:');         %*
xlabel('Time');
ylabel('Data');
title('A plot of some data');
axis([-1 6 -1 16]);                         %*

1) Explain what the two lines marked with asterices do.  Use MatLab's help facility if necessary.

Some other plots

Generate a new m-file and enter the following:

x=0:0.1:4;
p1=[0.475 -2.150  3.500];
p2=[0.285 -2.214  5.499 -5.071  3.500];
y1=polyval(p1,x);
y2=polyval(p2,x);
plot(x,y1,x,y2);
xlabel('x');
ylabel('Polynomials');
legend('p1','p2');

You should get a graph similar to the one shown.  If you attended class on Tuesday and managed to stay awake, you'll note that these are the two curves I used to demonstrate minimization.

2) Write short comments for each of the lines in the m-file.  Again, you may find it helpful to refer to the MatLab help facility.

Some 3-D Plotting

At the Matlab prompt (in the command window) type the command

>> [x,y]=meshgrid(1:0.5:3,1:4)

3) Explain what the meshgrid command does.

Start a new m-file and enter and run the following code:

[x,y]=meshgrid(0:0.1:4,0:0.1:4);
z=polyval(p1,x).*polyval(p2,y);
surf(x,y,z)
xlabel('x');
ylabel('y');

You should get a surface. 

4) Explain it's shape in terms of the two polynomials.  It may be useful to rotate the image (hit the Rotate-3D button and use the mouse to rotate) until you are looking at the x or y axes.

Try the following while observing the figure

>> colormap(hot)
>> colormap(bone)
>> colormap(jet)

To see more colormaps, try

>> doc colormap

Also try the following (one at a time)

>> shading interp

>> light

>> lighting phong

>> material metal

>> material dull

>> alpha(0.5)

Change your perspective (make sure you can see your figure when you run these commands).

>> for i=1:100, camorbit(1,0); pause(0.1); end

>> for i=1:100, camzoom(1+0.001*i); pause(0.1); end

The Least Squares Method (a line through the origin y=mx).

You have probably all fit a line to a set of noisy data.  This can be made more precise through a technique called called the least squares method.  The graph below represents a series of measured values of a variable y, taken at certain times t, and a crude attempt to fit a line through the points.  Note that the line goes through the origin so it is defined by a single parameter, the slope (m).

The data values are given as above

t=0:5; %Define a time vector.
y=[0.59 2.98 6.16 9.09 11.91 15.36]; %Define some data values.

To find the least squares line that best fits the data, we first define the function for a line  that goes through the origin with slope a1.

We want to find the value of the coefficient a1 (i.e., the slope) that minimizes the sum squared error (sse) between the line and the data:

The quantities ti and yi are the experimental values taken from the graph.  The quantity (y(ti)-yi) is just the error between the fitted line and the experimental value for the ith datum.  We square it so the positive errors and the negative errors both contribute to the sum squared error.  The summation goes to 6 because there were 6 data values. We can re-write the expression as:

5)  Write a script that calculates and plots sse on the vertical axis as the value of a1 is varied from 0 to 6, with a difference of 0.1 between successive values  (a1 should be the horizontal axis).  Obviously, by looking at the data, the minimum of sse should be when a1 is approximately 3.  You will probably need to write a for-loop to do this to iterate through the values of a1 -- you don't need a loop to calculate an individual value of a1.  See below for a hint on how you might calculate sse (you need to use MatLab's "sum" function)..

Now the task at hand is to find the value of a1 that minimize sse.  You can think of the sse as a function of a1.  You know that to minimize a function, you just take the derivative and set it to zero. 

To find the value of a1 that minimizes sse, set the derivative equal to zero.

or

The quantity a1 is easily calculated using MatLab's "sum" command.  (The "sum" command can also be used to calculate sse).

>> a1=sum(y.*t)/sum(t.^2)

6) Explain how the line given above correlates to the equation for a1.  In particular, explain the need for the two periods ("y.*t", "t.^2").  What is the value for a1?

7) Plot the original data (as red circles) and the least squares fit to the data (a dashed blue line).  Turn in the plot.

The Least Squares Method (y=mx+b).

Start a new m-file and enter and run the following code:

ti=[2.5 5.0 7.0 11.0 12.0];
yi=[1.2 2.2 4.0 10.0 13.8];
plot(ti,yi,'rx')
xlabel('t')
ylabel('y')
title('E5 Plot, t vs y')
axis([0 15 0 15])

You should get a graph like the one shown.

To find the least squares line that best fits the data, we first define the function for a line (with intercept a1 and slope a2).

We want to find the values of the coefficients a1 and a2 that minimizes the sum squared error (sse) between the line and the data:

The quantities ti and yi are the experimental values taken from the graph.  The quantity (y(ti)-yi) is just the error between the fitted line and the experimental value for the ith datum.  We square it so the positive errors and the negative errors both contribute to the sum squared error.  The summation goes to 5 because there were 5 data values. We can re-write the expression as:

Now the task at hand is to find the values of a1 and a2 that minimize sse.  You can think of the sse as a function of both a1 and a2.  You know that to minimize a function, you just take the derivative and set it to zero.  Because we have a function of two variables, a1 and a2, we need to use what are called partial derivatives.  A partial derivative for our purposes is the same as a standard derivative, you just assume the other variables are constant.  A partial derivative is denoted by the notation "", rather than the standard "d" for a standard derivative.  So we take the derivative of the sum squared error, sse, with respect to a1 and a2 and set them both to zero.

The first equation can be simplified as 

Similarly for the second equation

Since the ti and yi are known, this simplifies to two equations in two unknowns, a1 and a2.   If we let 

and 

We can reframe the problem as a matrix problem, with the vector a unknown.

The only difficulty remains calculating the sums.  

>> ti=[2.5 5.0 7.0 11.0 12.0];
>> yi=[1.2 2.2 4.0 10.0 13.8];
>> c11=length(ti);
>> c12=sum(ti);
>> c21=c12;
>> c22=sum(ti.*ti);
>> b1=sum(yi);
>> b2=sum(ti.*yi);
>> C=[c11 c12; c21 c22]
C =
    5.0000   37.5000
   37.5000  345.2500

>> b=[b1 b2]'
b =
   31.2000
  317.6000

>> a=inv(C)*b
a =
   -3.5569
    1.3063

So the best fit line is y(t)=-3.56+1.31t.

>> y=a(1)+a(2)*ti;
>> plot(ti,yi,'rx',ti,y,'b:')
>> axis([0 15 0 15])
 

Make a surface plot of sse 

with the two axes being a1 and a2.  Make sure the axes are labeled.  The MatLab command "surfc" also plots a contour plot and can make the minimum of the surface easier to find.

8) Turn in the plot and show that the minimum is where you expected it according to your calculations above.  Note: the minimum can be hard to see because the error gets large very quickly.  To accentuate the minimum, you can plot the logarithm of the surface.  This tends to de-emphasize large values relative to the smaller (minimum) value.

Extending the Least Squares Method

This part is optional, only do it if you have the time and inclination to do more.

The method described above can be extended to any polynomial.  For example, it could be used to find the three coefficients that define the parabola given by:

9) Extend the least squares method to a second order polynomial, and find the best fit values for a1, a2, and a3.  Show your work, and include a plot of the parabola plotted along with the original data.