Session 16: ODE-IVP and Numerical Integration 4

• Download English-US transcript (PDF)

The following content is provided under a Creative Commons license. Your support will help MIT OpenCourseWare continue to offer high quality educational resources for free. To make a donation or to view additional materials from hundreds of MIT courses, visit MIT OpenCourseWare at ocw.mit.edu.

PROFESSOR: All right, I think we're gathered. So let's get going. We were going to DAEs today, but we decided to postpone for one day and do a little catch up on different things that you should know. And I thought I'd first back up a minute and recall what you've learned so far. So we started out talking a lot about linear algebra and, particularly, the problems of this type where you're given a matrix and you're given a vector, and you want to figure out what the x is that makes this equation true.

And you've become extremely familiar with the convenient Matlab backslash function. You've probably used this a lot of times. And if you remember, this is one of the only problems where we often have a unique solution. We know a solution exists. And we can do a finite number of steps in order to get the solution. And actually, the method we have is really good. So it almost always gives us the solution and to machine accuracy. So this is a brilliant method. And so everything is based on this in our whole field, actually.

So we learned this first. We learned it takes n cubed operations in order to get the solution, where n is the size of the vectors. And when we do the solution, the x we get for backslash, invariably, it either is a solution, solves this equation to machine accuracy, in the sense that you multiply m times x, you get b to as many digits as you can read.

However, we talked about how m could be poorly conditioned so that there might be multiple x's that would solve this equation to machine accuracy. And so that's the problem with ill conditioning. So it's not that it doesn't solve the equation. It's that we're not really sure that's the right x. You could get a bunch of different x's that, basically, solve it as well as we can tell. All right, so do you guys remember this?

All right, and so that also implies that if you make a very small change in b, you might get a gigantic change in x. Because even for one b, you have a quite a range of x's that could work, all right. Then we try to solve problems like this. Right, systems of equations, not just linear ones, but nonlinear ones as well. And what we decided to do as one of our best methods is Newton-Raphson, where we keep on doing j backslash f [INAUDIBLE] negative f.

Right, so all we did in this method, a lot of times, is we just kept on calling this method over and over again. All right, so this took order of n cubes. This thing is going to take order of n, I don't know what you call it, steps or something. Maybe if you guys have a word for it. What do you call the iterations? Iterations. So it takes the number of iterations times n cubed.

Now, if you really can make Newton-Raphson work, and you have a really great initial guess, the number of iterations is really tiny. Because Newton-Raphson is such a fantastically great convergence. However, we're usually not such good guessers. And so it might take a lot of iterations. And in fact, we use the dogleg method to solve this as we've discussed about where you mix up different methods. And you have a trust region and all this baloney. And so it might actually be quite a lot of iterations to get there.

And then we went and tried to solve this problem. And what we decided to do for both of these was almost do the same thing. So I said, well, we can really solve this one by instead solving this. Right, that's really what we ended up doing.

That's what happens inside fsolve is it converts the problem of finding f of x equals 0 into a minimization problem to try to make the norm of f as small as possible. Because the answer will be when the norm is zero. And if f of x equals 0 has a solution, then the global minimum here is a solution, as long as your initial guess is good enough that you're in the valley that leads to the true minimum. Then this will work. And that's actually what fsolve does.

So these two guys get kind of coupled together, the minimization problem and the root-finding problem, or the problem of solving systems of non-linear equations. And in both of them, we're going to call the first method a million times. OK, so that's the way we're working so far. And actually, when you do fmincon, which was the next level, then you actually end up looking for saddle points, which, again, we look for as minimizing the norm of a gradient squared, which gets it back into this kind of form. And so they're all the same thing.

And you can, again, work out estimates of how much effort it takes. And it's something times n cubed, Typically, if most your work is this. Now, I should make a comment that, in practice, that's not always true. So in practice, sometimes, although the highest order term is n cubed, the linear algebra programs that do this are so great that the prefactor is tiny.

So in order to get the limit where the n cubed term was actually dominate the CPU time, you had to have a gigantic n. And in fact, the limit, it's often the effort to construct the matrix j. And so it really goes as n squared. So although this nominally is like n cubed, in practice, a lot of times, it's n squared. And I say it for these other methods.

