Lecture 61

In our last lecture, we have talked about trapezoidal rule for numerical integration. In numerical integration, we start from a set of discrete points x i, y i, where i goes from 0 to n, that means n plus 1 discrete points, and we try to find out the definite integral x equals x naught to x n, y dx. The general procedure for it, as we have discussed in last lecture, is that first we have to fit a polynomial y through these discrete points, and to do so, here I am showing in Lagrangian form, it would be a product of x minus x j divided by x i minus x j multiplied by f of x i. Once we fit the polynomial, we can now use that polynomial equation to derive that integral y dx. When I apply this x naught and x n to limits, I get an expression of the form of w i f of x i. These w i's are essentially the integration of this x minus x j, x i minus x j product of them, and once I apply the limit x naught to x n, it will become a some sort of a real number, it's a constant. So, essentially at the end, the integration is pretty much a linear combination of f of x i, where this f of x i is known here to this y i. This general idea that we have applied in trapezoidal rule is still applicable for Simpson's rule. The only difference here that you will see is that this weight w i will be different in Simpson's rule as compared to the trapezoidal rule. And in fact, in most numerical scheme, we will look for different types of weights so that we can reduce the error in numerical integration. So, let's see how we do so in Simpson's one third rule. Let us say I have a function here y equals f of x, and I have selected two points x naught y naught to x 1 y 1. And you remember that we can find this integral, this is a definite integral going from x equals x naught to x equals x 1. And I am integrating y, which is a function of x. And the general approach of course to fit a polynomial and then integrate. To fit a polynomial we use this Lagrangian form l i f x i. And then once we integrate, we get the form of the summation i equals 0 to n w i f x i as we have just discussed. So, the point now becomes that how do we fit this polynomial. And you have seen that the simplest case of fitting this polynomial is a linear fit. When we do a linear fit, the corresponding integration rule that we have developed in our last lecture is known as trapezoidal rule. We use this term trapezoidal because in linear fit, as I fit this using a straight line, the integration becomes area of this trapezium. Because of this trapezium shape, the rule is known as trapezoidal rule. You can see here clearly this integration y dx x naught to x 1 is x 1 minus x naught by 2 multiplied by y naught plus y 1. This is essentially area of this trapezium over here. This has an error of the order of h to the power 3 locally. But if I add many such terms and make a composite trapezoidal rule, then this error will be of the order of h to the power 2. So, that will be a second order accurate scheme. Locally, this is a third order accurate scheme. Now, if I make this composite, then I can go from x naught to x n and the equation becomes h by 2 multiplied by y naught plus 2 y 1 plus 2 y 2 up to plus 2 y n minus 1 plus y n and as I said that the order of error will be h square. Now, we have linear fit and there we are going from h 3 h 2 this sort of order of error. We hope that if we have a better fit and not linear, which tells us that if we have for instance, a quadratic fit, then hopefully, we will be able to have a much better error estimation. And when we have a quadratic fit, we say that this is Simpson's one third rule. Let us say how we do that. Okay, so let us say I have this y equals f x and I am going from x naught y naught to x 2 y 2 between this I want to do an integration that is y dx x naught to x. If I use trapezoidal rule, then I fit this straight line between these two points and then I try to find out the area of this trapezium here. That should be the area should be equal to our estimated integral. When we do so, of course, you can say that I miss out a lot of space here, this much I am missing out. So, I have lots of error now. I want to reduce the error. To reduce the error, I want to fit a higher degree polynomial. To fit a higher degree polynomial, we need more points. So, let us say that I know another point x 1 y 1 between this x naught y naught and x 2 y 2. And let us say this is my point x 1 y 1. If I fit these three points here, x naught y naught x 1 y 1 and x 2 y 2 using a quadratic function or degree two polynomial, and then I find this area over here. This is not a trapezium because I can now see that over the trapezium we are also capturing a little more. So, we can clearly we can see that this area will be closer to the original area of the curve, which means that hopefully we will go closer to the analytical value of the integration. Now, let us see that how we will do so. We once again start from the Lagrangian polynomial interpolation between three points. So, I have L i f of x i, i goes from 0 to n, and my n here is 2 because I am going from x equal to x naught, x equal to x 1, x equal to x 2, three points. So, my y will be L naught f x naught plus L 1 f x 1 plus L 2 f x 2. Recall that our L i is a product of j equal to 0 to n, x minus x j divided by x i minus x j, and j cannot be equal to i. With this, I see that n in this case is 2. So, I replace that by 2. The remaining expression remains just the same. So, let us now find L naught. So, in L naught, I have j equals 0 to 2, i equals 0. So, I put a 0 here, x naught. Since i equals 0, j cannot be 0. So, I have j equals 1 and j equals 2. So, I will have two terms in this product, x minus x j, x naught minus x j, j can be 1 and 2, which means x minus x 1, x naught minus x 1, x minus x 2, x naught minus x 2. This product is my L naught, which is a function of x. Similarly, L 1, here I have j going from 0 to 2, j is not equal to 1. I have x minus x j divided by x 1 minus x j. Since j is not equal to 1, j could be either 0 or 2. So, I have x minus x naught divided by x 1 minus x naught multiplied by x minus x 2 multiplied by x 1 divided by x 1 minus x 2. This gives you my L 1. And now finally, L 2. So, this is j equal to 0 to 2, j is not equal to i and i here is 2. So, j cannot be equal to 2. So, j could be either 0 or 1. And I have x minus x j, x 2 minus x j, which means I have x minus x naught, x 2 minus x naught and x 2 minus x 1 in the second term. So, this gives me my expression of L 2. So, now, I have the equation of y, all I have to do now is to integrate this equation. So, now, see that this particular equation, if I integrate once again, f x naught comes out of the integration sign, inside I will have L naught and dx. Similarly, in second term, f x 1, which is a constant comes out of integration sign, inside I have L 1 dx. Then in third term f of x 2 comes out of integration sign and inside the integral I have L 2 dx. Now, for simplicity, we assume that this point x 1 y 1, this x 1 is just in the middle of x naught and x 2. So, if it is so, then I can say x 2 minus x 1 is same as x 1 minus x naught and that is what is our step size now, h. So, it is not, although in this figure I am comparing like trapezoidal and Simpson, but in case of trapezoidal, the step size will be twice of h in Simpson, this is in between x 1 is just in between. So, this is simply half of trapezoidal. Now, once we take this, simplify this, because my L naught, L 1 and L 2 are slightly complicated expression, doing the integral is tedious, but once you simplify with these assumptions that x 2 minus x 1 equals x 1 minus x naught equals h, then this integral becomes slightly easier. If you conduct this integral, which the detail I am leaving for your exercise, you will see that the expression will come as h by 3 f of x naught plus 4 f x 1 plus f of x 2. Because of this 1 by 3 terms, we call this Simpson's 1 third rule. Similarly, to this 1 third rule, if we take like instead of 3 points, if we take 4 points, we can get higher order and we can get many different Simpson's rule like Simpson's 3 8 rule and many other higher order Simpson's rules. But usually it is not necessary for most cases to go beyond quadratic. If I now for a realistic integration process between x naught to x n, I divide the entire area in many, many different segments just like what we did in case of composite trapezoidal rule. If I do the same thing here, then as you can understand that I will have this f x naught plus 4 f x 1 plus f x 2, then the next integral will come h by 3 f x 2 plus 4 f x 3 plus f x 4 this type. So, this f x 2 will be added together and thus I will have 2 of this term. So, I have y naught plus 4 y 1 plus 2 y 2 plus 4 y 3 plus 2 y 4, etc, etc. This is how it will go. And thus at the end, I will have 4 of y n minus 1 plus y n. I am not going through the detail, but you should be able to show that Simpson's rule is locally, this will be a fifth order accurate scheme. And globally, this is a fourth order accurate scheme. So, recall our trapezoidal rule was locally, it was third and globally, it was second order. So, now you can see that the order of accuracy has increased. The order of accuracy only tells us that how much the error depends on the step size or grid size. If I change the grid size, how fast the error will change. It does not talk about the actual error. So, don't be confused over it. If we reduce the step size, then the error in Simpson's rule will drop much faster than the error in trapezoidal rule. Only this much we can say from this idea of order of accuracy. So, now if we compare the trapezoidal versus Simpson, and I am writing it in a composite manner, because this is the manner we actually use, we never use two points or three points, we use many points, sometimes hundreds or thousands of points. And I am writing this here, this is h t, that is the step size in trapezoidal rule. This is the step size in Simpson's rule. And as you can see that they may not be same. So, that is why I am just putting a subscript here. And in trapezoidal rule, when step sizes are same, if they are same, then you can see that in trapezoidal rule, the order of accuracy is of h square. Here in Simpson's rule, this is h to the power 4. So, this is this is supposed to be better. We can generalize this idea that given certain point xi f of xi, the integral of x naught to xn y dx, where y equals f of x is a linear combination of f of xi s and this w s are the weight. The general idea, general expression like this is known as Newton-cotes formula. And this Newton-cotes formula, just like we have seen two examples of this trapezoidal rule and Simpson's one third rule, there could be many, many, many, in fact, infinitely many different Newton-cotes formula. The idea is that given certain xi and f xi, I have to figure out certain weight so that I can get the value of the definite integral to my desired level of accuracy. Now, it is generally expected that if we go for higher order polynomial fitting, just like from trapezoidal to Simpson's, we'll probably get more accurate result. The second point is that the grid size or step size that we have taken, which we have taken uniform in our example, but it is not necessary for them to be uniform. It could be non-uniform depending on the situation, depending on the function, it could be made non-uniform. Finally, if you are applying Simpson's rule, you can see that in a composite Simpson's rule, I need first 1, 2, 3 points to fit my first quadratic, then I need 1, 2, 3 points to fit my second quadratic, then I need 1, 2, 3 points to fit my third quadratic. So, essentially I can have 1, 2, 3 points or 1, 2, 3, 4, 5 points or 1, 2, 3, 4, 5, 6, 7 points like that. So, anytime I want to use Simpson's, I need odd number of data points. If you have even number, then Simpson's one-third rule you cannot use. You must have odd numbers. Trapezoidal, of course, does not have any such restriction, any even or odd number, you can still use it. Let us now close this with a small exercise. You repeat the last problem that you did using Simpson's one-third rule for 0 to pi sin x dx. You repeat this exercise for same number of n, that is 10, 50, 100, and you check it out that which of these two methods, trapezoidal versus Simpson's, which one of them is giving less error and how those errors are changing in these two methods as I change the value of n. So, that's all for this lecture. Thank you and I'll see you in next class.