# E5 Matlab 2

In this lab you will again be using MatLab.

## To Turn in:

By Friday, 10/1/10 at 5 pm, post on your website the answers to the questions in pink, as with Matlab 1.  You are encouraged to work alone, with help from Wizards or classmates, but if you work with other, please make sure that all collaborator's names are on the posted document.

## Writing Programs (called m-files or scripts) in MatLab

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 or m-file start MatLab and then go to File->New->M-file

1) Begin by writing a program that prompts the user for the (x,y) coordinates of two points in the Cartesian plane, then prints out the distance between the points.

Save your file on your network user folder (usually the "H:" drive).  Note that the filename can have no blank spaces in it, and should not begin with a number.  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."

For the second program open a new M-file and 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');

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 or m-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]);                         %*```

2) 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.

3) 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)`

4) 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.

5) Explain its 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 fit a line to a set of noisy data.  This procedure can be made more precise through a technique called the least squares method (because it minimizes the sum of the squarea of differences between data and interpolation).  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:

6)  Write a script (m-file) 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 as shown in the following Matlab code line:   (The "sum" command can also be used to calculate sse).

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

7) Explain how this code line 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?

8) Plot the original data (as red circles) and the least squares fit to the data (a dashed blue line).  Post the plot on your website, along with the answers to the previous (pink) questions.

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

9) Post on your website 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:

10) 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.