And then a lot of them get into how good the initial guesses are. Because the number of iterations can get crazy. And in fact, maybe, it never converges if your initial guess is bad. So that can be a really critical thing. [INAUDIBLE] is better. Because you can't change n I mean, how many unknowns you've got is how many unknowns you've got. But you can change the number iterations by, say, more clever initial guessing. And that was the whole concept about continuation and homotopy and all this kind of stuff was trying to help reduce this number.

Last time, we talked about implicitly solving ODEs. And so we had an update formula. So we were solving an ODE. The initial value problem, y of-- Right, and we said that, first, we went through the explicit methods. And they were perfectly straightforward to do. And you don't have to solve any systems of nonlinear equations to do them. And so that was great.

But then we pointed out that, a lot of times, they're numerically unstable. And so then you have to instead use implicit methods. But the implicit methods generally work this way. They have your new value of y. What you're trying to compute is your current value of y plus delta t sums g. It depends on the size of delta t and y new. Maybe y old also.

And these g's are things like the Crank-Nicolson, which you're going to work on in the homework, and we show it in class too. And so there's different update formulas of people suggesting how to do this. The trick of this is that the y new is inside the function. So therefore, this is really another problem like we had before. You can really rewrite this that way by just putting this on the one side. Put them all on one side of zero. And now, I'm trying to figure out what y new is. And I'm going to solve that using this method. And I'm going to solve that by using this method.

All right, so now, how is the scaling going to go? At each time step, I have to solve one of these guys. How much effort did that take? We did that, right. Here it is, n [? after ?] times n squared. And then how many times do I have to do it?

AUDIENCE: It's h times 7 times [INAUDIBLE].

PROFESSOR: Right it depends on how many times steps I have [INAUDIBLE] my delta t. So it's something like t final minus t 0 over delta t. That's the number of times steps. But I have constant delta t. Times the number of iterations that it takes to do the nonlinear solve times n squared, which is a cost of forming the Jacobian matrix.

Now, this can be pretty large. So as we talked about, this might be as big 10 to the eighth in a bad case. Certainly be less than 10 to the eighth because otherwise, you'll have numerical problems. Maybe, it's really going to be 10 to the sixth for a real problem. So you have a million time steps. At each times step, you have to solve the number of iter-- the solve costs this much. Number of iterations might be what? How many iterations do you guys take to solve nonlinear equations? You solve a lot of them now. How many iterations does it take?

AUDIENCE: 10.

PROFESSOR: 10, does that 10 sound right?

AUDIENCE: Yeah, [INAUDIBLE].

PROFESSOR: You've got an initial guess 10. If you have a medium initial guess, maybe 100. It it's more than 100, you're not going to converge, right?

AUDIENCE: Yeah.

PROFESSOR: Yeah, OK. So this is order of 10 maybe. So the million 10 and then n squared just depends on how big your problem is. So in my group, the problems we typically do, n is approximately 200. So this is like 200 squared. And then you can see that the number of iterations, and this is maybe a little bit overestimated. Maybe I can get away with 10 to the fifth, something like that. But you can see it starts to get pretty big.

But for writing it, as you guys are writing software, it's not that big a deal. So I can assign it as part 1 of a three-part homework problem. It's the right ODE solver. It solves implicit equations. And you can write it. And the reason you can is because you've already written methods that solve this. Actually, Matlab did it for you. But you write it yourself. I assigned to the [? 10:10 ?] students. So you write the [? backsub ?] solution to solve the equations.

And we we've already wrestled this. And I think you guys wrote one of these yourself. And also, we can use fsolve, which is just another little bit better implementation of the same thing. And then now, you write your next level. So you can see you're doing a composite process, where you're doing something at the top level. It's calling something at a lower level. That's calling something at a lower level.

And then what you really have to be concerned about if you're worried about CPU time or whether the thing is going to solve before the homework is due, is how big does this get. How many operations are you really asking the computer to do? Now, remember the computer is fast. Right, so it can do a gigahertz or something. So you can do like 10 to the eighth somethings per second.

And what should I say about this? I have not included, but I should that it's-- well, I'm assuming here that the function evaluation is sort of order n of how many variables they have. OK, sometimes, function evaluations are much more difficult than that, primarily because, sometimes, they're actually including further nested stuff inside. So your f is actually coming from some really complicated calculation itself. And then you have another big giant factor of how much it costs to do f. But on this right here, I'm just assuming it goes as n.

