 Okay good morning. So I am going to get back to the demo on Laplace's equation. What we will do today is I will start with the same 5 by 5 grid that is we will have 3 interior grid points. Last time I do not know whether you realized it or not but we were doing Jacobi iteration just so that you realize that I was doing Jacobi iteration right. It was simultaneous relaxation, all the points were changed in one shot okay that is basically what was happening. We will also try to see whether we can do Gauss-Seidel okay. So we will do that on our 5 by 5 system and 5 by 5 grid basically remember gives you 3 by 3 interior grid points. So the actual system of equations that you are solving is a 9 by 9 system of equations. So we will do that with that 5 by 5 system and then maybe we will try something larger right and for the larger one I will basically use a package not package but a program that I have already written so that I do not sit down and we can do it manually I mean it does not matter I can write the program here it is not a big deal. Maybe even for the larger one I will just do it by do it by hand it is not it is just averaging right. Then if you have time today we will also do SOR I do not know whether you have tried looking at SOR right but I already have a program written for SOR for finding the optimal omega value it is not to aid in finding the optimal omega value right it does not actually hunt for the it does not do an optimization does not solve an optimization problem okay so that is as far as that goes. The other thing is if you have tried SOR right so that is remember it was over relaxation but it was still successive over relaxation it was for the successive iteration that is for the Gauss Seidel okay. My suggestion is try to find out what happens if you use a relaxation parameter with Jacobi iteration just for the fun of it just see what happens if you try right if you try using a relaxation parameter with the original Jacobi iteration algorithm that we came up with is that fine that is okay. So just try to find out what happens and we will revisit that question what I just suggested now we will revisit it a few weeks from now I will come back to that point a few weeks from now okay. So I am starting off I Python when you use it of course you do not have to do this no color I mean I am doing no color simply because the red that comes here does not come out well on the video okay that is the only reason why I am shutting off color I just wanted to show you that there is a difference between last class and this class just so that we do not run into any confusion okay so we will set it up quickly so let me just to so that is the that is the plotting utility and this package Maya v just so that you know is built on top of another package called visualization toolkit which is a C++ package so if you want to access any of these 3D plotting utilities functions right either in C or C++ you can easily write a program and tie it to visualization toolkit it is called VTK the library is called live VTK okay I will just type it out so that you know what it is I am not going to the library is called live VTK what Maya v does is it places a layer of an interface between Python and the VTK toolkit the visualization toolkit that is basically what is happening so all the visualization because I had a question in the last class about visualizing what tools we use your writing programs in C and C++ and so on right so fundamentally all the visualization programs are basically written in C C++ Fortran right and what you see here in Python is normally a Python interface to an existing program in Fortran so there are a lot of lot of packages out there that you can use okay so grace plot for instance is a front end to something called grace which has a history of its own either XM grace or XMGR whatever it is so you can you can you can try it interface directly to grace okay you can try a package that I like to use is called YGL I will just type it again because this package basically tries to emulate the silicon graphics library called GL right so which which is very easy to use so you can try I mean as I said that but there are lots of packages out there that you can do use and finally of course is a suggestion if you can put a dot on the screen and this is relevant to what we are doing in this course right so the screen that you are seeing right up front is what I am using on my laptop is 1024 x 768 dots if you can put a dot on the screen you can draw a line on the screen you understand and if you can draw a line on the screen you can plot okay is that fine okay so there is always always ways by which you can plot right there is always ways by which you can plot maybe at a maybe in a later class I will show you do a few demos of various plotting packages that you can use okay the other part as far as programming goes before I get into the demo itself because of the nature of this course you are writing all your programs in right I encourage you to write your programs in C C++ of Fortran and so on but when you eventually get to code development right so the way I do my code development is very often I start my initial code development in Python is that fine right and then as in when so if I have my program functionally correct the program is functionally correct you understand what I am saying it functions properly then I start through a process of writing either C routines Fortran routines or C++ routines to replace elements of that Python program with something that runs faster as required so as a consequence I am able to sit when I am in this process I am always able to interact through Python interact with my program always interact and if it is batch then I use my Python again to write a script to run the batch program that is the idea so in a sense using Python is a programming language to steer large programs that is the that is what you are looking for that is the ability that you want to get to write at some point in the future but right now my suggestion is if you want to try to do that do that but my suggestion is keep writing programs in C C++ Fortran till you make sure that you are totally comfortable doing okay let us get on with it. So I am going to import nump because I need an array right so I will just like last time I will create a 5x5 grid going from 0.5j I think by now the syntax is most probably I hope familiar to you 0.2 1.2 and that should create for me a 5x5 mesh okay that should create for me a 5x5 mesh those are the two things that we looked at in the last class and of course as a consequence just like I did last time I will say p is x star x – y star y okay right and we will we will set to 0 okay and of course at this font size it looks a bit funny but you can see that you do the right thing phi of – 1, 1, 1 – 1, oh sometimes that does not make sense we will do the right thing right it is what I is telling you this is what I is telling you the reason why I do this is yes I make mistakes we all make mistakes right so that I am not trying to make it I am not trying to make it perfect now you can see that all the interior points is 0 it is fine okay let me just type up because this helps me at a later point so phi of 1, – 1, 1, – 1 equals what is it today I will get a little lazy 0.25 times what we have 0, – 2 that shifted to the left plus shifted to the right shifted up or shifted down plus that shifted up fine so that in fact will give me that in fact will give me the first iteration okay is that fine everyone so we can and we can go through this iteration what if I wanted to do Gauss-Seidel instead so you can go through this iteration what I will do is I will quickly I will just quickly redo this so I reset my phi I set it equal to 0 okay because I have I need to store I need to store my error error shall I say Jacobi equals that okay fine how many iterations shall I do 10 iterations 15 iterations or let us do 100 iterations just to be I in range 100 okay I hope I have got everything right so I just copy that out see if that works whatever has to be deleted error is x squared – y star y – phi phi I that fine I will just square it I will take the sum so err.append fine possibly even take the square root of that in p.square root sum of the squares is that fine this should give me what I want and I make a mistake err.j when you get clever that is what happens err.j let me set up my phi again please bear with me err okay that is fine okay so err.j seems to have lot of zeros and way up there somewhere so obviously the suggestion for 10 or 15 iterations was pretty good pretty good suggestion and the interesting thing is that the error afterwards is 0 is identically 0 okay so it is actually possible that 2 iterates get the same it converges to the same value you need to think about what that means in this context so what about Gauss-Seidel I want to do Gauss-Seidel would give me I set phi as x squared – y squared so let me create err.j for Gauss-Seidel I set phi as x squared – y squared I set the make sure that the initial guess is right okay and then Gauss-Seidel what changes something here changes so in Gauss-Seidel I will have to I cannot do that phi I you know I cannot do that one whatever blah blah blah so I have to let me let me I have to actually iterate through the whole thing right I have to actually iterate through the whole thing so what I will do is what I will do is let me just type it up for i in range 100 yeah what do you want to do now for j in range 1, 4 for k in range, 4 you have already done this this is most probably a hassle but anyway it is okay phi of i, j equals 0.25 times phi of okay as I said this is most probably the last time I will bother you by writing this j, k let us say wake up Ramakrishna j plus 1, k plus phi of j minus 1, k yeah you will have to plus phi of j, k plus 1 plus phi of j, k minus 1 is that okay everyone and it is updated I am using the same array so it is updated immediately okay so that is as far as that loop is concerned do I need to do anything else there I do not need to do anything else there so that is for the first for loop that is for the second for loop so err equals x star x minus y star y minus phi and this copying from above as I said just bear with me I do not want to make a mistake here in p dot square root of what is that fine okay so now we can plot these two and see what they look like if you want to just see errg a lot of zeros okay now we can plot these two so from grace plot import my grace plot I create the plotting utility what do I want to plot plot two things very similar as I have always said use a log scale so go to a log scale apply I do not know why it went to 10 power minus 4 1 e minus 16 even 10 power minus 16 it seems to just okay fine and as I indicated earlier as I indicated earlier Gauss Seidel converges faster than right Gauss Seidel converges faster than Jacobi iteration and it sort of if you squint at it it does look as though it is twice as fast on a log scale the slope does look as though it is twice as fast is that fine is that okay right so I am not going to see how far down it goes it is going to suddenly drop to 0 right that was very clear so one of the things of course what I normally do is I do not just allow the program to run in the background I usually plot this graph while the program is running so that you only look at the graphical output rather than what I will do is why do not why do not we run 101 by 101 or something of that sort we will run a big one run a big one so I will recreate this so I will just basically say x, y is 101 what do you want 101 is fine 101 101 so that gives me 101 by 101 grid okay and everything else should actually work as it is because there is nothing in there that was specific to the number of grid points except for the Gauss Seidel we have to be a bit careful so what I will do is I will make this e or r j equals I will just overwrite it it does not matter or maybe we can make it 101 either way go back to the Gauss Seidel I am sorry Jacobi iteration so that should do Jacobi iteration for us is that fine j should be 101 that fine thank you okay any questions something they basically if you think about the x, y coordinate system at any given point I have an x, y right so that m grid is essentially generating for me the array of x's and y's that represent the mesh underlying grid does that make sense I can go to the board maybe it is fine we will see if it shows yeah so essentially what we are doing is no no it is fine it is okay essentially what essentially what we are saying essentially what we are saying is so this is 0, 0 this is what I am generating right so this is 0, up in this case 0.3313 0, 0.66 approximately right so the x matrix consists of all the x's and the y matrix consists of all the y's right so if you are somewhere in between so this would be 0.66, 0.33 and so on so the x matrix consists of all the x entries and the y matrix consists of all the y entries okay that is basically what it does is that fine yeah sorry about that anyway so let us get back to this I think this code is fine as it stands so if I run this now it should give me okay and we did this 100 times is there any chance anything would have happened 100 times do you want to run it longer okay why do not we repeat the process I will just make this I will just make this G repeat the process for Gauss-Seidel and Gauss-Seidel you remember I have to change quite a few things I have to make this E number or G 101 and I have to change the limits to 100 is that fine okay and that takes a little longer because now we are sort of looping through I am not using the inherent iteration capabilities of numpy or the matrix abilities of numpy but doing it myself manually in Python okay so G.plot errj101 it has non-positive values what happened let us try errg maybe I made a mistake somewhere yeah both have non-positive values what is x, x is a 101 by 101 matrix you just wanted to make sure that was right how is that possible how is it possible for square root to get a non-positive value it does not make sense let us try that one more time make sure I have not made a mistake anywhere I initialize 5 I set it equal to 0 we will do a 1000 iterations but of course that should make a difference just in case there is some issue there I will just do that as you can see even a 1000 iterations if you use numpy it is quite fast and let me plot that just to make sure that I am not going to I am going to get the same thing then we will just have to go on data has non-positive values still do not get it but anyway it is okay it does not matter the point that I want to the real really the critical point that I want to make out here right is first the blue line which is Gauss-Seidel does seem to converge faster still then Jacobi iteration one thing good but the bad thing is which we have which is what we expect that the rate of convergence is quite slow right look at what 5 by 5 is doing down here and look at what my Gauss-Seidel is doing up a go Jacobi is doing up there it is going to take forever for it to get down to a small enough value that is 10 power-6 or 10 power-16 is what we did in the other one am I making sense is that okay right so the convergence is actually extremely slow so just as a let me let me do then a I will just restart this just to be clear that there are no threads running and so on let us just now do SOR or something of that something little different okay right and we will try to see whether we can animate the code while it runs so from maybe I could have done it with the earlier thing but it does not matter okay just to be fine right okay we will do maybe we will do one other thing what we will do is we will see how this how these iterations go I just wanted to show you that in the last class I use Maya V to draw contour plots I just wanted to show you that you can actually do surface plots okay so why do not I why do not I just do that so I am going to do redo the 101 by 101 grid I created the 101 by 101 grid right P is x star x-y star y I do not know how many of you have tried using Maya V so far but it is relatively straight forward to use what I do is I say V equals m lab dot surface x,y,phi this is to show you what the initial solution looks like and that should generate something that looks like that like what is that one of the things that I like about one of the things that I like about Maya V is that you can interact with it okay so let me maximize it okay I am going to add I am going to add a legend to it so legend so that we can see what it is so it goes from minus 1 to plus 1 and the coordinate system is at the bottom left hand corner if you cannot see the coordinate system I will make that larger okay the coordinate system is at the bottom left hand corner right fine so this is the solution that we are expecting right so if you are saying wait a minute what is this what is this okay let me put an outline list like I did last time and as I said the thing that I like about it is that you can interact with it right so you can see from where the coordinate system or from in fact from the fact that the solution is 0 0 at this far corner tells you that is x equals 0 y equals 0 right so and along this diagonal along this diagonal x equals y so the solution is 0 so here it is a positive value right and here it is a negative value fine it goes from plus 1 to minus 1 so what happens when I set my solution when I initialize the solution that is the final solution that is what we are supposed to get okay if I just so that you see what is happening here most probably you already know what is happening I set it to 0.0 that is our initial condition okay our initial condition except at the boundary where it is x squared minus y squared everywhere else it is 0 that is our initial condition is that fine everyone okay so now that we know that Gauss-Seidel and Jacobi we can actually see the solution evolve so what we will do is I will say let us do Gauss-Seidel let us do Jacobi okay for I in Jacobi is faster right Gauss-Seidel is slower which for the sake of a demo is better but anyway how many iterations do you want me to do at the rate at which it was going most probably 5000 iterations or something is most probably what I mean 5000 iterations okay fine let us see if we can find the Jacobi control R PHI I will just search back okay that is Jacobi iteration is that fine everyone okay so and I am going to plot only the PHI so we have computed the PHI V dot M lab scalars at PHI let us run it and you can see the solution evolve and my computer is having problem is not going to rescale till I get to the last version okay that is fine so all you can make out so this 5000 is quite a lot all you can make out is that flat portion seems to be diminishing at a rather slow rate and it is really going to take for almost ever in fact on the video itself I do not know whether this will be visible this is going to take forever to converge so it looks like even 5000 may not quite make it right you can see that little gray spot in the middle you can see that little gray spot in the middle diminishing and going towards the solution remember this is Jacobi right go sidle will basically run twice as fast okay and that is basically it so essentially the thing that you can do is that you can visualize you can visualize as I said maybe I should not have minimized it but you can visualize this code running as the code evolves it is actually possible for you to look at the solution right and remember that just because on the screen it looks close to the solution does not mean that it is close to the solution screen resolution is very small is that fine okay so this is these are some skills that possibly I am not I am only showing you the tools being used these are some skills that I think you should try to acquire that you are able to visualize a that you are able to visualize your data be that you are able to do it while your code is actually running right there are sometimes especially if you are going to do it with a scripting language like python or something of that sort there is a certain advantage to being able to interact with the code that you allow it to finish executing certain number of time steps or certain number of iterations or whatever it is and then come back and interact with the code okay fine so that patch is almost gone so actually at 5000 it seems to be quite close maybe 5000 was on the large side let me let see fine okay meanwhile are there any questions let me see if I maximize this whether it will I give it time it will so even if I try to interact with it right now because my because it is completely tied you think about it 101 by 101 is a 10000 by it is 10000 unknowns it is solving a 10000 by 10000 system we are doing 5000 5000 iteration that is a lot of calculation so it is also a good idea to do the computation as to how many operations we are performing beforehand before getting into it like I just did right now maybe I should have just said 1000 instead of 5000 anyway it is okay it does not matter it should be close to completion so the other point of course is that if you try to do this 5000 iterations directly in C if you have already done this you know most probably seems a lot faster than this okay so that is the big advantage of doing it in C or C plus plus of OTRAN so you get you pay a price for the ability to interact with the code I am tempted to just kill it maybe I will just kill it okay fine let us get back we will start python again now we will do we have just enough time to do SOR so I have a Laplace equation solver that I have already written and SOR basically what it does is it is going to run SOR for me for that Laplace equation solver right now it is hardwired for a 41 by 41 okay hardwired for a 41 by 41 grid so I will create something that will run SOR for me okay and it is very it has very simple it has very simple way to do it so what does this do let me let me what does run SOR do what am I doing what run SOR does is it is going to run SOR for the values going from 0.0 to 2 point right and I am taking 20 values in between it is going to run for 20 values it is going to take 10 iterations is that fine I think I have really worn you guys out with this anyway okay let us see let me just run it and see what we get and this S has an error which is the error versus omega right I have S has an omega those are the omega values that we ran and it has a set of errors is that fine okay right so I will import just to plot it just quickly do this okay and plot there we go right and anytime we are plotting error we really need to this looks relatively flat okay this is look that encouraging but then we go to the power of the log scale anytime we plot error we always look at it on the log scale then look that flat now right so it looks like whatever we are looking for this is important by the way so you have to look at it in the log scale whatever we are looking for the optimal omega value that we are looking for is somewhere there is everybody with me the optimal omega value that we are looking for is somewhere there so let me hold that graph because I want to plot on top of it okay and I will create another SOR solver for me okay but okay let us go back to that graph now I want to hunt I want to try out another set okay so what value shall I take what value is omega where shall I hunt I look at this graph where shall I hunt 1.5 well 1.2 you know is not going to do it 1.5 to 1.9 may be right is conservative bear in mind when you are doing this that the convergence that you have is non-uniform convergence that is what we want okay so this limit is going to drift right does that make sense does it seem or it is just going to continue down 1.9 to so I say omega I will say W NP dot I will create a vector the syntax is very similar to M grid from 1.5 to 1.9 and maybe we will take 10 samples instead of 20 okay 10 samples okay and I will create another SOR solver and pass it this omega so S1 should run 10 iterations S1 should run 10 iterations for and S1 also has an error so I will plot as my plot S1 dot omega versus S1 dot error well 10 iterations what did you expect I mean it is going to be the same thing the only point is that the curve is smoother because we have more points that is all that is happened okay so I fortunately I can go more right so I have a we will go for a little more time I have a little go more which will go for a little more time for 10 more iterations so and if I do that okay then I can plot that again and what do I get aha so it is better and it is a little sharper so it seems to be somewhere there okay it seems to be somewhere there is that okay in fact what I will do I will get rid of some of this part or let it is okay I will get rid of some of this part it does not look as sharp because I have removed that part but right so it seems to be around what you say well yeah I will be conservative still 1.65 okay that is 1.7 1.65 to 1.85 why why why risk it okay so now I go back and I set this from 1.65 to 1.85 I am doing this because this is a demo I mean what you would do is you would hunt for it maybe use bisection or something you would hunt for it in a more systematic fashion the other ways by which we can do this so now I will create a S2 so that does only 10 that is no that is no good so we will do an S2 go more and we will plot S2 as you would expect it is the same you have done it twice it is the same it is a smoother curve we will do it some more okay so you can see it is getting tighter and it seems to be around between 1.75 and 1.8 so it is possible for you to actually hunt for an optimal omega right and trust me this curve here is not going to move that much okay if I were to plot if I were to do just for the fun of it if I were to run S remember our friend S there are 20 of them so it will take a little more time maybe I will do it twice because I did the other one twice okay if I were to plot this but you do not know what happened okay that is what you get so it is very clear that when you when I say optimal omega it is truly optimal omega and look at this look what is happening here you picked a value because you took 20 samples you picked a value there you picked a value here the optimal is actually somewhere in between optimal is actually somewhere in between okay so you can look at this and get be under the notion that it is at 1.6 actually it is not at 1.6 because of the way that you are taking samples okay there is an idea see you are looking at it and with your your eye you are actually doing some kind of finite difference right you are looking at it is increasing here decreasing there you are looking at slopes you are doing derivatives in your mind and you get the impression that the minima by setting the derivative to 0 is somewhere here but actually it is somewhere to the right and because of non-uniform convergence it looks as though it is drifting to that it is drifting to one side okay so the hunt for the optimal value of omega has to be done a little careful right but if you are going to make production runs I am going to make lots of runs I am going to solve this Laplace equation for different boundary conditions you know I am going to do 10,000 different boundary conditions then spending time finding that optimal omega makes a big difference because where here you are still at 1. something this is at 0.01 something right in fact it is 0.001 something am I making sense this is this minimum here just to give you an idea is at 0.003 at omega equals 1 it is at 0.3 please remember it is a log scale at omega equals 1 the residue is at 0.3 at omega equals 1.8 or 1.75 the residue is at 0.003 it is running 10 times faster it is converging right now it looks like 10 times faster is that fine okay and I casually made a suggestion that if you are changing boundary conditions that you can find the optimal omega and find out and solve for the various boundary conditions right the question that you should pop to me immediately is how do you know that the omega does not depend on the boundary condition do you think that the omega will depend on the right hand side optimal omega will depend on the right hand side or does it just depend on the matrix A solving A x equals B using SOR the question is does the optimal omega depend on the right hand side or does it just depend on the matrix A think back to the graph that I drew right and the steepest descent idea that we had which I mistakenly said Newton method at that steepest descent idea that okay fine okay. So I guess we will meet in tomorrow's class I would suggest that you try out some of these things for yourself okay we will meet in tomorrow's class thank you.