 Last class we looked at linear first order one dimensional wave equation right and how to solve it we thought we had a way to solve it, we tried forward time central space okay and forward time central space turned out on analysis on doing stability analysis, forward time central space turned out to be unconditionally unstable okay. So we are now going to try something else. So the equation that we are solving is this with the appropriate initial and boundary condition we have had a discussion earlier about initial and boundary condition right with the appropriate initial and boundary condition we decided that this is just a recap that we would use FTCS which means that the differential equation is represented at the point pq, q is an index in time p is an index in space right p-1q p-1q and pq-1 and we ended up with a scheme which was unconditionally unstable okay. So we decided the argument remain now we are we tried right we just tried central differences it worked in Laplace's equation with the same idea we just tried doing the discretization here it does it looks like when we did an analysis that it is not going to work okay we are doing forward difference in time the proposal last class was why not do forward difference in space and see whether that helps or not. So this is FTCS forward time centered space now we will try FTFS which is forward time forward space so we will use the point again we will represent the differential equation at the point pq right do not forget that we are representing the differential equation at the point pq right p-1 we will involve the points p-1q and pq-1 and this would be the scheme FTFS forward time forward space forward difference in both. So what does that give me so the dou u dou t term the dou u dou t term would give me a upq-1-upq divided by delta t the dou u dou x term would give me up plus 1q-up-upq divided by delta x right that is a forward difference in space that is a forward difference in space so if we solve again you will notice that except for one term which is in terms of at the time level q plus 1 all the others are at time level q. So we will retain the q plus 1 term on the left hand side and take everything else to the right hand side that is we can write the q plus 1 term explicitly in terms of the remaining terms okay so this is also an explicit scheme in that sense upq plus 1 is upq when I take this over to the right hand side this of course equals 0 I am sorry about that equals 0 when I take this over to the right hand side I get a negative sign so minus a delta t by delta x times up plus 1q-upq okay again we have an automaton again we have something you can just given the time given u at q you can go on to q plus 1 given the value of u at q you can go on to q plus 1 as I indicated in the last class we are going to see the so often that we give it a symbol that is called sigma okay it is called sigma and it looks like we have a similar simple scheme just like we had for FTCS right so but again there is a concern is it going to work or not okay so let us try to do the stability analysis for this now remember when we are doing the stability analysis we are involving only the points that are that are involved in this computation we are not talking about the whole computational domain we are not talking about the whole region we are not talking about the whole problem right we are not talking about the stability as in a global sense right so this is locally at this point if I go from point p point pq to pq plus 1 what is the gain if I go from pq to pq plus 1 what is the gain that is the question that we are asking is that fine everyone okay so what do I get just to remind you up plus 1 q is e power i theta upq I do not need the brackets where e power i theta was actually exponential i 2 pi n by l i 2 pi n by l i is square root of minus 1 fine that is only one that we need so upq plus 1 is upq minus sigma times I already know I am going to factor out the upq into e power i theta minus 1 is that okay just like we did last time we will compute the gain going from q to q plus 1 I can divide through by upq so the gain g is upq plus 1 divided by upq which equals 1 minus sigma times using there is a delta x and oh I am sorry thank you yeah yeah yeah yeah yeah n delta x that is very important delta x that n delta x okay theta has an n delta x defined in its definition because you were shifting by a delta x okay thank you right so what do we have here this is cos i theta cos theta plus i sin theta that is Euler's formula for e power i theta minus 1 we want the modulus of this to be less than 1 we are saying less than 1 the square is less than 1 gives us the same thing so mod g square is g times g bar g times its conjugate which equals we could have left it that way but I am going to expand it out expand this out so it will get a little messy we keep track of what is happening so that I do not make mistakes okay so what do I have maybe before I do this let me let me let me before I do this let us just expand this out so that I really do not make a mistake so I have a sigma here that is a 1 plus 1 plus sigma minus a sigma cos theta minus i sigma sin theta is that fine okay and its conjugate is the conjugate of course is 1 plus sigma plus sigma cos theta plus i sigma sin theta I am making the letters small but it is basically the sign flips on that there is a minus and minus that is a plus and the minus sigma cos theta and minus and minus that is a plus and minus sigma that is cos fine so this is 1 plus sigma plus sigma cos theta there I am sure each of you has a way by which you can do it I am just going to do it by brute force plus sigma squared sin squared theta standard expansion this gives me 1 plus sigma squared cos squared theta plus twice sigma plus minus it looks like that sign I am going to keep missing okay as long as you keep me honest it is fine okay 2 sigma minus 2 sigma cos theta minus 2 sigma squared cos theta plus sigma squared sin theta sin squared theta that is that last term so I have a sin squared theta plus cos squared theta that is going to give me a 1 is that fine everyone okay so what does this give me this gives me 1 plus 2 sigma squared plus 2 sigma minus 2 sigma cos theta minus 2 sigma squared cos theta and this has to be less than 1 a stability condition requires that this should be less than 1 right so far I have not I have just been evaluating the left hand side the stability condition requires that mod g squared less than 1 I am saying mod g squared less than 1 because we want mod g less than 1 okay fine so let me knock out that one this makes life a lot easier because now this quantity has to be less than 0 so in fact I can divide by 2 okay so I get sigma squared plus 2 sigma minus 2 sigma cos theta minus sigma squared plus sigma sigma squared cos theta less than 0 okay fine I can add I can subtract 1 from both sides fine okay okay where are we so now we can you can obviously see that sigma squared can be factored out sigma can be factored out so this looks like sigma squared plus sigma into 1 minus cos theta less than 0 theta equals 0 is the worst case scenario where it is where it is 0 right I want strictly less than anyway theta equals 0 again corresponds to the DC component right so except for that so if theta is not 0 I can just divide through by 1 minus cos theta right which says that for any other value of theta this quantity is going to be positive for any other value of theta theta greater than theta greater than 0 1 minus cos theta is greater than 0 so I can divide through by 1 minus cos theta fine so this gives me sigma into 1 plus sigma less than 0 so we actually have a stability condition okay when does this work no no it works sigma is let us just figure out what is the quantity let us see what it gets if sigma is less than 0 this will become negative right if sigma is less than 0 this will become negative whether you agree with sigma less than 0 or not sigma less than 0 there is a condition under which it will work supposedly right sigma is less than 0 but then we want to make sure that this stays positive right so sigma has to be greater than minus 1 so and sigma has to be greater than minus 1 what it basically says is minus 1 less than sigma less than 0 okay we have the condition now we can look at it and complain about right so what is the problem what is the problem maybe we go back there what is the problem a is greater than 0 sigma is a delta t by delta x and I want that to be negative right so I do not know I look at it and say we can infer from that that the condition will be stable right that the scheme will be stable if delta x is negative we do not want to go back in time we do not want to go backward in time that is not our intent right clearly if we go backward in time the scheme will work we do not want to do that the other possibility is it seems to tell us delta x is negative in my mind I would say well we are doing a forward difference is that telling me to do a backward difference right or a could be negative we have chosen we have chosen the problem a greater than 0 which means that if a is negative it would work right a is negative it would work so what I will say is well you know I get an inkling that may be instead of forward space I should use backward space the condition seems to tell me to go in that direction okay this condition seems to tell me that instead of using ftfs that I should use ftbs so instead of using forward time forward space I will use a stencil that looks like this so this is pq p-1q pq-1 and we suspect from the analysis that we have done right now we suspect that this should work okay we are now more hopeful that this will work so what does that give me by now you should be able to write it out okay right so upq plus 1 equals you all have the same form upq-sigma times upq-up-1q okay just look at the ftfs and you will see that this is what you get the sigma is got the delta t and delta x in it from the derivative is that fine everyone okay so as a consequence we will just quickly redo go through this process I do not want to get into trouble with signs I will get a little lazy here and make that p-1 that is e power-i theta okay so e power i theta of course is still defined that way that is fine so what do we have upq-1 equals upq-sigma times upq into 1-e power-i theta is that fine okay so all 3 of these schemes that we have got are explicit and here we have again upq-1 in this fashion okay divide through the gain g is upq-1 by upq which equals 1-sigma e power-i theta because this is I had written it last time as I just expanded it out-cos theta plus i sin theta is that right that we careful a lot with all the negative signs and therefore mod g squared equals is going to be a 1-sigma plus sigma cos theta okay 1-sigma plus sigma cos theta I am sort of skipping a few steps here I think you can keep track to make sure that I am not making any mistakes minus sigma that gives me a sigma squared sin squared theta is that fine everyone so that equals 1 plus sigma squared plus sigma squared cos squared theta looks very similar here is where it changes minus 2 sigma at flip sign plus 2 sigma cos theta what else plus 2 minus 2 sigma squared cos theta plus a sigma squared sin squared theta the sin squared and that cosine squared will again combine so this will give me something very similar 1 plus 2 sigma squared minus 2 sigma plus 2 sigma cos theta minus 2 sigma squared cos theta less than 1 okay subtracting 1 from both sides that gets knocked out right so right hand side becomes 0 what do we get here we can of course divide through by the 2 also okay so what do we get now so you have the sigma squared sigma squared which is a sigma squared into 1-cos theta and the minus sigma into 1-cos theta it will be sigma squared minus sigma into 1-cos sin theta less than 0 is that fine and this should do it for us so using the same argument divide through by 1-cos sin theta so you get sigma into 1-sigma less than 0 have I got that right and here because it is 1-sigma you saved so this tells us that sigma is greater than 0 and sigma is less than 1 okay so 0 less than sigma less than is that fine this is called the I make a mistake yeah sigma is greater than 0 and less than 1 did I make a mistake sigma into and I do not see it oh I am sorry okay this is greater than 0 I flipped it around okay sorry it should either be sigma into sigma minus 1 or I flipped I flipped the I made it 1-sigma actually that is true I do not know how I made that jump okay that is fine yeah it should be greater than 0 yeah it should be greater than 0 fine otherwise it would not have worked is that fine sigma is positive 1-sigma is positive so the product is positive okay sigma-1 of course then it had to be less than 0 that is fine thank you okay so what does this tell us the sigma value is positive so it works okay so this is called the remember I remind you again this is either called the Courant number Courant number or depending on who you talk to whom you talk Courant Levy-Fedrich number okay very often just refer to as the CFL named after the 3 very often just called the CFL okay so in computational fluid dynamics this number figures so often right so in the usual jargon will be so if you say oh my code I am having some difficulty people may ask you what Courant number or Uroni or what CFL or Uroni right what then they are not talking about this condition which itself is called the CFL condition Courant-Fedrich Levy condition so this condition itself is called the CFL condition the condition itself is called the CFL condition right condition itself is called the CFL condition but very often we do not know see we have just done this we have got this for what we have got this for forward time backward space for linear first order one dimensional wave equation with a greater than 0 propagation speed greater than 0 right this is very very restricted just it is not like FTBS this is the condition right for that problem so we do not know what conditions you would get otherwise right we do not know what are the other possible stability condition so this condition is called the CFL condition and this number is called the Courant number is the CFL number okay so it is normal for people to ask the question what CFL are you running or what Courant number are you running okay right so at a later date if some assignment you have some problem or whatever it is that may be a question I ask what Courant number are you running is that fine so since this number seems to be so important let us just sort of take a look at this number see what it is all about and see why FTBS works right there are two things that we need to do first that we have this number I get back here okay maybe I leave this I leave this figure here just so that you can see what is happening we will just look at what we have got here so A is greater than 0 sigma is A delta t by delta x right so this is the Courant number or the CFL what are its dimensions physical units it is non-dimensional okay so first of all we see that it is non-dimension it has no dimensions so A is a propagation speed which has units meters per second so it seems that delta x by delta t is also some kind of a propagation speed right I mean I am dividing A divided by some propagation speed this is how we would normally think about it in fluid mechanics so delta x by delta t comes from our grid so every time we take a time step it will we are every time I take a time step I am sort of propagating something that is at a distance delta x across a distance delta x when I go through a time step delta t I am propagating something that is at p-1 whatever property I have at p-1 I am causing its influence to propagate to p over a distance delta x so this is the speed delta x by delta t is the speed at which my grid is propagating the properties that I am right that my equation is trying to propagate equation is trying to propagate you this is what we found out this differential equation is propagating the property u at the speed A that is what this differential equation that is what we saw when we looked at characteristics what the scheme does is the scheme propagates the same property u right at a speed basically delta x by delta t effectively delta x by delta t that is what it looks like okay so the ratio A divided by delta x by delta t is the CFL number and this is called the grid speed I give this to you only as interest because very often this is not going to crop up that often in conversation right the CFL will occur more often in conversation this is by way of just information so delta x by delta t we will call the grid speed okay so it is non-dimensional that is fine so that looks like right now we will settle with that right we will come back to this we will come back to this propagation speed and what this means what it means a little later there is only one other thing that I want to tie with this so central differences did not work forward differences did not work a backward difference in space worked okay so if you look at again the fact that I am talking about the ratio of propagation speed to grid speed what is the speed with which it is propagating it is not just a speed right that there is also a direction involved in what direction is it being propagated so and that seems to be the problem here so in this case in this case we were propagating with this grid in both directions right so whereas our characteristics and I will just draw I will simply draw the characteristics just for just for our reference a characteristics are for this problem with a positive slope basically that means that the characteristics are I will use a different color chalk a characteristics are oriented in this fashion a characteristics are oriented in that fashion okay so the physics wants to propagate along this direction the physics wants to propagate along that direction in this case the numerics was also propagating in the opposite direction which is not what we want to do the numerics again wants to propagate in the opposite direction that is what we are trying to do here when we do forward space whereas here the numerics and the physics are propagating the same direction this idea that the numerics captures right the idea that the numerics captures what the underlying physics is doing right is basically called referred to as up winding so this is given a name up winding or if you are more interested in that water stream flowing that I talked about earlier or up streaming okay wind as in air flowing up winding okay up winding or up streaming fine so there are we will get back I have something to say about this at a later date but right now I just want you to be aware that the fact that we went through this process I want to know what is the meaning of the sigma and where it comes from this is the relationship right the fact that the sign change the sign change a change of sign the slope changes if a change of sign if a were negative then you would get a characteristic that is like this okay in which case forward space would work that is what this tells us that is what that was the clue that we got and we use that clue to get to the get to this point is that okay right so before I go on there are some remarks I want to make so the first thing that you have to see is that just because you know Laplace equation essentially it means that Laplace equation we sort of got a little lucky right we took Laplace's equation we did central differences I mean it was not a coincidence I obviously picked Laplace's equation as a problem because you remember this came under the category of simple problems and I obviously very deliberately picked Laplace's equation because it would work right but if you were trying it out you would you may you may try Laplace's equation and say oh it works this is very easy I know how to do finite difference methods right I am set so it does not work that way so you could take an equation you could use you could go through right you could go through a legal process saying I will do central differences I will do forward difference we went through the discretization and we came up with a scheme that would work if you have not tried to implement FTCS yet please implement FTCS right that is what I am trying to tell you is just because we know that FTCS does not work we know because we have done some analysis where first question is do you trust that analysis so check it out to see whether it works or not right if it is going to fail you have to have an idea as to how it fails I want you to see your program not functioning that is some because you will encounter this error see you may go you may discretize an equation at a later date right and I want you to be familiar with how these programs fail when they fail it is not enough do not just implement FTBS saying that oh I know this works the other 2 do not work I ignore the other 2 implement all 3 of them try all 3 of them and see how it works right and there is a little issue about boundary conditions as to how to apply boundary conditions so I want you unless you implement it you do not face that question okay so you have to try it you have to try it I want you to try it right I want you to try it so that you will see that you can try this this and this all 3 stencils okay you can try all 3 stencils. Now as I said normally when you go through this is what is going to happen the conclusion is just because we have a discretization that does not mean it works okay so now we have a scheme now we have a scheme that actually works maybe I will say something about the boundary condition and you know boundary conditions and so we are doing up winding we are doing up streaming I have suggested that you try FTCS but is there an issue with FTCS as I said normally what I would do is I would wait for you to try it out and come back and complain to me that there is a problem right but because we are in different setting let me go ahead and tell you what is going to happen what is going to happen let us try FTCS let us say let us say let us say you do FTCS okay so this is your x axis remember this is a pipe this pipe was of length L we can take without loss of generality we can take L is 1 right so the length of the pipe is 1 you go back to the problem that we pick and we can pick grid points along that length this of course is a time level Q and that Q plus 1 you would have similar set of grid points that is Q plus 1 and what does FTCS do FTCS uses these 3 the 4 but the middle one does not actually play a role in right because of the central space first derivative the middle one does not show but it does not matter it uses these 3 so this is the stencil it uses these 3 here and 1 there in order to march forward to step forward in time fine where have we prescribed boundary conditions remember for a greater than 0 what was the problem where have we prescribed boundary conditions we have prescribed boundary conditions at t equals 0 rather the initial conditions at t equals 0 right initial conditions at t equals 0 and at x equals 0 we have prescribed boundary conditions I will just say boundary condition is that fine everyone it is okay and we gave a reason we said look the characteristics are from left to right we gave a reason the characteristics are like this and therefore the boundary conditions are prescribed on that side right we gave a reason for it does it work with our notion of what differential equations are all about it does because you have a first derivative in time that tells me that I need an initial condition and you have first derivative in space right so in my mind when I say I am going to integrate with respect to time I will get one constant of integration that needs to be determined and if I am going to integrate with space I am going to get one constant in space that is going to be determined right so I need a condition along this and a condition along that we are satisfied mathematics is satisfied the physics is satisfied right what about the numerics look at this if you are doing forward time central space what is going to happen if you are try to implement FTCS what is going to happen so when you are calculating this grid point yes that works this is known this is known this is known what happens at this end when you want to calculate this for this grid point you need this you need the value on the right hand grid point am I making sense so the mathematics we have different we have boundary conditions that came from the differential equations from the mathematics we have boundary conditions which we have rationalized we agree that this is what the physics is right we gave this reason saying over the stream is flowing down I put chalk dust here once I put the chalk dust you are sort of committed you cannot insist assert at a later date right that the chalk dust flowing down should be something other than what you put right once you put the ink dye in or chalk dust or whatever it is tracer element into the stream you cannot then say no no no I do not want any any chalk dust at the other end it is there it is in you have no choice you cannot once it is in it is in right we are assuming that you are not able to grab it from in between there is no there is no source of sink in between okay so therefore the physical condition the physics of it are also satisfied but here the numerical scheme the numerical scheme requires a boundary condition the numerics requires a boundary condition in order to implement FTCS I need a boundary condition I will admit the fault is mine if FTCS is not you know it is not stable but if FTCS was stable I would use it I mean it is central difference second order accurate in space why would I not use it right and if I want to the fault is mine that I am going to use FTCS and once I decided to use it I need a boundary condition right I want to emphasize that neither the mathematics nor the physics insist on so we have to generate this boundary condition okay so how can we do it what do we do what are possible ways by which we can do it so one possibility is that you extrapolate right my project you mean extrapolate so you could take you can take whatever value here when it is the initial condition it does not matter but that is why I strategically called it Q I wanted to make sure it is not necessarily Q equals 0 right okay so you can you can so we do not know that there is no condition prescribed here so what you can do is you can copy you can copy maybe I will use you can copy the value at if this is the grid point n the last grid point from n-1 you can copy it to n then you will have one there right why am I copying it from n-1 to n because that is the direction in which my equation is propagating okay right there is another possibility this problem of using requiring a condition here will clearly be there even an FTFS but it will not be there for FTBS right it will not be there for FTBS in fact in FTBS we can solve the differential equation right here at the last point okay so in FTBS I am going to draw I am going to draw only this portion in FTBS these are the points that would be these are the points that are that part take that take part in this right so this is n this is n-1 this is Q this is Q plus 1 and in the case of FTBS you can actually solve the differential equation at this point okay at this point right and you can actually get the value here from the values at time level Q right so one possibility is that you just copy the other possibility is just for this one point you use FTBS is that fine just for this one point you use FTBS is that clear so there are different ways see now you can see why do we have these options we have these options because this is a problem that we have created by picking the particular scheme however it is our need that we need the value of U at that point okay is that fine I cannot emphasize that enough so I want you to be aware of this there are boundary conditions that your physics that are available from your physics right so when you solve the problem when you are I would not say solve the problem when you pose the problem you study the physical problem there are certain boundary conditions that are available okay then you go through a process of modeling and come up with a mathematical model the mathematical model will require boundary conditions right there are problems in which the mathematical model may require more boundary conditions than the physics provides am I clear the mathematical model may require more boundary conditions in the physics provides okay I can give you instances of such problems. Third possibility is given that you have the mathematical model and the physical model when you get to the computational model the computational model may require even more boundary conditions than these two require okay so in both cases the physics is what is available to you right because our physical intuition is what tells us what we need to do so when you get to the mathematical problem requiring more then you have to figure out a way by which you come up with those boundary conditions okay that is a little more difficult because they have to be compatible with the differential equations and so on in a similar fashion when you have the numerics right you need these boundary conditions you have to apply those boundary conditions in some fashion you have to generate those boundary conditions in some fashion fine okay are there any questions so now we have one possibility what we will do is we have done we have done forward time central space forward time backward space and forward time forward space and backward space okay we have done all three since forward time and backward space work you know see this is you have to look at the process that I am going through right remember yesterday's class I said well we are sort of grouping around trying to make sure that we come up with a scheme that works we tried FTCS it did not work that is somewhat of a shock then we are left with let us try let us figure out what to do if you are doing this for an equation for the first time that is how it would typically work I mean you what you try first off may not work so I say well backward space seem to do it why not backward time right I mean after all I tried forward space because I was doing forward time right and it got us somewhere so I asked the question why not backward time when I take a backward difference in time and if I am going to do a backward difference in time 3D fellow why not go back to central differences in space right unless I want to do second I want to do that right higher order scheme and unless unless the unless it fails unless it fails see this is I am just talking about pure right just from I want try to do the best that I can whether it is indeed the best that we can or not we have to see whether it is the best thing to do we have not judged we have not even computed any of this stuff we do not know how good these schemes are right now we are just as I said exploring or grouping around to find out what we are able to do so if I do backward time central space if central space does not work then we will try forward space backward space most likely we will try backward space right it is obvious right but as I said we will get a little greedy and see what happens we may pay the price for it so what does that mean that means the differential equation now is that P well we have a choice we can either represent the differential equation at PQ or represent the differential equation at Q plus 1 maybe we will do it at Q plus 1 for a change PQ plus 1 in the sense this index is our choice right so we will call the new grid point always Q plus 1 P minus 1 Q plus 1 P plus 1 Q plus 1 PQ and we are representing or approximating or representing the differential equation at that point at that point at that future point which we do not know the answer okay. So what do we get this is BTCS backward time central space what is do you do T approximately I am writing it out explicitly because I want to tell you where I am doing UPQ plus 1 minus UPQ divided by delta T right I mean it is obviously the same thing this is what we are saying earlier when we are talking about finite differential okay and do you do X at PQ plus 1 which is where we are approximating it UP plus 1 Q plus 1 and UP minus 1 Q plus 1 divided by 2 delta X substituting it to the differential equation going through the same process keeping all the terms that are in terms of Q plus 1 on the left hand side we will get UP plus 1 Q plus 1 minus UP minus 1 Q plus 1 by 2 delta X plus UPQ plus 1 times A delta T equals UPQ is that fine and again UP plus 1 Q plus 1 is UPQ plus 1 into E power I theta right and so on. So this becomes sigma by 2 E power I theta minus E power minus I theta times UPQ plus 1 plus UPQ plus 1 equals UPQ is that okay and what is this term 2 I sin theta. So if I find G now G equals 1 divided by 1 plus I sigma sin theta what is G it is UPQ plus 1 divided by UPQ is that fine I sort of quietly skipped a step there and therefore mod G has been less than 1 right or 1 is less than 1 plus sigma squared sin squared theta and we are very happy yes we finally have a scheme that is unconditionally stable. So in the next class we will see what this means see what we have done is that fine okay thank you.