So actually, let's just multiply this out. So this is what, 10 to the fourth, 10 to the 10th. So this is 10 to the 10th. 10 to the 10th kind of flops. But your computer does 4 times 10 to the ninth per second. Is that right? Is that [? what they ?] told us, four gigahertz? Yeah, so is this really not bad with just a few seconds on your computer. So that's why you can do these homework problems, and you still get them done. And if you have even [? porous ?] limitations, it takes 1,000 times longer. Still, 1,000 [INAUDIBLE] take 20 minutes on your computer.

So you're OK, but you can see we're starting to get up there. And it's not that hard to really make it a big problem if you do something wrong. And just a notation thing, a lot of these things we're doing are writing software that takes as input not just numbers or vectors. So those are things that we call functions. Right, things that take a number of vectors and input and give you, say, a function, a number of vectors and output.

But oftentimes, we actually want to take the names of functions as inputs. So for example, the root-finding problem, it takes as input what is f. Right, you have to tell it the name of a function for fsolve to be able to call it. So fsolve is an object that takes as input the name of a function. And things that do that are called functionals

So for example, this one is x is equal to this root-finding functional, fsolve. And I'm going to write with square brackets. It depends on f. That's it, right? [INAUDIBLE] f. So you give it the function f. And if fsolve's good, it should be able to tell you what the x is. Now, in reality, we help it out. So we give it x 0 [? 2 ?] just give it a better chance.

I write the square brackets just indicate that one of the things inside is a function, not just numbers. We do this a lot. So you guys have done this since you were kids. So derivatives Yeah, ddx of f. So we have, I don't know, for example, grad of f. That's a functional, the grad. The symbol means the operator, the functional, it operates on functions. You could write a general gradient function that takes any scalar valued function and computes the gradients with respect to the inputs.

This week, we give you one that computes a Jacobian. Give it any f of x. It gives you the Jacobian of f with respect to x. So that's a functional. So anyway, we do a lot with functionals. You've been using them all the time. I'm just trying get you to think about abstractly. And then there's a whole math about functionals. And there's a whole thing called calculus of variations. It's all about, if the objects you're worrying about are functionals, not functions, then you have another whole bunch of theorems and stuff you can do with that, which you'll probably end up doing before you get out of this place. But I'm not going to do that right now.

I would comment that one kind of functional that's particularly important to a lot of people in your research is the fact that the energy of any molecule is a functional of the electron density. And this is the basis of what's called density functional theory.

And so this is why a lot of people, you and many of you included, are going to end up using computer programs that are trying to find what the electron density shape is, things like molecular orbitals and stuff. You've probably heard these before. And then there's the functional that gives the energy of the molecule. And what you care about, mostly, is energy, because that connects to all your thermo. You're doing [? 1040, ?] it's all about the energies. And once you know all the energies, you can actually compute the entropies and all kinds of stuff too.

And so this is one application that's really important. It was a little surprising this was true and able to be proven to be true. And so the guy who did that got the Nobel Prize-- his name was Cohen-- a few years ago.

Actually, in liquid phase stuff, you work with Professor Blankschtein. They also have a version of this for properties of solution. And then where it's not the electron density, but it's actually the density of the material in the solution phase. Maybe Professor [? Swan ?] too. You guys do that? Yeah, so if you work with Professor [? Swan, ?] you might learn about that too.

So that's one where the actual functional is like a big focus of the entire project. And the key is that this is a function. The density is a function of the position. And then the functional converts that into a number. That's all page one here.

Now, I'm going to sort of change topics, but it's related. It's about the connections between numerical integration, interpolation, and basis functions. These things are all intimately tied up together. So I'll try to show you. So we try to solve a numerical integral.

And then let's look at what we're really doing. When we're doing this, we typically only do a few function evaluations. So we have a few points where we know the numbers f of tn. We've picked a few t's. We evaluated f. And from those few points, that small number of function evaluations, we want to get a good estimate of this whole integral.

Now, it's a little goofy, actually, if you think about what we're trying to do. So suppose you're trying to do an integral from negative 1 to 1. We have some function that has a point here, a point here, a point here. And just for concreteness, let's say this is 1 over 26. And this is 1. And this is 1 over 26. And I have a particular function in mind that has those three points.

Now, the first thing to keep in mind is there's an infinite number of functions that go through those three points. Now, what we've been doing is we've been fitting some function to these points and then integrating that function that we fitted. So when you're in high school, you did the rectangle rule and what the function was used was this and this. And you knew the integral of this function and the integral of this function. And you added them up. And you decided from the rectangle rule that the integral was equal to 1 and 1/26.

