In this lab you will again be using MatLab.

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.

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.

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.

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 colormapAlso 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

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 a

_{1}.We want to find the value of the coefficient a

_{1}(i.e., the slope) that minimizes the sum squared error (sse) between the line and the data:

The quantities t

_{i}and y_{i }are the experimental values taken from the graph. The quantity (y(t_{i})-y_{i}) is just the error between the fitted line and the experimental value for the i^{th}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 a_{1} is varied from 0 to 6, with a difference of
0.1 between successive values (a_{1 }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 a_{1} -- you don't
need a loop to calculate an individual value of a_{1}. 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 a

_{1}that minimizesse. You can think of thesseas a function of a_{1}. You know that to minimize a function, you just take the derivative and set it to zero.

To find the value of a

_{1}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.

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 a

_{1}and slope a_{2}).

We want to find the values of the coefficients a

_{1}and a_{2}that minimizes the sum squared error (sse) between the line and the data:

The quantities t

_{i}and y_{i }are the experimental values taken from the graph. The quantity (y(t_{i})-y_{i}) is just the error between the fitted line and the experimental value for the i^{th}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 a

_{1}and a_{2}that minimizesse. You can think of thesseas a function of both a_{1}and a_{2}. 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, a_{1}and a_{2}, 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 a_{1}and a_{2}and set them both to zero.

The first equation can be simplified as

Similarly for the second equation

Since the t

_{i}and y_{i}are known, this simplifies to two equations in two unknowns, a_{1}and a_{2}. If we let

and

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

aunknown.

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.3063So 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.

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 a_{1}, a_{2}, and a_{3}.
Show your work, and include a plot of the parabola plotted along
with the original data.