 So welcome back everyone and we are at the last topic in the first part of our workshop. So far what we have done is we have hopefully built the necessary background that is required to get on with the numerical solution techniques. The idea was so far that some basic fluid mechanics at least minimally as is required will be covered and especially the governing equations which are normally used in typical flow situations, single phase laminar of course should be familiar to everyone. Again let me take a minute to point out that if you read our objectives for this particular workshop carefully we had decided that we want to give out information to put together a first level CFD course and therefore many of the complications that normally we see in real life situations we cannot really handle. One of the examples is for example turbulence. So we are not really talking about turbulence modeling at all and normally turbulence modeling is taken up in a advanced CFD class after having understood the basic CFD methodologies. Then we get on to other activities like turbulence modeling or several questions have already been asked in the last few days on for example two phase flows. So these are situations which are far more complicated from simply the physics involved in them and unless a clear understanding of the physics involved in them is not there, there is no point really getting into the CFD side of let us say turbulence or something like two phase flow modeling or even a multi scale modeling nowadays and simulation as they do now. Our objective is that way very simple. A first level CFD course which you can actually offer in your colleges and institutions typically at the final year undergraduate level or at the first year post graduate level and that is why we have restricted ourselves to only the standard single phase fluid flow equations and this is what I wanted to emphasize again before proceeding. So getting on with our first type of numerical solution method what we will talk about in the session today now and the first session tomorrow is what we call one way of solving these governing equations using the computer and this is what we will call the finite difference method or a finite difference methodology. The slide number one that is projected right now is a general slide in the sense that it has nothing to do with finite differencing as such. This talks about the general numerical solution methodology. I think we have talked about this a little bit but then since we are now starting to get into the numerical solutions let us look at it again. Any numerical solution process begins with the governing equations, the choice of governing equations. Now the governing equations are the ones that we have already seen plus there can be more if required the choice is really up to the modular as to what sort of set of governing equations to choose in order to simulate a certain situation. As we have already seen the governing equations are partial differential equations they are coupled they are non-linear and clearly the computers do not solve differential equations in that sense. So we need to actually convert those differential equations into something that the computer handles and solves which is essentially algebraic equations. And the process of converting the governing equations into algebraic equations is what is called as discretization and this discretization is usually carried on a computational mesh or a computational grid. There are different methods of discretization meaning different methods by which you can achieve algebraic equations from the governing partial differential equations. And the first method which is a very classic and old method in some sense called the finite differencing is what we are going to discuss briefly in the session today and tomorrow. Please keep in mind that the focus of the course later starting from tomorrow late morning session will be on the finite volume technique of numerical solution. So the discretization that we are going to really focus on is the finite volume discretization. However to at least get going into this numerical solution methods we will begin with the finite differencing. It is quite likely that many of you are already familiar with this finite differencing. So for them this will serve as a very brief and fast introduction. So anyway once you obtain the algebraic equations using this discretization process you need to solve them and then there is usually some sort of a solution algorithm available to solve these algebraic equations after which you obtain the so called a numerical solution rather which is essentially a set of numbers each set of numbers corresponding to one discrete point let us say within the domain. So the domain of interest in which you are doing this simulation is divided into a set of discrete points and the numerical solution is obtained in terms of numbers pertaining to various these points. So in order to make sense of these numbers usually what we do is what we call post processing and I was mentioning a couple of times that one of the post processing activities is perhaps plotting the streamlines for a flow solution. So that it gives you an idea of how the flow is looking, how the flow is behaving. So that is one standard post processing technique. Other post processing techniques will involve looking at contour plots looking at velocity vector plots and then maybe you can derive secondary fields such as a vorticity field and plot that vorticity field in terms of vorticity contours and so on. So most of these details you will go through the next week when formally the finite volume technique will be discussed. In the session today and tomorrow what we are going to do is we are going to just about introduce this finite differencing method and that too not in the general sense what we are going to do is we will take our two last problems that we worked out in our exact solutions. If you remember one was this unsteady flow between two parallel plates and the other one was a two dimensional steady state conduction. So with respect to those two problems we will simply outline what a finite differencing approach does so that it is a fairly restrictive discussion and again the reason for this restrictive discussion on finite differencing is because finite volume technique is what we are going to focus on. One of the reasons is that most of the commercial codes that are available today employ the finite volume techniques so that makes a case for the finite volume techniques to be incorporated in the course. So let us get on with it. Let me first take the first problem which was our transient flow between infinite parallel plates problem which was the last but one problem in the exact solutions discussion and here if you remember through some initial simplifications we had obtained a governing differential equation which looked as it is projected right now it is the partial derivative of u with respect to time is equal to the kinematic viscosity times the second derivative of u with respect to y plus a source term and that source term if you remember involved the pressure gradient so I think it was minus 1 over rho dp dx and we had said that this is going to be a constant as far as our solution is concerned. So the governing equation for our finite difference technique discussion to begin with at least is going to be this unsteady diffusion type equation with a source term. If you remember the boundary conditions are the u velocity at y equal to 0 which is which was the bottom plate for all times was equal to 0 and the u velocity at y equal to h which is the top plate at all times was equal to capital U. The initial condition was that the u velocity everywhere that is for all y's at time equal to 0 is equal to 0 so we begin with everything stationary and then at time equal to 0 set the top plate in motion in the x direction with the velocity of capital U and then we are basically going to figure out how the velocity field evolves in time. The analytical solution of this problem is already available with us from our discussion in the previous session. So if you look at some of the standard terminologies that CFD people use this is what can be considered as a time marching problem meaning that at time equal to 0 an initial condition is provided and the solution evolves in time going from time 0 to a later time. And as I said this present governing equation that we are seeing is what is called in the applied mathematicians terms a diffusion equation. If the source term is not there it is really the diffusion equation but you can say that it is a diffusion equation with a source term and the source term is a constant. So now what do we do in a finite difference solution? If you recall the procedure we convert governing equations which are in partial differential equation form into an algebraic equation. So here what we first do is we divide the domain which in this case is only a one dimensional domain in the y direction. So we divide the y direction extent from y equal to 0 to y equal to h into a number of node points or grid points as they are called. So in this case I am showing the grid points as 1 at y equal to 0, 2, 3 etc. with capital N as the last grid point or the last node point at y is equal to h. As I mentioned earlier also the node points or the grid points are the points in the domain where the numerical solution will be obtained. So in that sense there will be no information available for the solution between 1 and 2, between 2 and 3 and so on. We will simply have a representative value in some sense at 0.1, 0.2, 0.3, 0.4 and so on. The first point here is a boundary node at y equal to 0, capital N is another boundary node at the other boundary and at present these grid points have been placed with a uniform spacing between them. There is no need really to have a uniform spacing but we want to keep things as simple as they can at least when we begin and that is the reason we are keeping the grid points separated from each other at a uniform spacing. So this delta y is what we are going to call grid spacing and is uniform as I said in the present case. So that is the grid setting for our problem at hand. It is very straightforward because it is a one dimensional problem. So what we do in case of finite differencing approach is let us just consider one representative node point in the interior of the domain which I am going to denote as I here. It has a right neighbor which is I plus 1 at a distance of delta y which has a left neighbor as I minus 1 again at a distance of delta y and we will try to come up with appropriate algebraic expressions which we will replace the original partial differentiations or partial differential in the original PDE with and then somehow solve the algebraic equations. So the technique to achieve these algebraic equations in the case of a finite difference method is to employ a Taylor series expansion which we are quite familiar with now to express the value of variable at a neighboring point with respect to the value of the variable at let us say the point under consideration and therefore if you want to look at the present problem u is our variable. So I want to express u at the I plus first point which is the right neighbor in some sense for I and that is expressed as a standard Taylor series expansion where we write u at I plus 1 is equal to u at I plus the derivative of u with respect to y at I times the spacing between these two points which is delta y plus the second derivative of u with respect to y at I times delta y squared over 2 factorial and so on. So it is up to you as to how many terms you want to add in this Taylor series expansion. We are going to restrict ourselves to add the most two terms primarily because you will typically see that these are the standard expressions that people will employ for finite differencing and we want to keep the computation time to a somewhat manageable time so that we do not want to unnecessarily keep on adding additional terms in the Taylor series expansion but the choices completely up to the user as to how many terms are to be maintained in a Taylor series expansion such as this. So equation number one here is expression of variable value at I plus 1 in terms of variable value at I. Equation number two likewise is variable value at I minus 1 which is the left side neighbor in terms of the variable value at I and so in that sense we are moving by a distance of negative delta y. So therefore, u at I minus 1 is u of I minus du dy at I times delta y because there is a square of the distance involved it is plus d 2 u dy squared at I times delta y squared over 2 factorial and so on. So with these two we are in a position to proceed further. So if we add just equation 1 and 2 what you will see is that on the left hand side you will generate u of I plus 1 plus u of I minus 1. First term on the right hand side will give you u i plus u i that is 2 times u i. The second term on the right hand side is actually going to cancel it because it is the same with the sin inversion and the third term on the right hand side will be d 2 u dy squared times delta y squared because 2 factorial is 2. So half plus half will give you 1 and that is what we are going to see here. Adding 1 and 2 you get u i plus 1 plus u of I minus 1 is 2 times u i plus the second derivative of u with respect to y at I multiplied by the square of the grid spacing and if you rearrange this you will write an approximation in some sense for the second derivative of u with respect to y at the grid point I equal to u i plus 1 minus 2 u i. So I am bringing this 2 u i on the left hand side plus u i minus 1 which already existed there whole thing divided by delta y squared and I am adding something which actually I am not going to talk about because again going back to the first day when we outlined the objectives we said that we want to actually make sure that we focus on the numerical methods and not the numerical analysis. So we are really not going to talk about much on the order of accuracy of these kinds of approximations and so on. But those who are interested should always take up reading in some books and depending on how many terms you have kept in your Taylor series expansion you can come up with what is called as an order of accuracy for the approximation that you are generating. Now remember that we are generating an approximation here for the second derivative of u with respect to y at the grid point and that approximation is given by the expression on the right hand side plus what is called as the truncation error which simply represents the terms that you have discarded from the Taylor series expansion. So it is a very crude way of saying what is truncation error but I will leave with that with the understanding that we are focusing more on the methods and not on the numerical analysis where all these errors and stability criteria and so on will be very rigorously worked out. So I had been at the beginning of this workshop a list of books and in that list of books there are two very well written basic CFD books which involve this finite differencing approach. One was by John Anderson and the other one was Hoffman and Cheyenne. Either of those will actually provide lot of these truncation error type details and those who are interested in pursuing finite differencing later are requested to look into those books because for the purpose of this workshop we are not going to focus on those terms. Just for the purpose of information I have included this truncation error to signify that we have discarded terms of the order higher than second order in the Taylor series expansion. So that is as far as the second derivative of u with respect to y at a given grid point is concerned we have generated an approximate expression for that. As I said I am going to stick to this particular equation with which we are discussing everything. So you see the first term on the right hand side is where this second derivative of u with respect to y occurs and that is what we have been talking about so far. The only other partial derivative that occurs in this equation is the time derivative of u on the left hand side. So we need to come up with an approximation for the partial derivative of u with respect to time. To do that we first realize and understand what is meant by this time marching problem. So what I mean by this is let me just quickly draw a simple sketch on the whiteboard. So this is at time equal to 0 let us say we have all these grid points in the domain. From time equal to 0 we will move to time equal to delta t let us say we move by one time step and thereby at all these grid points 1, 2, 3, 4, etcetera the solution will be updated from the old time step to a new time step. So that you can think about a derivative with respect to time as a new value of the variable at the grid point minus the old value of the variable at the grid point divided by the time step or the time difference between these two values which is delta t and that can be considered as one expression for generating an approximation for the time derivative. So let us go back and see what we are doing here now. So on the left here what I have is u suffix i at time t plus delta t which is essentially the value of the variable at a given grid point at the new time step clearly is expressed in terms of the value of the variable at the same grid point at t which is the old time step plus du dt which is the derivative of u with respect to time at the grid point of interest multiplied by delta t which is our first order Taylor series expansion and that is it. I have chosen to use only one expression or one term in the Taylor series expansion. If you rearrange this I will obtain an approximation for the time derivative of u at the grid point i as u at t plus delta t minus u at t divided by delta t which is what was shown on the sketch a minute ago and rather than now writing the time steps as t and t plus delta t etcetera I will start following the perhaps standard nomenclature of using a superscript n. So n n would imply that we are talking about a given time level n plus 1 is the next time level. So going back to my sketch in general you can say that this is nth time level and the n plus first time level which is same as new and old in our earlier terminology. Just like we have gone ahead in time and formed the approximation for a time derivative at a given grid point we can decide or choose to go backward in time. So with respect to a given level I can choose to go backward by a step of delta t and then write u at the ith grid point for the time step t minus delta t is equal to u at the same grid point at the given time level t minus the time derivative of u at the grid point multiplied by delta t again choosing to use only the first term in the Taylor series expansion and if you rearrange this utilizing our superscript notation of n I will have u i at n minus u i at n minus 1 which is the previous time step over delta t and the terminologies that will be employed here are a forward time difference and a backward time difference. So this is what I will want to want to write. Actually you can point out here that these are exactly the same if you just replace n by n minus 1 you will get the same thing is what you can say and actually that is right in that sense there is no difference between these two unless you club these along with the other term you do not really have any meaning by themselves to call them either a forward difference or backward difference. So let us see what that means having done this we have our two expressions either we go forward in time or we go backward in time to obtain the time difference. Remember that we have partial derivative of u with respect to time equal to nu times second derivative with respect to y of u plus s the source term we have obtained some approximations for the time derivative we have obtained an approximation for the spatial derivative now we will put it together and see what happens. We adjust to complete that discussion the spatial derivative that we obtained earlier for the second derivative of u with respect to y at the grid point is popularly called as a central difference. Obviously you can utilize only the two grid points ui and ui plus 1 or ui minus 1 and ui and come up with different type of representation but here the most standard one that we are using is what is called as the central difference. So now we will put together the equation in the sense that we can use the forward time approximation with the central difference or the backward time approximation with the central difference in order to replace the original partial differential equation into what we will call as a finite difference equation. So far I did not talk about at what time level we actually evaluate this spatial derivative and now is the time to talk about that so what we will say is that if we evaluate the time derivative sorry if we evaluate the spatial derivative at the given time level n and utilize this with the time difference we generate the first equation here which is what is called as a forward difference or forward time central space scheme. Why central space? Because the second derivative in space has been obtained with what is called as a central difference approximation whereas I am utilizing a forward time difference. So that is the way to look at it. So now if you see this equation here on the bottom the first equation on the bottom you will see that I have ui at n plus 1 minus ui at n over delta t is equal to nu multiplied by the expression for the central difference the approximation for the second derivative or the expression for the central difference evaluated at n the present time level plus our source term or I can choose to combine the backward difference and the central difference in the sense that the time derivative will be evaluated as the backward difference with the spatial derivative getting evaluated at the current time step. So you can see how the difference what is the difference between these two guys are in both cases the spatial derivative is getting evaluated at the current time step whereas in the top equation here we are marching forward in time in the bottom equation we have gone backward in time. So these two approaches for solving the same equation will be called as forward time central space and backward time central space. So now we will rewrite the equations in slightly different form before that though I think let us what I have done is I have taken this forward time central space equation and we will try to discuss that completely. So this is what the forward time central space form is and if you rearrange this to read the value of variable at a given grid point at the new time step which is n plus 1 superscript that will be given as the value of the variable at the same grid point at the previous time step or the current time step I should say plus nu multiplied by this delta t over delta y squared that is simply in the rearrangement whole thing multiplying whatever you see inside the bracket here ui plus 1 at the nth level minus 2 ui at the nth level plus ui minus 1 at the nth level plus the source term which existed here and it will get multiplied by delta t. So this is what we are going to call as an explicit method this forward time central space form is also going to be called as the explicit method in our case and we will see why we are calling this as explicit. So let us look at how things are happening. So let me go to my white board. So this is our expression for the explicit method or the forward time central space method. So let us try to draw the picture for two different time levels. This is the nth time level and this is the nth plus first time level and let me look at a grid point and it is two neighbors. So this is i, i plus 1, i minus 1 and this is the same ith grid point at a later time step. If you go back to our equation at the top the value of the variable at the ith grid point for the n plus first time step is expressed in terms of from the right hand side the value of the variable at these three grid points for the previous time step. In other words the information from these three grid points the ith grid point, i minus first grid point and i plus 1 grid point at the previous time step level is getting used to previous meaning I should say the current time step level is getting used to calculate the value of the variable at the ith grid point in the next time level. So that is what we had shown here on the sketch going from n to n plus 1 going back again to my white board. The reason this is called as an explicit method is because this is considered to be a time marching process as we said. Let us say that we start from time equal to 0. So let this nth level be specially time equal to 0 and then let this be time equal to delta t the first time step from time equal to 0. At time equal to 0 all these values of the variable will be available through the ith initial condition specification. So that the value of variable at ith grid point i minus 1 i plus 1 and all such grid points on the domain or in the domain will be known at time equal to 0. Knowing that you visit a given grid point let us say it is our ith grid point and you want to update the value of the variable at this ith grid point you have the expression available to you which is completely known in terms of the information. Do we know the value of u at i at the nth which is our 0 level? Yes we do. New delta t and delta y will be known because new is a physical parameter the kinematic viscosity which will be provided as part of the problem statement this is new. Delta t is a time step that we choose similarly delta y is the grid spacing that we choose. So these are the user inputs in some sense going in here all values of the variable at the nth time step which is our 0 time step right now or 0 time level all values inside the bracket here are known. The source is a known source which means that in this case it is going to be provided as part of the problem statement delta t again is the time step which will be user chosen and therefore this entire right hand side at a given time step will be entirely known and knowing that you can go and update the value of the variable at the grid point i and similarly you can do the same thing for the next grid point and the next grid point and the next grid point and so on. So we visit each grid point one by one so you can think about a loop that is operating within your program which will visit each grid point one by one and will update the value of the variable at each of these grid points using the expression here where it will use the values of the variables at the same grid point and the two neighbors on one on either side along with the source information etcetera and will calculate the updated value of this variable at the grid point of interest. So since the entire right hand side is supposed to be explicitly known at the given time step we call this approach as an explicit method. So you can imagine that if you have a if you want to think about this as a programming exercise we may want to define two arrays in our program if you want to keep at least track of the old time value and the new time value. So one array will be essentially called let us say you old and the other array you can call a new new and with these two arrays you can essentially go ahead in time or time march as they say starting from time equal to 0 and going in steps of delta t. So once the entire solution is updated for the domain at time equal to delta t we will then go to the level time equal to 2 delta t let us say in a similar fashion. Now that the entire solution is available at time equal to delta t using this known solution at time equal to delta t you will find out one by one again the solution at time equal to 2 delta t and so on in this manner the solution process proceeds and we calculate our data time evolution of the variable u within the domain. If you want to output the value at any particular time you can provide that information within your loop as you are marching in time. So when the time value is reached at some predetermined value you can instruct your program that write the values of u in my domain in an output file which you can later use to plot the value of the variable at that particular time step and so on. That is how the codes that have been generated in Sylab are written so that you want to know how the value of the variable u is evolving with time within the domain so that you need to store the values of u at various times starting from 0 and then you can plot those to see how the variation is evolving in time within the domain. So this is what we call an explicit method remember that this boxed expression is essentially the one line that you need to know as the algorithm for this case everything else is additional business as long as you utilize this correctly the explicit method will be coded that is absolutely not a problem. Other thing is that the boundary conditions will be separately stored. So the boxed expression is essentially for the interior grid points as we say or the interior node points in the particular case that we are talking about right now our problem of interest the boundary node 1 at y equal to 0 and the boundary node capital N at y equal to h will be assigned the value of variable depending on our boundary conditions. So that is what has been mentioned here that the boundary conditions are that essentially u 1 at all time steps is equal to 0. So it does not have to be necessarily n plus 1 it is essentially all time steps. So it is always 0 no matter what and similarly the value of variable at the right boundary node is equal to capital U. So that is always there and it can be hard coded right at the beginning because these are time independent boundary conditions so they will always be there or you can always in the loop also keep on including them every time as explicitly added information in the time loop so to say. Actually the explicit method is very very simple it is very easy to code there is hardly anything required as an expertise to code this and it is one of the most simple and elegant methods that one can think of. The only drawback is that unless you choose the time step in a problem such as this or in general when using an explicit method in a certain manner the solution will diverge or will explode and the method does not work. So that feature with which you figure out under what conditions the explicit method is going to provide you correct solution without exploding or blowing is given by what is called as stability analysis and that stability analysis is again part of the numerical analysis that we employ with these numerical methods. We are not going to get into the stability analysis per say in detail I have just mentioned the condition that the stability analysis gives for the present problem of interest the diffusion problem that we are talking about if one wants to use explicit method and it turns out that the condition for the explicit method to work without exploding or blowing is the kinematic viscosity times the delta t the time step divided by delta y the whole bracket square meaning delta y squared I should say should be less than one half. If this ratio turns out to be higher than one half then you will see that the explicit method does not work it simply crashes. The formal treatment of stability analysis is very elegant actually but we are not going to discuss that since our focus is not really on numerical analysis those who are interested in knowing the details of stability analysis can look into that Andersons or Hoffman's books and they can look at the details if they want but as far as the present course is concerned we are not going to get into the stability analysis of an explicit method. So now let us see what happens in the case in the other case that is where we had chosen to evaluate the spatial difference at the current time step and we are using the time difference in the backward fashion. So you are going from the current time step to the previous time step. So this is the last equation on the slide is what we will choose and we start looking at it. So this is again the same equation repeated at the top of slide number 8 right now and what I will do is just to bring it in a so called standard form I will replace n by n plus 1 so that n minus 1 will become n and n here will become n plus 1. So this is usually the way it is represented so that is the same there is no difference between the way it is represented here and it is represented at the top is just that this time index the superscript n is essentially a dummy index if you want or a dummy superscript and I am simply changing it from n to n plus 1 that is not really going to change anything. Now let us see what happens here if you rearrange this now we are going to see that the value of the variable at the ith grid point for the n plus first time step which is our next time step from the current time step is going to use the information on the grid point value of the variable at the previous time step but as far as the neighbors are concerned it is actually going to utilize the information on the neighbors at the time step where we are trying to evaluate the value itself. So basically the value of u i plus 1 here and i minus 1 at the n plus first time step is not really known and therefore the right hand side here in the equation that is presently written on the board and box is not explicitly known and therefore this method which is called as the backward time central space is also called as an implicit method because the right hand side is not known explicitly or completely let us say we have to evaluate the right hand side in some sense as part of the solution that is the way one can interpret these words implicit. So let us see how it can be done then we rewrite the backward time central space form or the implicit method form again this is exactly the same equation that is repeated and we rearrange this as some constant with a minus sign appearing multiplying the value of the variable at i minus first grid point at the next time step plus 1 plus some constant so that is another constant times u of i at n plus 1 time step minus the same constant that was there here times u at i plus first grid point at the n plus 1 or the next time step equal to the value of the variable at the ith grid point at previous time step or the current time step I should say sorry plus s which is our source times delta t the only thing that is known here is the right hand side which will be known at a given time step either if the given time step as is a initial time equal to 0 part of initial condition or at any particular time step we have already arrived and found the solution and now from there we want to go ahead so that the current time step information on u in the entire domain would be considered to be known and therefore the right hand side since the source is also known the entire right hand side will be known. What about these coefficients nu multiplied by delta t over delta y squared 1 plus 2 nu delta t over delta y squared etcetera these are also known because as we discussed earlier nu which is the kinematic viscosity will be provided as the problem statement delta t which is the time step which is chosen by the user so is the grid spacing delta y so that this entire constant you can treat it as known and you can evaluate it as 1 value a let us say. So, I am replacing this minus nu delta t over delta y squared with the symbol a similarly 1 plus 2 times nu times delta t over delta y squared I am using the symbol b and again this will come out as a same as what we had in here and so in a and of course the right hand side is the right hand side which is a known quantity at a given time step is given a symbol c just for the sake of brevity. So let us see now what happens in this case let me go to the whiteboard I should probably write this here as minus is better and this as plus that is a better way of doing it. So now let us see what this means so if I choose i equal to 2 because i is equal to 1 is our left boundary grid point and the left boundary grid point is not to be considered in this expression because the boundary condition at i equal to 1 is anyway known so the solution is essentially known at i equal to 1 so that we need not bother about. So if I use i equal to 2 I have one equation which let me write so i sorry a times i is equal to 1 means this is u 1 at n plus 1 plus b u 2 n plus 1 plus a u 2 n plus 1 plus b u 3 n plus 1 equal to c 2 n this really is known because this is the boundary condition u at 1 for whatever time step that you are talking about. So if you want since this is known you can transpose it to the other side and let us do that so b times u 2 n plus 1 plus a times u 3 n plus 1 equal to c 2 at n minus a times u 1 at n plus 1 which is known from the boundary condition so this is equation number 1. If I choose i equal to 3 I will generate a times u 2 n plus 1 plus b times u 3 n plus 1 equal to c 2 n plus 1 equal to c 2 n plus 1 plus a times u 4 at n plus 1 equal to c 3 at n this is 2. One more equation I will write a times u 3 n plus 1 plus b times u 4 n plus 1 plus a times u 5 n plus 1 equal to c 4 n 3 and so on. So you can generate one such equation for each of the interior grid points moving one by one let us say and you can always imagine that finally you will reach the right hand boundary node point capital N which is what we have called and for capital N the solution will be known again as the boundary condition so we do not have to bother about that. But let us look at what has happened if you look at this carefully you are basically trying to solve for let me try to underline those variables u 2 u 3 u 4 u 5 etcetera all the way up to u n minus 1 at the n plus first level and you are trying to solve this at the given or at the next time step n plus 1 you are trying to solve this at the next time step n plus 1 simultaneously because you do not know u 2 you do not know u 3 you do not know u 4 you do not know u 5 but you are generating all these equations which can be solved simultaneously at the time step n plus 1 and if you are able to solve these simultaneously you will basically obtain a solution of the value of variables u in all interior grid points starting from 2 to n minus 1 capital N minus 1 because u 1 and u capital N are known as boundary conditions. So having discussed this let us go and read the line at the bottom of the slide we have essentially explicitly worked out a few of these equations and realized that the boxed equation on the slide number 9 here will essentially imply that we are dealing with a system of linear algebraic equations for the unknown interior grid point values i equal to 2 to capital N minus 1 of the variable u at the n plus first time level. So what this implicit method or backward time central space method requires is a solution algorithm which will handle and solve this system of linear algebraic equations and that is what I have tried to show in a matrix form so to say what I have tried to show here is I am not now using the superscripts n plus 1 and n the coefficient matrix will contain these coefficients a, b and c solution vector as it is called is going to be a column vector which is going to involve u 1 to u capital N and looks like I have included here the boundary nodes as well for which the solution is actually known but that is fine. The right hand side will involve this constant c which are typically known at every time step. So the right hand side is known the solution vector is going to be determined as we solve this linear system at every time step. So for every time step the solution here in the case of this implicit method will require a solution of this entire linear system of equations from time step to time step. So let us just look at this matrix and make sure that the matrix representation is exactly as what we had written in a few minutes before when we were writing out those equations. So if you look at the matrix equation it is 1 multiplied by u 1 equal to 0 which is the boundary condition. So the first row multiplying the first column will give you only this 0 which is 1 multiplied by u 1 equal to 0 that is our boundary condition at the left node that is fine. Similarly immediately if you see the last row multiplying the column vector u will give us these are all 0 so nothing will come from there 1 multiplied by u suffix capital N is equal to u which is capital U that is the right side boundary condition. So the first row and the last row in the matrix equation is giving us the boundary condition so that is known and we are not going to talk about that. Let us look at the second row for example. So the second row multiplied by the column vector will give us a times let me go to my whiteboard from which I can write that b times u 2 plus a times u 3 equal to c 2 minus a times u 1 and this is exactly what we had written down a few minutes earlier as our equation number one as we had written earlier. Only thing is that in the equation number one when we were writing we were explicitly adding the superscript N plus 1 also. Now here it is understood that we are seeking the solution for the u variable at N plus first level whereas whatever is on the right hand side in these equations is at Nth level so that is the understanding other than of course this boundary condition which is known for all time steps. Let us look at one more equation from the matrix. This we looked as a multiplied by u 1 plus b multiplied by u 2 plus a multiplied by u 3 equal to c 2. Now 0 will multiply u 1 so nothing a times u 2 plus b times u 3 plus a times u 4 is going to give me c 3. Let me just write that which was our earlier equation 2 where we had of course the superscript N plus 1 and so on. So the idea is that a statement written as a times u i minus 1 plus b times u i plus a times u i plus 1 equal to c i appropriate superscript is essentially as i goes from 1 to N if you want or 2 to N minus 1 whichever way you want to select it is going to be a system of linear algebraic equations and that system of linear algebraic equations has been expressed as a coefficient matrix multiplied by the column vector of the unknowns which are u 2 u 3 etcetera all the way up to u n minus 1 equal to a right hand side where the right hand side is essentially a known right hand side at the given time step. So at every time step going from N to N plus 1 knowing the solution at the nth time step you solve this linear system to obtain the solution at the N plus 1 first time step and this process continues from time step to time step. So you can imagine that in the case of an implicit method such as what we are describing right now the amount of computation will be more as compared to the explicit method because we are solving this linear system from time step to time step in order to obtain the solution. So in that sense the computation is indeed more the advantage of this implicit method or the backward time central space method is that it is unconditionally stable which means that there is no restriction on the time step that we had in the case of explicit method. So this means that no restriction on delta t you can choose any value of the time step and the method will not explore the method will actually give you a solution how accurate that solution is a different matter but at least the method will not explore and blow up it will actually give you a solution. Let us just stop here for a minute and see what we have done. So we are talking about this finite difference approach the finite difference approach involves using these Taylor series expansions to replace the partial derivatives in an original partial derivative. So this is a differential equation with appropriate finite difference approximations as we have done here. So what you see right now here is a finite difference approximation for a second order sorry for a second order derivative d2u dy squared at a given grid point. Similarly we have a finite difference representation for the time derivative. We have put together the finite difference representation and in doing so we convert the partial differential equation into a finite difference equation and it is the finite difference equation that actually gets solved in the computer. In the case of this explicit method which is also called as the forward time central space method it is what we need to do is we just essentially visit one by one grid point and the value of the variable at the next step time step that is is updated based on the value of the variable at the previous time step which is known from our previously evaluated solution or it is the initial condition itself in an explicit manner meaning that the entire right hand side here is completely known at time equal to n and you can evaluate immediately the solution at n plus 1. On the other hand when it goes to the implicit method we realize that we end up generating a system of linear algebraic equations and the solution at the time step n plus 1 is obtained using the solution of this system of linear algebraic equations where the right hand side here as is shown in this matrix equation is known from the previous time step. So from time step to time step implicit method will involve repeated solutions of this linear algebraic equations whereas explicit method does not. Therefore the implicit method is more computationally extensive or expensive but it offers the advantage that there is no restriction on the time step that you can choose it will always work no matter what the time step is the accuracy of the solution as a time progresses is something that may not be good if you take very large time steps but otherwise the method will not crash as would the explicit method. What you realize is a special situation that is occurring in this particular matrix representation what you realize is that the only non-zero components in this matrix are aligned with the main diagonal and one diagonal on either side of the main diagonal and this is what is popularly called as a tri diagonal matrix, tri diagonal coefficient matrix and there is a very efficient solution algorithm for solving a linear system such as what we see on the board where the coefficient matrix is a tri diagonal matrix it is what is called as the Thomas algorithm or the tri diagonal matrix algorithm which is TDMA. So with that I will stop today. Thank you very much.