But then some other people in high school are a little smarter. And they said it's not so smart [INAUDIBLE]. This thing was symmetrical to begin with. And now, I'm fitting it with a function that's wildly unsymmetrical. So let's instead not evaluate those points. Let's instead evaluate the midpoints. And maybe we find out that the midpoint values are here and here. And so I would assume that the function is just a flat line. That looks not so good either, because it doesn't get through these points.

But anyway, for the two points I knew, if I only used those two points, that would be the midpoint method. And I would get the I midpoint method. And I'd get something. And I did it for this particular function I have in mind. And it was 8 over 29.

And then some other people say, well, let's use the trapezoid rule. So let's approximate the function is this over the interval, which looks kind of sensible. And if you do that, it turns out you get the same numbers as you do with the rectangle rule. So I trapezoid.

Now, let me just as a warning, if you just do the rectangle rule and then the trapezoid rule, and you got the same number both times, you might be tempted to conclude you have converged. And you have the true solution, because you did two methods of really different orders, and you got exactly the same number. And you're super happy, and you're ready to publish your paper.

But of course, the function might not look like this. So the particular function that I always had in mind was f of t is equal to 1 over 1 plus 25 t squared. OK, that was just the function I had in mind. And the way that function looks like is like that. And the real integral of that function, which you guys if you remember how do it from high school, you can actually evaluate that integral. It's nothing to do with either 8/29 or 1 in 126. It's something completely different.

Now, those of you who went to really advanced high schools, you also learned Simpson's rule too. And Simpson's rule says let's fit it to a parabola. And so that would give you something like this. And if you do Simpson's rule, I ended up getting this if I did the arithmetic right. 1 and 1/3 plus 2/78. Who knew? OK, which also has little relation to the true value of the integral.

Now, if you're smart, and you did trapezoid and then I it to Simpson's, then you say, oh my goodness, these are not really that similar. And so I'm not converged, and I'd better do something else. So I just want you to be aware when we do these rules, what we're actually doing is we're fitting our points to some continuous function that we know how to integrate. And then we do the analytical integral of that function, which is a fine thing to do.

But just be aware that it's like extrapolation. It's like we're imagining what the function is doing between these points, when we actually haven't tested it. Now, we could. We could just call more points. We could calculate. But that would cost us more function evals. And we're trying to save our computer time. So we want to finish off before the TG so we can go buy some beer. And so we have to take some choices about what to do. And we can't do an infinite number of points for sure.

All right, so in all these methods, what we're doing is we're approximating our true f with some approximate f, which is equal to a sum of some basis functions like this. So we pick some basis functions. We know how to integrate. We fit, somehow determine what these c's are. And then we decide that the integral is equal to the integral of f tilde, not the actual original f.

And then that's actually going to be the sum of ck times the integral. And we carefully chose basis functions that we don't how to do this integral; analytically. So we just plug that in, and we just have to figure out that the c's are to get this to work. OK, so that's what you really did when you're doing this trapezoid rule and Simpson's rule, and all that kind of stuff. And it's just different choices of what the basis functions were that you used.

Now, there's a lot of choices about how to do that kind of approximation to make f tildes that are sort of like f's with some limited amount of information. And it generally has to do with how many k's we have. So if k is less than n, then it's really an approximation. Because there's no way that we would expect that with just a few parameters, we could get a lot f of tn. So we would expect that f tilde of tn is not equal to f of tn, at least not at all the points. Because if n is-- if we have a lot of points and we only have a few parameters, we don't expect to be able to fit it perfectly.

So we know this is what we expect. It might happen if, by some luck, we chose a basis function which was actually the true function f, then it would actually fit perfectly. But most of the time, not. By Murphy's Law, in fact, it's going to be completely wrong. But this is not such a bad idea. We may come back to this later.

Another possibility is to choose k equals to n. And then usually, we can force f tilde of tn to equal f of tn. So at the points that we have, if we have enough parameters equal to the number of points, then often, we can get an exact fit. We can force the function to go through all the points. And when we do this, this is called interpolation.

Now, when should you do the one or the other? When should I make k as big as n? Or when should I use k less than n? A really crucial thing is whether you really believe your points. If you really believe you know f of tn to a very, very high number of significant figures, then it make sense. You know that's true, so you should force your fitting function to actually match that.

