1 00:00:08,730 --> 00:00:10,320 PROFESSOR: OK. 2 00:00:10,320 --> 00:00:11,330 Here we are. 3 00:00:11,330 --> 00:00:12,800 We've got our data. 4 00:00:12,800 --> 00:00:15,540 We've managed to plot it. 5 00:00:15,540 --> 00:00:21,960 Now let's go about constructing a matrix which 6 00:00:21,960 --> 00:00:27,670 will enable us to fit a polynomial to this data. 7 00:00:27,670 --> 00:00:30,160 First of all, I need a matrix. 8 00:00:30,160 --> 00:00:32,759 Here is a way to construct a matrix. 9 00:00:32,759 --> 00:00:41,760 We put s equal to-- The function zeros is a function which 10 00:00:41,760 --> 00:00:45,040 produces a matrix with all zeros in it, 11 00:00:45,040 --> 00:00:50,190 and I want this matrix to have dimensions endpoints 12 00:00:50,190 --> 00:00:52,330 by endpoints. 13 00:00:52,330 --> 00:00:56,400 So this command should produce such a matrix, 14 00:00:56,400 --> 00:00:59,240 and lo and behold, it does. 15 00:00:59,240 --> 00:01:07,490 Let me just clean up this display a little bit. 16 00:01:07,490 --> 00:01:09,230 At the moment, it's a little bit annoying 17 00:01:09,230 --> 00:01:14,340 that I'm plotting out absolutely everything that I have. 18 00:01:14,340 --> 00:01:17,162 This won't be too much trouble when I only have six points, 19 00:01:17,162 --> 00:01:18,620 but it might be in a lot of trouble 20 00:01:18,620 --> 00:01:20,820 if I had a lot more points. 21 00:01:20,820 --> 00:01:25,860 Here's how to stop our MATLAB or Octave from printing out 22 00:01:25,860 --> 00:01:28,190 the results of a command. 23 00:01:28,190 --> 00:01:31,480 All you do is you put a semicolon 24 00:01:31,480 --> 00:01:33,320 at the end of the line, like that. 25 00:01:33,320 --> 00:01:35,180 Or, like that. 26 00:01:35,180 --> 00:01:39,780 If I make that change, so you hit and run again, 27 00:01:39,780 --> 00:01:42,080 it no longer prints out x and y. 28 00:01:42,080 --> 00:01:49,626 It's still doing the plotting, and it's printing out s 29 00:01:49,626 --> 00:01:52,630 as being zeros of endpoints endpoints. 30 00:01:52,630 --> 00:01:54,370 That's not terribly interesting either, 31 00:01:54,370 --> 00:01:57,880 so I'm going to put a code on at the end of that line 32 00:01:57,880 --> 00:02:01,570 and prevent s from being printed out. 33 00:02:01,570 --> 00:02:03,590 Now, zeros is not what I want. 34 00:02:03,590 --> 00:02:10,152 What I want is a matrix which has the appropriate values 35 00:02:10,152 --> 00:02:14,950 to enable me to fit polynomials. 36 00:02:14,950 --> 00:02:20,260 And those values I'm going to generate by using a loop. 37 00:02:20,260 --> 00:02:24,030 In general, it's not a great idea within MATLAB or Octave 38 00:02:24,030 --> 00:02:28,330 to use loops, but it's perfectly possible to do so. 39 00:02:28,330 --> 00:02:33,370 Loops are generated in the following way. 40 00:02:33,370 --> 00:02:37,150 Supposing I choose i to be the index of my loop, 41 00:02:37,150 --> 00:02:41,660 I can let i range from one to endpoints 42 00:02:41,660 --> 00:02:48,370 by writing 4i equals 1 Coulomb endpoints. 43 00:02:48,370 --> 00:02:58,380 The end of that for loop is indicated by writing end. 44 00:02:58,380 --> 00:03:02,640 Actually, in Octave, it's endfor, but in MATLAB it's end, 45 00:03:02,640 --> 00:03:05,340 and so I'm just going to put end and that works in Octave 46 00:03:05,340 --> 00:03:06,770 as well. 47 00:03:06,770 --> 00:03:09,660 So that would loop over one index, 48 00:03:09,660 --> 00:03:13,320 but actually I want two indices. 49 00:03:13,320 --> 00:03:17,125 So I'm going to put another loop inside of the first loop 50 00:03:17,125 --> 00:03:21,320 that I had and I'm going to make its index j 51 00:03:21,320 --> 00:03:26,590 and I'm going to let it also run from one to endpoints. 52 00:03:26,590 --> 00:03:31,410 And I need to end that loop, too. 53 00:03:31,410 --> 00:03:34,860 Let's do my indentation neatly so I 54 00:03:34,860 --> 00:03:41,790 can understand this loop when I come back and read it later, 55 00:03:41,790 --> 00:03:46,625 and now that's a loop which will loop over all 56 00:03:46,625 --> 00:03:51,270 of the indices of the matrix s. 57 00:03:51,270 --> 00:03:55,530 To refer to a particular element of the matrix s, 58 00:03:55,530 --> 00:04:01,760 I write s of i and j, where i and j are indices. 59 00:04:01,760 --> 00:04:06,020 So I'm going to set the ij-th element of s 60 00:04:06,020 --> 00:04:13,510 equal to an appropriate value to use for polynomial fitting, 61 00:04:13,510 --> 00:04:19,820 and that value is x-- actually x I have to refer to 62 00:04:19,820 --> 00:04:24,390 by its index, so that's xi-- and I'm going to raise it 63 00:04:24,390 --> 00:04:28,780 to the power j minus 1. 64 00:04:28,780 --> 00:04:31,280 So if j is 1, which is the lowest element, 65 00:04:31,280 --> 00:04:35,180 this would raise it to the 0 power, which would just give me 66 00:04:35,180 --> 00:04:38,420 the result of being x to the 0, which is 1, 67 00:04:38,420 --> 00:04:46,550 and for higher powers, it would go on appropriately. 68 00:04:46,550 --> 00:04:50,120 I don't want to print this out at every time 69 00:04:50,120 --> 00:04:51,520 that I calculate it, so I'm going 70 00:04:51,520 --> 00:04:54,670 to end the line with a semicolon. 71 00:04:54,670 --> 00:04:57,690 And then at the end, perhaps I do want to see s, 72 00:04:57,690 --> 00:05:02,320 so let me just write s, and that will cause s to be printed out. 73 00:05:02,320 --> 00:05:04,480 So here I am. 74 00:05:04,480 --> 00:05:10,220 I run it again, and now s consists of this 6 75 00:05:10,220 --> 00:05:16,680 by 6 matrix whose values, if I look at the top line, 76 00:05:16,680 --> 00:05:19,840 you can see that-- Well, first of all, 77 00:05:19,840 --> 00:05:22,190 let's look at the first column. 78 00:05:22,190 --> 00:05:26,590 The first column is all the values of j equals 1, 79 00:05:26,590 --> 00:05:28,060 and they're all x to the 0. 80 00:05:30,730 --> 00:05:35,270 And then this next column is x to the power 1, 81 00:05:35,270 --> 00:05:40,990 so you can see that it ranges from a 6 to 1. 82 00:05:40,990 --> 00:05:45,380 The next column is x to the power 2, 83 00:05:45,380 --> 00:05:47,800 and so it runs from a smaller value up to 1. 84 00:05:47,800 --> 00:05:53,070 So of course all the values along the bottom row are 1, 85 00:05:53,070 --> 00:05:58,580 and the values along the top row are 1/6 to the power 0, 86 00:05:58,580 --> 00:06:01,420 1, 2, 3, and so forth. 87 00:06:01,420 --> 00:06:08,870 So that's my s and that's what I'm going to use as my matrix. 88 00:06:08,870 --> 00:06:15,740 Now, I can take the inverse of s by writing 89 00:06:15,740 --> 00:06:21,160 some new variable name. 90 00:06:21,160 --> 00:06:24,140 I'm going to call it sinv to remind me 91 00:06:24,140 --> 00:06:26,180 that it's the inverse of s, and put 92 00:06:26,180 --> 00:06:29,300 it equal to a function which is called inv, which 93 00:06:29,300 --> 00:06:33,620 stands for inverse, called s. 94 00:06:33,620 --> 00:06:38,285 So this command finds me the inverse of the matrix s. 95 00:06:38,285 --> 00:06:42,160 And when I run that, I find that s inverse 96 00:06:42,160 --> 00:06:49,055 is given by this matrix here, whose values will 97 00:06:49,055 --> 00:06:52,560 be very hard to guess. 98 00:06:52,560 --> 00:06:57,360 Now, the way to find the coefficients 99 00:06:57,360 --> 00:07:04,030 of a polynomial fit to my original data 100 00:07:04,030 --> 00:07:08,770 are to put those coefficients equal to the matrix product 101 00:07:08,770 --> 00:07:16,060 of the inverse of s, which is sinv times the y values, which 102 00:07:16,060 --> 00:07:19,190 we learned about in the lecture. 103 00:07:19,190 --> 00:07:23,250 So that would give me the coefficients. 104 00:07:23,250 --> 00:07:27,420 When I run that, I find out what the coefficients, 105 00:07:27,420 --> 00:07:31,550 or the different powers, of x are. 106 00:07:31,550 --> 00:07:37,630 So in other words, these coefficients-- 107 00:07:37,630 --> 00:07:42,590 the first line is the coefficient of the power x 108 00:07:42,590 --> 00:07:48,640 to the power 0, the next one is x to the power 1, and so on. 109 00:07:48,640 --> 00:07:51,920 So I've actually found the coefficients of the fit. 110 00:07:51,920 --> 00:07:58,260 For the moment, it's not obvious whether I've done it correctly. 111 00:07:58,260 --> 00:08:03,980 I could certainly plot out the values 112 00:08:03,980 --> 00:08:08,060 of this fit at various places, and that's 113 00:08:08,060 --> 00:08:09,585 what I'm going to do. 114 00:08:09,585 --> 00:08:14,070 I'm going to plot it actually with more points 115 00:08:14,070 --> 00:08:15,340 than I started with. 116 00:08:15,340 --> 00:08:23,100 So let me define a new parameter, np2. 117 00:08:23,100 --> 00:08:24,825 I might make it equal to, let's say, 40. 118 00:08:27,890 --> 00:08:35,559 I could then have a second variable, x2, 119 00:08:35,559 --> 00:08:44,930 equal to something running from 1 to np2 120 00:08:44,930 --> 00:08:48,540 transposed, divided by np2. 121 00:08:48,540 --> 00:08:52,710 So that will be an x array running from a small value 122 00:08:52,710 --> 00:08:57,420 to 1 with, in this case, 40 steps. 123 00:08:57,420 --> 00:09:02,580 I'll set, for now, y2 equal to x2. 124 00:09:02,580 --> 00:09:07,720 That's just created an array that 125 00:09:07,720 --> 00:09:13,990 is of the same length as x, and is actually equal to x, 126 00:09:13,990 --> 00:09:16,710 but I'm actually going to end up setting y 127 00:09:16,710 --> 00:09:18,990 to be something else in a minute. 128 00:09:22,900 --> 00:09:32,812 And let me, for now, put fj, a new variable, equal to the c 129 00:09:32,812 --> 00:09:33,895 that I've just calculated. 130 00:09:38,542 --> 00:09:40,000 This won't do anything interesting, 131 00:09:40,000 --> 00:09:42,470 because I've commented those out, and so when I run it, 132 00:09:42,470 --> 00:09:44,402 I just get the same result that I had before. 133 00:09:44,402 --> 00:09:46,350 So I haven't actually done anything. 134 00:09:46,350 --> 00:09:51,860 Now I'm going to have another couple of nested for loops. 135 00:09:51,860 --> 00:10:03,110 I'm going to say for j equals 1,np2, and for i equals 1:np2. 136 00:10:09,450 --> 00:10:14,400 Let's f set the value of fj, which 137 00:10:14,400 --> 00:10:20,040 is a matrix of the vector, strictly speaking, column 138 00:10:20,040 --> 00:10:23,870 vector of the same length as c. 139 00:10:23,870 --> 00:10:33,050 Let's set that the i-th element of that equal to the x value, 140 00:10:33,050 --> 00:10:38,460 which is x2, because this is the second array that I'm creating, 141 00:10:38,460 --> 00:10:47,340 evaluated at parameter j to the power j minus 1. 142 00:10:47,340 --> 00:10:56,868 So that's actually evaluating x2 to the appropriate power. 143 00:10:56,868 --> 00:10:58,367 So it shouldn't have been j minus 1. 144 00:10:58,367 --> 00:11:00,140 It should have been i minus 1. 145 00:11:00,140 --> 00:11:03,140 So this is x2 j to the i minus 1. 146 00:11:06,824 --> 00:11:08,740 Again, I don't want to print all of those out, 147 00:11:08,740 --> 00:11:12,240 so I'm going to put a semicolon at the end. 148 00:11:12,240 --> 00:11:22,726 So that's calculated by fj value, 149 00:11:22,726 --> 00:11:32,015 and now I'm going to put y2 of j equal to-- Well, what 150 00:11:32,015 --> 00:11:38,090 I basically want to do is form the dot product between c, 151 00:11:38,090 --> 00:11:41,825 the coefficients, and the value of fj 152 00:11:41,825 --> 00:11:44,020 that I've just calculated. 153 00:11:44,020 --> 00:11:52,920 So I want c1 times fj1 plus c2 times fj2 and so forth. 154 00:11:52,920 --> 00:11:54,752 One way of doing that-- there are 155 00:11:54,752 --> 00:11:56,290 various different ways of doing it-- 156 00:11:56,290 --> 00:11:58,070 is to take the transpose of c. 157 00:11:58,070 --> 00:12:01,390 C remember, is a column vector. 158 00:12:01,390 --> 00:12:02,500 Let's transpose it. 159 00:12:02,500 --> 00:12:04,770 That's going to give me a row vector. 160 00:12:04,770 --> 00:12:07,980 I can then take the product of that row vector 161 00:12:07,980 --> 00:12:12,970 with the column vector fj. 162 00:12:15,890 --> 00:12:20,470 And so taking to the matrix product of a column 163 00:12:20,470 --> 00:12:23,640 vector times a row vector is equivalent to taking 164 00:12:23,640 --> 00:12:27,870 the dot product of two vectors. 165 00:12:27,870 --> 00:12:32,380 So that has produced a value y2 of j, 166 00:12:32,380 --> 00:12:37,480 which is equal to the sum of the coefficients times the values 167 00:12:37,480 --> 00:12:40,820 of x to the appropriate power. 168 00:12:40,820 --> 00:12:48,100 So that's what I want to form in that matrix. 169 00:12:48,100 --> 00:12:50,550 Again, I'm-- 170 00:12:53,085 --> 00:12:54,710 It seems as though I've made a mistake. 171 00:13:05,690 --> 00:13:10,410 And I'm not quite sure what I did wrong, 172 00:13:10,410 --> 00:13:16,630 but I need to find out. 173 00:13:16,630 --> 00:13:18,345 So let's have a look. 174 00:13:21,500 --> 00:13:26,770 It's saying that I made my mistake at line 29 of column 8. 175 00:13:26,770 --> 00:13:36,080 29 is this, and it's saying that I have a 1 by 6, 176 00:13:36,080 --> 00:13:41,290 and op2 is 40 by 1. 177 00:13:41,290 --> 00:13:46,870 I don't know why that is 40 by 1. 178 00:13:46,870 --> 00:13:55,120 I'm a bit surprised by it, because-- Ah. 179 00:13:55,120 --> 00:13:57,080 I see what I did wrong. 180 00:13:57,080 --> 00:13:59,550 This should have been endpoints. 181 00:13:59,550 --> 00:14:05,040 So what I did was, this inner for loop 182 00:14:05,040 --> 00:14:07,420 shouldn't have been over all of the j values. 183 00:14:07,420 --> 00:14:13,330 The j values refer to the length of the new vectors 184 00:14:13,330 --> 00:14:14,660 that I produced. 185 00:14:14,660 --> 00:14:18,330 It should've been over the shorter version. 186 00:14:18,330 --> 00:14:19,360 So let's see. 187 00:14:19,360 --> 00:14:22,270 Let's save that and hope that's corrected my error. 188 00:14:22,270 --> 00:14:22,770 Yes. 189 00:14:22,770 --> 00:14:24,250 Lo and behold, it has. 190 00:14:24,250 --> 00:14:27,270 I haven't actually done anything to show what the result is, 191 00:14:27,270 --> 00:14:33,030 but at least now I've got the result. 192 00:14:33,030 --> 00:14:40,330 Now what I want to do is, I want to plot these two new matrices, 193 00:14:40,330 --> 00:14:50,810 x2 and y2, and actually I'd like to overplot them 194 00:14:50,810 --> 00:14:52,530 on my old plot. 195 00:14:52,530 --> 00:14:57,935 If I just say plot like this, save this and plot it again, 196 00:14:57,935 --> 00:14:59,160 that isn't what happens. 197 00:14:59,160 --> 00:15:02,490 Instead, what happens is just the curve that I've just 198 00:15:02,490 --> 00:15:09,325 calculated is plotted, and my previous plot is wiped out. 199 00:15:12,370 --> 00:15:19,210 What I can do to prevent that happening is to say, hold on. 200 00:15:19,210 --> 00:15:27,230 Only MATLAB slash Octave would have something like hold 201 00:15:27,230 --> 00:15:28,950 on as a command. 202 00:15:28,950 --> 00:15:32,980 But anyway, hold on basically says, retain the data 203 00:15:32,980 --> 00:15:36,620 that you've already got in this plot and add some more data on. 204 00:15:36,620 --> 00:15:39,680 So let's try this. 205 00:15:39,680 --> 00:15:43,600 So now what we see is the data that I've plotted out 206 00:15:43,600 --> 00:15:48,670 in my first plot, which was up here, is held 207 00:15:48,670 --> 00:15:51,790 and the second plot is plotted. 208 00:15:51,790 --> 00:15:56,910 Notice that this fit does indeed go through the six points 209 00:15:56,910 --> 00:16:07,080 that I originally retained, and yet, 210 00:16:07,080 --> 00:16:09,070 that line actually does some funny things. 211 00:16:09,070 --> 00:16:13,780 Particularly at the end here, it goes negative. 212 00:16:13,780 --> 00:16:17,670 There's another trick that one can do in plotting, 213 00:16:17,670 --> 00:16:23,988 and one can say things like axis manual, 214 00:16:23,988 --> 00:16:29,120 and I think that, if I'm lucky, that will fix this plot. 215 00:16:29,120 --> 00:16:30,800 No, it won't. 216 00:16:30,800 --> 00:16:32,880 Never mind. 217 00:16:32,880 --> 00:16:36,780 Anyway, there's a way-- and I'll show it later-- 218 00:16:36,780 --> 00:16:41,530 to prevent the plot rescaling and the required command 219 00:16:41,530 --> 00:16:42,760 is axis manual. 220 00:16:42,760 --> 00:16:49,010 But at the moment, I think, before I do it, 221 00:16:49,010 --> 00:16:52,960 I've got to say hold off. 222 00:16:52,960 --> 00:16:57,740 And then if I re-run it, I will find that lo and behold, 223 00:16:57,740 --> 00:17:01,940 I've got a plot of the six points and the line 224 00:17:01,940 --> 00:17:03,940 that I fitted, which is a polynomial fitted 225 00:17:03,940 --> 00:17:05,839 to those lines. 226 00:17:05,839 --> 00:17:09,760 And lo and behold, it fits.