 Now we want to explore, see we have seen, we do a demo, we can see that as the number of grid points increases that the number of iterations increases. So for given CPU time per grid point per iteration, the number of iterations increases in a very rapid fashion because the spectral norm is so large, it is getting closer and closer to 1, right, it is going to take forever to essentially come, remember we got a geometric series, right. So if you say that xn plus 1 is rho j times xn or en plus 1 is rho j times en, right and you want e to go to 0, the fixed point is 0. If this is very close to 1, if it is 1 of course it is not going to get there, if it is very close to 1 it is going to take forever for the error to go to 0, is that fine, okay. So it is a problem, that is a problem, right. So we have to think of ways by which are there ways by which we can make the code or the iterations go faster. So given that the cost the CPU time per grid point per iteration is the same, you cannot change that. You may not mind doing a little more work that is increasing that CPU time per grid point per iteration if we are able to converge to the solution faster, okay. So this is only Jacobi, I am not really going to sit down and work through this for Gauss-Seidel, but we did Gauss-Seidel and it does turn out that Gauss-Seidel is marginally faster than Jacobi, right. But the problem is that if rho Jacobi is 0.99 and rho Gauss-Seidel is 0.98, right, okay. You in fact this is twice as fast as that, right, it is ironical 0.98, 0.99 into 0.99 gives me 0.98 essentially, right. So this is in fact twice as fast as that, but in an absolute sense it is no use to me, right. It is only 0.98. I want something that is much smaller, right, I want something that is much smaller. So we have to think of, we have to come up with ways by which we can do this. What are possible ways by which we can make the program run faster, give me suggestions. What are possible ways by which we can make the program run faster or iterations converge faster? What are the things that we have done? We have an initial guess, so you could get an improved initial guess, right, we get an improved initial guess. Clearly, see if you manage to guess the solution, you get the answer in one iteration. So clearly an improved initial guess is something that is worthwhile, higher grade point. So one possibility is that you somehow use, you use the initial, make an initial guess which is a poor guess, 0 everywhere on a course grid, what we will call it, you call it a course grid, right. And use that course, use that as an initial guess for, so you use, you take the 10 by 10, you take the solution on the 10 by 10 as an initial guess for the 20 by 20, right. In fact, you can imagine that you use 10 by 10 as an initial guess for 20 by 20 and then use 20 by 20 as an initial guess for 40 by 40, right. So you could have a hierarchy of grids, we look at this at the later in the semester, there is a scheme called multigrid method which sort of exploits this, but as far as you are concerned, a good initial guess would get us to the answer faster, okay. So there seems a possibility there, but as I said we will do that a little later, what else, anything else that we can do, okay. So one thing that we, one thing that we will try, see all we have, all we have, if you think about it, for Gauss-Seidel, all we have is Vn plus 1 PQ is 0.25 Vn P plus 1 Q, you have seen me write this so many times now, Vn P minus 1 Qn plus 1 plus Vp Q plus 1 N plus 1, right. So this is Gauss-Seidel, this is Gauss-Seidel. Now this is all we have, we have phi, the latest phi and we have the prior phi's and the improved we are using, you know, what more can we do. So one of the things that we can try, right. So this is sort of like just try it out because this is all we have to play with is instead of calling this Vn plus 1, I will call this phi star and I will write Vn plus 1 at PQ, in fact is Vn PQ plus phi star PQ and I want to pick, I want to take a linear combination or what is called a convex combination of these two in order to form that, right. So I will pick a parameter omega and when omega equals 1, I want to get back Gauss-Seidel, okay. When omega equals 1, I want to get back Gauss-Seidel. So omega equals 1 means that this should be omega, that should be 1 minus omega, omega equals 1 will get me back the original equation, okay. And omega equals 0, there is no progress, okay. So now we have done a linear combination. You can say well what is the big deal, what have we gained by this, right. What have we gained by this process? It turns out that there is an optimal value of omega for which convergence is greatly improved, okay. There is an omega value for which it is greatly improved. There is an omega equals omega opt for which the convergence of this iteration scheme, right is faster than the original Gauss-Seidel, okay. It is faster than the original Gauss-Seidel, is that fine. Remember that this was called successive relaxation, Gauss-Seidel by itself is called successive relaxation. You can try various values of omega. For Laplace's equation it turns out that omega believe it or not is greater than 1, optimal value of omega is greater than 1, okay. And for that reason it is called successive over relaxation or SOR, successive over relaxation. So instead of successive relaxation which is what Gauss-Seidel was, right. Because omega is greater than 1 it is like over relaxation, right successive over relaxation, is that fine. So we have introduced this omega. Right now you have to take my word for it that it actually improves. It actually improves the convergence rate. The, as I said we are not going to go through this whole analysis for it to show that it actually improves the convergence rate. We will do it in a more empirical fashion. You can actually write a program and try it out. And the question would be, right, if you do not have an expression for omega for the optimal omega. As it turns out for Laplace's equation on a unit square with a uniform grid you can actually come up with an expression for an optimal omega, okay. Right? There are enough books out there that will give you, look at anything on matrix iterative analysis by RS Varga comes to mind. There are lots of books out there that will give you expressions that they have already derived expressions for an optimal omega, fine. But these things work for a unit square uniform distribution. What do we do when we are going to actually solve problems from fluid flow equations or something that is quite complex? There is no expression. So what you have to do is you have to get the, you have to, you have to find the optimal omega by experimentation. You have to hunt for it. Is that fine? The advantage is once you have found it, once you have found it, you are able to do production runs so to speak. You can make multiple runs with that optimal omega, right. So how do we find this optimal omega? How do we do this? There are two possible ways, two possible iterations. One of them may be better than the other. We will see that. There are two possible ways that you can do it. So pick values of omega and right now without an explanation I will say pick values of omega between 0 and 2.0. We look at where this 0 and 2.0 comes from. Pick values of omega between 0 and 2, right. Say at equivalent rolls 0.1, 0.2, 0.3, 0.4. Pick values of omega between 0 and 2. And then you can do SOR, perform SOR iterations, right. Perform SOR iterations. Now comes the critical thing. Perform SOR iterations till, right. There are two possibilities. Perform SOR iterations till En less than say 0.1. That is one possibility. Some predetermined value, okay. And then you can ask the question, so this N will be a function of omega. That is what we are saying. That N will be a function of omega. You will say allow it to drop by 0.1 and this N will be a function of omega. And what we hope is that if I get omega versus N omega versus omega, we hope that we get something like this and we know where there is a minimum, okay. There is a danger here though. I am always wary of this setting of value 0.1. What is the danger? Yeah, what if there is a value of omega for which it does not converge? You keep on running and it is going and it is going. You will never get to 0.1. It never drops by 0.1. See, there is always that danger. There is always that danger. Yes, I know for Laplace's equation, we just showed that the iterations are going to converge, right. But for SOR, we have not really done that. So there is always that danger. If you are solving the Navier-Stokes equations or whatever, there is always the danger that you try to get that 0.1 drop and it just does not happen, that the program runs forever, right. So in all of this kind of exploratory work, I am citing this example. I am citing this because this is a very important point that I want to make. For all exploratory work, you always bound the amount of computation that you are willing to perform. You always limit the amount of computation you are willing to perform, right. So you do not say, we do not say till en is less than 0.1. That is not the way to do it, right. What you do is you say we do say 10 iterations or 100 iterations. Then you have limited the number of iterations that you are basically going to say, I am willing to do this much and no more. In which case you can then plot what would you get? What would you plot? Number of iterations is fixed, right. So this would be en versus omega, right. And now you may get, again we hope that you get something of that sort. We do not know. We hope that you get something of that sort. We do not know. Is that fine? Okay. So we can run, we can run multiple. So you can take 10 values. For example, what I would do is I would take 10 values. I would then turn around and say that let me explore for the optimal omega. Having got this, right, having computed this graph, let me explore for the optimal omega in this gap, right. And it may turn out that I actually get something like this. I deliberately shown it sort of shifting to one side. Why have I shown it shifting towards one side? We will see whether this actually happens. This is those situations that they talk about math classes but you have not really encountered that often. This is a situation where you have non-uniform convergence. In fact, we want non-uniform convergence. We want for different omegas the convergence rate to be different, right. We are looking for a situation where the convergence is non-uniform. Specifically, we are going out and looking for it. Am I making sense? We have picked a parameter and we are hoping that I really, really, really hope that the convergence rate is non-uniform. If it is uniform, then we have not gained anything, right, okay. So we will hunt and we will try to narrow down. Now it is up to you how much effort you are willing to put in to narrow down that value. Am I making sense? And if the drop is sufficiently fast, you should be able to get an order of magnitude speed up or something like that, an order of magnitude speed up or more, fine. So that we will sort of look in a look at a demo, right and see what happens. We will try to figure out what happens. Now we have to address one issue. We have SOR. We have some mechanism and algorithm by which. So I do not want you to do this. Well, I want you to do this just for fun, right. If I say I do not want you to do this, then you should go ahead and try and do it just to see what is Ramakrishna talking about why should I do it. Just try to do this just for fun. But once the classroom, once the class, the learning process is over, you do not usually do this. What you do is you do this. You bound, you limit the amount of computation that you do in exploratory, right. You make sure that it does not, it cannot become unbounded, okay. So right now I will remove that. So we have an algorithm. You can find the optimal omega. Why is the hunt limited in this case, in this particular case? Why is it limited to 0 to 2.0? What is the reason? So let us find out, okay. Remember the matrix A, what is this matrix A? This is Af equals B. B had the boundary conditions. This is the discrete version of Laplace's equation using central differences. Nabla square phi equals 0 plus boundary conditions turned out to be right using central differences, the matrix of that form. And what was the structure of the matrix A? A was symmetric, right. A was symmetric. Remember that A was symmetric meaning A equals A transpose, okay. Now consider a function that I just casually introduced. A, well I call it x equal to phi. Q of phi is 1 half phi transpose A phi minus phi transpose B. You write this in the standard form because otherwise I will keep making this mistake of replacing writing x instead of phi. We just write this in the standard form Ax equals B, right, where x and B are vectors. So Q of x is 1 half. You could put a negative sign in front of it if you want but it does not matter. x transpose Ax minus x transpose B. This x has a little tilde in front of it to indicate it is not our coordinate. So what is Q? Given an x or given a phi, given an x or given a phi, it gives me a number, okay. What does Q? Given an x or given a phi, depending on which equation you are looking at, it gives me a number. Is that okay? Fine. It is mapping this vector into one number, into the real line, okay. It is quadratic. How do I find its minimum? How do you find the minimum of a function? Take the gradient and set it equal to 0. So what we want is grad Q and set it equal to 0, okay. We want to take grad Q and set it equal to 0. I think maybe this is better if we use index notation to do the gradient so that you do not get confused. Q of xi, i is just a subscript is one-half xj for xi aij xj minus xi. This has a summation over i and j and this has a summation over i. So the repeated index indicates that there is a summation. If the subscript is repeated, that indicates there is a summation. So I am not going to. So what is the gradient? Doh Q, doh xi is one-half aij, doh Q, doh xk, I am sorry, one-half doh xi, doh xk, aij xj plus one-half xi, aij, doh xj, doh xk minus, doh xi, doh xk bi. I hope all of you are familiar with your index notation. Doh xi, doh xk would be the matrix representation would be a unit vector. It would be an identity matrix but in index notation would be delta ik, delta ik equals 1 if i equals k equals 0 otherwise, okay. So doh x1, doh x1 is 1, doh x2 is 0, okay. Doh x1, doh x3 is 0, right. Doh x2, doh x1 is 0, doh x2, doh x2 is 0. The derivative of the thing with itself is 1. The derivative of the thing with something else is 0, okay. That is basically what it says. By itself is 1, if something else is 0, okay. And remember that repeated index means summation. So what does that mean? You can check this out. I will let you verify this for yourself. So this gives me one-half akj xj plus one-half xi ai k minus bk. X i, ai k is the same as xj, aj, I mean, you are summing over, you are going to sum over the i, right. Whether you sum over i or sum over j, it does not make a difference. So this in fact turns out to be one-half akj xj plus one-half xj ajk minus bk and you want to set this equal to 0. But because a is symmetric, akj is the same as ajk. Akj is the same as ajk because a is symmetric. That is why that symmetry is very important, right. Symmetry is very important. Is that fine? So what does that result in then? Those two-halves add up and you get Ax equals b. Writing it back in vector notation, you get Ax equals b. So solving Ax equals b is equivalent to minimizing q of x, right. When a is symmetric, solving Ax equals b is equivalent to minimizing q of x. Is that fine? Everyone? So we look at the problem in one dimension. It is easier to understand. We can draw figures, right. So we look at a problem in one dimension. For one-dimensional equivalent of this, so the tildas will go away. We will just use x because it is a scalar. So qx in fact will be one-half ax squared plus bx. This is q as a function of x. And what is dou q dou x or dq dx in this case because it is only dependent on one variable, Ax plus b equals 0 telling us that the minimum is that corresponds to minus b by a. From your understanding of quadratic equations, you already know that to be a fact. Let us graph this. Now in general, ax squared plus bx is not symmetric about the, figure is not symmetric about this. But I will allow you, you go back and look at your quadratic equations and you will see that what I am saying is valid in general, right. That is if I were to draw a quadratic, the minimum is here. That is minus b by a, right. Minimum is there minus b by a. Of course if I were to translate the coordinate system by some c, I would get two roots about that minimum. And independent of the orientation of the quadratic, those two roots will be symmetric about that minimum. Is that fine? I want you to go back and check to make sure that what I am saying that you agree with. So if this is your xn, if that is your xn, because this is a scalar equation, life is very easy. After all that is what you are doing. You are taking a derivative and setting it equal to 0 and you will get this point immediately. That is what you are doing. You are solving for that point. This is xn plus 1. What I am saying is I call it x star, okay. I am calling it x star. And I am going to say that xn plus 1 is 1 minus omega times xn plus omega times x star. Is that fine? If I subtract out xn from here from this equation, I get xn plus 1 minus xn, which is delta xn, which is the correction that I need to add to my xn to get the new iteration, is in fact omega times delta x star, okay. So you just look at this. This is delta x star. That is delta x star. When omega is 1, this is the value that you get. When omega is 0, delta x star is 0 and what is the value of q that you get? This is q. What is the value of q that you get? You get the original value. When omega is 2, you swing over to the other side, symmetrically over to the other side and the q value is the same. For any omega between the values 0 and 2, you will get q values. If you choose omega depending on what omega value that you choose, you will get a value in between of x and the q will decrease. So these iterations for omega values between 0 and 2 will generate a sequence of qs, q1, q2, q3, so on. It will generate a sequence of qs which are decreasing and they are bounded from below. So you know that there is a limit. You know there is a limit and you are going to reach that limit. Does that make sense? As long as you pick omega values between 0 and 2, if you pick omega 2.1, it is going to take you there. The q is increasing. If you pick omega less than 0, you are going to go here which is also increasing. If you pick omega 0, you are not going to shift from this point and if you pick omega 2, you are going to oscillate from one point to the other point. So I would suggest that you try this out. You will see how this goes. Are there any questions? Thank you.