However, if you get the f of tn, for example, from an experiment or from a stochastic simulation that has some noise in it, then you might not really want to do this. Because you'll be fitting the noise in the experiment as well. And then you might want to use a smaller number of parameters to make the fits. OK, so you can do either one. You're in charge. Right, you get the problem. If somebody tells you, please give me this integral, they don't care how you do it. They just want you to come back and give them the true value of the integral to some significant figures. So it's up to you to decide what to do. So you can decide what to do.

What a lot of people do in the math world is this one, interpolation, where you're assuming that the f of tn's are known exactly. And so your function evaluating function is great. And in that case, you can write down this nice linear equation. So we're trying to make f of tn with a true function to be equal to ck phi n of t. Sorry, phi k of t. Which n in [INAUDIBLE] k is equal to n.

And I can rewrite this this way. This is like a set of numbers. So I can write as a vector. I'll call it y. And this thing has two indices on it, so it looks like a matrix. So I'll call that m. And c is my vector of unknown parameters. And wow, I know how to solve this one. So c is equal to m backslash y.

Now, for this to work and have a unique solution, I need m not to be singular. And what that means is that the columns of m have to be linearly independent. So what are the columns of m? They are phi k of t1, phi k of t2, phi k of tn. And I want that this column is not a sum of the other columns.

So this is kind of requirement for j not equal [? to k ?] This is kind of a requirement if I'm going to do this kind of procedure is I can't choose a set of basis functions that's going to have one of them be a linear combination of the other ones. Because then I'm going to get a singular matrix. My whole process is going to have a big problem right from the beginning. So I have to ensure that this is not true. Yeah, that this is true. These are not equal to some of them.

Another way to write that is to choose the the columns that phi to be orthogonal to each other. And so I can write it this way. OK, so I want them to be orthogonal to each other.

And then you can see how this kind of generalizes. If I had a lot of t's, this would sort of look like the integral of phi k of t, phi j of t is equal to 0. And so this gives the idea of orthogonal functions. So this is orthogonal functions, and my discrete point's tn. And this is the general concept of orthogonal functions.

And then when I define it this way, I need to actually decide my own a's and b's that I want to use. And if you use different ones, you actually have different rules about what's orthogonal to each other. And also, sometimes, people are fishy, and they [? spit ?] a little w of t in here just to do, because they have some special problem where this w keeps showing up or something, and they want to do that. But that's up to them. All right, we won't specify that.

But this is the general idea of making basis functions orthogonal to each other. And it's very analogous to making vectors orthogonal to each other, which is what we did here. These are actually vectors. Right, it's just as the vector gets longer and longer, it gets more and more like the continuous function, and we have more t points. And eventually, it looks like that integral as you go to infinity.

So anyway, so a lot of times, people will choose functions which are designed to be orthogonal. And one way to do it is just to go with a continuous one to start. Now, we can choose any basis functions we want. Over here, we already chose a whole bunch. You guys have done this already. I don't know if you knew you were doing it. But you chose a lot of different basis functions completely different from each other. You can choose the different ones if you want.

And one popular choice is polynomials. But it's not the only choice. And in fact, a little bit later, we're going to, in the homework problem six, I think we're going to have a problem where you have to use a different kind of basis functions. You can use any basis functions you want. And there's kind of ones that match naturally to the problem.

However, people know a lot about polynomials. And so they've decided to use polynomials a lot. And so that's the most common choice for numerical integration is polynomials. One thing is polynomials are super easy to integrate, because that was the first thing you learned, right, how to do the intervals of polynomials. And also, they're easy to evaluate. It only takes n operations to evaluate a polynomial of nth order. And there's a ton of theorems, and things are known about polynomials. We have all kinds of special properties. And we have all kinds of convenient tools for handling them. So let's go with polynomials for now.

Now, I want to warn you, though, that how you order polynomials are really not very good fitting functions or interpolating functions for a lot of problems. And in particular, for this one I showed here where this is the function. Here, that function, you can try to fit that with higher and higher order polynomials. So I can have my function. The function value at 0 is 1. So the function looks kind of reasonable. Like that, symmetrical, better than I drew it.

And you can put some number of points in here. And then you can do an interpolating polynomial that goes through all the points. And if you do it for this one with five points, you'll get something that looks like this. It goes negative, comes up like this, goes negative again. So that's one you get if you use five points across this interval. If you use 10 points, then you get one that looks like this where the peak value here is two where as the max of the original function's just one. And the original function's always positive, but they both go negative.

