 So, in our last lecture we looked at the way of using Taylor series approximations to discretize OD boundary value problem. So we are able to convert a differential equation and the two boundary conditions into a set of algebraic equations. So the problem got transformed from the original infinite dimensional space to a problem on a finite dimensional spaces which was set of linear or non-linear algebraic equations which have to be solved. So today let us look at some examples and let us see whether we can extend this concept to solving partial differential equations. So I wrote a generic boundary value problem, discretization of, so we had this original problem was or I would say y is equal to T of x, this original problem was, so you have this differential equation which holds over entire domain 0 to 1, this is the differential equation and then you have two boundary conditions, one is that, so you have boundary condition 1 which is f1 and you have second boundary condition which is f2, we have two boundary conditions. Then what the way we proceeded next was we discretize the domain, we discretize the domain, we mark these grid points such that using equidistant grid points we discretize the domain, we mark the grid points and then we had this notation that a dependent variable u, ui is equal to u of, ui was the value that this dependent variable takes at z equal to zi that we denoted as ui, note that this u here is the approximate solution, the real solution of this would be some u star z which would be a continuous function, twice differentiable continuous function and we are not able to in general find out this solution in many situations where the operator is not linear, we are not able to find this solution, so we want to discretize and solve this problem numerically and not analytically. So we have created this grid points, we have marked the dependent variable value at grid points as ui and then we converted this boundary value problem, the differential equation first we converted into a set of algebraic equations, so that was done by psi, so we approximated the second derivative using Taylor series, we approximated the first derivative using Taylor series, so at all the internal grid points, see here these are internal grid points z2, z3, z4 up to zn these are all internal grid points and then z1 and zn plus 1 are two boundary points, z1 and zn plus 1 are two boundary points whereas, so at all the internal grid points we enforce this equation, please note that enforcing this equation does not mean we are solving the exact problem, we are solving an approximate problem and for the approximate solution we enforce this equal to 0, there is an alternate way of solving this problem is to create two hypothetical grid points, one on the either side of the end points and then discretize, either way is fine if delta z is small, so for the time being let us write the equation in which we use one sided derivative approximations, so this what I get here is the approximation, this approximation is nothing but in the notation that we used earlier it is nothing but y cap is equal to or y tilde is equal to t cap x tilde, so we started with some problem we discretized and we got this discretize problem, so we started with a differential equation and two boundary conditions, we got set of algebraic equations, so the transform problem is completely different from the original problem, we hope that if you choose delta z is small then the Taylor series approximation is a good approximation and then this transformation will give you a solution which is close to the true solution, this will not this way of transforming will not recover the original solution, it will only recover approximate solution, if you reduce delta z then you will get better and better solution but the price that you have to pay is that smaller the delta z more number of equations and then you will realize that these equations are coupled for ith equation requires u value of i plus 1 and i minus 1, so these equations are coupled, so if you let us say take 100 internal points delta z is very very small if I take 100 internal points then I will have 100 equations here and then 2 additional equations, 102 equations in 102 unknowns to be solved simultaneously, how do you solve this Newton Raphson method or if these turn out to be linear equations solve them using you know ax equal to b simple Gaussian elimination, so you are transforming the problem and then solving it let us look at some examples today, so this is the background and we catch up from here, so my first example is going to be a tubular reactor example, so my first example is going to be T R AM or tubular reactor with axial mixing, I am not going to derive the equations I am going to just state the equations and then we will see how to discretize that is the key thing here, so this is a simple reaction where a goes to b, a very very simple reaction where a goes to b and the ODE BVP that is given to you is 1 by Peclet number into d 2 c by d z square by d c by d z minus this is Danckworth number into c square is equal to 0, so this should be obeyed between 0 z 1, okay this is my differential equation and then I have two boundary conditions, so my boundary conditions here, so my first boundary condition at z is equal to 0 the derivative at z is equal to 0 is Peclet number into c concentration at 0 minus 1 and my second boundary condition at, so this is my second boundary condition I am just going to apply the method of finite difference, this method that we are developing is called as finite difference method because we have developed approximation of first and second order derivatives using Taylor series approximation we have developed finite difference approximations okay, so now my first task what is this operator psi here, this operator psi here is nothing but this differential equation okay, so this second order derivative I am going to replace by its approximation at any ith point okay, same thing is true about this derivative I am going to replace approximation of this derivative at any ith point and so I am going to force the residual the term that you get is called as residual the residual is forced to 0 at all the internal grid points at the two boundary points I am going to use the two boundary conditions, so this equation will now be replaced this equation will now be replaced by 1 by Peclet number, so this will be ci plus 1 minus 2 ci plus ci minus 1, now what are these ci, c1, c2, c3, c1, c2, c3 are concentrations at the grid points, so at this point the concentration is 1 at this point is 2, c3 at in general at this point is ck, ck plus 1 and so on cn and cn plus 1, so these are the dependent variable values at the grid points okay, at any ith location at any ith location I am discretizing the differential equation, so the first derivative at ith location the first term that is second derivative is replaced by this then minus ci plus 1 minus ci minus 1 divided by delta z divided by 2 delta z then what is the next term, da ci square is equal to 0, this algebraic equation this converted algebraic equation should hold at all the internal grid points i is equal to 2, i is equal to 3 up to i is equal to n okay, so suppose I create let us say total grid points is 100, so 98 internal grid points 2 boundary points okay then I will get 98 equations like this okay each equation is related to the neighboring values of each equation is related to the neighboring values of c okay, how do I discretize this, what is my first boundary condition, my first boundary condition is this will be c1 this will be c2 minus c1 by delta z is equal to Peclet number into c1 minus 1, so this is one more equation this is one more equation I have just I have just converted this into a difference equation finite difference method alright I have converted the first boundary condition into a difference equation now I will convert the second boundary condition into a difference equation, so this will be cn plus 1 minus cn by delta z is equal to 0, how many equations I have now not 3 n plus 1 equations how many unknowns I have n plus 1 unknowns, if I solve this n plus 1 equations in n plus 1 unknowns I will reconstruct a solution an approximate solution of the boundary value problem okay, difficult to solve analytically because of this ci square and in general I have taken very very simple reaction okay it could be 2a going to so 2a going to be something like that, so I have taken a very simple reaction this could be a more complex equation here okay, so this n plus 1 equation in n plus 1 unknowns I have to solve them simultaneously to now I want to point out something else here in addition to talking about discretization I want to talk about some nice structure that appears okay I am going to rearrange this set of equations and then write them in a specific form to give you inside about something that will be talking about little later I am going to talk about a sparse system, now this particular equation is a very nice equation if you look at if you look at it there are n plus 1 variables but only 3 variables appear in one equation only 3 variables appear okay and when you have large number of equations you can actually make use of this fact that only 3 variables are appearing in each equation these kind of equations will give rise to what are called as sparse systems okay so we will side by side as a side note I will also introduce this idea of a sparse system and then later on we will be looking at special algorithms to deal with the sparse systems but right now it just to give you an idea which is so this first equation and this together I am going to rewrite in a matrix form which will give me a sparse equation okay so if I if I collect all the terms of you know ci, ci plus 1 and ci minus 1 together then I can write this equation as alpha ci plus 1 plus sorry minus beta ci plus gamma ci minus 1 is equal to da ci square okay where alpha is 1 minus 1 by delta g square okay I get equation like this n is equal to 2 3 up to n okay I will just group the terms together I will just group the terms together okay and then we have two more equations we have two more equations coming from the boundary conditions I am going to combine this and I am going to write this as a matrix equation okay so if I include the two additional boundary conditions okay what are the boundary conditions one boundary condition will give me right one boundary condition will give me this equation the other boundary condition will give me ci n plus 1 is equal minus cn is equal to 0 so these I have n plus 1 equations in n plus 1 unknowns and I am going to write them in a specific form I am going to write them in this particular form so this particular equation will be set of n non-linear equations in n unknowns you might be wondering initially when I started talking about Newton Raphson method in multiple multivariable domain where do you get this where do you get these equations classic example okay one boundary value problem is giving me if I take 100 grid points I will get 100 equations okay so large number of equations this equal to see this is like this is like matrix A this is the matrix A operating on vector x which is f of x you have to solve ax equal to f of x what is x here x here is this vector c1 c2 cn plus 1 in abstract form this equation this equation one matrix equation is nothing but matrix A operating on x gives me f of x I have to solve this I have to solve this problem or in other words I have to solve ax minus fx equal to 0 ax minus fx equal to 0 you can solve this by Newton Raphson you can solve this by simple iterative methods for example I could I could solve this by a simple method which is iterative method of this form ax k plus 1 is equal to f of x k I start with a guess x not I start with a guess x not so I start with a guess of concentration profile what is x don't forget that this is concentration profile this is where your input as a chemical engineer will come into picture you have to give a sensible concentration profile if you suppose give here negative numbers it is not going to work so giving initial guess somebody had asked is it very important it is very very important how do you give initial guess okay so this is my concentration profile discretize concentration profile and I can give a guess if I put this guess x not here I can generate the new guess by solving ax equal to fx not then x1 I can substitute get x2 x2 I can substitute get x3 and so on and then I can see whether you know difference between this and this is going to 0 or not by checking norm of this and see whether algorithm converges okay this is one way of solving right now that is not important what is important is that if I solve it like this this a I have to solve ax equal to b see because once I put x0 here okay this f of x is a known value okay I have to solve for x1 by solving the linear system ax1 is equal to okay then I have to solve another linear system ax2 is equal to fx1 and so on so I have to solve large number of linear systems linear algebraic equations now if I come back here if I come back here okay and then if I look at this matrix there is something special about this matrix this is a sum of three diagonal but the upper one is not strictly diagonal so this is a this is a matrix this is called as a tri diagonal matrix it's a banded matrix there are a lot of zeros how many elements this matrix has let's say if I take if I take 100 grid points how many elements this will have 100 plus 100 if I take 1000 grid points this will have 1000 plus 1000 and as we said more than number of grid points better is the approximation so I am forced to take smaller delta z if I take smaller delta z more number of equations matrix will have large number of elements okay if I develop some special methods that can deal with this banded matrices then this iterative calculations will become very easy okay so later on we are going to look at this methods for sparse systems this is this is a sparse so we will develop special methods for dealing with sparse linear systems okay which will reduce the computation because suppose you need 200 iterations 200 times if you have to invert a 1000 cross 1000 matrix not a nice thing to do okay instead if you develop a method which takes into account that there are lots of zeros okay you will be able to do calculations very very fast that's where all these sparse matrix method will come into that is a side note we will need this very often and when we do the discretization using finite difference you will see that almost every time you hit into this sparse matrices sparse matrix okay some form of the other okay so let's move on to some other example but it is not necessary that this finite difference method can be used to solve only boundary value problems it can be used to construct approximate solutions of partial differential equations as well okay I will take the same tran problem and show that how you can convert this problem a partial differential equation right now I took the steady state problem I will look at the problem which is time varying now in this case we got algebraic equations in the partial differential equation we might end up with something else okay so what is it that will end up with so my second example is a PD say same unsteady state same tubular reactor with axial mixing same tram okay with except that I am going to consider the unsteady state condition not the steady state condition okay so here so let's say this is a reaction in which 2A goes to B and then we have this dou C by dou t that is rate of change of concentration inside the reactor okay this is 1 by Peclet number now unlike the previous case where I was looking at the steady state behavior I am looking at a transient behavior I am looking at a unsteady state behavior so dou C by dou t is not 0 okay the concentration is a function of time and space not just space okay earlier we looked at the boundary value problem that came that arises when you look at a steady state of behavior of this particular system okay I still have the same two boundary conditions but I will also have an initial condition now okay there are two boundary conditions at z equal to 0 z equal to 1 and is also initial condition in time okay so at t equal to 0 so at t equal to 0 I have some concentration profile I have some initial concentration profile okay inside this is given by fz okay this is my initial condition there is some concentration profile inside and then I have two boundary conditions sorry this is not 1 this is t so this boundary condition is at z is equal to 1 this boundary condition is at z is equal to 0 so the solution now solution now is actually function of two independent variables one is time one is space another is time at z equal to 0 for all time this boundary condition holds at z is equal to 1 at z is equal to 1 at all time this boundary condition holds okay I want to now discretize this particular problem using finite difference method okay using finite difference method what I am going to do is I am going to discretize in space and leave the time untouched for the time being how will I solve the resulting problem something that we will see later but we will convert into a standard form which can be solved using a standard tool okay so in this case this partial differential equation I am going to convert into a set of ordinary differential equation okay I am going to discretize in space not discretize in time okay so what I will get is a set of differential algebraic equations what I will get is a set of differential algebraic equations so how will I discretize this okay the same trick that we did earlier in space we are going to you know we are going to denote this grid points z1 z2 z3 and then you have this c1 okay again I have discretized I have discretized my domain except now the dependent variable c that is concentration at the grid points is not only function of space is also function of time okay I will discretize only in space so c1 so the discretization in space appears through these indices 1 2 3 4 5 up to n plus 1 each one of them is a continuous each one of them is a continuous function of time so I am going to only discretize in space okay so now what will be the residuals I write I enforce the partial differential equation only at the grid points actually the original partial differential equation holds at every point inside the domain at every time we are not able to do that because we will get large number of equations okay if you start infinite number of equations if you start forcing it at every point we have to discretize and okay so now if I discretize this I will get dci by dt see the differential equation is forced at the ith grid point okay what will be this 1 by Peclet number into ci plus 1 t minus 2 ci function of time well minus ci plus 1 t minus ci minus 1 t by 2 delta z what I got here what I got here is set of ordinary differential equations I started with the partial differential equations I started with one partial differential equations one partial differential equation got converted into large number of ordinary differential equations okay so in this case the original operator t is partial differential equation the transformed operator is set of algebraic and differential equations why I am saying algebraic because I have two more conditions here okay one of them well in this particular case you will get two more differential equations what are the two additional differential equations the two additional differential equations that you get here are well this is a derivative in space so you will get algebraic equations here not the differential equations so what will be here c 2 t minus c 1 t by delta z I am taking forward difference approximation here and is equal to Peclet number into c 1 t minus 1 so this is one algebraic equation okay what is the second algebraic equation so what I get here is a differential algebraic system okay I got differential equations in time and I got two algebraic equations actually in this is a simple problem in this particular problem I can eliminate two variables and I can convert this problem into n minus one ordinary differential equations I can treat it like that or alternatively I can treat those ordinary differential equations and this two algebraic equations together into a single DAE differential algebraic system at towards the end of the course will be looking at how to solve differential algebraic systems so there are two ways to go from here one is to use these two constraint eliminate c 1 and eliminate c n plus 1 you get equations only in terms of only in terms of c 2 to c n n minus one differential equations in n minus one unknowns and then well you cannot stop here you also have to discretize initial condition this is a differential equation it will need initial condition so you cannot just stop at here discretizing the boundary conditions what about the initial condition what will be the initial condition so initial condition will be c 1 t c 1 0 is equal to f z 1 c 2 0 is equal to f z 2 and c n plus 1 0 is equal to f so these are the initial conditions these two algebraic constraints coming out of the boundary conditions and then n minus one ordinary differential equations have to be solved simultaneously okay have to be solved simultaneously so we have transformed the problem which was a partial differential equation which is difficult to solve analytically because you have c a square term appearing okay analytical solution could be difficult in this particular case in general if you have a more complex reaction rate equation it will be very difficult to solve this problem analytically you have to solve it numerically okay so you have to convert this problem into set of ordinary differential equation initial value problem okay or differential algebraic equation with initial value specified okay so these equations have to be solved simultaneously the transform problem is different from the original problem so remember this let us look at some other partial differential equations okay is the idea clear what is happening here here we transform a partial differential equation into a set of ordinary differential equations earlier a boundary value problem was transformed into set of algebraic equations which were non-linear which have to be solved simultaneously remember this set of ordinary differential equations is not a linear ordinary differential equations this has to be solved numerically using some method like Euler method or whatever we will be developing these methods later on right now I am just worried about problem transformation okay I am only worried about problem transformation yeah no I have n plus 1 variables and I have n plus 1 equations out of this n minus 1 equations are differential to our algebraic okay so I can eliminate two variables from this ordinary differential equations using the two algebraic constraints and convert into n minus 1 ODE initial value problem that is possible in this particular case it is possible if your boundary conditions had some non-linearity it is not possible to convert into very easily into ordinary differential equation then it will be a DAE system it has to be solved as a DAE system okay my third example is going to be a partial differential equation in two dimensions okay so in two spatial dimensions so this is a model of a furnace okay temperature distribution in the furnace the three walls of the furnace are insulated and at a constant temperature idealization okay and then there is a you know convective heat transfer from one of the faces I mean if you ask me how do I model this particular room temperature distribution in this room I would use equation something like this well why why is this like a furnace equation well you can see that this phase is not insulated okay we could say that from this phase there is convective heat transfer to outside okay these three walls let's assume as an idealization are insulated perfectly at the constant temperature where is the heat being generated well each one of you is like a bulb of 40 watts so there are 40 students into 40 watts so there is so much heat being generated inside this room okay so the temperature inside this room is function of two variables x and y okay the way the heat is distributed the heat sources are distributed all of you are sitting along different places okay this is like a so what is this model so this is the famous Laplace equation okay let's initially consider the steady state problem all of you are perfectly generating the same amount of heat okay at all time there is perfect steady state okay now there will be four boundary conditions here at four different boundaries right so we said that the three boundary conditions are so this equation should hold between 0 x 1 I have normalized the length okay between 0 and 1 so this partial differential equation should hold at every point at every point inside this room okay and then I have four conditions my four conditions are x is equal to 0 my t is equal to t star insulated boundary okay then x is equal to 1 t is equal to t star okay when I say t here which means t along the boundary okay I'm not I'm taking shortcuts and not writing the full then at y equal to 0 I have t is equal to so I should write t x 0 is equal to t star so accordingly I should write even these equations okay t at all x so three boundaries have a constant temperature and the fourth boundary so I have convective heat transfer at the fourth boundary at y equal to 1 I have k so I have four boundary conditions here I have four boundary conditions here okay now how am I going to discretize this I'm going to use Taylor series approximation in one dimension at a time okay in one dimension at a time so I'm going to I'm going to discretize this domain okay by constructing grid points or my domain is something like this I'm going to construct grid points here so this is my one two three four and so on this is my okay so my finally I'm going to get nx plus one grid points along x direction and let's say n y plus one grid points along y direction one to n y plus one one to nx plus one so so if you draw these lines here parallel to the so in general in general I'm concerned about forcing the differential equation at some i z point okay I'm constructing a grid here I'm discretizing see earlier we had only one special dimension we discretize in only one special dimension now I'm going to discretize in two special dimensions x and y okay and then these partial derivatives are going to be approximated in two special dimensions so my partial differential equation where should this where should this partial differential equation hold everywhere actually it should hold at all the points okay all the points even for one dimension we created hundred grid points suppose you create hundred like this and hundred like this how many internal grid points will be there hundred cross hundred okay we'll have hundred cross hundred internal grid points and you want to force the partial differential equation at every grid point at every grid point right so in general at i z point I can write this equation as well we are going to use this notation ti j is equal to t x i y j so temperature this is to simplify the discretization process I'm going to use t i j is temperature at x i y j okay so this is ith vertical line and this is zh parallel line okay at this point so I'm going to write this equation and so what will be what will be the the first one will be ti plus 1 j minus 2 ti j plus ti minus 1 j by delta x square plus ti j plus 1 see I'm taking partial derivatives in x direction partial derivative in y direction 2 ti j plus 2 okay how many such equations I'll get i going from 2 3 2 nx j going from 2 3 2 ny how many algebraic equations I get from a one partial differential equation I'll get nx minus 1 into ny minus 1 right okay where are the additional equations coming from boundary conditions so I should discretize the boundary conditions and then take all these equations together okay how many equations and how many unknowns nx plus 1 cross nx plus 1 100 cross 100 if you take 100 grid points very modest requirement just imagine this room 100 grid points not not too many if you want a better solution probably you should go 1000 cross 1000 but the number of variables will be very very large I'm not I'm not getting into the physics of the problem given this problem we want to solve it let's so now we will require additional equations right at the boundaries so what are the boundary equations no no let's not let's not get into the physics you do that in other course transferred course ask that question I'm just worried right now about problem discretization okay so now say this this will give additional equations so I will get t star at x is equal to 0 x is equal to 0 I will get all these points right along y so I will get t 1 j is equal to t star j going from 1 to ny plus 1 then the second boundary condition will be t nx plus 1 j is equal to t star j going from 1 to ny plus 1 okay then you discretize the third boundary condition okay and then this fourth boundary condition you discretize using forward difference okay sorry backward difference we discretize this using backward difference okay at all the grid points at all the grid points along x okay yeah mixed derivative you don't get second order derivative of the boundary so why don't we develop an approximation a difference equation can we develop a difference equation so you take first derivative in x then you take derivative in y so you will get large number of equations you'll get large number of equations which are coupled in general and they have to be this will also give rise to a sparse system you can see that any equation here how many equations how many variables appear only neighboring four variables appear if you take t ij ijth point there are only five variables appearing how many variables we have we have large number of variables nx plus 1 cross ny plus 1 but in one equation there are only five variables you can expect that if you try to solve this numerically you will get sparse system okay you get a sparse system and then of course we will be looking at solving sparse linear systems or sparse non-linear systems separately but the is the idea clear how do you do discretization okay so in this case the original problem is a partial differential equation discretize problem is a set of non-linear or linear algebraic equations which have to be solved simultaneously which have to be solved simultaneously so the transform problem is different the solution that you get from the transform problem is an approximate solution okay not the original solution so in the next class we'll see one more example one more way of discretizing this and then we'll start with some other methods of discretizing