Home » Courses » Nuclear Science and Engineering » Essential Numerical Methods » Tutorial Videos » Octave/MATLAB® for Beginners, Part 2: Fitting Data and Plotting
Flash and JavaScript are required for this feature.
Download the video from iTunes U or the Internet Archive.
Description: Plot the line that results from our fit on the data. Construct a matrix with looping, inverse, transposition functions. Use 'hold on' & 'hold off' features. Manual axis configuration. Using semicolon to clean up command line display. Debugging errors.
Instructor: Prof. Ian Hutchinson
Octave/MATLAB® for Beginne...
PROFESSOR: OK. Here we are. We've got our data. We've managed to plot it. Now let's go about constructing a matrix which will enable us to fit a polynomial to this data. First of all, I need a matrix. Here is a way to construct a matrix. We put s equal to-- The function zeros is a function which produces a matrix with all zeros in it, and I want this matrix to have dimensions endpoints by endpoints. So this command should produce such a matrix, and lo and behold, it does.
Let me just clean up this display a little bit. At the moment, it's a little bit annoying that I'm plotting out absolutely everything that I have. This won't be too much trouble when I only have six points, but it might be in a lot of trouble if I had a lot more points.
Here's how to stop our MATLAB or Octave from printing out the results of a command. All you do is you put a semicolon at the end of the line, like that. Or, like that. If I make that change, so you hit and run again, it no longer prints out x and y. It's still doing the plotting, and it's printing out s as being zeros of endpoints endpoints.
That's not terribly interesting either, so I'm going to put a code on at the end of that line and prevent s from being printed out. Now, zeros is not what I want. What I want is a matrix which has the appropriate values to enable me to fit polynomials. And those values I'm going to generate by using a loop.
In general, it's not a great idea within MATLAB or Octave to use loops, but it's perfectly possible to do so. Loops are generated in the following way. Supposing I choose i to be the index of my loop, I can let i range from one to endpoints by writing 4i equals 1 Coulomb endpoints.
The end of that for loop is indicated by writing end. Actually, in Octave, it's endfor, but in MATLAB it's end, and so I'm just going to put end and that works in Octave as well. So that would loop over one index, but actually I want two indices. So I'm going to put another loop inside of the first loop that I had and I'm going to make its index j and I'm going to let it also run from one to endpoints.
And I need to end that loop, too. Let's do my indentation neatly so I can understand this loop when I come back and read it later, and now that's a loop which will loop over all of the indices of the matrix s. To refer to a particular element of the matrix s, I write s of i and j, where i and j are indices.
So I'm going to set the ij-th element of s equal to an appropriate value to use for polynomial fitting, and that value is x-- actually x I have to refer to by its index, so that's xi-- and I'm going to raise it to the power j minus 1.
So if j is 1, which is the lowest element, this would raise it to the 0 power, which would just give me the result of being x to the 0, which is 1, and for higher powers, it would go on appropriately. I don't want to print this out at every time that I calculate it, so I'm going to end the line with a semicolon. And then at the end, perhaps I do want to see s, so let me just write s, and that will cause s to be printed out.
So here I am. I run it again, and now s consists of this 6 by 6 matrix whose values, if I look at the top line, you can see that-- Well, first of all, let's look at the first column. The first column is all the values of j equals 1, and they're all x to the 0. And then this next column is x to the power 1, so you can see that it ranges from a 6 to 1.
The next column is x to the power 2, and so it runs from a smaller value up to 1. So of course all the values along the bottom row are 1, and the values along the top row are 1/6 to the power 0, 1, 2, 3, and so forth. So that's my s and that's what I'm going to use as my matrix.
Now, I can take the inverse of s by writing some new variable name. I'm going to call it sinv to remind me that it's the inverse of s, and put it equal to a function which is called inv, which stands for inverse, called s. So this command finds me the inverse of the matrix s. And when I run that, I find that s inverse is given by this matrix here, whose values will be very hard to guess.
Now, the way to find the coefficients of a polynomial fit to my original data are to put those coefficients equal to the matrix product of the inverse of s, which is sinv times the y values, which we learned about in the lecture. So that would give me the coefficients. When I run that, I find out what the coefficients, or the different powers, of x are.
So in other words, these coefficients-- the first line is the coefficient of the power x to the power 0, the next one is x to the power 1, and so on. So I've actually found the coefficients of the fit. For the moment, it's not obvious whether I've done it correctly. I could certainly plot out the values of this fit at various places, and that's what I'm going to do.
I'm going to plot it actually with more points than I started with. So let me define a new parameter, np2. I might make it equal to, let's say, 40. I could then have a second variable, x2, equal to something running from 1 to np2 transposed, divided by np2. So that will be an x array running from a small value to 1 with, in this case, 40 steps.
I'll set, for now, y2 equal to x2. That's just created an array that is of the same length as x, and is actually equal to x, but I'm actually going to end up setting y to be something else in a minute. And let me, for now, put fj, a new variable, equal to the c that I've just calculated.
This won't do anything interesting, because I've commented those out, and so when I run it, I just get the same result that I had before. So I haven't actually done anything. Now I'm going to have another couple of nested for loops. I'm going to say for j equals 1,np2, and for i equals 1:np2.
Let's f set the value of fj, which is a matrix of the vector, strictly speaking, column vector of the same length as c. Let's set that the i-th element of that equal to the x value, which is x2, because this is the second array that I'm creating, evaluated at parameter j to the power j minus 1. So that's actually evaluating x2 to the appropriate power.
So it shouldn't have been j minus 1. It should have been i minus 1. So this is x2 j to the i minus 1. Again, I don't want to print all of those out, so I'm going to put a semicolon at the end. So that's calculated by fj value, and now I'm going to put y2 of j equal to-- Well, what I basically want to do is form the dot product between c, the coefficients, and the value of fj that I've just calculated. So I want c1 times fj1 plus c2 times fj2 and so forth.
One way of doing that-- there are various different ways of doing it-- is to take the transpose of c. C remember, is a column vector. Let's transpose it. That's going to give me a row vector. I can then take the product of that row vector with the column vector fj. And so taking to the matrix product of a column vector times a row vector is equivalent to taking the dot product of two vectors.
So that has produced a value y2 of j, which is equal to the sum of the coefficients times the values of x to the appropriate power. So that's what I want to form in that matrix. Again, I'm--
It seems as though I've made a mistake. And I'm not quite sure what I did wrong, but I need to find out. So let's have a look. It's saying that I made my mistake at line 29 of column 8. 29 is this, and it's saying that I have a 1 by 6, and op2 is 40 by 1. I don't know why that is 40 by 1. I'm a bit surprised by it, because-- Ah. I see what I did wrong. This should have been endpoints.
So what I did was, this inner for loop shouldn't have been over all of the j values. The j values refer to the length of the new vectors that I produced. It should've been over the shorter version. So let's see. Let's save that and hope that's corrected my error. Yes. Lo and behold, it has. I haven't actually done anything to show what the result is, but at least now I've got the result.
Now what I want to do is, I want to plot these two new matrices, x2 and y2, and actually I'd like to overplot them on my old plot. If I just say plot like this, save this and plot it again, that isn't what happens. Instead, what happens is just the curve that I've just calculated is plotted, and my previous plot is wiped out.
What I can do to prevent that happening is to say, hold on. Only MATLAB slash Octave would have something like hold on as a command. But anyway, hold on basically says, retain the data that you've already got in this plot and add some more data on. So let's try this.
So now what we see is the data that I've plotted out in my first plot, which was up here, is held and the second plot is plotted. Notice that this fit does indeed go through the six points that I originally retained, and yet, that line actually does some funny things. Particularly at the end here, it goes negative.
There's another trick that one can do in plotting, and one can say things like axis manual, and I think that, if I'm lucky, that will fix this plot. No, it won't. Never mind. Anyway, there's a way-- and I'll show it later-- to prevent the plot rescaling and the required command is axis manual.
But at the moment, I think, before I do it, I've got to say hold off. And then if I re-run it, I will find that lo and behold, I've got a plot of the six points and the line that I fitted, which is a polynomial fitted to those lines. And lo and behold, it fits.
MIT OpenCourseWare makes the materials used in the teaching of almost all of MIT's subjects available on the Web, free of charge. With more than 2,200 courses available, OCW is delivering on the promise of open sharing of knowledge. Learn more »
© 2001–2015
Massachusetts Institute of Technology
Your use of the MIT OpenCourseWare site and materials is subject to our Creative Commons License and other terms of use.