So even within the range of the interpolation, the polynomials can have wild swings. And so this first one was for a fifth order polynomial, or maybe fourth order, five points. And this one is for the next 10th order or 9th order polynomial. I'm not sure which one. Yes.

AUDIENCE: What points are these polynomials trying to fit? Or does it seem like [INAUDIBLE]?

PROFESSOR: So I didn't draw very well. When you guys go home, do this. Right, on Matlab it will take you a minute to solve this equation where you put in the function values here. And put in the polynomial values that [INAUDIBLE] here for the monomials, t to the first power, to the second power, to the third power, like that. And just do it, and you can get the [? pot ?] yourself.

And just to convince yourself about the problem, and this is very important to keep in mind. Because we're so used to using polynomials that think they're the solution to everything. And Excel does them great. But it's not this solution to everything. Nonetheless, we're going to plow on and keep using them for a few minutes.

So one possibility when I'm choosing this basis set is I could choose what are called monomials, which is phi k is equal to x to the k minus 1. And this one is super simple. And so a lot of times it satisfies this orthogonality thing, which I showed you before. But it's a really bad choice usually.

And what it is is that the matrix you have to solve for this interpolation, if you have evenly-spaced points and you have x to the k, it turns out that this matrix usually has a condition number that grows exponentially with the number of terms in your basis, the number of how many monomials you include. So condition numbers that grow exponentially is really bad idea. So this is not what people use. So you might be tempted to do this, and maybe you did it once in your life. I won't blame you if you did. But it's not a really smart idea.

So instead, a lot of the mathematicians of the 1700s and 1800s spent a lot of time worrying about this. And they'd end up deciding to use other polynomials that have people's names on them. So a guy named Lagrange, the same guy who did the multiplier. He suggested some polynomials. They're kind of complicated. They're given in the book.

But what they really are is I'm taking the monomial basis set, and I'm making a new basis set. So I want phi l to be equal to the sum dlk of phi k, where this phi k it the monomials. So really, you can choose any d you want. You're just taking linear combinations of the original monomial basis set to make up polynomials. And if you cleverly choose the d, you can choose it to have special properties.

So Lagrange gave a certain choice of the d that is brilliant in the sense that, when you solve the problem m times c is equal to y to try to do your interpolation, it turns out that the c's you want, ck, is equal to f of tk. So you don't have to do any work to solve what the c's are. It's just the f values directly. No chance for condition numbers, no problems, that's it. So that's the Lagrange polynomials.

Now, Newton had a different idea. So this is for Lagrange polynomials. If you choose a Lagrangian [? polynomial ?] basis, then this is it. And you do a square system. That's it. If you choose the Newton, even he got into this business too. And he made some polynomials. And he made ones that have a special property that if you solve it-- so he chose a different d. You still have to do some solve. You get the solution.

And then you decide that I don't like that interpretation. I want to put some more points in. I'm going to do a few more function evaluations. I'm going to add some more. So from the first solve, I have the values of c1, c2, c3, up to c the number of points I had initially.

And now, I want to add one more data point, or add one more f of tn, one more function evaluation. I want to figure out with c of n plus 1 is, what the next one is. It turns out that it's his polynomials. The values these c's don't change when I make this matrix a little bit bigger by adding one more y. And I add one more parameter. And I make one more row in the matrix, one more column in the matrix.

When I do that, it's a unique property that none of these other c's change. The solution is going to be the same, except for the n plus one [? level. ?] So then I can write a neat, cute little equation just for the n plus 1 one. I don't have to mess with the previous ones. And so I'm just adding a basis function that sort of polishes up my original fit. So that's a nice property too.

Whereas with the Lagrange ones, you have to recompute everything. And all the c's will change. Actually, the c's won't change. I take that back. You add points. The c's won't change, but the actual polynomials will change. The d will change. Because the d depends on the choice of the time points in a complicated way in Lagrange.

So now, you start to think, oh my goodness, I have all these choices, I can use Lagrange. I could use Newton. I could use monomials. Maybe, I should go back to bed. And you think, well, could someone just tell me what the best one is? And I'll just do that one. I don't need to learn about all this historical stuff about development of polynomial mathematics. I'd just like you to just tell me what the answer is.

