 Okay, good morning. So what we were doing last time is we were looking at the stability analysis for heat equation right, we were doing the linearized stability analysis for heat equation. To remind you as to why we were doing that, I had suggested that FTCS was unstable because the second derivative term that showed up in the modified equation for FTCS at the wrong sign was negative right, the coefficient for the second derivative was negative right. And I had suggested that we could add artificial dissipation of our own right, we could just add the second derivative term and try to knock out that negative term possibly make it even positive okay, that was the suggestion. And then I had sort of casually mentioned that you can add as much dissipation as you want, you can add as much artificial dissipation as you want okay, that is fine. So right, so what we are looking at now is what is the effect of this artificial dissipation. So if you add this, if you can add as much as you want right, then we had concluded that the equation with that much artificial dissipation was almost like heat equation, if you added an enormous amount of artificial dissipation that was like heat equation and was there a consequence to that. So we are just looking at the stability analysis of heat equation, we had already done it right, I will just quickly repeat it. So what did we have, the heat equation, 1D heat equation looks like that, if you use FTCS, if you apply FTCS to this, again the problem by itself is not complete right, I mean you have to apply the boundary, I have to tell you what are the boundary conditions, what is the actual physical problem that we are solving. But since this stability analysis, von Neumann stability analysis is only at a grid point right, we can just discretize this and ask the answer that question. So FTCS applied to this is UPQ plus 1 is UPQ plus, I will say lambda times UP plus 1Q minus 2UPQ plus UP minus 1Q, where lambda is kappa delta T by delta x square okay. And we went, when we went through the same process, when we went through the same process of writing UP plus 1Q as E power i n delta x UPQ, incidentally in some books you will see this called the shift operator because it shifts us by one grid point okay. So some books, you may see this called the shift operator okay. So doing this, what did we get? We got UP plus 1, UPQ plus 1 divided by UPQ which is G, the gain, what did this work out to be? 1 plus 2 lambda times cos theta minus 1. And we wanted the modulus of this to be less than 1 and therefore we got the condition that 0 less than lambda less than 1 half, is that fine? That is 0 less than kappa delta T by delta x square less than 1 half. So if kappa is very large, so the question is what dominates, what was the CFL condition, the regular condition that we called the CFL condition, that was 0 less than A delta T by delta x less than 1 okay. Depending on the relative magnitudes of A and kappa, you will end up applying the appropriate stability condition. So if you make your artificial dissipation, the heat equation it is kappa but if you add it as artificial dissipation, it would be mu 2. If you make the artificial dissipation too large, so large that mu 2 dominates, for all practical purposes it becomes like heat equation, then this is really the stability condition that you need to be worried about okay. So if A is of the order of 1 and kappa is of the order of 1 million, then this is really the stability condition that you have to be bothered about okay. Of course in general you could say that it seems that the general scheme, general explicit scheme right using three points seems to be of the form UPQ plus 1 is some A times UPQ UP plus 1Q plus B times UPQ plus C times UP minus 1Q. Am I making sense? In general, an explicit scheme involving the points P plus 1, P, P minus 1 would look something like this okay. Maybe what you can do is you can try out, just try out see how you do the stability analysis for this. Just substitute, go through the same process, just go through the same process and you will see that you can come up with the stability condition in terms of A, B and C okay. You can make a general statement, you can come up with the stability condition in terms of A, B and C. You can try it out, if you have any difficulties get back to me, is that right okay. So this is as far as the 1D heat equation is concerned, so you can add artificial dissipation, you have to be very careful how much artificial dissipation that you add right, you cannot add too much of it. The other thing that you want to remember is the artificial dissipation term that you are going to add, you are going to add to the equation that you are solving. That is if the equation way back when I told you that U is a perturbation but we will go back to situation where U is the variable that you are solving for. So the equation that you are solving A dou U dou t plus A dou U dou x and you want this to be equal to 0, you want this to be equal to 0 but instead of 0 you are going to solve this. You are clearly not solving the equation that you set out to solve okay. You are clearly not so it is not as though you are doing FTBS setting sigma equals 1 right, so the solution is contaminated. So from your point of view you want to keep mu 2 as small as possible, you want to keep mu 2 as small as possible is that fine okay. So and we rationalize this, we justify this saying anyway I am solving the modified equation right, so I am just fixing the modified equation to my satisfaction. So this term possibly can be picked exactly from the modified equations. When I do a demonstration I will actually show you how to pick this right, how we would pick this term. We will try out a few things and see what it does to the modified equation okay. So in a later class I will actually do a demo of deriving the modified equation but in a automated fashion and then we will see whether we can add different terms and see what that does okay right. In this conversation I am going to do a little aside now, we are going to take a step aside right because I am with heat equation, I am going to continue with heat equation just for a brief period right and then we will come back to this conversation right. I want to look at 2D heat equation simply because of the stability condition that I have got right simply because of the nature of the stability condition that I have got I am going to make a point, I am going to use it to make a point, I want to look at 2D heat equation. 2D heat equation is dou u dou t equals kappa dou squared u dou x squared plus dou squared u dou y squared, it is an isotropic material so kappa there is only one thermal conductivity right thermal conductivity is not changing based on orientation okay. So if I did FTCS for this, what would I get? So now I will use the superscript for time because otherwise I will have too many subscripts in the denominator, I will use p and r for x and y and that is q plus 1 okay q is in time, this is along x that is along y uprq fine plus lambda x up plus 1rq minus 2 uprq plus up minus 1rq plus lambda sub y upr plus 1q minus 2 uprq plus upr minus 1 where lambda x is kappa delta t by delta x squared lambda y is kappa delta t by delta y squared. So to keep life simple I am going to make delta x equals delta y equals h okay so I can then combine these I can then combine these. So what does this give me? uprq plus 1 equals all of this I am going to now do the same thing that I did for in the 2D case what I did for the 1 in the 1D case right. So I do not know whether I can if I just write it out whether it is going to be a hassle for you. So this is uprq maybe I do not skip a step plus lambda times what is this going to be remember now delta x and delta y are the same right. So i n delta x right you will have i n delta x and i m delta y okay i n h and i m y. So that is going to give you 2 sets of thetas okay maybe I will just write it out and then I will explain what I am doing. So this will give me a e power i theta minus 2 plus e power minus i theta times uprq and the other one will give me a lambda e power i alpha minus 2 plus e power minus i alpha uprq where theta is i theta n delta x and alpha is m delta y but we have decided delta x equals delta y but the wave numbers are still different okay. However much we simplify the x and y coordinates keep on falling apart it does not matter okay. So what does it give me therefore the gain when you take one time step what is the gain 1 plus is going to be very similar to what we had earlier 2 lambda cos theta minus 1 plus 2 lambda cos alpha minus 1 okay and when is mod g less than 1 when is mod g less than 1. So you can combine these two basically right when because you basically have to look at when these values are the largest. So that will correspond I will let you work that out that will correspond to 0 less than lambda less than 1 by 4 okay. Just like we did in last class you can just work out mod g less than 1 and you will see this is the condition that we get right and we can sort of guess that if it were 3D Laplace equation it will be 0 less than lambda less than 1 6th okay you will just get one extra term for the z coordinate okay right. But this is really what I was interested in this condition when I saw the one half I said okay let us just let us take an aside and look at. So lambda equals 1 4th what happens here when lambda is 1 4th what happens to this equation when lambda is 1 4th what happens to that equation when lambda is 1 4th when lambda is 1 4th you get a very familiar looking equation upr q plus 1 is 1 4th upr q plus oh sorry p plus 1r q up minus 1r q plus minus 4 z right that gets cancelled plus upr plus 1q plus upr minus 1q right which was our iterative solution to Laplace's equation fine this is a iterative solution to Laplace's equation we will revisit this I could not resist taking this aside we will revisit this later but what this basically says is marching heat equation in time is the same as sweeping Laplace's equation in space okay. So now there are two different ways by which we can get solutions to Laplace's equation either you take Laplace's equation add a time derivative and solve the unsteady equation and allow this solution to evolve in time or you do an iterative method Gauss Seidel or whatever it is or Gauss Jordan and if you are doing Jacobi iteration not Gauss Jordan if you are doing Jacobi iteration it looks like all you are doing is solving heat equation what if lambda instead of being 1 4th or 1 8th what would you get can you guess lambda instead of being 1 4th if I take lambda equals 1 8th what would you get that is like picking some kind of a relaxation parameter that is like picking some kind of relaxation parameter anyway if you want we can have this as I said if you want we can have this conversation later but that is like picking the relaxation parameter. So you see you have UPR Q plus 1 if you picked it as 1 8th right you would get a linear combination of UPR at Q right omega times that and 1 minus omega times that whichever way you want you would get a linear combination okay you would get a linear combination the only constraint that you have is the only constraint that you have is that you have a stability condition here okay you have a stability condition maybe I can work it out instead of since I have gone down this path maybe I can just work it out what happens when lambda is 1 8 what do we get from that equation can you just tell me UPR Q plus 1 equals UPR Q lambda is 1 8th maybe you can just read it out that is plus 1 by 8 UP plus 1 R Q minus 2 UPR Q plus UP minus 1 R Q plus 1 8 UP R plus 1 Q minus 2 UPR Q plus UP R minus 1 Q okay there are 4 of these therefore it becomes half so that becomes 1 half UPR Q plus 1 half UP plus 1 R Q plus UP minus 1 R Q plus UP R plus 1 Q plus UP R minus 1 Q divided by 4 so the 1 8th I have written it as an average plus 1 half this is like taking omega equals half this is like taking omega equals half okay I do not know how many of you tried when I said why do not you try SOR with Jacobi I do not know if you have tried it you try SOR with Jacobi anyone well if you have tried it if you have tried it you would have found that for omega greater than 1 it would not have worked right it would not work but because for Jacobi iteration you have a stability condition that says lambda has to be less than 1 fourth that is omega can only be at the most 1 fine Gauss-Seidel it can go up to 2 Jacobi it cannot Jacobi there is a stability condition that says that lambda is less than 1 fourth or which corresponds to omega equals 1 okay if lambda is greater than 1 fourth you would not get the average right and it would not work fine okay right so this sort of connects I wanted to connect what we were doing with Laplace's equation right with all this evolving in time I want you to understand that so marching in time to a steady state solution is the same as sweeping in space so there is no sense getting worked up saying that oh you are going evolving in time I am just doing sweeping in space I am not you have an extra coordinate I do not have that coordinate no they are both basically the same right but what it does is it gives you a different perspective the same algorithm it gives you a different way of looking at it so as long as you keep that in mind that whether I marching in time or sweeping in space okay that there is a there is an equivalence okay that we are fine is that okay right so that is the that is the end of the aside that we will get back to we will get back to where we were that is the end of that aside we will get back to where we were what we are talking about now is how much dissipation can we add we saw that if you add too much dissipation you add way too much dissipation the stability condition changes right and the stability condition gets actually worse the time steps that you have to take will get much smaller so you have to be very careful whether how much dissipation you may be tempted to add a lot of artificial dissipation that does not quite work the second thing is if you add a lot of if you add a lot of artificial the artificial dissipation that you add right now the way we are adding it you could say we are adding it explicitly right we are adding it at the current time level the way you are adding the artificial dissipation actually contaminates the solution right so if you say hey wait a minute I am not actually solving the equation that I set out to solve anyway I am solving the modified equation I am going to add artificial dissipation that is one argument the other thing is look you know you are solving the modified equation instead of improving it so that the modified equation gets closer to the actual equation why are you making it worse right that is a counter argument that you can get right so we have we have to be a bit careful how you handle this but you have to have an awareness that whether you like it or not when you solve the problem there is numerical dissipation that is showing up because the modified equation is not the same as the original equation you have to have an awareness as to what it is right the second thing is if you are going to add something to it you have to add it carefully right okay let me try to let me try to come to this problem in a different from a different direction let me try to come to this problem from a different direction okay what if you had dou u dou t plus a dou u dou x equals 0 and I did not add I have been writing this and I did not do that or what if you had instead of a dou u dou x you had either a as a function or you had you had something else you had some u right the kind of equation we are used to influence mechanics would be something like this where u is the solution and therefore you do not know what it is right you do not know what it is beforehand a priori you do not know what is the value of u okay so what if you had the situation how do you do how do you ensure that you are doing FTBS how do you ensure that you are doing FTBS right we will look at this equation in greater detail a little later but right now I am using it just for motivation how do you ensure that you are doing FTBS or I am sorry how do you make sure that you are using up winding right it is not so much FTBS but how do you make sure that you are doing up winding right if you do not know what the sign is you have a scheme you do not know what the sign is how would we make sure that we are doing FTBS so there are different ways by which you could do it of course you could have you could have a if then kind of a discretization right so the algorithm then becomes it is a true algorithm you turn around and say if a is greater than 0 use backward space if a is less than 0 right use forward space that is possibility a equals 0 I do not know right a equals 0 of course does something to this particular equation so life becomes easier a equals 0 can be a headache right again we will see what the what that headache is at a later time so the other possibility is that you can ask the question is there a way for us to create a switch an automatic mechanism so that I do not have this conditional state right possibility so the question is what is what is that quantity and what is this quantity so if I divide this by 2 right so if if a is positive if a is positive this would be 0 right and that will be that will be a if a is negative this is going to be 0 and this would be would be a right would be a right so now we have something that looks like a switch something that is 0 1 so you could then turn around and say that UPQ plus 1 the objective here is I want you to see where we started off with we started off Laplace equation central differences we tried it out we tried out various things some things work some things did not work how do you develop these algorithms what is the way by which we are we are grouping around but as we get along you get better at it right and there are lots of little little tools that you can use to construct algorithms you have to get an idea as to how these things happen so this is UPQ I seem to have automatically shifted to a superscript it does not matter right minus a plus mod a by 2 delta t by delta x into what do you want here UPQ minus UP minus 1Q thank you minus a minus mod a UP plus 1Q minus UPQ is that fine so this will do it this will automatically switch one way to do it so automatically switch you do not have the conditional statements but then you are doing a lot of work right you do not have the conditional statements you eliminated the conditional statements but you are doing a lot of work you are going to evaluate these terms independent of whether that is 0 or not right you are going to add them all up and throw them away and you are going to add them all up and it may turn out that what you have fortunately here you are not going to end up with round off error because they are identical right but these kinds of algorithms you have to be a bit careful so you are going to you are going to calculate all of this and just chuck it because that happens to be 0 what is the other possibility what is the other possibility in the last class we saw something what is the other possibility what was the difference between central difference and forward difference or backward difference you remember I want to add something to this what was the difference between central difference and backward space forward central space and backward space you remember second derivative but the coefficient is very important okay so it is something like I will add mod a because we do not know the sign away I will add mod a right delta x by 2 dou square that is I hope I do not want the continuous thing I want the up plus 1 q minus 2 up q you can go check this from your last class plus up minus 1 q divided by delta x is that okay that is fine why did I do mod a here I do not know the sign I want you to check that both f tcs and f tbs what is the difference between f tcs f tcs minus f tbs f tcs minus f tfs okay and see what you get so if you add this quantity what is it going to do it is going to convert the central difference to backward space you understand and by taking this sign out of the game I have essentially made it made sure that whether it is forward difference or whether this is positive or negative that is going to cancel out you can just try it out and see okay and it will always be upwind that is another way to do it you add the right kind of artificial dissipation right you add the right kind of artificial dissipation of course from a stability point of view what we are talking about earlier I cannot afford to have this a to be negative that is pretty obvious that is mu 2 negative right that is mu 2 negative does that make sense that is mu 2 b if a is negative that means mu 2 is negative it is going to diverge right so it is pretty clear that it has to be modulus of a let us see there are different clues there are different clues that we have as to why we are doing what we do either you you can do it from here to see what is the correction term that I have to do in order to change the central difference scheme to a upwinding scheme or you can look at it from the modified equation that we have got and say over the coefficient has to be positive is that fine everyone so there is a way there are clear cut ways by which you can determine what happens when you add a specific term but the addition of this term does not eliminate the second derivative term the addition of this term only converts a central difference scheme to a one sided difference one sided first order scheme so we lost the order you lost the order of the scheme okay are there any questions so while we are at it while we are at this so this is so we have seen that you can do FT so we have FTCS FTBS FTFS right we have to use either FTBS or FTFS depending on which way the stream is flowing so it is better to do FTCS a centered scheme see this is one way to look at it one argument and just add something to it so that it becomes upwind biased or the other thing is to say that you do upwinding directly right you can get so whether you are doing upwind I would say if you are doing this you are also doing upwinding it is very clear if you are doing this you are also doing upwinding right so there is no sense getting into an argument as to whether you are doing central differences is just a matter of in order for the error to decay I am going to have dissipation and in order to get the dissipation I need to do something okay so there is no sense getting into an argument as to whether you are doing central differences or but the minute you say it is FTCS plus a correction term then we can ask the question can I eliminate the higher order can I eliminate the higher order terms and get a more accurate scheme that is a different story right we will look at that as I said when I do the demo okay now that we have done this I know in the beginning of the class I said I am not going to do a survey of a lot of schemes and so on but I am going to do a few schemes now just to show you just to go on with this philosophy of how do these schemes evolve how do we develop these schemes right and I will get you to a point where you should be if you really want right if you go out look at all the schemes that are out there and you say no I do not like this I have a better idea you should be able to come out with something on your own is that fine okay so we come here what we will do is and all of these as I said we have clues for these things from what we have done so far we have clues for these things from what we have done so far so earlier when we derived modified equation what had I suggested when we did modified equation what did we do we expanded we expanded using Taylor series and we substituted for individual terms in Taylor series okay then you can ask the question why you do not develop a scheme using Taylor series this is a classical technique using Taylor series to solve differential equation that is if you have UPQ plus 1 I guess maybe I will skip to the subscript is UPQ plus let us do Taylor series delta t dou u dou t at the point pq plus delta t squared by 2 factorial dou squared u dou t squared at the point pq and so on fine and we did the modified equation of course we wrote the right hand side but here I am going to stop at this point I say wait a minute why go there at all you already have a trick what is dou u dou t dou u dou t is minus a dou u dou x again I am back to a greater than 0 that situation where a is greater than 0 and what is dou squared u dou t squared a squared dou squared u dou x squared substitute back UPQ plus 1 is UPQ plus delta t a dou u dou x is UP plus 1q minus UP minus 1q by 2 delta x minus sign I do not know I always forget that minus sign okay fine plus a squared delta t squared by 2 delta x squared I do the open bracket here but I am going to write it here what is it or maybe I will write it in the bottom a squared delta t squared by 2 delta x squared up plus 1q minus 2upq plus up minus 1q how is that we just hooked up a scheme this is called the lax when drop scheme lax when drop scheme right all you just basically do is you just go through do a Taylor series expansion classically it is solution to differential equations using Taylor series it was normally done with ODEs is this application to PDEs we get the credit for it okay so here you have it FTCS there FTCS there so this would be sigma and that would be sigma squared by 2 and the scheme comes with its dissipation add right you can try out you can go through do the stability analysis for this I am not going to do this anymore right these I leave as exercise you can try out the stability analysis for it you can find out the modified equation and see what is the nature of the modified equation right you can try to find out what is the nature of the modified equation is that fine everyone okay what else I am just going to I am just going to know in a freewheeling fashion connect all the bits and pieces that we have done to see whether we can come up with other schemes that is all right I am just going to do a few of these before we sort of end this whole thing of wave linear wave equation right okay so I am going to just sort of try around something called a two step Laxman draw method they do the following I will draw the stencil in this case simply because I am actually drawing grid lines simply because it is easier to understand with the grid lines so this is p plus 1 p p minus 1 p plus 1 at time level q that is q plus 1 I seem to have drawn a line in between I will also draw another line in between here lots of extra lines right so that is this is a point pq plus 1 basically that is what we want I am going to do this in two steps I am going to do this two different ways right I will tell you what is the two step scheme final two step scheme so what we can do is you have the value here you have the value here you had this in between value you could do FTCS and find that so humor me I will call that p plus half right I know I mean it is a we were counting it is integers but we will call it p plus half p minus half good even with counting I have to make sure I get the sign right today is my day for negative sign anyway p minus half so then this would be p minus half value I would expect that would be q plus half right so you can see that if there was a way by which I could use FTCS to find that and a way by which I could use FTCS to find that then I could repeat the process and use FTCS between these two to find that is that fine okay two step process the process what I just described now is not the two step lax vendor of process it is a two step process so let me write that over first and then I will tell you what is the improvement so you can say up minus half maybe I will go back to super scripts here q plus half equals up minus half q minus a delta t by delta x up minus up minus 1 these are all at q is that fine everyone how do I find that up minus half q give me a suggestion so I can take the average up minus half q is up plus up plus up minus 1 divided by 2 so in a similar fashion up plus half q well it would be the same I mean so obviously you are going to find this value using these two right for each of these intervals you do that so I will go directly to up q plus 1 so here is the first first suggestion first suggestion up q plus half minus a delta t by delta x we are assuming delta x is equal everywhere up plus half q minus up minus half q does that make sense and any time to find this intermediate value if you do not know it take the average of the neighbors okay everyone it is fine so do you expect this to be stable using FTCS you expect it to be stable it should be unstable before you start the stability analysis you first predict what you expected to be and see whether you get right and in order to do the stability analysis you have up q plus 1 remember you have to get this in terms of right hand side has to be all upqs which basically means that you substitute for the p minus half in terms of p is and p plus 1 you understand get eliminate all the half that is the easiest way to do it just make the substitutions eliminate all the halves so that you finally get upq upq q plus 1 in terms of upqs okay take the ratio and you can do the stability analysis fine this is a two step method the two step laxman draw method you stop at this point and you say wait a minute there is something here why bother with taking FTCS here again I have the value here if I take a time derivative across this I am taking a central difference at this point I do FTCS to get that FTCS to get that I am sorry FTCS to get that I do FTCS to get this point FTCS to get that point here I can do central time I can actually do central time I can get a higher order accuracy in time I can do a central time is that fine I can do CTCS that means the second step instead of using this instead of using this in the two step laxman draw method you would write upq plus 1 is upq minus the rest of the stuff central difference in time central difference in space okay one of the reasons why I do not do this is it can get dreary right it can get tiresome I am just putting up indices out there right I can already see it on your faces it can get it can get it can be it can vary out okay this is the two step laxman draw method it is done in two steps first you get to q plus half and then you do the full thing full delta t fine okay I did this so that you are aware that you can you do not have to do it it does not have to be done in one step you can do it in multiple steps the time integration can be done in multiple steps in the whole class of rangekuta schemes and so on where the time integration is actually done in multiple steps so there are multi-step methods right just so that you are aware of it a third thing that we did right back at the beginning when we were doing finite differences what have we done so far amongst the various things we have done FTCS we have done BTCS and of course we have done lot of other stuff but I am more interested in FTCS and BTCS before we derived finite differences using Taylor series we did it for forward difference and backward difference do you remember the first time I introduced central difference what I did you remember what I did I took the average of the forward difference and backward difference right I looked at the truncation error term for forward difference and backward difference and said wait a minute these are the same magnitude but opposite sign why do not I take the average I have forward time central space backward time central space why do not I take the average why do not I take a combination of these two you understand right or better still I can take a weighted average not why just average now we are into things like SOR and you know things so why not just why not a weighted average why not a theta times FTCS plus a 1-theta times BTCS that sounds vague right theta equals half would be an average theta equals half would be an average how do we do this theta times UPQ plus 1 equals theta times UPQ minus theta sigma by 2 UP plus 1 minus UP minus 1 what is BTCS 1-theta times what was BTCS sigma by 2 UP plus 1 Q plus 1 minus UP minus 1 Q plus 1 plus UPQ plus 1 equals UPQ is everybody with me is that fine okay so add them up we will see what we get you add them up what do you get there is a 1-theta UPQ plus 1 just to point out this always happens a theta UPQ plus 1 so the theta will go away right so you get a from here a minus sigma by 2 UP minus 1 Q plus 1 a plus UPQ plus 1 plus sigma by 2 UP plus 1 Q plus 1 equals again as I pointed out the theta and minus theta will cancel UPQ minus theta sigma by 2 UP plus 1 Q minus UP minus 1 Q I lost some 1-theta somewhere it should be a hope that is not too messy fine fine and theta equals half theta equals half you get a very famous technique called the Crank Nicholson technique but here you would have to solve a system of equations but clearly you can take various values of theta right clearly you can take various values of theta I would suggest that you try to do the stability analysis for this you try to do the stability analysis for the two step Lax-Vendorff method all of these all of these schemes P plus 1 P P plus 1 P minus 1 Q Q plus 1 both of these schemes represent the differential equation approximate the differential equation at the midpoint both of these schemes represent or approximate the differential equation at the midpoint fine okay okay so I will see you in the next class thank you.