1 00:00:01,540 --> 00:00:03,910 The following content is provided under a Creative 2 00:00:03,910 --> 00:00:05,300 Commons license. 3 00:00:05,300 --> 00:00:07,510 Your support will help MIT OpenCourseWare 4 00:00:07,510 --> 00:00:11,600 continue to offer high quality educational resources for free. 5 00:00:11,600 --> 00:00:14,140 To make a donation or to view additional materials 6 00:00:14,140 --> 00:00:18,100 from hundreds of MIT courses, visit MIT OpenCourseWare 7 00:00:18,100 --> 00:00:18,980 at ocw.mit.edu. 8 00:00:23,865 --> 00:00:25,740 PROFESSOR: All right, I think we're gathered. 9 00:00:25,740 --> 00:00:28,270 So let's get going. 10 00:00:28,270 --> 00:00:30,510 We were going to DAEs today, but we 11 00:00:30,510 --> 00:00:33,812 decided to postpone for one day and do a little catch up 12 00:00:33,812 --> 00:00:35,520 on different things that you should know. 13 00:00:38,850 --> 00:00:43,110 And I thought I'd first back up a minute 14 00:00:43,110 --> 00:00:47,050 and recall what you've learned so far. 15 00:00:47,050 --> 00:00:50,970 So we started out talking a lot about linear algebra 16 00:00:50,970 --> 00:00:54,300 and, particularly, the problems of this type 17 00:00:54,300 --> 00:00:56,640 where you're given a matrix and you're given a vector, 18 00:00:56,640 --> 00:00:58,306 and you want to figure out what the x is 19 00:00:58,306 --> 00:01:00,330 that makes this equation true. 20 00:01:00,330 --> 00:01:02,010 And you've become extremely familiar 21 00:01:02,010 --> 00:01:06,340 with the convenient Matlab backslash function. 22 00:01:06,340 --> 00:01:08,330 You've probably used this a lot of times. 23 00:01:08,330 --> 00:01:11,430 And if you remember, this is one of the only problems where we 24 00:01:11,430 --> 00:01:14,310 often have a unique solution. 25 00:01:14,310 --> 00:01:16,320 We know a solution exists. 26 00:01:16,320 --> 00:01:18,630 And we can do a finite number of steps 27 00:01:18,630 --> 00:01:20,310 in order to get the solution. 28 00:01:20,310 --> 00:01:22,870 And actually, the method we have is really good. 29 00:01:22,870 --> 00:01:25,272 So it almost always gives us the solution 30 00:01:25,272 --> 00:01:27,580 and to machine accuracy. 31 00:01:27,580 --> 00:01:30,030 So this is a brilliant method. 32 00:01:30,030 --> 00:01:32,640 And so everything is based on this in our whole field, 33 00:01:32,640 --> 00:01:35,140 actually. 34 00:01:35,140 --> 00:01:36,750 So we learned this first. 35 00:01:36,750 --> 00:01:43,200 We learned it takes n cubed operations 36 00:01:43,200 --> 00:01:46,012 in order to get the solution, where 37 00:01:46,012 --> 00:01:47,220 n is the size of the vectors. 38 00:01:51,140 --> 00:01:59,470 And when we do the solution, the x we get for backslash, 39 00:01:59,470 --> 00:02:01,690 invariably, it either is a solution, 40 00:02:01,690 --> 00:02:03,537 solves this equation to machine accuracy, 41 00:02:03,537 --> 00:02:05,245 in the sense that you multiply m times x, 42 00:02:05,245 --> 00:02:08,900 you get b to as many digits as you can read. 43 00:02:08,900 --> 00:02:12,220 However, we talked about how m could be poorly conditioned 44 00:02:12,220 --> 00:02:13,870 so that there might be multiple x's 45 00:02:13,870 --> 00:02:16,880 that would solve this equation to machine accuracy. 46 00:02:16,880 --> 00:02:19,800 And so that's the problem with ill conditioning. 47 00:02:19,800 --> 00:02:21,770 So it's not that it doesn't solve the equation. 48 00:02:21,770 --> 00:02:24,077 It's that we're not really sure that's the right x. 49 00:02:24,077 --> 00:02:25,910 You could get a bunch of different x's that, 50 00:02:25,910 --> 00:02:29,890 basically, solve it as well as we can tell. 51 00:02:29,890 --> 00:02:31,986 All right, so do you guys remember this? 52 00:02:31,986 --> 00:02:34,450 All right, and so that also implies 53 00:02:34,450 --> 00:02:36,550 that if you make a very small change in b, 54 00:02:36,550 --> 00:02:39,496 you might get a gigantic change in x. 55 00:02:39,496 --> 00:02:41,995 Because even for one b, you have a quite a range of x's that 56 00:02:41,995 --> 00:02:45,930 could work, all right. 57 00:02:45,930 --> 00:02:48,570 Then we try to solve problems like this. 58 00:02:51,880 --> 00:02:55,470 Right, systems of equations, not just linear ones, 59 00:02:55,470 --> 00:02:57,400 but nonlinear ones as well. 60 00:02:57,400 --> 00:03:00,890 And what we decided to do as one of our best methods 61 00:03:00,890 --> 00:03:06,115 is Newton-Raphson, where we keep on doing j 62 00:03:06,115 --> 00:03:10,680 backslash f [INAUDIBLE] negative f. 63 00:03:10,680 --> 00:03:15,200 Right, so all we did in this method, a lot of times, 64 00:03:15,200 --> 00:03:19,410 is we just kept on calling this method over and over again. 65 00:03:19,410 --> 00:03:21,520 All right, so this took order of n cubes. 66 00:03:21,520 --> 00:03:24,730 This thing is going to take order of n, 67 00:03:24,730 --> 00:03:27,589 I don't know what you call it, steps or something. 68 00:03:27,589 --> 00:03:29,130 Maybe if you guys have a word for it. 69 00:03:29,130 --> 00:03:31,560 What do you call the iterations? 70 00:03:31,560 --> 00:03:33,770 Iterations. 71 00:03:33,770 --> 00:03:38,746 So it takes the number of iterations times n cubed. 72 00:03:38,746 --> 00:03:42,110 Now, if you really can make Newton-Raphson work, 73 00:03:42,110 --> 00:03:43,910 and you have a really great initial guess, 74 00:03:43,910 --> 00:03:45,440 the number of iterations is really tiny. 75 00:03:45,440 --> 00:03:47,606 Because Newton-Raphson is such a fantastically great 76 00:03:47,606 --> 00:03:48,440 convergence. 77 00:03:48,440 --> 00:03:51,410 However, we're usually not such good guessers. 78 00:03:51,410 --> 00:03:53,720 And so it might take a lot of iterations. 79 00:03:53,720 --> 00:03:56,120 And in fact, we use the dogleg method 80 00:03:56,120 --> 00:03:59,060 to solve this as we've discussed about where you mix up 81 00:03:59,060 --> 00:03:59,840 different methods. 82 00:03:59,840 --> 00:04:02,445 And you have a trust region and all this baloney. 83 00:04:02,445 --> 00:04:04,070 And so it might actually be quite a lot 84 00:04:04,070 --> 00:04:05,194 of iterations to get there. 85 00:04:09,170 --> 00:04:12,140 And then we went and tried to solve this problem. 86 00:04:19,459 --> 00:04:22,145 And what we decided to do for both of these was almost 87 00:04:22,145 --> 00:04:23,550 do the same thing. 88 00:04:23,550 --> 00:04:28,810 So I said, well, we can really solve this one 89 00:04:28,810 --> 00:04:30,190 by instead solving this. 90 00:04:41,510 --> 00:04:44,820 Right, that's really what we ended up doing. 91 00:04:44,820 --> 00:04:47,310 That's what happens inside fsolve is 92 00:04:47,310 --> 00:04:49,215 it converts the problem of finding 93 00:04:49,215 --> 00:04:53,490 f of x equals 0 into a minimization problem 94 00:04:53,490 --> 00:04:56,040 to try to make the norm of f as small as possible. 95 00:04:56,040 --> 00:04:59,320 Because the answer will be when the norm is zero. 96 00:04:59,320 --> 00:05:01,660 And if f of x equals 0 has a solution, 97 00:05:01,660 --> 00:05:04,690 then the global minimum here is a solution, 98 00:05:04,690 --> 00:05:06,850 as long as your initial guess is good enough 99 00:05:06,850 --> 00:05:11,320 that you're in the valley that leads to the true minimum. 100 00:05:11,320 --> 00:05:12,590 Then this will work. 101 00:05:12,590 --> 00:05:14,620 And that's actually what fsolve does. 102 00:05:14,620 --> 00:05:18,490 So these two guys get kind of coupled together, 103 00:05:18,490 --> 00:05:23,805 the minimization problem and the root-finding problem, 104 00:05:23,805 --> 00:05:25,180 or the problem of solving systems 105 00:05:25,180 --> 00:05:26,830 of non-linear equations. 106 00:05:26,830 --> 00:05:29,740 And in both of them, we're going to call the first method 107 00:05:29,740 --> 00:05:31,400 a million times. 108 00:05:31,400 --> 00:05:36,940 OK, so that's the way we're working so far. 109 00:05:36,940 --> 00:05:39,220 And actually, when you do fmincon, 110 00:05:39,220 --> 00:05:41,890 which was the next level, then you actually 111 00:05:41,890 --> 00:05:43,390 end up looking for saddle points, 112 00:05:43,390 --> 00:05:46,540 which, again, we look for as minimizing 113 00:05:46,540 --> 00:05:51,340 the norm of a gradient squared, which gets it back 114 00:05:51,340 --> 00:05:53,254 into this kind of form. 115 00:05:53,254 --> 00:05:54,670 And so they're all the same thing. 116 00:05:57,470 --> 00:06:01,270 And you can, again, work out estimates 117 00:06:01,270 --> 00:06:03,160 of how much effort it takes. 118 00:06:03,160 --> 00:06:06,430 And it's something times n cubed, Typically, 119 00:06:06,430 --> 00:06:09,500 if most your work is this. 120 00:06:09,500 --> 00:06:12,880 Now, I should make a comment that, in practice, that's 121 00:06:12,880 --> 00:06:14,060 not always true. 122 00:06:14,060 --> 00:06:17,320 So in practice, sometimes, although the highest order term 123 00:06:17,320 --> 00:06:20,820 is n cubed, the linear algebra programs that do this 124 00:06:20,820 --> 00:06:24,370 are so great that the prefactor is tiny. 125 00:06:24,370 --> 00:06:26,910 So in order to get the limit where 126 00:06:26,910 --> 00:06:29,280 the n cubed term was actually dominate the CPU time, 127 00:06:29,280 --> 00:06:31,230 you had to have a gigantic n. 128 00:06:31,230 --> 00:06:34,170 And in fact, the limit, it's often the effort 129 00:06:34,170 --> 00:06:35,770 to construct the matrix j. 130 00:06:35,770 --> 00:06:37,920 And so it really goes as n squared. 131 00:06:37,920 --> 00:06:40,589 So although this nominally is like n cubed, in practice, 132 00:06:40,589 --> 00:06:41,880 a lot of times, it's n squared. 133 00:06:44,610 --> 00:06:46,680 And I say it for these other methods. 134 00:06:46,680 --> 00:06:48,540 And then a lot of them get into how good 135 00:06:48,540 --> 00:06:49,620 the initial guesses are. 136 00:06:49,620 --> 00:06:53,467 Because the number of iterations can get crazy. 137 00:06:53,467 --> 00:06:55,050 And in fact, maybe, it never converges 138 00:06:55,050 --> 00:06:56,820 if your initial guess is bad. 139 00:06:56,820 --> 00:06:58,904 So that can be a really critical thing. 140 00:06:58,904 --> 00:06:59,820 [INAUDIBLE] is better. 141 00:06:59,820 --> 00:07:01,380 Because you can't change n I mean, 142 00:07:01,380 --> 00:07:04,469 how many unknowns you've got is how many unknowns you've got. 143 00:07:04,469 --> 00:07:06,510 But you can change the number iterations by, say, 144 00:07:06,510 --> 00:07:08,410 more clever initial guessing. 145 00:07:08,410 --> 00:07:11,550 And that was the whole concept about continuation and homotopy 146 00:07:11,550 --> 00:07:13,110 and all this kind of stuff was trying 147 00:07:13,110 --> 00:07:15,210 to help reduce this number. 148 00:07:18,440 --> 00:07:23,490 Last time, we talked about implicitly solving ODEs. 149 00:07:23,490 --> 00:07:26,740 And so we had an update formula. 150 00:07:26,740 --> 00:07:28,320 So we were solving an ODE. 151 00:07:34,890 --> 00:07:37,122 The initial value problem, y of-- 152 00:07:41,010 --> 00:07:44,580 Right, and we said that, first, we went 153 00:07:44,580 --> 00:07:46,440 through the explicit methods. 154 00:07:46,440 --> 00:07:48,687 And they were perfectly straightforward to do. 155 00:07:48,687 --> 00:07:51,270 And you don't have to solve any systems of nonlinear equations 156 00:07:51,270 --> 00:07:52,020 to do them. 157 00:07:52,020 --> 00:07:53,500 And so that was great. 158 00:07:53,500 --> 00:07:55,375 But then we pointed out that, a lot of times, 159 00:07:55,375 --> 00:07:56,980 they're numerically unstable. 160 00:07:56,980 --> 00:08:01,050 And so then you have to instead use implicit methods. 161 00:08:01,050 --> 00:08:03,960 But the implicit methods generally work this way. 162 00:08:03,960 --> 00:08:05,210 They have your new value of y. 163 00:08:05,210 --> 00:08:09,660 What you're trying to compute is your current value of y 164 00:08:09,660 --> 00:08:13,930 plus delta t sums g. 165 00:08:13,930 --> 00:08:18,305 It depends on the size of delta t and y new. 166 00:08:18,305 --> 00:08:19,518 Maybe y old also. 167 00:08:24,000 --> 00:08:26,870 And these g's are things like the Crank-Nicolson, which 168 00:08:26,870 --> 00:08:28,590 you're going to work on in the homework, 169 00:08:28,590 --> 00:08:31,364 and we show it in class too. 170 00:08:31,364 --> 00:08:33,030 And so there's different update formulas 171 00:08:33,030 --> 00:08:36,059 of people suggesting how to do this. 172 00:08:36,059 --> 00:08:39,870 The trick of this is that the y new is inside the function. 173 00:08:39,870 --> 00:08:45,470 So therefore, this is really another problem 174 00:08:45,470 --> 00:08:48,160 like we had before. 175 00:08:48,160 --> 00:08:49,820 You can really rewrite this that way 176 00:08:49,820 --> 00:08:51,860 by just putting this on the one side. 177 00:08:51,860 --> 00:08:53,455 Put them all on one side of zero. 178 00:08:53,455 --> 00:08:55,455 And now, I'm trying to figure out what y new is. 179 00:08:55,455 --> 00:08:58,350 And I'm going to solve that using this method. 180 00:08:58,350 --> 00:09:01,800 And I'm going to solve that by using this method. 181 00:09:01,800 --> 00:09:05,960 All right, so now, how is the scaling going to go? 182 00:09:05,960 --> 00:09:12,140 At each time step, I have to solve one of these guys. 183 00:09:12,140 --> 00:09:15,199 How much effort did that take? 184 00:09:15,199 --> 00:09:15,990 We did that, right. 185 00:09:15,990 --> 00:09:19,020 Here it is, n [? after ?] times n squared. 186 00:09:19,020 --> 00:09:21,050 And then how many times do I have to do it? 187 00:09:23,894 --> 00:09:25,728 AUDIENCE: It's h times 7 times [INAUDIBLE].. 188 00:09:25,728 --> 00:09:27,852 PROFESSOR: Right it depends on how many times steps 189 00:09:27,852 --> 00:09:29,320 I have [INAUDIBLE] my delta t. 190 00:09:29,320 --> 00:09:34,976 So it's something like t final minus t 0 over delta t. 191 00:09:34,976 --> 00:09:36,350 That's the number of times steps. 192 00:09:36,350 --> 00:09:38,230 But I have constant delta t. 193 00:09:38,230 --> 00:09:40,960 Times the number of iterations that it 194 00:09:40,960 --> 00:09:44,920 takes to do the nonlinear solve times n 195 00:09:44,920 --> 00:09:50,940 squared, which is a cost of forming the Jacobian matrix. 196 00:09:50,940 --> 00:09:53,084 Now, this can be pretty large. 197 00:09:53,084 --> 00:09:55,000 So as we talked about, this might be as big 10 198 00:09:55,000 --> 00:09:59,914 to the eighth in a bad case. 199 00:09:59,914 --> 00:10:02,330 Certainly be less than 10 to the eighth because otherwise, 200 00:10:02,330 --> 00:10:03,621 you'll have numerical problems. 201 00:10:03,621 --> 00:10:05,600 Maybe, it's really going to be 10 202 00:10:05,600 --> 00:10:08,910 to the sixth for a real problem. 203 00:10:08,910 --> 00:10:10,430 So you have a million time steps. 204 00:10:10,430 --> 00:10:14,740 At each times step, you have to solve the number of iter-- 205 00:10:14,740 --> 00:10:16,310 the solve costs this much. 206 00:10:16,310 --> 00:10:17,780 Number of iterations might be what? 207 00:10:17,780 --> 00:10:19,280 How many iterations do you guys take 208 00:10:19,280 --> 00:10:21,634 to solve nonlinear equations? 209 00:10:21,634 --> 00:10:22,800 You solve a lot of them now. 210 00:10:22,800 --> 00:10:26,245 How many iterations does it take? 211 00:10:26,245 --> 00:10:27,134 AUDIENCE: 10. 212 00:10:27,134 --> 00:10:28,800 PROFESSOR: 10, does that 10 sound right? 213 00:10:32,084 --> 00:10:33,250 AUDIENCE: Yeah, [INAUDIBLE]. 214 00:10:33,250 --> 00:10:35,260 PROFESSOR: You've got an initial guess 10. 215 00:10:35,260 --> 00:10:37,177 If you have a medium initial guess, maybe 100. 216 00:10:37,177 --> 00:10:39,635 It it's more than 100, you're not going to converge, right? 217 00:10:39,635 --> 00:10:40,300 AUDIENCE: Yeah. 218 00:10:40,300 --> 00:10:41,133 PROFESSOR: Yeah, OK. 219 00:10:41,133 --> 00:10:43,025 So this is order of 10 maybe. 220 00:10:43,025 --> 00:10:45,700 So the million 10 and then n squared 221 00:10:45,700 --> 00:10:47,710 just depends on how big your problem is. 222 00:10:47,710 --> 00:10:52,060 So in my group, the problems we typically do, 223 00:10:52,060 --> 00:10:54,070 n is approximately 200. 224 00:10:54,070 --> 00:10:57,610 So this is like 200 squared. 225 00:10:57,610 --> 00:11:00,345 And then you can see that the number of iterations, and this 226 00:11:00,345 --> 00:11:01,845 is maybe a little bit overestimated. 227 00:11:01,845 --> 00:11:05,330 Maybe I can get away with 10 to the fifth, something like that. 228 00:11:05,330 --> 00:11:08,350 But you can see it starts to get pretty big. 229 00:11:08,350 --> 00:11:11,980 But for writing it, as you guys are writing software, 230 00:11:11,980 --> 00:11:13,090 it's not that big a deal. 231 00:11:13,090 --> 00:11:17,180 So I can assign it as part 1 of a three-part homework problem. 232 00:11:17,180 --> 00:11:18,270 It's the right ODE solver. 233 00:11:18,270 --> 00:11:20,346 It solves implicit equations. 234 00:11:20,346 --> 00:11:21,220 And you can write it. 235 00:11:21,220 --> 00:11:24,330 And the reason you can is because you've already written 236 00:11:24,330 --> 00:11:25,330 methods that solve this. 237 00:11:25,330 --> 00:11:26,620 Actually, Matlab did it for you. 238 00:11:26,620 --> 00:11:27,703 But you write it yourself. 239 00:11:27,703 --> 00:11:29,380 I assigned to the [? 10:10 ?] students. 240 00:11:29,380 --> 00:11:32,890 So you write the [? backsub ?] solution 241 00:11:32,890 --> 00:11:34,780 to solve the equations. 242 00:11:34,780 --> 00:11:38,340 And we we've already wrestled this. 243 00:11:38,340 --> 00:11:40,724 And I think you guys wrote one of these yourself. 244 00:11:40,724 --> 00:11:42,140 And also, we can use fsolve, which 245 00:11:42,140 --> 00:11:44,139 is just another little bit better implementation 246 00:11:44,139 --> 00:11:45,040 of the same thing. 247 00:11:45,040 --> 00:11:49,000 And then now, you write your next level. 248 00:11:49,000 --> 00:11:52,180 So you can see you're doing a composite process, where 249 00:11:52,180 --> 00:11:54,234 you're doing something at the top level. 250 00:11:54,234 --> 00:11:55,900 It's calling something at a lower level. 251 00:11:55,900 --> 00:11:57,850 That's calling something at a lower level. 252 00:11:57,850 --> 00:11:59,200 And then what you really have to be concerned 253 00:11:59,200 --> 00:12:00,520 about if you're worried about CPU time 254 00:12:00,520 --> 00:12:03,061 or whether the thing is going to solve before the homework is 255 00:12:03,061 --> 00:12:07,062 due, is how big does this get. 256 00:12:07,062 --> 00:12:08,770 How many operations are you really asking 257 00:12:08,770 --> 00:12:09,780 the computer to do? 258 00:12:09,780 --> 00:12:11,815 Now, remember the computer is fast. 259 00:12:11,815 --> 00:12:13,690 Right, so it can do a gigahertz or something. 260 00:12:13,690 --> 00:12:16,390 So you can do like 10 to the eighth somethings per second. 261 00:12:19,320 --> 00:12:21,050 And what should I say about this? 262 00:12:21,050 --> 00:12:27,290 I have not included, but I should that it's-- 263 00:12:27,290 --> 00:12:30,415 well, I'm assuming here that the function evaluation 264 00:12:30,415 --> 00:12:34,080 is sort of order n of how many variables they have. 265 00:12:34,080 --> 00:12:35,980 OK, sometimes, function evaluations 266 00:12:35,980 --> 00:12:38,177 are much more difficult than that, 267 00:12:38,177 --> 00:12:40,510 primarily because, sometimes, they're actually including 268 00:12:40,510 --> 00:12:42,125 further nested stuff inside. 269 00:12:42,125 --> 00:12:44,500 So your f is actually coming from some really complicated 270 00:12:44,500 --> 00:12:46,090 calculation itself. 271 00:12:46,090 --> 00:12:48,190 And then you have another big giant factor 272 00:12:48,190 --> 00:12:50,949 of how much it costs to do f. 273 00:12:50,949 --> 00:12:53,240 But on this right here, I'm just assuming it goes as n. 274 00:12:57,900 --> 00:12:59,680 So actually, let's just multiply this out. 275 00:12:59,680 --> 00:13:02,330 So this is what, 10 to the fourth, 10 to the 10th. 276 00:13:02,330 --> 00:13:04,320 So this is 10 to the 10th. 277 00:13:04,320 --> 00:13:08,550 10 to the 10th kind of flops. 278 00:13:08,550 --> 00:13:12,630 But your computer does 4 times 10 to the ninth per second. 279 00:13:12,630 --> 00:13:14,231 Is that right? 280 00:13:14,231 --> 00:13:16,240 Is that [? what they ?] told us, four gigahertz? 281 00:13:16,240 --> 00:13:19,570 Yeah, so is this really not bad with just a few seconds 282 00:13:19,570 --> 00:13:20,500 on your computer. 283 00:13:20,500 --> 00:13:22,570 So that's why you can do these homework problems, 284 00:13:22,570 --> 00:13:23,960 and you still get them done. 285 00:13:23,960 --> 00:13:25,330 And if you have even [? porous ?] limitations, 286 00:13:25,330 --> 00:13:26,440 it takes 1,000 times longer. 287 00:13:26,440 --> 00:13:28,856 Still, 1,000 [INAUDIBLE] take 20 minutes on your computer. 288 00:13:31,340 --> 00:13:34,640 So you're OK, but you can see we're starting to get up there. 289 00:13:34,640 --> 00:13:37,580 And it's not that hard to really make it a big problem 290 00:13:37,580 --> 00:13:39,320 if you do something wrong. 291 00:13:44,010 --> 00:13:51,740 And just a notation thing, a lot of these things we're doing 292 00:13:51,740 --> 00:13:58,220 are writing software that takes as input not just numbers 293 00:13:58,220 --> 00:13:59,889 or vectors. 294 00:13:59,889 --> 00:14:01,680 So those are things that we call functions. 295 00:14:01,680 --> 00:14:03,890 Right, things that take a number of vectors and input 296 00:14:03,890 --> 00:14:07,235 and give you, say, a function, a number of vectors and output. 297 00:14:07,235 --> 00:14:09,860 But oftentimes, we actually want to take the names of functions 298 00:14:09,860 --> 00:14:11,150 as inputs. 299 00:14:11,150 --> 00:14:13,880 So for example, the root-finding problem, 300 00:14:13,880 --> 00:14:16,330 it takes as input what is f. 301 00:14:16,330 --> 00:14:19,570 Right, you have to tell it the name of a function for fsolve 302 00:14:19,570 --> 00:14:20,890 to be able to call it. 303 00:14:20,890 --> 00:14:23,830 So fsolve is an object that takes as input 304 00:14:23,830 --> 00:14:26,290 the name of a function. 305 00:14:26,290 --> 00:14:35,600 And things that do that are called functionals 306 00:14:35,600 --> 00:14:43,870 So for example, this one is x is equal to this root-finding 307 00:14:43,870 --> 00:14:48,295 functional, fsolve. 308 00:14:48,295 --> 00:14:50,380 And I'm going to write with square brackets. 309 00:14:50,380 --> 00:14:51,742 It depends on f. 310 00:14:51,742 --> 00:14:52,450 That's it, right? 311 00:14:52,450 --> 00:14:55,630 [INAUDIBLE] f. 312 00:14:55,630 --> 00:14:57,550 So you give it the function f. 313 00:14:57,550 --> 00:14:59,060 And if fsolve's good, it should be 314 00:14:59,060 --> 00:15:02,080 able to tell you what the x is. 315 00:15:02,080 --> 00:15:04,450 Now, in reality, we help it out. 316 00:15:04,450 --> 00:15:07,970 So we give it x 0 [? 2 ?] just give it a better chance. 317 00:15:11,034 --> 00:15:13,700 I write the square brackets just indicate that one of the things 318 00:15:13,700 --> 00:15:18,110 inside is a function, not just numbers. 319 00:15:18,110 --> 00:15:20,010 We do this a lot. 320 00:15:20,010 --> 00:15:24,050 So you guys have done this since you were kids. 321 00:15:24,050 --> 00:15:29,880 So derivatives Yeah, ddx of f. 322 00:15:33,720 --> 00:15:38,224 So we have, I don't know, for example, grad of f. 323 00:15:41,160 --> 00:15:44,510 That's a functional, the grad. 324 00:15:44,510 --> 00:15:47,740 The symbol means the operator, the functional, 325 00:15:47,740 --> 00:15:49,090 it operates on functions. 326 00:15:49,090 --> 00:15:50,980 You could write a general gradient function 327 00:15:50,980 --> 00:15:54,472 that takes any scalar valued function and computes 328 00:15:54,472 --> 00:15:56,180 the gradients with respect to the inputs. 329 00:15:58,654 --> 00:16:00,820 This week, we give you one that computes a Jacobian. 330 00:16:00,820 --> 00:16:02,230 Give it any f of x. 331 00:16:02,230 --> 00:16:05,710 It gives you the Jacobian of f with respect to x. 332 00:16:05,710 --> 00:16:07,270 So that's a functional. 333 00:16:07,270 --> 00:16:09,001 So anyway, we do a lot with functionals. 334 00:16:09,001 --> 00:16:10,500 You've been using them all the time. 335 00:16:10,500 --> 00:16:12,730 I'm just trying get you to think about abstractly. 336 00:16:12,730 --> 00:16:15,577 And then there's a whole math about functionals. 337 00:16:15,577 --> 00:16:17,910 And there's a whole thing called calculus of variations. 338 00:16:17,910 --> 00:16:20,080 It's all about, if the objects you're worrying about 339 00:16:20,080 --> 00:16:22,300 are functionals, not functions, then you 340 00:16:22,300 --> 00:16:24,739 have another whole bunch of theorems and stuff 341 00:16:24,739 --> 00:16:26,530 you can do with that, which you'll probably 342 00:16:26,530 --> 00:16:28,570 end up doing before you get out of this place. 343 00:16:28,570 --> 00:16:30,195 But I'm not going to do that right now. 344 00:16:32,740 --> 00:16:35,500 I would comment that one kind of functional that's 345 00:16:35,500 --> 00:16:41,800 particularly important to a lot of people in your research 346 00:16:41,800 --> 00:16:49,120 is the fact that the energy of any molecule 347 00:16:49,120 --> 00:16:51,640 is a functional of the electron density. 348 00:16:56,127 --> 00:16:57,710 And this is the basis of what's called 349 00:16:57,710 --> 00:16:58,793 density functional theory. 350 00:17:09,160 --> 00:17:11,888 And so this is why a lot of people, you and many of you 351 00:17:11,888 --> 00:17:14,304 included, are going to end up using computer programs that 352 00:17:14,304 --> 00:17:17,752 are trying to find what the electron density shape is, 353 00:17:17,752 --> 00:17:19,460 things like molecular orbitals and stuff. 354 00:17:19,460 --> 00:17:21,490 You've probably heard these before. 355 00:17:21,490 --> 00:17:23,560 And then there's the functional that gives 356 00:17:23,560 --> 00:17:24,215 the energy of the molecule. 357 00:17:24,215 --> 00:17:26,006 And what you care about, mostly, is energy, 358 00:17:26,006 --> 00:17:28,510 because that connects to all your thermo. 359 00:17:28,510 --> 00:17:31,402 You're doing [? 1040, ?] it's all about the energies. 360 00:17:31,402 --> 00:17:32,860 And once you know all the energies, 361 00:17:32,860 --> 00:17:34,443 you can actually compute the entropies 362 00:17:34,443 --> 00:17:35,730 and all kinds of stuff too. 363 00:17:35,730 --> 00:17:40,972 And so this is one application that's really important. 364 00:17:40,972 --> 00:17:42,430 It was a little surprising this was 365 00:17:42,430 --> 00:17:44,900 true and able to be proven to be true. 366 00:17:44,900 --> 00:17:47,400 And so the guy who did that got the Nobel Prize-- his name 367 00:17:47,400 --> 00:17:48,850 was Cohen-- a few years ago. 368 00:17:51,540 --> 00:17:53,365 Actually, in liquid phase stuff, you 369 00:17:53,365 --> 00:17:55,180 work with Professor Blankschtein. 370 00:17:55,180 --> 00:17:58,105 They also have a version of this for properties of solution. 371 00:17:58,105 --> 00:17:59,980 And then where it's not the electron density, 372 00:17:59,980 --> 00:18:03,100 but it's actually the density of the material in the solution 373 00:18:03,100 --> 00:18:03,600 phase. 374 00:18:03,600 --> 00:18:04,891 Maybe Professor [? Swan ?] too. 375 00:18:04,891 --> 00:18:05,620 You guys do that? 376 00:18:05,620 --> 00:18:07,090 Yeah, so if you work with Professor [? Swan, ?] you 377 00:18:07,090 --> 00:18:08,390 might learn about that too. 378 00:18:08,390 --> 00:18:10,110 So that's one where the actual functional 379 00:18:10,110 --> 00:18:14,486 is like a big focus of the entire project. 380 00:18:14,486 --> 00:18:16,110 And the key is that this is a function. 381 00:18:16,110 --> 00:18:18,100 The density is a function of the position. 382 00:18:18,100 --> 00:18:20,650 And then the functional converts that into a number. 383 00:18:26,025 --> 00:18:27,310 That's all page one here. 384 00:18:29,850 --> 00:18:33,280 Now, I'm going to sort of change topics, but it's related. 385 00:18:33,280 --> 00:18:37,360 It's about the connections between numerical integration, 386 00:18:37,360 --> 00:18:39,969 interpolation, and basis functions. 387 00:18:39,969 --> 00:18:42,010 These things are all intimately tied up together. 388 00:18:42,010 --> 00:18:43,540 So I'll try to show you. 389 00:18:43,540 --> 00:18:46,240 So we try to solve a numerical integral. 390 00:19:00,320 --> 00:19:03,170 And then let's look at what we're really doing. 391 00:19:03,170 --> 00:19:04,970 When we're doing this, we typically 392 00:19:04,970 --> 00:19:07,190 only do a few function evaluations. 393 00:19:07,190 --> 00:19:12,260 So we have a few points where we know the numbers f of tn. 394 00:19:12,260 --> 00:19:14,570 We've picked a few t's. 395 00:19:14,570 --> 00:19:16,190 We evaluated f. 396 00:19:16,190 --> 00:19:18,909 And from those few points, that small number 397 00:19:18,909 --> 00:19:21,200 of function evaluations, we want to get a good estimate 398 00:19:21,200 --> 00:19:22,712 of this whole integral. 399 00:19:22,712 --> 00:19:24,170 Now, it's a little goofy, actually, 400 00:19:24,170 --> 00:19:25,970 if you think about what we're trying to do. 401 00:19:25,970 --> 00:19:28,386 So suppose you're trying to do an integral from negative 1 402 00:19:28,386 --> 00:19:29,134 to 1. 403 00:19:29,134 --> 00:19:31,550 We have some function that has a point here, a point here, 404 00:19:31,550 --> 00:19:34,070 a point here. 405 00:19:34,070 --> 00:19:37,010 And just for concreteness, let's say this is 1 over 26. 406 00:19:37,010 --> 00:19:38,000 And this is 1. 407 00:19:38,000 --> 00:19:39,724 And this is 1 over 26. 408 00:19:39,724 --> 00:19:41,390 And I have a particular function in mind 409 00:19:41,390 --> 00:19:43,700 that has those three points. 410 00:19:43,700 --> 00:19:45,620 Now, the first thing to keep in mind 411 00:19:45,620 --> 00:19:48,211 is there's an infinite number of functions that 412 00:19:48,211 --> 00:19:49,460 go through those three points. 413 00:19:51,990 --> 00:19:55,850 Now, what we've been doing is we've 414 00:19:55,850 --> 00:20:01,840 been fitting some function to these points 415 00:20:01,840 --> 00:20:05,365 and then integrating that function that we fitted. 416 00:20:05,365 --> 00:20:07,740 So when you're in high school, you did the rectangle rule 417 00:20:07,740 --> 00:20:15,720 and what the function was used was this and this. 418 00:20:15,720 --> 00:20:17,820 And you knew the integral of this function 419 00:20:17,820 --> 00:20:20,220 and the integral of this function. 420 00:20:20,220 --> 00:20:21,760 And you added them up. 421 00:20:21,760 --> 00:20:24,940 And you decided from the rectangle rule 422 00:20:24,940 --> 00:20:28,019 that the integral was equal to 1 and 1/26. 423 00:20:33,052 --> 00:20:34,760 But then some other people in high school 424 00:20:34,760 --> 00:20:35,880 are a little smarter. 425 00:20:35,880 --> 00:20:38,262 And they said it's not so smart [INAUDIBLE].. 426 00:20:38,262 --> 00:20:39,970 This thing was symmetrical to begin with. 427 00:20:39,970 --> 00:20:41,595 And now, I'm fitting it with a function 428 00:20:41,595 --> 00:20:42,890 that's wildly unsymmetrical. 429 00:20:42,890 --> 00:20:45,729 So let's instead not evaluate those points. 430 00:20:45,729 --> 00:20:47,270 Let's instead evaluate the midpoints. 431 00:20:47,270 --> 00:20:49,186 And maybe we find out that the midpoint values 432 00:20:49,186 --> 00:20:51,060 are here and here. 433 00:20:51,060 --> 00:20:58,569 And so I would assume that the function is just a flat line. 434 00:20:58,569 --> 00:21:00,610 That looks not so good either, because it doesn't 435 00:21:00,610 --> 00:21:01,190 get through these points. 436 00:21:01,190 --> 00:21:02,690 But anyway, for the two points I knew, 437 00:21:02,690 --> 00:21:03,770 if I only used those two points, that 438 00:21:03,770 --> 00:21:05,240 would be the midpoint method. 439 00:21:05,240 --> 00:21:07,800 And I would get the I midpoint method. 440 00:21:07,800 --> 00:21:08,870 And I'd get something. 441 00:21:08,870 --> 00:21:11,390 And I did it for this particular function I have in mind. 442 00:21:11,390 --> 00:21:12,590 And it was 8 over 29. 443 00:21:15,149 --> 00:21:16,690 And then some other people say, well, 444 00:21:16,690 --> 00:21:18,030 let's use the trapezoid rule. 445 00:21:18,030 --> 00:21:19,440 So let's approximate the function 446 00:21:19,440 --> 00:21:25,310 is this over the interval, which looks kind of sensible. 447 00:21:25,310 --> 00:21:26,810 And if you do that, it turns out you 448 00:21:26,810 --> 00:21:29,101 get the same numbers as you do with the rectangle rule. 449 00:21:29,101 --> 00:21:30,447 So I trapezoid. 450 00:21:34,830 --> 00:21:38,915 Now, let me just as a warning, if you just do the rectangle 451 00:21:38,915 --> 00:21:40,290 rule and then the trapezoid rule, 452 00:21:40,290 --> 00:21:41,915 and you got the same number both times, 453 00:21:41,915 --> 00:21:44,082 you might be tempted to conclude you have converged. 454 00:21:44,082 --> 00:21:45,873 And you have the true solution, because you 455 00:21:45,873 --> 00:21:47,730 did two methods of really different orders, 456 00:21:47,730 --> 00:21:49,524 and you got exactly the same number. 457 00:21:49,524 --> 00:21:50,940 And you're super happy, and you're 458 00:21:50,940 --> 00:21:53,622 ready to publish your paper. 459 00:21:53,622 --> 00:21:55,830 But of course, the function might not look like this. 460 00:21:55,830 --> 00:22:00,510 So the particular function that I always had in mind was f of t 461 00:22:00,510 --> 00:22:05,520 is equal to 1 over 1 plus 25 t squared. 462 00:22:05,520 --> 00:22:07,460 OK, that was just the function I had in mind. 463 00:22:07,460 --> 00:22:16,616 And the way that function looks like is like that. 464 00:22:16,616 --> 00:22:18,240 And the real integral of that function, 465 00:22:18,240 --> 00:22:20,160 which you guys if you remember how do it from high school, 466 00:22:20,160 --> 00:22:21,950 you can actually evaluate that integral. 467 00:22:21,950 --> 00:22:26,520 It's nothing to do with either 8/29 or 1 in 126. 468 00:22:26,520 --> 00:22:28,020 It's something completely different. 469 00:22:28,020 --> 00:22:30,530 Now, those of you who went to really advanced high schools, 470 00:22:30,530 --> 00:22:32,135 you also learned Simpson's rule too. 471 00:22:32,135 --> 00:22:35,150 And Simpson's rule says let's fit it to a parabola. 472 00:22:35,150 --> 00:22:38,930 And so that would give you something like this. 473 00:22:38,930 --> 00:22:41,870 And if you do Simpson's rule, I ended up getting this 474 00:22:41,870 --> 00:22:43,330 if I did the arithmetic right. 475 00:22:46,390 --> 00:22:51,420 1 and 1/3 plus 2/78. 476 00:22:51,420 --> 00:22:52,520 Who knew? 477 00:22:52,520 --> 00:22:55,070 OK, which also has little relation 478 00:22:55,070 --> 00:22:57,750 to the true value of the integral. 479 00:22:57,750 --> 00:23:00,000 Now, if you're smart, and you did trapezoid and then I 480 00:23:00,000 --> 00:23:02,790 it to Simpson's, then you say, oh my goodness, these 481 00:23:02,790 --> 00:23:03,990 are not really that similar. 482 00:23:03,990 --> 00:23:06,448 And so I'm not converged, and I'd better do something else. 483 00:23:08,324 --> 00:23:11,800 So I just want you to be aware when 484 00:23:11,800 --> 00:23:13,900 we do these rules, what we're actually doing 485 00:23:13,900 --> 00:23:19,360 is we're fitting our points to some continuous function 486 00:23:19,360 --> 00:23:21,040 that we know how to integrate. 487 00:23:21,040 --> 00:23:22,623 And then we do the analytical integral 488 00:23:22,623 --> 00:23:25,570 of that function, which is a fine thing to do. 489 00:23:25,570 --> 00:23:28,814 But just be aware that it's like extrapolation. 490 00:23:28,814 --> 00:23:30,730 It's like we're imagining what the function is 491 00:23:30,730 --> 00:23:34,610 doing between these points, when we actually haven't tested it. 492 00:23:34,610 --> 00:23:35,305 Now, we could. 493 00:23:35,305 --> 00:23:36,596 We could just call more points. 494 00:23:36,596 --> 00:23:37,388 We could calculate. 495 00:23:37,388 --> 00:23:39,179 But that would cost us more function evals. 496 00:23:39,179 --> 00:23:41,127 And we're trying to save our computer time. 497 00:23:41,127 --> 00:23:42,710 So we want to finish off before the TG 498 00:23:42,710 --> 00:23:43,900 so we can go buy some beer. 499 00:23:43,900 --> 00:23:51,370 And so we have to take some choices about what to do. 500 00:23:51,370 --> 00:23:53,914 And we can't do an infinite number of points for sure. 501 00:23:53,914 --> 00:23:56,080 All right, so in all these methods, what we're doing 502 00:23:56,080 --> 00:23:59,110 is we're approximating our true f 503 00:23:59,110 --> 00:24:04,300 with some approximate f, which is equal to a sum of some basis 504 00:24:04,300 --> 00:24:12,657 functions like this. 505 00:24:12,657 --> 00:24:13,990 So we pick some basis functions. 506 00:24:13,990 --> 00:24:15,820 We know how to integrate. 507 00:24:15,820 --> 00:24:18,730 We fit, somehow determine what these c's are. 508 00:24:18,730 --> 00:24:20,890 And then we decide that the integral 509 00:24:20,890 --> 00:24:27,260 is equal to the integral of f tilde, not 510 00:24:27,260 --> 00:24:29,360 the actual original f. 511 00:24:29,360 --> 00:24:31,690 And then that's actually going to be the sum 512 00:24:31,690 --> 00:24:33,320 of ck times the integral. 513 00:24:37,047 --> 00:24:38,630 And we carefully chose basis functions 514 00:24:38,630 --> 00:24:42,260 that we don't how to do this integral; analytically. 515 00:24:42,260 --> 00:24:44,030 So we just plug that in, and we just 516 00:24:44,030 --> 00:24:48,040 have to figure out that the c's are to get this to work. 517 00:24:48,040 --> 00:24:49,750 OK, so that's what you really did 518 00:24:49,750 --> 00:24:51,290 when you're doing this trapezoid rule and Simpson's rule, 519 00:24:51,290 --> 00:24:52,630 and all that kind of stuff. 520 00:24:52,630 --> 00:24:54,671 And it's just different choices of what the basis 521 00:24:54,671 --> 00:24:56,200 functions were that you used. 522 00:24:59,495 --> 00:25:01,120 Now, there's a lot of choices about how 523 00:25:01,120 --> 00:25:04,720 to do that kind of approximation to make f tildes that 524 00:25:04,720 --> 00:25:09,320 are sort of like f's with some limited amount of information. 525 00:25:09,320 --> 00:25:12,140 And it generally has to do with how many k's we have. 526 00:25:12,140 --> 00:25:15,845 So if k is less than n, then it's really an approximation. 527 00:25:15,845 --> 00:25:17,380 Because there's no way that we would 528 00:25:17,380 --> 00:25:20,690 expect that with just a few parameters, 529 00:25:20,690 --> 00:25:22,780 we could get a lot f of tn. 530 00:25:22,780 --> 00:25:30,870 So we would expect that f tilde of tn is not equal to f of tn, 531 00:25:30,870 --> 00:25:32,710 at least not at all the points. 532 00:25:32,710 --> 00:25:34,030 Because if n is-- 533 00:25:34,030 --> 00:25:36,670 if we have a lot of points and we only have a few parameters, 534 00:25:36,670 --> 00:25:40,900 we don't expect to be able to fit it perfectly. 535 00:25:40,900 --> 00:25:43,410 So we know this is what we expect. 536 00:25:43,410 --> 00:25:45,577 It might happen if, by some luck, 537 00:25:45,577 --> 00:25:47,410 we chose a basis function which was actually 538 00:25:47,410 --> 00:25:50,510 the true function f, then it would actually fit perfectly. 539 00:25:50,510 --> 00:25:54,100 But most of the time, not. 540 00:25:54,100 --> 00:25:57,740 By Murphy's Law, in fact, it's going to be completely wrong. 541 00:25:57,740 --> 00:26:00,440 But this is not such a bad idea. 542 00:26:00,440 --> 00:26:03,020 We may come back to this later. 543 00:26:03,020 --> 00:26:06,890 Another possibility is to choose k equals to n. 544 00:26:06,890 --> 00:26:21,390 And then usually, we can force f tilde of tn to equal f of tn. 545 00:26:21,390 --> 00:26:24,310 So at the points that we have, if we have enough parameters 546 00:26:24,310 --> 00:26:26,860 equal to the number of points, then often, 547 00:26:26,860 --> 00:26:27,880 we can get an exact fit. 548 00:26:27,880 --> 00:26:31,510 We can force the function to go through all the points. 549 00:26:31,510 --> 00:26:35,771 And when we do this, this is called interpolation. 550 00:26:43,480 --> 00:26:46,620 Now, when should you do the one or the other? 551 00:26:46,620 --> 00:26:48,310 When should I make k as big as n? 552 00:26:48,310 --> 00:26:50,170 Or when should I use k less than n? 553 00:26:50,170 --> 00:26:52,030 A really crucial thing is whether you really 554 00:26:52,030 --> 00:26:53,820 believe your points. 555 00:26:53,820 --> 00:26:56,050 If you really believe you know f of tn 556 00:26:56,050 --> 00:26:58,780 to a very, very high number of significant figures, 557 00:26:58,780 --> 00:26:59,872 then it make sense. 558 00:26:59,872 --> 00:27:01,330 You know that's true, so you should 559 00:27:01,330 --> 00:27:05,170 force your fitting function to actually match that. 560 00:27:05,170 --> 00:27:08,650 However, if you get the f of tn, for example, from an experiment 561 00:27:08,650 --> 00:27:11,890 or from a stochastic simulation that has some noise in it, 562 00:27:11,890 --> 00:27:14,410 then you might not really want to do this. 563 00:27:14,410 --> 00:27:16,630 Because you'll be fitting the noise 564 00:27:16,630 --> 00:27:18,530 in the experiment as well. 565 00:27:18,530 --> 00:27:24,280 And then you might want to use a smaller number of parameters 566 00:27:24,280 --> 00:27:25,642 to make the fits. 567 00:27:25,642 --> 00:27:26,850 OK, so you can do either one. 568 00:27:26,850 --> 00:27:28,240 You're in charge. 569 00:27:28,240 --> 00:27:29,674 Right, you get the problem. 570 00:27:29,674 --> 00:27:31,840 If somebody tells you, please give me this integral, 571 00:27:31,840 --> 00:27:33,610 they don't care how you do it. 572 00:27:33,610 --> 00:27:35,485 They just want you to come back and give them 573 00:27:35,485 --> 00:27:38,986 the true value of the integral to some significant figures. 574 00:27:38,986 --> 00:27:40,610 So it's up to you to decide what to do. 575 00:27:40,610 --> 00:27:42,252 So you can decide what to do. 576 00:27:42,252 --> 00:27:43,960 What a lot of people do in the math world 577 00:27:43,960 --> 00:27:46,630 is this one, interpolation, where 578 00:27:46,630 --> 00:27:50,490 you're assuming that the f of tn's are known exactly. 579 00:27:50,490 --> 00:27:53,560 And so your function evaluating function is great. 580 00:27:53,560 --> 00:27:57,920 And in that case, you can write down this nice linear equation. 581 00:27:57,920 --> 00:28:01,440 So we're trying to make f of tn with a true function 582 00:28:01,440 --> 00:28:06,006 to be equal to ck phi n of t. 583 00:28:06,006 --> 00:28:07,469 Sorry, phi k of t. 584 00:28:12,260 --> 00:28:14,790 Which n in [INAUDIBLE] k is equal to n. 585 00:28:14,790 --> 00:28:17,900 And I can rewrite this this way. 586 00:28:17,900 --> 00:28:20,446 This is like a set of numbers. 587 00:28:20,446 --> 00:28:21,570 So I can write as a vector. 588 00:28:21,570 --> 00:28:23,540 I'll call it y. 589 00:28:23,540 --> 00:28:25,520 And this thing has two indices on it, 590 00:28:25,520 --> 00:28:27,300 so it looks like a matrix. 591 00:28:27,300 --> 00:28:30,384 So I'll call that m. 592 00:28:30,384 --> 00:28:34,620 And c is my vector of unknown parameters. 593 00:28:34,620 --> 00:28:37,560 And wow, I know how to solve this one. 594 00:28:37,560 --> 00:28:40,730 So c is equal to m backslash y. 595 00:28:46,200 --> 00:28:49,770 Now, for this to work and have a unique solution, 596 00:28:49,770 --> 00:28:52,200 I need m not to be singular. 597 00:28:52,200 --> 00:28:54,090 And what that means is that the columns of m 598 00:28:54,090 --> 00:28:55,530 have to be linearly independent. 599 00:28:55,530 --> 00:28:57,120 So what are the columns of m? 600 00:28:57,120 --> 00:29:10,330 They are phi k of t1, phi k of t2, phi k of tn. 601 00:29:14,160 --> 00:29:21,315 And I want that this column is not a sum of the other columns. 602 00:29:36,030 --> 00:29:38,100 So this is kind of requirement for j 603 00:29:38,100 --> 00:29:40,542 not equal [? to k ?] This is kind of a requirement 604 00:29:40,542 --> 00:29:42,250 if I'm going to do this kind of procedure 605 00:29:42,250 --> 00:29:46,390 is I can't choose a set of basis functions 606 00:29:46,390 --> 00:29:48,850 that's going to have one of them be a linear combination 607 00:29:48,850 --> 00:29:49,810 of the other ones. 608 00:29:49,810 --> 00:29:50,980 Because then I'm going to get a singular matrix. 609 00:29:50,980 --> 00:29:52,210 My whole process is going to have a big problem 610 00:29:52,210 --> 00:29:53,540 right from the beginning. 611 00:29:53,540 --> 00:29:56,290 So I have to ensure that this is not true. 612 00:29:56,290 --> 00:29:57,550 Yeah, that this is true. 613 00:29:57,550 --> 00:29:59,290 These are not equal to some of them. 614 00:30:02,400 --> 00:30:14,570 Another way to write that is to choose the the columns that phi 615 00:30:14,570 --> 00:30:16,869 to be orthogonal to each other. 616 00:30:16,869 --> 00:30:18,160 And so I can write it this way. 617 00:30:34,760 --> 00:30:37,502 OK, so I want them to be orthogonal to each other. 618 00:30:37,502 --> 00:30:39,585 And then you can see how this kind of generalizes. 619 00:30:39,585 --> 00:30:42,280 If I had a lot of t's, this would sort of 620 00:30:42,280 --> 00:30:49,432 look like the integral of phi k of t, phi j of t is equal to 0. 621 00:30:55,770 --> 00:30:58,770 And so this gives the idea of orthogonal functions. 622 00:30:58,770 --> 00:31:03,490 So this is orthogonal functions, and my discrete point's tn. 623 00:31:03,490 --> 00:31:07,410 And this is the general concept of orthogonal functions. 624 00:31:07,410 --> 00:31:08,880 And then when I define it this way, 625 00:31:08,880 --> 00:31:12,849 I need to actually decide my own a's and b's that I want to use. 626 00:31:12,849 --> 00:31:14,640 And if you use different ones, you actually 627 00:31:14,640 --> 00:31:16,015 have different rules about what's 628 00:31:16,015 --> 00:31:17,280 orthogonal to each other. 629 00:31:17,280 --> 00:31:20,040 And also, sometimes, people are fishy, 630 00:31:20,040 --> 00:31:22,030 and they [? spit ?] a little w of t 631 00:31:22,030 --> 00:31:25,249 in here just to do, because they have some special problem where 632 00:31:25,249 --> 00:31:26,790 this w keeps showing up or something, 633 00:31:26,790 --> 00:31:27,831 and they want to do that. 634 00:31:27,831 --> 00:31:29,425 But that's up to them. 635 00:31:29,425 --> 00:31:30,800 All right, we won't specify that. 636 00:31:30,800 --> 00:31:33,200 But this is the general idea of making basis functions 637 00:31:33,200 --> 00:31:35,540 orthogonal to each other. 638 00:31:35,540 --> 00:31:36,980 And it's very analogous to making 639 00:31:36,980 --> 00:31:38,300 vectors orthogonal to each other, which 640 00:31:38,300 --> 00:31:39,133 is what we did here. 641 00:31:39,133 --> 00:31:41,140 These are actually vectors. 642 00:31:41,140 --> 00:31:43,430 Right, it's just as the vector gets longer and longer, 643 00:31:43,430 --> 00:31:44,990 it gets more and more like the continuous function, 644 00:31:44,990 --> 00:31:46,580 and we have more t points. 645 00:31:46,580 --> 00:31:48,420 And eventually, it looks like that integral 646 00:31:48,420 --> 00:31:49,660 as you go to infinity. 647 00:31:52,040 --> 00:31:53,540 So anyway, so a lot of times, people 648 00:31:53,540 --> 00:31:57,290 will choose functions which are designed to be orthogonal. 649 00:31:57,290 --> 00:31:59,720 And one way to do it is just to go with a continuous one 650 00:31:59,720 --> 00:32:00,219 to start. 651 00:32:02,690 --> 00:32:05,570 Now, we can choose any basis functions we want. 652 00:32:05,570 --> 00:32:07,430 Over here, we already chose a whole bunch. 653 00:32:07,430 --> 00:32:09,180 You guys have done this already. 654 00:32:09,180 --> 00:32:10,070 I don't know if you knew you were doing it. 655 00:32:10,070 --> 00:32:11,570 But you chose a lot of different basis functions completely 656 00:32:11,570 --> 00:32:13,260 different from each other. 657 00:32:13,260 --> 00:32:16,170 You can choose the different ones if you want. 658 00:32:16,170 --> 00:32:21,680 And one popular choice is polynomials. 659 00:32:21,680 --> 00:32:23,180 But it's not the only choice. 660 00:32:23,180 --> 00:32:25,755 And in fact, a little bit later, we're going to, 661 00:32:25,755 --> 00:32:27,410 in the homework problem six, I think 662 00:32:27,410 --> 00:32:28,640 we're going to have a problem where 663 00:32:28,640 --> 00:32:30,530 you have to use a different kind of basis functions. 664 00:32:30,530 --> 00:32:32,238 You can use any basis functions you want. 665 00:32:32,238 --> 00:32:35,110 And there's kind of ones that match naturally to the problem. 666 00:32:35,110 --> 00:32:38,930 However, people know a lot about polynomials. 667 00:32:38,930 --> 00:32:42,770 And so they've decided to use polynomials a lot. 668 00:32:42,770 --> 00:32:44,530 And so that's the most common choice 669 00:32:44,530 --> 00:32:46,500 for numerical integration is polynomials. 670 00:32:46,500 --> 00:32:48,870 One thing is polynomials are super easy to integrate, 671 00:32:48,870 --> 00:32:51,036 because that was the first thing you learned, right, 672 00:32:51,036 --> 00:32:52,700 how to do the intervals of polynomials. 673 00:32:52,700 --> 00:32:54,507 And also, they're easy to evaluate. 674 00:32:54,507 --> 00:32:56,090 It only takes n operations to evaluate 675 00:32:56,090 --> 00:32:58,340 a polynomial of nth order. 676 00:32:58,340 --> 00:33:01,304 And there's a ton of theorems, and things 677 00:33:01,304 --> 00:33:02,470 are known about polynomials. 678 00:33:02,470 --> 00:33:03,800 We have all kinds of special properties. 679 00:33:03,800 --> 00:33:04,850 And we have all kinds of convenient tools 680 00:33:04,850 --> 00:33:05,600 for handling them. 681 00:33:05,600 --> 00:33:08,710 So let's go with polynomials for now. 682 00:33:08,710 --> 00:33:12,490 Now, I want to warn you, though, that how you order polynomials 683 00:33:12,490 --> 00:33:17,230 are really not very good fitting functions 684 00:33:17,230 --> 00:33:20,020 or interpolating functions for a lot of problems. 685 00:33:20,020 --> 00:33:22,420 And in particular, for this one I 686 00:33:22,420 --> 00:33:25,620 showed here where this is the function. 687 00:33:25,620 --> 00:33:31,682 Here, that function, you can try to fit that with higher 688 00:33:31,682 --> 00:33:32,890 and higher order polynomials. 689 00:33:32,890 --> 00:33:33,973 So I can have my function. 690 00:33:38,670 --> 00:33:41,700 The function value at 0 is 1. 691 00:33:41,700 --> 00:33:44,094 So the function looks kind of reasonable. 692 00:33:47,060 --> 00:33:50,650 Like that, symmetrical, better than I drew it. 693 00:33:50,650 --> 00:33:56,300 And you can put some number of points in here. 694 00:33:56,300 --> 00:33:58,550 And then you can do an interpolating polynomial that 695 00:33:58,550 --> 00:34:01,640 goes through all the points. 696 00:34:01,640 --> 00:34:07,400 And if you do it for this one with five points, 697 00:34:07,400 --> 00:34:11,409 you'll get something that looks like this. 698 00:34:11,409 --> 00:34:18,150 It goes negative, comes up like this, goes negative again. 699 00:34:18,150 --> 00:34:21,570 So that's one you get if you use five 700 00:34:21,570 --> 00:34:23,969 points across this interval. 701 00:34:23,969 --> 00:34:26,340 If you use 10 points, then you get 702 00:34:26,340 --> 00:34:41,960 one that looks like this where the peak value here 703 00:34:41,960 --> 00:34:45,280 is two where as the max of the original function's just one. 704 00:34:45,280 --> 00:34:47,170 And the original function's always positive, 705 00:34:47,170 --> 00:34:49,520 but they both go negative. 706 00:34:49,520 --> 00:34:54,130 So even within the range of the interpolation, 707 00:34:54,130 --> 00:34:56,590 the polynomials can have wild swings. 708 00:34:56,590 --> 00:35:00,740 And so this first one was for a fifth order polynomial, 709 00:35:00,740 --> 00:35:03,280 or maybe fourth order, five points. 710 00:35:03,280 --> 00:35:06,630 And this one is for the next 10th order or 9th order 711 00:35:06,630 --> 00:35:07,350 polynomial. 712 00:35:07,350 --> 00:35:09,760 I'm not sure which one. 713 00:35:09,760 --> 00:35:10,260 Yes. 714 00:35:10,260 --> 00:35:12,760 AUDIENCE: What points are these polynomials trying to fit? 715 00:35:12,760 --> 00:35:15,430 Or does it seem like [INAUDIBLE]?? 716 00:35:15,430 --> 00:35:18,220 PROFESSOR: So I didn't draw very well. 717 00:35:18,220 --> 00:35:19,960 When you guys go home, do this. 718 00:35:19,960 --> 00:35:23,100 Right, on Matlab it will take you a minute 719 00:35:23,100 --> 00:35:26,720 to solve this equation where you put in the function values 720 00:35:26,720 --> 00:35:27,870 here. 721 00:35:27,870 --> 00:35:33,670 And put in the polynomial values that [INAUDIBLE] 722 00:35:33,670 --> 00:35:37,215 here for the monomials, t to the first power, 723 00:35:37,215 --> 00:35:39,340 to the second power, to the third power, like that. 724 00:35:39,340 --> 00:35:43,230 And just do it, and you can get the [? pot ?] yourself. 725 00:35:43,230 --> 00:35:45,230 And just to convince yourself about the problem, 726 00:35:45,230 --> 00:35:47,042 and this is very important to keep in mind. 727 00:35:47,042 --> 00:35:49,250 Because we're so used to using polynomials that think 728 00:35:49,250 --> 00:35:50,708 they're the solution to everything. 729 00:35:50,708 --> 00:35:51,830 And Excel does them great. 730 00:35:51,830 --> 00:35:55,010 But it's not this solution to everything. 731 00:35:55,010 --> 00:35:56,890 Nonetheless, we're going to plow on and keep 732 00:35:56,890 --> 00:35:58,098 using them for a few minutes. 733 00:36:02,790 --> 00:36:08,200 So one possibility when I'm choosing this basis set 734 00:36:08,200 --> 00:36:10,280 is I could choose what are called 735 00:36:10,280 --> 00:36:17,260 monomials, which is phi k is equal to x to the k minus 1. 736 00:36:21,370 --> 00:36:23,575 And this one is super simple. 737 00:36:26,200 --> 00:36:28,390 And so a lot of times it satisfies 738 00:36:28,390 --> 00:36:32,170 this orthogonality thing, which I showed you before. 739 00:36:32,170 --> 00:36:36,280 But it's a really bad choice usually. 740 00:36:36,280 --> 00:36:38,480 And what it is is that the matrix 741 00:36:38,480 --> 00:36:41,170 you have to solve for this interpolation, 742 00:36:41,170 --> 00:36:46,960 if you have evenly-spaced points and you have x to the k, 743 00:36:46,960 --> 00:36:49,330 it turns out that this matrix usually 744 00:36:49,330 --> 00:36:52,120 has a condition number that grows exponentially 745 00:36:52,120 --> 00:36:55,430 with the number of terms in your basis, 746 00:36:55,430 --> 00:36:57,880 the number of how many monomials you include. 747 00:36:57,880 --> 00:37:00,010 So condition numbers that grow exponentially 748 00:37:00,010 --> 00:37:03,400 is really bad idea. 749 00:37:03,400 --> 00:37:05,252 So this is not what people use. 750 00:37:05,252 --> 00:37:06,710 So you might be tempted to do this, 751 00:37:06,710 --> 00:37:08,540 and maybe you did it once in your life. 752 00:37:08,540 --> 00:37:10,180 I won't blame you if you did. 753 00:37:10,180 --> 00:37:11,980 But it's not a really smart idea. 754 00:37:11,980 --> 00:37:14,350 So instead, a lot of the mathematicians 755 00:37:14,350 --> 00:37:16,270 of the 1700s and 1800s spent a lot 756 00:37:16,270 --> 00:37:19,040 of time worrying about this. 757 00:37:19,040 --> 00:37:24,250 And they'd end up deciding to use other polynomials that 758 00:37:24,250 --> 00:37:26,050 have people's names on them. 759 00:37:26,050 --> 00:37:35,020 So a guy named Lagrange, the same guy 760 00:37:35,020 --> 00:37:38,340 who did the multiplier. 761 00:37:38,340 --> 00:37:40,594 He suggested some polynomials. 762 00:37:40,594 --> 00:37:41,760 They're kind of complicated. 763 00:37:41,760 --> 00:37:43,140 They're given in the book. 764 00:37:43,140 --> 00:37:47,440 But what they really are is I'm taking the monomial basis set, 765 00:37:47,440 --> 00:37:49,650 and I'm making a new basis set. 766 00:37:49,650 --> 00:37:58,430 So I want phi l to be equal to the sum dlk of phi k, 767 00:37:58,430 --> 00:38:02,640 where this phi k it the monomials. 768 00:38:02,640 --> 00:38:04,927 So really, you can choose any d you want. 769 00:38:04,927 --> 00:38:06,510 You're just taking linear combinations 770 00:38:06,510 --> 00:38:10,920 of the original monomial basis set to make up polynomials. 771 00:38:10,920 --> 00:38:12,690 And if you cleverly choose the d, 772 00:38:12,690 --> 00:38:14,610 you can choose it to have special properties. 773 00:38:14,610 --> 00:38:16,860 So Lagrange gave a certain choice 774 00:38:16,860 --> 00:38:20,220 of the d that is brilliant in the sense 775 00:38:20,220 --> 00:38:23,130 that, when you solve the problem m 776 00:38:23,130 --> 00:38:27,720 times c is equal to y to try to do your interpolation, it turns 777 00:38:27,720 --> 00:38:34,440 out that the c's you want, ck, is equal to f of tk. 778 00:38:34,440 --> 00:38:37,310 So you don't have to do any work to solve what the c's are. 779 00:38:37,310 --> 00:38:39,930 It's just the f values directly. 780 00:38:39,930 --> 00:38:44,530 No chance for condition numbers, no problems, that's it. 781 00:38:44,530 --> 00:38:47,790 So that's the Lagrange polynomials. 782 00:38:47,790 --> 00:38:51,500 Now, Newton had a different idea. 783 00:38:51,500 --> 00:38:55,105 So this is for Lagrange polynomials. 784 00:39:00,352 --> 00:39:04,510 If you choose a Lagrangian [? polynomial ?] basis, 785 00:39:04,510 --> 00:39:05,200 then this is it. 786 00:39:05,200 --> 00:39:06,650 And you do a square system. 787 00:39:06,650 --> 00:39:07,950 That's it. 788 00:39:07,950 --> 00:39:12,065 If you choose the Newton, even he got into this business too. 789 00:39:12,065 --> 00:39:13,770 And he made some polynomials. 790 00:39:13,770 --> 00:39:16,800 And he made ones that have a special property 791 00:39:16,800 --> 00:39:19,830 that if you solve it-- 792 00:39:19,830 --> 00:39:21,270 so he chose a different d. 793 00:39:21,270 --> 00:39:23,460 You still have to do some solve. 794 00:39:23,460 --> 00:39:25,620 You get the solution. 795 00:39:25,620 --> 00:39:28,500 And then you decide that I don't like that interpretation. 796 00:39:28,500 --> 00:39:30,060 I want to put some more points in. 797 00:39:30,060 --> 00:39:32,130 I'm going to do a few more function evaluations. 798 00:39:32,130 --> 00:39:33,420 I'm going to add some more. 799 00:39:33,420 --> 00:39:38,010 So from the first solve, I have the values of c1, c2, c3, 800 00:39:38,010 --> 00:39:41,550 up to c the number of points I had initially. 801 00:39:41,550 --> 00:39:44,510 And now, I want to add one more data point, 802 00:39:44,510 --> 00:39:47,811 or add one more f of tn, one more function evaluation. 803 00:39:47,811 --> 00:39:49,435 I want to figure out with c of n plus 1 804 00:39:49,435 --> 00:39:51,330 is, what the next one is. 805 00:39:51,330 --> 00:39:54,420 It turns out that it's his polynomials. 806 00:39:54,420 --> 00:39:56,310 The values these c's don't change 807 00:39:56,310 --> 00:39:59,140 when I make this matrix a little bit bigger by adding one more 808 00:39:59,140 --> 00:40:01,060 y. 809 00:40:01,060 --> 00:40:04,060 And I add one more parameter. 810 00:40:04,060 --> 00:40:05,830 And I make one more row in the matrix, 811 00:40:05,830 --> 00:40:08,070 one more column in the matrix. 812 00:40:08,070 --> 00:40:10,360 When I do that, it's a unique property 813 00:40:10,360 --> 00:40:12,675 that none of these other c's change. 814 00:40:12,675 --> 00:40:15,300 The solution is going to be the same, except for the n plus one 815 00:40:15,300 --> 00:40:17,530 [? level. ?] So then I can write a neat, cute little equation 816 00:40:17,530 --> 00:40:18,613 just for the n plus 1 one. 817 00:40:18,613 --> 00:40:20,920 I don't have to mess with the previous ones. 818 00:40:20,920 --> 00:40:23,461 And so I'm just adding a basis function that sort of polishes 819 00:40:23,461 --> 00:40:25,350 up my original fit. 820 00:40:25,350 --> 00:40:28,180 So that's a nice property too. 821 00:40:28,180 --> 00:40:30,550 Whereas with the Lagrange ones, you 822 00:40:30,550 --> 00:40:32,110 have to recompute everything. 823 00:40:32,110 --> 00:40:33,276 And all the c's will change. 824 00:40:38,739 --> 00:40:40,030 Actually, the c's won't change. 825 00:40:40,030 --> 00:40:40,738 I take that back. 826 00:40:40,738 --> 00:40:41,440 You add points. 827 00:40:41,440 --> 00:40:44,010 The c's won't change, but the actual polynomials will change. 828 00:40:44,010 --> 00:40:45,720 The d will change. 829 00:40:45,720 --> 00:40:49,200 Because the d depends on the choice of the time 830 00:40:49,200 --> 00:40:51,750 points in a complicated way in Lagrange. 831 00:40:56,046 --> 00:40:57,920 So now, you start to think, oh my goodness, I 832 00:40:57,920 --> 00:40:59,450 have all these choices, I can use Lagrange. 833 00:40:59,450 --> 00:41:00,241 I could use Newton. 834 00:41:00,241 --> 00:41:01,550 I could use monomials. 835 00:41:01,550 --> 00:41:03,740 Maybe, I should go back to bed. 836 00:41:03,740 --> 00:41:06,912 And you think, well, could someone just 837 00:41:06,912 --> 00:41:08,120 tell me what the best one is? 838 00:41:08,120 --> 00:41:08,810 And I'll just do that one. 839 00:41:08,810 --> 00:41:11,018 I don't need to learn about all this historical stuff 840 00:41:11,018 --> 00:41:13,220 about development of polynomial mathematics. 841 00:41:13,220 --> 00:41:15,450 I'd just like you to just tell me what the answer is. 842 00:41:15,450 --> 00:41:20,940 So Gauss, typically, figured out what the best on is. 843 00:41:20,940 --> 00:41:26,050 And he said, you know, if we're going to do integrals, 844 00:41:26,050 --> 00:41:28,340 we'd really like to just write our integrals this way. 845 00:41:28,340 --> 00:41:32,600 So our approximate value to be equal to some weight 846 00:41:32,600 --> 00:41:38,340 functions times the function evaluations at some times. 847 00:41:38,340 --> 00:41:40,170 Right, that's easy to evaluate. 848 00:41:40,170 --> 00:41:41,322 That would be good. 849 00:41:41,322 --> 00:41:43,780 If I pre-computed the weight functions [? and ?] the times, 850 00:41:43,780 --> 00:41:45,560 then I put any function f in here, 851 00:41:45,560 --> 00:41:47,530 I could just do this like zip. 852 00:41:47,530 --> 00:41:51,486 It's one dot product of the f of tn's I'll have the solution. 853 00:41:51,486 --> 00:41:52,360 So I want to do that. 854 00:41:52,360 --> 00:41:55,390 So what's the best possible choices of the w's and the t's 855 00:41:55,390 --> 00:41:57,660 to make this really be close? 856 00:41:57,660 --> 00:41:59,370 Well, the problem is that there's 857 00:41:59,370 --> 00:42:00,420 an infinite number of functions that somebody 858 00:42:00,420 --> 00:42:01,860 might want to integrate. 859 00:42:01,860 --> 00:42:02,910 I can't do them all. 860 00:42:02,910 --> 00:42:05,280 So he said, let's just find the one that 861 00:42:05,280 --> 00:42:12,220 works for as many monomials as possible to the highest order, 862 00:42:12,220 --> 00:42:14,790 as far out as I can go. 863 00:42:14,790 --> 00:42:20,310 And it turns out I have nw's. 864 00:42:20,310 --> 00:42:21,750 I have nt's. 865 00:42:21,750 --> 00:42:25,280 So I have two n parameters that I can fiddle around with. 866 00:42:25,280 --> 00:42:28,270 And so I should get something like two n's in my polynomials, 867 00:42:28,270 --> 00:42:30,560 maybe 2n minus 1. 868 00:42:30,560 --> 00:42:34,630 So what I want is that I want to choose this 869 00:42:34,630 --> 00:42:39,060 so that for the monomial 1, and the monomial t, 870 00:42:39,060 --> 00:42:42,880 and the monomial t squared, that this i is actually 871 00:42:42,880 --> 00:42:45,270 going to be exact. 872 00:42:45,270 --> 00:42:49,190 So what I want is for the monomial 1, 873 00:42:49,190 --> 00:42:52,970 I want the summation of wn. 874 00:42:52,970 --> 00:42:55,760 My f for the monomial 1 is always 1. 875 00:42:55,760 --> 00:42:57,620 So it's times 1. 876 00:42:57,620 --> 00:42:59,600 I want this to equal the integral from a 877 00:42:59,600 --> 00:43:06,330 to b of 1dt, which is b minus a. 878 00:43:06,330 --> 00:43:09,730 OK, so that gives me one equation for the w's. 879 00:43:09,730 --> 00:43:15,310 Then for the next one, I want the summation wn 880 00:43:15,310 --> 00:43:21,478 times tn to equal the integral of a to b of t. 881 00:43:21,478 --> 00:43:28,570 dt, which is equal to b squared over a number 2 minus a squared 882 00:43:28,570 --> 00:43:32,390 over 2 if I did my calculus correctly. 883 00:43:32,390 --> 00:43:33,796 And so there's another question. 884 00:43:33,796 --> 00:43:35,420 I have another question for these guys. 885 00:43:35,420 --> 00:43:37,919 And I keep going until I get enough equations that determine 886 00:43:37,919 --> 00:43:39,395 all the w's and c's. 887 00:43:39,395 --> 00:43:40,895 All right, so I do it for t squared. 888 00:43:40,895 --> 00:43:49,410 So I can wn tn squared is equal to [INAUDIBLE] b t squared dt, 889 00:43:49,410 --> 00:43:52,260 and so on, OK. 890 00:43:52,260 --> 00:43:54,510 And so if you do this process, you 891 00:43:54,510 --> 00:43:57,840 get the weights and the points to use 892 00:43:57,840 --> 00:44:00,950 that will give you what's called the Gaussian quadrature rule. 893 00:44:00,950 --> 00:44:05,470 Now, the brilliant thing about this is that by this setup, 894 00:44:05,470 --> 00:44:08,770 it's going to be exact for integrating any polynomial up 895 00:44:08,770 --> 00:44:10,641 to order 2n minus 1. 896 00:44:13,347 --> 00:44:17,990 So if I choose three, here I'll be 897 00:44:17,990 --> 00:44:20,570 able to determine-- if I choose three points, 898 00:44:20,570 --> 00:44:22,550 I'll have six parameters. 899 00:44:22,550 --> 00:44:27,898 So I'll be able determine up to the 11-- 900 00:44:27,898 --> 00:44:29,814 11th order polynomial would integrate exactly. 901 00:44:33,610 --> 00:44:36,482 No, that's wrong. 902 00:44:36,482 --> 00:44:38,480 n's only 3. 903 00:44:38,480 --> 00:44:41,370 I have six parameters. 904 00:44:41,370 --> 00:44:42,420 I can't do this. 905 00:44:42,420 --> 00:44:44,601 2 times 3 minus 1, what's that? 906 00:44:44,601 --> 00:44:46,800 5, there we go. 907 00:44:46,800 --> 00:44:50,180 To the fifth order polynomial integrate exactly, 908 00:44:50,180 --> 00:44:52,580 I'll have an [? error ?] for sixth order polynomial. 909 00:44:52,580 --> 00:44:56,280 So this is a rule. 910 00:44:56,280 --> 00:45:02,160 So before, if I gave you three points like here, most of you 911 00:45:02,160 --> 00:45:03,720 probably would use Simpson's rule. 912 00:45:03,720 --> 00:45:05,295 Fit this to a parabola. 913 00:45:05,295 --> 00:45:09,030 It has three parameters, a, ax squared plus bx plus c, 914 00:45:09,030 --> 00:45:10,955 right, for a second order polynomial. 915 00:45:10,955 --> 00:45:11,580 I would fit it. 916 00:45:11,580 --> 00:45:12,330 Then I'd integrate that. 917 00:45:12,330 --> 00:45:14,710 That would be the normal thing that a lot of people would do. 918 00:45:14,710 --> 00:45:15,510 I don't know if you would do that. 919 00:45:15,510 --> 00:45:17,402 But a lot of people would do that. 920 00:45:17,402 --> 00:45:19,110 So that would be that what you would get. 921 00:45:19,110 --> 00:45:20,100 And what [? error ?] would you get? 922 00:45:20,100 --> 00:45:21,766 You'd have error that would be the scale 923 00:45:21,766 --> 00:45:25,110 of the width of the intervals to the fourth power. 924 00:45:25,110 --> 00:45:28,430 Right, because Simpson's rule exact up to the cubic. 925 00:45:28,430 --> 00:45:31,035 You guys remember this from high school? 926 00:45:31,035 --> 00:45:33,160 I don't know if you remember delta t to the fourth. 927 00:45:33,160 --> 00:45:34,170 Yeah, OK. 928 00:45:34,170 --> 00:45:41,260 So Simpson's rule has an order of delta t to the fourth. 929 00:45:41,260 --> 00:45:48,530 But Gauss's method has order of delta t to the sixth. 930 00:45:48,530 --> 00:45:51,810 So this is for three points. 931 00:45:51,810 --> 00:45:54,060 In Gauss's method, you can extend it to as many points 932 00:45:54,060 --> 00:45:56,550 as you want. 933 00:45:56,550 --> 00:45:58,220 All right, so that's kind of nice. 934 00:45:58,220 --> 00:45:59,890 You only do three points, you get 935 00:45:59,890 --> 00:46:02,160 [? error ?] of the sixth order. 936 00:46:02,160 --> 00:46:04,030 So that seems pretty efficient. 937 00:46:04,030 --> 00:46:05,305 So people kept going. 938 00:46:05,305 --> 00:46:08,490 And so the very, very popular one 939 00:46:08,490 --> 00:46:10,860 use, actually, seven points for Gauss. 940 00:46:14,940 --> 00:46:20,270 So you'd be able to get to the n minus 1. 941 00:46:20,270 --> 00:46:23,220 So you'd be exact for 13. 942 00:46:23,220 --> 00:46:26,832 So if you have delta t to the 14th power as your error, 943 00:46:26,832 --> 00:46:28,450 that's pretty good. 944 00:46:28,450 --> 00:46:36,760 And then what people often do is they'll add eight more points. 945 00:46:36,760 --> 00:46:39,220 And they make another method that uses all these points. 946 00:46:39,220 --> 00:46:44,894 It's called [? Conrad's ?] method. 947 00:46:44,894 --> 00:46:46,810 He made another method that gives you the best 948 00:46:46,810 --> 00:46:49,476 possible solution if you had the 15 points, given that 7 of them 949 00:46:49,476 --> 00:46:52,090 are already fixed by the Gauss procedure. 950 00:46:52,090 --> 00:46:56,340 And so you can get the integral using Gauss's method 951 00:46:56,340 --> 00:46:57,740 and using [? Conrad's ?] method. 952 00:47:07,480 --> 00:47:11,590 So you can compute the i using Gauss's method, the seven 953 00:47:11,590 --> 00:47:14,122 points, the i with the [? Conrad ?] method. 954 00:47:14,122 --> 00:47:15,580 I don't know how to spell his name, 955 00:47:15,580 --> 00:47:18,260 Actually I think it's like this. 956 00:47:18,260 --> 00:47:19,930 15 points. 957 00:47:19,930 --> 00:47:22,580 And then you can take the difference between these two. 958 00:47:22,580 --> 00:47:24,080 And if these two are really similar, 959 00:47:24,080 --> 00:47:26,121 then you might think your integral's pretty good. 960 00:47:26,121 --> 00:47:27,896 Because they're pretty different methods, 961 00:47:27,896 --> 00:47:29,770 and if they get the same number, you're good. 962 00:47:29,770 --> 00:47:31,690 And people have done a lot of numerical studies 963 00:47:31,690 --> 00:47:33,231 in this particular case, because it's 964 00:47:33,231 --> 00:47:35,350 used a lot in actual practice. 965 00:47:35,350 --> 00:47:39,750 And it turns out that the [? error ?] 966 00:47:39,750 --> 00:47:45,760 is almost always less than this quantity to the 3 Gauss power. 967 00:47:45,760 --> 00:47:48,480 So you have an error estimate. 968 00:47:48,480 --> 00:47:51,760 So you do 15 function evaluations. 969 00:47:51,760 --> 00:47:54,770 You get a pretty accurate value for your integral, 970 00:47:54,770 --> 00:47:59,480 maybe, if you can be fit well with whatever order polynomial, 971 00:47:59,480 --> 00:48:02,120 13th order polynomial. 972 00:48:02,120 --> 00:48:04,850 13th order polynomial, yeah. 973 00:48:04,850 --> 00:48:08,190 And so any function that can be fit pretty 974 00:48:08,190 --> 00:48:09,690 well with the 13th order polynomial, 975 00:48:09,690 --> 00:48:11,700 you should get pretty good value in the integral. 976 00:48:11,700 --> 00:48:14,210 And you can check by comparing to the [? Conrad ?] estimate, 977 00:48:14,210 --> 00:48:15,550 which uses a lot more points. 978 00:48:15,550 --> 00:48:17,950 And if two of these guys give the same number 979 00:48:17,950 --> 00:48:19,550 or pretty close, then you're pretty 980 00:48:19,550 --> 00:48:21,027 confident that you're good. 981 00:48:21,027 --> 00:48:23,360 And there's even a math proof about this particular one. 982 00:48:23,360 --> 00:48:26,810 So this is like really popular, and it's 983 00:48:26,810 --> 00:48:28,310 a lot better than Simpson's rule, 984 00:48:28,310 --> 00:48:29,600 but it's kind of complicated. 985 00:48:29,600 --> 00:48:31,100 And part of it you have to recompute 986 00:48:31,100 --> 00:48:35,000 these tn's and wn's for the kind of problem you have. 987 00:48:35,000 --> 00:48:35,666 Yeah. 988 00:48:35,666 --> 00:48:41,360 AUDIENCE: If you're using 15 points already, 989 00:48:41,360 --> 00:48:45,440 why not just split your step size smaller? 990 00:48:45,440 --> 00:48:48,185 And if you did that, how would the error [INAUDIBLE] 991 00:48:48,185 --> 00:48:49,075 the simplest? 992 00:48:49,075 --> 00:48:50,117 Use a trapezoid. 993 00:48:50,117 --> 00:48:51,700 PROFESSOR: Trapezoid or Simpson's rule 994 00:48:51,700 --> 00:48:52,370 of something like that. 995 00:48:52,370 --> 00:48:52,590 Yes. 996 00:48:52,590 --> 00:48:54,381 AUDIENCE: How would the error with those 15 997 00:48:54,381 --> 00:48:59,450 steps using trapezoid compare to [? Conrad ?] and Gauss? 998 00:48:59,450 --> 00:49:02,190 PROFESSOR: So we should probably do the math. 999 00:49:02,190 --> 00:49:06,410 So say we did Simpson's rule, and we have delta t, 1000 00:49:06,410 --> 00:49:11,960 and we divide it by 15 or something. 1001 00:49:11,960 --> 00:49:14,120 Then we can figure out what this is, 1002 00:49:14,120 --> 00:49:17,919 but it's still going to have the scaling to the fourth power. 1003 00:49:17,919 --> 00:49:19,460 So the prefactor will be really tiny, 1004 00:49:19,460 --> 00:49:22,320 because it will be 1/15 to the fourth power. 1005 00:49:22,320 --> 00:49:24,239 But the scaling is bad. 1006 00:49:24,239 --> 00:49:26,530 So that if I decide that my integral's not good enough, 1007 00:49:26,530 --> 00:49:29,780 and I need to do more points as I double the number of points 1008 00:49:29,780 --> 00:49:33,180 to get better accuracy, my thing won't improve that much. 1009 00:49:33,180 --> 00:49:35,430 Because it's still only going to use the fourth power, 1010 00:49:35,430 --> 00:49:36,590 whereas the Gauss' one say it's going 1011 00:49:36,590 --> 00:49:38,330 use the sixth power if I only use the three points. 1012 00:49:38,330 --> 00:49:40,204 Where we use 15 points, is going to get scale 1013 00:49:40,204 --> 00:49:43,110 to some 16th power or 15th power or some really high power. 1014 00:49:43,110 --> 00:49:47,880 So for doing sequences, the Gauss' one is a lot better. 1015 00:49:47,880 --> 00:49:51,110 However, what people do is they usually 1016 00:49:51,110 --> 00:49:52,730 pre-solve the system of equations. 1017 00:49:52,730 --> 00:49:55,064 Because this system of equation does not depend. 1018 00:49:55,064 --> 00:49:56,480 When you do it with the monomials, 1019 00:49:56,480 --> 00:49:58,310 it doesn't depend what f is. 1020 00:49:58,310 --> 00:50:01,527 You can figure out the optimal t's and w's ahead of time. 1021 00:50:01,527 --> 00:50:03,110 And so they [INAUDIBLE] tables of them 1022 00:50:03,110 --> 00:50:05,900 depending on where you want to go. 1023 00:50:05,900 --> 00:50:08,510 People recomputed them, and then you can get them that way. 1024 00:50:11,217 --> 00:50:12,050 Time is running out. 1025 00:50:12,050 --> 00:50:14,060 I just want to jump up. 1026 00:50:16,820 --> 00:50:19,020 Now, for single integrals, I didn't have to tell you 1027 00:50:19,020 --> 00:50:19,520 and of this. 1028 00:50:19,520 --> 00:50:22,960 Because you could have just solved it with your ODE solver. 1029 00:50:22,960 --> 00:50:26,870 So you can just go to ode15s or ode45, 1030 00:50:26,870 --> 00:50:28,400 and you'll get good enough. 1031 00:50:28,400 --> 00:50:30,080 It won't be as efficient as Gaussian. 1032 00:50:30,080 --> 00:50:31,044 But who cares. 1033 00:50:31,044 --> 00:50:32,210 You still get a good answer. 1034 00:50:32,210 --> 00:50:33,440 You can report to your boss. 1035 00:50:33,440 --> 00:50:35,810 It'll take you three minutes to call 1036 00:50:35,810 --> 00:50:38,099 it e45 to integrate the one-dimensional integral. 1037 00:50:38,099 --> 00:50:40,640 The real challenge comes with the multidimensional integrals. 1038 00:50:40,640 --> 00:50:44,580 So in a lot of problems, we have more than one variable 1039 00:50:44,580 --> 00:50:47,600 we want to integrate over. 1040 00:50:47,600 --> 00:50:51,380 And so if you have an integral, say, of a 1041 00:50:51,380 --> 00:50:58,620 to b from the lower bound of x to upper bound of x. 1042 00:50:58,620 --> 00:51:03,209 So this is dx, dy, f of xy. 1043 00:51:03,209 --> 00:51:04,750 This is the kind of integral that you 1044 00:51:04,750 --> 00:51:06,489 might need to do sometimes. 1045 00:51:06,489 --> 00:51:08,530 OK, because you might have two variables that you 1046 00:51:08,530 --> 00:51:10,180 want to integrate over. 1047 00:51:10,180 --> 00:51:15,070 And so how we do this, the best way to think about it, I think, 1048 00:51:15,070 --> 00:51:22,362 is that f of x is equal to the integral of l 1049 00:51:22,362 --> 00:51:27,695 of x u of x dy f of xy. 1050 00:51:30,600 --> 00:51:32,524 And then, if I knew what that was, 1051 00:51:32,524 --> 00:51:34,690 I would evaluate that say with a Gaussian quadrature 1052 00:51:34,690 --> 00:51:35,231 or something. 1053 00:51:35,231 --> 00:51:40,070 So I'd say that this i 2d is approximately equal to the sum 1054 00:51:40,070 --> 00:51:41,540 of some weights-- 1055 00:51:41,540 --> 00:51:44,150 I'll label them x just to remind you where they came from-- 1056 00:51:46,680 --> 00:51:52,190 times f of xn. 1057 00:51:52,190 --> 00:51:53,910 And I would use, say, Gaussian quadrature 1058 00:51:53,910 --> 00:51:56,243 to figure out what the best value of the x's and the y's 1059 00:51:56,243 --> 00:51:57,110 are to use. 1060 00:51:57,110 --> 00:51:59,780 And then I would actually evaluate this guy 1061 00:51:59,780 --> 00:52:01,730 with a quadrature two. 1062 00:52:01,730 --> 00:52:08,800 So I'd say f of xk or xn is equal to the integral from l 1063 00:52:08,800 --> 00:52:17,100 of xn to u of xn dy f of xn y. 1064 00:52:17,100 --> 00:52:20,260 And so this itself, I would give it another quadrature. 1065 00:52:20,260 --> 00:52:24,390 This is going to be equal some other weights. 1066 00:52:24,390 --> 00:52:32,495 nj f of xn yj. 1067 00:52:32,495 --> 00:52:34,910 Is that all right? 1068 00:52:34,910 --> 00:52:37,130 And so this is nested inside this. 1069 00:52:37,130 --> 00:52:40,350 So now, I have a double loop. 1070 00:52:40,350 --> 00:52:43,240 So if I took me, I don't know, 15 points, 1071 00:52:43,240 --> 00:52:45,740 15 function evaluations to get a good number for my integral 1072 00:52:45,740 --> 00:52:47,930 in one dimension, then it might take 1073 00:52:47,930 --> 00:52:53,240 me 15 squared function evaluations to get 1074 00:52:53,240 --> 00:52:55,240 a pretty good integral for two. 1075 00:52:55,240 --> 00:52:57,260 15 might be a little optimistic, so maybe 100 1076 00:52:57,260 --> 00:52:59,210 is more like if you really want a good number. 1077 00:52:59,210 --> 00:53:02,577 So I need 101 dimension, 100 in the second dimension squared. 1078 00:53:02,577 --> 00:53:04,910 Because I subdivided the integral into six equal pieces, 1079 00:53:04,910 --> 00:53:06,020 say. 1080 00:53:06,020 --> 00:53:08,030 And then now, so this is OK. 1081 00:53:08,030 --> 00:53:08,980 So you can do this. 1082 00:53:08,980 --> 00:53:10,640 It's only 10,000 function evaluations. 1083 00:53:10,640 --> 00:53:13,630 No problem, you've got a good gigahertz computer. 1084 00:53:13,630 --> 00:53:15,260 Now, in reality, in a lot of them, 1085 00:53:15,260 --> 00:53:17,100 we do have a lot more dimensions than that. 1086 00:53:17,100 --> 00:53:20,420 So in my life, I do a lot of electronic structure integrals. 1087 00:53:20,420 --> 00:53:25,937 We have integrals that are dx 1d y1 d z1 for the first electron 1088 00:53:25,937 --> 00:53:27,020 interacting with a second. 1089 00:53:30,670 --> 00:53:34,820 And then something that I'm trying to integrate. 1090 00:53:34,820 --> 00:53:37,920 It depends on, say, the distance between electron one 1091 00:53:37,920 --> 00:53:38,730 and electron two. 1092 00:53:41,250 --> 00:53:46,386 OK, so this is actually six integrals. 1093 00:53:46,386 --> 00:53:48,720 And so instead of being a 2 in the exponent, 1094 00:53:48,720 --> 00:53:51,700 this is going to be like a 6 in the exponent. 1095 00:53:51,700 --> 00:53:54,350 So now, this is 10 to the 12th. 1096 00:53:54,350 --> 00:53:57,270 That means that just evaluating one integral this way 1097 00:53:57,270 --> 00:53:59,740 might take me 100 seconds if I can evaluate 1098 00:53:59,740 --> 00:54:02,026 the functions like that. 1099 00:54:02,026 --> 00:54:04,150 But in fact, the functions I want to evaluate here, 1100 00:54:04,150 --> 00:54:06,100 these guys will take multiple operations 1101 00:54:06,100 --> 00:54:07,910 to evaluate the function. 1102 00:54:07,910 --> 00:54:09,790 And so this turns out I can't do. 1103 00:54:09,790 --> 00:54:12,980 I can't even do one integral this way with its six integrals 1104 00:54:12,980 --> 00:54:14,270 deep. 1105 00:54:14,270 --> 00:54:18,490 And so this is the general problem 1106 00:54:18,490 --> 00:54:19,834 of multidimensional integration. 1107 00:54:19,834 --> 00:54:21,500 It's called the curse of dimensionality. 1108 00:54:21,500 --> 00:54:24,680 Each time you add one more integral, you're nested there. 1109 00:54:24,680 --> 00:54:28,220 All of a sudden, the effort goes up tremendously and more, 1110 00:54:28,220 --> 00:54:30,410 a couple of orders of magnitude more. 1111 00:54:30,410 --> 00:54:33,050 And so we'll come back later in the class 1112 00:54:33,050 --> 00:54:35,150 about how to deal with cases like this 1113 00:54:35,150 --> 00:54:37,070 when you have high dimensionality cases to try 1114 00:54:37,070 --> 00:54:38,390 to get better scaling. 1115 00:54:38,390 --> 00:54:40,240 Instead of having it being exponential, 1116 00:54:40,240 --> 00:54:47,060 it was really 100 to the n right now, 100 to the n dimensions. 1117 00:54:47,060 --> 00:54:50,070 [? Can we ?] figure out some better way to do it. 1118 00:54:50,070 --> 00:54:52,330 All right, thank you.