So Gauss, typically, figured out what the best on is. And he said, you know, if we're going to do integrals, we'd really like to just write our integrals this way. So our approximate value to be equal to some weight functions times the function evaluations at some times. Right, that's easy to evaluate. That would be good. If I pre-computed the weight functions [? and ?] the times, then I put any function f in here, I could just do this like zip. It's one dot product of the f of tn's I'll have the solution. So I want to do that.

So what's the best possible choices of the w's and the t's to make this really be close? Well, the problem is that there's an infinite number of functions that somebody might want to integrate. I can't do them all. So he said, let's just find the one that works for as many monomials as possible to the highest order, as far out as I can go. And it turns out I have nw's. I have nt's. So I have two n parameters that I can fiddle around with. And so I should get something like two n's in my polynomials, maybe 2n minus 1.

So what I want is that I want to choose this so that for the monomial 1, and the monomial t, and the monomial t squared, that this i is actually going to be exact. So what I want is for the monomial 1, I want the summation of wn. My f for the monomial 1 is always 1. So it's times 1. I want this to equal the integral from a to b of 1dt, which is b minus a. OK, so that gives me one equation for the w's.

Then for the next one, I want the summation wn times tn to equal the integral of a to b of t. dt, which is equal to b squared over a number 2 minus a squared over 2 if I did my calculus correctly. And so there's another question. I have another question for these guys. And I keep going until I get enough equations that determine all the w's and c's. All right, so I do it for t squared. So I can wn tn squared is equal to [INAUDIBLE] b t squared dt, and so on, OK.

And so if you do this process, you get the weights and the points to use that will give you what's called the Gaussian quadrature rule. Now, the brilliant thing about this is that by this setup, it's going to be exact for integrating any polynomial up to order 2n minus 1. So if I choose three, here I'll be able to determine-- if I choose three points, I'll have six parameters. So I'll be able determine up to the 11-- 11th order polynomial would integrate exactly.

No, that's wrong. n's only 3. I have six parameters. I can't do this. 2 times 3 minus 1, what's that? 5, there we go. To the fifth order polynomial integrate exactly, I'll have an [? error ?] for sixth order polynomial.

So this is a rule. So before, if I gave you three points like here, most of you probably would use Simpson's rule. Fit this to a parabola. It has three parameters, a, ax squared plus bx plus c, right, for a second order polynomial. I would fit it. Then I'd integrate that. That would be the normal thing that a lot of people would do. I don't know if you would do that. But a lot of people would do that. So that would be that what you would get.

And what [? error ?] would you get? You'd have error that would be the scale of the width of the intervals to the fourth power. Right, because Simpson's rule exact up to the cubic. You guys remember this from high school? I don't know if you remember delta t to the fourth. Yeah, OK.

So Simpson's rule has an order of delta t to the fourth. But Gauss's method has order of delta t to the sixth. So this is for three points. In Gauss's method, you can extend it to as many points as you want. All right, so that's kind of nice. You only do three points, you get [? error ?] of the sixth order. So that seems pretty efficient. So people kept going.

And so the very, very popular one use, actually, seven points for Gauss. So you'd be able to get to the n minus 1. So you'd be exact for 13. So if you have delta t to the 14th power as your error, that's pretty good. And then what people often do is they'll add eight more points.

And they make another method that uses all these points. It's called [? Conrad's ?] method. He made another method that gives you the best possible solution if you had the 15 points, given that 7 of them are already fixed by the Gauss procedure. And so you can get the integral using Gauss's method and using [? Conrad's ?] method.

So you can compute the i using Gauss's method, the seven points, the i with the [? Conrad ?] method. I don't know how to spell his name, Actually I think it's like this. 15 points. And then you can take the difference between these two. And if these two are really similar, then you might think your integral's pretty good. Because they're pretty different methods, and if they get the same number, you're good. And people have done a lot of numerical studies in this particular case, because it's used a lot in actual practice.

And it turns out that the [? error ?] is almost always less than this quantity to the 3 Gauss power. So you have an error estimate. So you do 15 function evaluations. You get a pretty accurate value for your integral, maybe, if you can be fit well with whatever order polynomial, 13th order polynomial. 13th order polynomial, yeah.

And so any function that can be fit pretty well with the 13th order polynomial, you should get pretty good value in the integral. And you can check by comparing to the [? Conrad ?] estimate, which uses a lot more points. And if two of these guys give the same number or pretty close, then you're pretty confident that you're good. And there's even a math proof about this particular one. So this is like really popular, and it's a lot better than Simpson's rule, but it's kind of complicated.

And part of it you have to recompute these tn's and wn's for the kind of problem you have. Yeah.

AUDIENCE: If you're using 15 points already, why not just split your step size smaller? And if you did that, how would the error [INAUDIBLE] the simplest? Use a trapezoid.

PROFESSOR: Trapezoid or Simpson's rule of something like that. Yes.

AUDIENCE: How would the error with those 15 steps using trapezoid compare to [? Conrad ?] and Gauss?

PROFESSOR: So we should probably do the math. So say we did Simpson's rule, and we have delta t, and we divide it by 15 or something. Then we can figure out what this is, but it's still going to have the scaling to the fourth power. So the prefactor will be really tiny, because it will be 1/15 to the fourth power.

But the scaling is bad. So that if I decide that my integral's not good enough, and I need to do more points as I double the number of points to get better accuracy, my thing won't improve that much. Because it's still only going to use the fourth power, whereas the Gauss' one say it's going use the sixth power if I only use the three points. Where we use 15 points, is going to get scale to some 16th power or 15th power or some really high power.

So for doing sequences, the Gauss' one is a lot better. However, what people do is they usually pre-solve the system of equations. Because this system of equation does not depend. When you do it with the monomials, it doesn't depend what f is. You can figure out the optimal t's and w's ahead of time. And so they [INAUDIBLE] tables of them depending on where you want to go. People recomputed them, and then you can get them that way.

Time is running out. I just want to jump up. Now, for single integrals, I didn't have to tell you and of this. Because you could have just solved it with your ODE solver. So you can just go to ode15s or ode45, and you'll get good enough. It won't be as efficient as Gaussian. But who cares. You still get a good answer. You can report to your boss. It'll take you three minutes to call it e45 to integrate the one-dimensional integral.

The real challenge comes with the multidimensional integrals. So in a lot of problems, we have more than one variable we want to integrate over. And so if you have an integral, say, of a to b from the lower bound of x to upper bound of x. So this is dx, dy, f of xy. This is the kind of integral that you might need to do sometimes. OK, because you might have two variables that you want to integrate over.

And so how we do this, the best way to think about it, I think, is that f of x is equal to the integral of l of x u of x dy f of xy. And then, if I knew what that was, I would evaluate that say with a Gaussian quadrature or something. So I'd say that this i 2d is approximately equal to the sum of some weights-- I'll label them x just to remind you where they came from-- times f of xn. And I would use, say, Gaussian quadrature to figure out what the best value of the x's and the y's are to use.

And then I would actually evaluate this guy with a quadrature two. So I'd say f of xk or xn is equal to the integral from l of xn to u of xn dy f of xn y. And so this itself, I would give it another quadrature. This is going to be equal some other weights. nj f of xn yj. Is that all right?

And so this is nested inside this. So now, I have a double loop. So if I took me, I don't know, 15 points, 15 function evaluations to get a good number for my integral in one dimension, then it might take me 15 squared function evaluations to get a pretty good integral for two. 15 might be a little optimistic, so maybe 100 is more like if you really want a good number. So I need 101 dimension, 100 in the second dimension squared. Because I subdivided the integral into six equal pieces, say. And then now, so this is OK. So you can do this. It's only 10,000 function evaluations. No problem, you've got a good gigahertz computer.

Now, in reality, in a lot of them, we do have a lot more dimensions than that. So in my life, I do a lot of electronic structure integrals. We have integrals that are dx 1d y1 d z1 for the first electron interacting with a second. And then something that I'm trying to integrate. It depends on, say, the distance between electron one and electron two. OK, so this is actually six integrals.

And so instead of being a 2 in the exponent, this is going to be like a 6 in the exponent. So now, this is 10 to the 12th. That means that just evaluating one integral this way might take me 100 seconds if I can evaluate the functions like that. But in fact, the functions I want to evaluate here, these guys will take multiple operations to evaluate the function. And so this turns out I can't do. I can't even do one integral this way with its six integrals deep.

And so this is the general problem of multidimensional integration. It's called the curse of dimensionality. Each time you add one more integral, you're nested there. All of a sudden, the effort goes up tremendously and more, a couple of orders of magnitude more. And so we'll come back later in the class about how to deal with cases like this when you have high dimensionality cases to try to get better scaling. Instead of having it being exponential, it was really 100 to the n right now, 100 to the n dimensions. [? Can we ?] figure out some better way to do it. All right, thank you.