 So, we have been looking at finite difference method that is used for discretizing boundary value problems as well as partial differential equations and in particular in my last lecture I was looking at a partial differential equation Laplace equation and then we approximated the partial derivatives in spatial directions and then what results is a set of linear algebraic equations. So, we start with the partial differential equation and then process of discretization resulted in set of linear algebraic equations. The problem was so we had this partial differential equation so x and y are normalized spatial coordinates and then we have this partial differential equation that holds on the interior points and then I said this was a example of a furnace and then we use finite difference method to discretize this. What you get here is a set of linear algebraic equations what I want to show again with the same example is that it depends upon the method you choose to discretize okay it is not that this particular problem and finite difference method will always yield set of linear algebraic equations I am going to do one more method called method of lines which will yield ordinary differential equations boundary value problem. So, depending upon how you choose to discretize you will get different types of approximations from the same problem that is very important okay and using the same method same approach that is finite difference okay of course when we use some other method which will be discussing for example starting from towards the end of this lecture and next lecture we will talk about orthogonal collocation so that is another way of discretizing and of course will yield different way of formulating or formulating the approximate problem. So here we had these boundary conditions at x is equal to 0 we have condition that T x T 0 y is equal to T star then at x equal to 1 we have this condition T 1 y is equal to T star the third boundary condition is again so it is like saying that three walls are insulated then there is a convective heat transfer from 1 at y equal to 1 at y equal to 0 we have T temperature at x 0 is equal to T star so three walls are insulated and the fourth one we have convective at y equal to 1 we have this boundary condition so we have these four boundary conditions and when we use finite difference method what we got was a set of linear algebraic equations what we had done was we discretize the domain we discretize the domain so this is 1, 2, 3 and so on so we discretize the domain we constructed we demarcated grid points along y axis and along x axis and then we imposed the partial differential equation at the grid points or in the terms of finite difference method we force the residual to 0 okay we force the residual to 0 at the grid points and this forcing residual to 0 at the grid points gave rise to large number of algebraic equations and same thing is true about boundary conditions we use finite difference at the boundary and then we get additional equations at the boundary so total number of equations that we have finally is how many equations you have how many variables you have nx plus 1 cross ny plus 1 equations and you have same number of variables okay now instead of writing those equations again I want to just point out something if you actually try to rearrange this what you get is set of linear algebraic equations and if you rearrange them into a matrix form the standard you know ax equal to b form then you get a very special kind of matrix I do not want to write the matrix I just want to show how it looks graphically okay it is a it is a special diagonal matrix it is a sparse matrix what you get here and see this matrix will have nonzero elements along first three diagonals there will be lot of zeros there will be a nonzero diagonal again there will be a zero here zero diagonals and zero diagonals so you get a sparse matrix here if you rearrange all the equations so this is you know I am t11 or I think you have to eliminate the boundaries so you get t22 to tnx and y if you collect all the if you eliminate the boundary conditions and then collect all the unknowns together and write this as one equation you will get this you will get a vector b here which is coming out of f of x y the values of f of x at different points f of x y at different points and then so this is see this here in this case this is the source term and as I explained to you that the sources if you are modeling this particular room okay then the sources would be each one of you know each one of you is a heat source like a 40 watt bulb and the distribution of these bulbs can be given by f x y you know at which points these bulbs are located and then that would lead to a particular solution a particular temperature distribution so the solution here actually the true solution is a surface with one temperature value attached to each point each xy in this we are not able to solve this equation at every point we are able to solve it only at finite number of points you will get a good solution if you create more number of grid points okay the solution will improve with more and more number of grid points but the computations increase with more and more number of grid points okay and this matrix which you see here which you will get will be a banded matrix and this will be a sparse matrix this is a side note that we will be looking at sparse matrices little later but I just want to point out that I just want to provide a motivation for looking at sparse matrices why do we have to look at special methods for matrices which are filled with large number of zeros because you know you can save computation time when you have large matrices suppose this n x and n y you know this happens to be the other number of variables okay and this matrix will be n x cross n x see the other number of variables so this matrix will be n x cross n y cross n x cross n y this will be a huge matrix just remember total number of variables are these so this if you eliminate boundary conditions you will get something like n x cross n y variables in this vector and this matrix is n x cross n y cross n x cross n y so if you take 100 100 grid points just imagine what is the size of this matrix inverting this kind of a matrix is possible but you know you are wasting lot of computer time because there are lot of zeros and those will probably give zeros when you invert you do not have to you know waste your energy in finding a zeros on the right hand side if you can devise some intelligent algorithm you can save computation time okay now is this the only way of discretizing we discretize both in x and y okay we discretize both in x and y why first of all why do I get this banded structure because when I discretize okay each equation will have only five variables associated with it so I will have this variable this variable this variable this variable so this is i this is this point is i, j this point will be i plus 1, j this will be i minus 1, j and so on so 1, 2, 3, 4 points neighboring points are appearing in each equation and such you have large number of coupled equations you have to solve them together simultaneously to arrive at the solution okay the solution surface is being approximated by solving the equation at discrete number of points since it is an approximation what you get is an approximate solution cannot be the true solution okay for this particular problem in many cases you will be able to solve the problem analytically that is a different story okay but if there is a non-linearity then you are in trouble for example I can make this problem non-linear just by saying that this thermal diffusivity alpha is a function of temperature right now I assume constant temperature you know probably if you are doing a modeling in this room the temperature variation is not too much it is a good assumption that alpha is constant but if there is too much variation of temperature in the region where you are modeling then it is worth modeling alpha as a function of temperature and then immediately the equations which are linear algebraic equations will become non-linear algebraic equations okay leave that aside but let us not talk about non-linearity let us look at this problem again through another method called method of lines okay so what I am going to do now this is the method which belongs to the same finite difference class or it could reappear when you do orthogonal collocations instead of discretizing in both dimensions x and y I could choose to discretize only in one dimension say x okay and I want to retain the differential operator in y direction okay I will discretize only in x direction I want to retain differential operator in y direction this is another way of discretizing so this will give me method of lines so what I am going to do here so what I am going to do here is instead of drawing grid points in both x and y direction I am not going to mark I am not going to mark the I am going to draw parallel lines I am going to mark grid points along along y okay and I am going to draw parallel lines here well the draw the lines in my drawing are not looking parallel but I intend to draw parallel lines okay so I am going to draw parallel lines here okay I want to discretize only in one direction and not in the other direction so I have made a mistake here well I want to draw lines parallel to x axis so I want to draw lines parallel to x axis well it is not really important whether you discretize you keep you discretize only in x direction or whether you discretize in y direction it is not that important in this particular problem because of asymmetric boundary conditions it is good to discretize along x direction so what I am going to do now I am going to mark this variable now this is my notation okay my notation is T i y okay T i y is a variation of temperature along y direction in along the ith line parallel to y axis I have drawn lines parallel to y axis okay and T i y is the variation of temperature okay along ith line parallel to the y axis okay I am just fixing x i so this is this is my x equal to x i this is ith point and T i y is nothing but okay so now the way I want to discretize this equation now is okay now look at this equation the way I have discretized this I have discretized only in x special coordinate x I have chosen not to discretize in y okay now I have converted the partial differential operator along y to total differential because now x has been taken care because of discretization okay so this becomes an ordinary differential equation this becomes an ordinary differential equation now we have to use boundary conditions we have to use boundary conditions right and boundary conditions because of boundary conditions you will get is this one ordinary differential equation how many you get at all the internal gate points remember at all the internal gate points so 2, 3 up to nx and this will be the last point is nx plus 1 so till nx till last parallel line we get these differential equations so i goes from 2, 3 up to nx so you get ordinary differential equations you get ordinary differential equations not one ordinary differential equation you get coupled with a set of ordinary differential equations right ith value ith value or ty ty is related to ti plus 1y and ti minus 1y right so the temperature on the temperature on ith line is a function of temperature on this line and this line is a function of so what you get here is a set of coupled ordinary differential equations not one ordinary differential equation okay you get large number of coupled ordinary differential equations suppose we discretize with 100 grid points internal grid points 100 coupled ordinary differential equations need to be solved simultaneously okay because neighboring variables appear in each one of them and now we can use the boundary condition to complete the problem so we have this two boundary conditions t1y is equal to t star and tnx plus 1y is equal to t star you can use these two conditions to eliminate some variables okay you can use this boundary condition too or you may have to solve it as a differential algebraic system okay so you have these conditions in addition plus you will have to have two conditions which will now become boundary conditions for this set of ordinary differential equations so those will arise because of this and this so if I just write them in well we do not have here okay I will write it here this is one of the boundary conditions okay and the second one will come because of discretization of this okay and that would be Ki1 so discretization of this boundary condition will give you this set of equations okay now what you get here what you get here is a set of coupled second order ordinary differential equations which are these subject to boundary conditions at y equal to 0 and at y equal to 1 okay so same problem finite difference method but instead of choosing to discretize in one in both the directions if I discretize only in one direction I get ODE boundary value problem same original PD I get ODE boundary value problem so it depends upon how you choose to discretize okay this problem could yield set of linear algebraic equations if you choose to discretize in one particular way this problem would yield set of ordinary differential equations boundary value problem if you choose to discretize another way so it depends upon how you choose to discretize the problem okay now I am just going to write one more variation of the same problem we have been looking at the steady state equation what if I decide to look at in addition time variation this is a steady state Laplace equation what about time variation if I include dou t by you know then I will have to then discretization can yield well then there is a problem in using method of lines because if you discretize only in one dimension you will get up another partial differential equation so which you again have to discretize so it depends upon the problem at hand so suppose I have here a modified version of this problem so the modified version I would just create here itself if I add if I add here okay and my right hand side is instead of f of x y I will make it time dependent I will make it time dependent okay and all these conditions at boundaries will hold for all time okay all these conditions at the boundaries will hold at all time and then I will also have initial conditions so I will have to give a initial distribution of temperature in the room so what I want to know is how the temperature surface is changing as a function of time okay when I am including time derivative here okay I am looking at the spatial distribution of temperature in this room or in a furnace or whatever whatever the condition is and then I am interested in time evolution of this temperatures okay not just the steady state okay so this problem will also have additional initial conditions coming up so the boundary conditions have to hold for all t so this has to hold for all time t including this boundary condition and then you will have t is equal to 0 you will have temperature at x y x y 0 is equal to some function say h x y so this is the initial condition so this is the initial distribution of temperature in the room okay and then I want to solve this problem what should I do here can you suggest something what if I decide to discretize into spatial dimensions I discretize in x and y I leave temperature as it is what will I get I will get ordinary differential equation initial value problem I will get ordinary differential equations initial value problem large number of ordinary differential equations which are coupled okay so here if I you know if I now take I am not going to completely write the solution you can work it out I am just going to hint at the solution so I am going to call now t i j x i or I am going to call t i j t is equal to temperature at x i y j and t so like in the previous case I have created a grid I have created grid here and so on I have created a grid and I am at i j th point the time variation I am going to write like this and then well I will discretize this problem as d t i j t by d t plus alpha times now I will approximate this as t i plus 1 comma j minus 2 t i j plus t i minus 1 comma j divided by delta x square plus so if I discretize into dimensions I get large number of equations okay large number of ordinary differential equations in time but in this case in this case I will not get a boundary value problem just think about it in this case I will get a ordinary differential equations initial value problem okay so like now we are not talking about how to solve these problems we are just talking about transformations we want to take a problem transform it into a solvable or computable form how to actually compute that solution will come to that later okay I want to separate these two things I just want to give you a view point that well if I take a partial differential equations or ordinary differential equations I know from we are Starr's theorem that any continuous function can be approximated using polynomials so I use Taylor series approximation construct local approximations of the differential operator and force the so-called residuals to zero that gives rise to an approximate problem okay approximate problem would look completely different from the original problem a PD is giving rise to set of algebraic equations linear algebraic equations PD giving rise to set of ordinary differential equations boundary value problem or a PD giving rise to set of ordinary differential equations initial value problem. So we bring all these into some standard forms and then we apply the standard tools to solve these problems okay that is what you should realize there is no unique way of coming up with the discretization you could have your preferences solve the problem at hand so in this place you can complete the initial conditions you can replace all you can remove the so this could become a DAE differential algebraic system with initial value problem initial condition specified or if you eliminate it could be just ordinary differential equations initial value problem okay so till now we were looking at Taylor series approximation okay and as I said this is only one of the methods that are used for discretization other two methods are interpolation polynomial interpolation and then least square approximation so now we will stop here with Taylor series approximation we will now move on to polynomial interpolation and let us see how interpolation can be used to discretize the same set of problems what I have deliberately done is in the lecture notes is that I have taken same set of problems and discretize them using different approaches okay so that you will get a better insight into what is really happening at the discretization stage okay so now let us move on to interpolation interpolation is something which probably you are introduced in your undergraduate curriculum can anyone tell me what is interpolation what is an interpolating polynomial what is an interpolating polynomial when do you do interpolation or interpolating functions you are given a set of points and values so you are given some set of so you actually want to develop an approximation for a function continuous function let us take a continuous function say u or to be very precise let me see to begin with polynomial interpolation well you can also do functional interpolation in general functional interpolation but we are interested right now in polynomial interpolation because we will use the polynomial interpolation ideas to come up with approximations or discretization of some operators okay so now polynomial interpolation is we are given some values of say dependent variable uz and independent variable z so we know some values of u okay just before I move on to this let me clarify one point in finite difference method I kept on using equi-space grid points okay it is not necessary that you should have equi-space grid points that was only for the convenience of development understanding actually you can use non equi-space grid points okay now for example if you are now when to use non equi-space grid points and then how to space them will depend upon your understanding of the physical problem okay see for example if you have a pfr okay if you have a fluxor reactor and let us say most of the action is in the initial part of the pfr and later you know the concentration defiles becomes steady or it does not change too much it becomes flagged okay it is worth putting more grid points in the initial part which are at closer which are closer to each other and then putting sparse grid points later because there is not much thing happening at the end of the reactor okay so this depends upon your understanding of the problem in some cases you may want to place close grid points in the at the two ends and you know sparse grid points in the middle it depends upon the problem okay so there is nothing sacrosanct about equi-space grid points equi-space grid points only make development on the board simple but in general when you and as a beginner you might start with equi-space grid point but then if you know about some problem where there is more variation of a particular variation in particular region you can put closer grid points and you can have sparse grid points at that depends upon your understanding of the system okay so let us move on to this polynomial interpolation now what I want to do here is I have this function I have this function uz okay as a function of independent variable z and I know the function values at these discrete points okay at these discrete points in the domain so I know so this is my point number one this is point number two this is in general ith point and this is my n plus one point okay from z1 to zn I know value of a function okay now what I want to do is I want to construct a polynomial which passes through each one of them okay I know this is a continuous function of independent variable z okay by we are star's theorem I can approximate it as a polynomial function okay we are star's theorem allows me to approximate it as a polynomial function I want to construct a polynomial function that passes through all the points okay this is different from least square approximations which you do in your experimental work you know where you fit a line so the line may not pass through every point that is least square approximation we will we will look at least square approximations little later but right now we are concerned about constructing a polynomial that goes through every point okay that goes through every point so my problem is to construct a polynomial pz I am going to construct a polynomial approximation pz now this is alpha 0 alpha 1 if I have n plus 1 points I can fit a polynomial of order n that exactly passes through all these points I have to I need a polynomial of order n because the polynomial must pass through every point so this is the constraint okay so what does this constraint mean so I have these points z z1 z2 z n plus 1 I have these points okay what I want is that value of this polynomial at each of these points should be exactly equal to ui that is an interpolating polynomial okay so I have this set of constraints so alpha 0 plus alpha 1 z1 up to alpha n z1 raise to n is equal to u1 alpha 0 plus alpha 1 z2 and so on so I get n plus 1 equations I get these n plus 1 equations in n plus 1 unknowns what are the unknowns here alpha 0 alpha 1 alpha 2 up to alpha n there are n plus 1 unknowns okay and there are n plus 1 equations how do you solve this so this can be simply transformed into 1 z1 z1 square z1 raise to n 1 z2 z2 square z2 raise to n so this is simply solving ax equal to b okay well unfortunately life is not so simple because it turns out that if you take if you want to develop a higher order interpolation polynomial nth order interpolation polynomial in general above 3 or 4 this matrix here it looks very simple right ax equal to b this is known to you these values are known to you function values are known to you right now okay what you do not know is alpha 0 to alpha n and you know this matrix right you know the point at which so you know z1 you know z2 so you can create this matrix we will see in the course of looking at ill conditioned matrices we will see that this is one of the highly ill conditioned matrix okay highly ill conditioned matrix difficult to invert you can get so higher order interpolation polynomials if you just go by this brute force approach can lead you to problems well one way to get out of the situation is to use orthogonal polynomials this polynomials are not orthogonal this polynomials are not orthogonal so that is why you get into trouble if you use orthogonal polynomials you can come out of this problem and then this matrix is very nicely behaved and so on but right now let us not worry about non the that part right now let us assume that this matrix is invertible so how do I get alpha 0 to alpha n so if this is my a matrix okay this is my a matrix and let us say this is my capital U vector then and if I call this simply as theta parameter vector then you know at least one paper I can write theta is equal to a inverse U right a matrix is known to me U vector I know the function values at these points theta I can get what I get is a interpolating polynomial this is called interpolating polynomial well I am going to use the idea of interpolating polynomial to discretize boundary value problems to discretize partial differential equations what I will be getting there or what I will the method that is used to do that is called as method of orthogonal collocations why in this word orthogonal comes into picture we will look at those details but the basis is this that if you know a value at certain number of points then theta can be expressed as a inverse U okay and this polynomial which you get here if you are able to calculate all these alpha 0 to alpha n correctly this polynomial will pass exactly through all these points that is very important okay this polynomial will pass through each one of these points it is a nth order polynomial that passes through each one of these points now the trouble I told you is that constructing a very high order polynomial through large number of points is can give you to some ill conditioning and when you are discretizing see we saw in finite difference method what was the what was the message take home message that you know if you have more number of grid points better the solution okay so when you are discretizing this is true for even for if you use interpolation approach to discretize so you would need large number of grid points okay but that would mean you have to fit a large order polynomial and then you can get into ill conditioning problem now when I am going to do the development of orthogonal collocation method I am not going to bother too much about this ill conditioning I am just going to develop a method to and convey the approach but what you really need is spline interpolation do you know what is spline interpolation cubic spline interpolation well it's it's I have discussed it in detail in the lecture notes I will just give you a brief idea what is cubic spline interpolation the exact equations you can see here in section 4.2 but let me explain the idea the basic idea is this that instead of fitting one giant polynomial of high order which passes through all the equation we could choose to we could choose to construct piecewise polynomial approximation okay so I could construct one approximation between these two points I could construct another approximation between these two points third approximation between these two points what I have to make sure is that these neighboring approximations are continuous okay neighboring approximations are continuous typically what is done is that we fit a cubic polynomial between two neighboring points okay and make sure that they are smoothly joined okay so how do you make sure that they are smoothly joined well we make sure that the first derivative and the second derivative is match for the neighboring polynomials okay so this is called as piecewise polynomial approximation and this is many times required because so what is what is the idea here is that you still have the same number of points okay you still have u i that is u z i for z i or i going from 1 to up to n plus 1 okay you still have these grid points okay now what I am going to do here in piecewise polynomial approximation well let me be little careful about the terminology piecewise polynomial interpolation approximation will be using in some other context okay so piecewise polynomial interpolation so what I do here is that I fit a polynomial between two neighboring points and then I make sure that there is a continuity between two neighboring polynomials okay so I will have a polynomial p1z which is alpha 01 plus alpha 1 comma 1 now there are many polynomials so I have to have you know a notation which is you know which has two coefficients the coefficients will have alpha 01 so this is the first polynomial 0 coefficient okay then z minus z1 plus alpha 21 z minus z1 square plus alpha 31 z minus z1 cube and this holds between z1 z z2 p2z is alpha 02 plus alpha 12 z minus z2 plus alpha 22 z minus z2 square plus alpha 32 z minus z2 cube so this polynomial holds between z2 z and z3 and likewise you know between each segment I fit a polynomial in this approach this is called piecewise polynomial approximation okay this is called piecewise polynomial approximation now how do I make sure there is continuity I have to put some conditions here so one condition is of course one condition is of course you know these values one condition come from these values okay the other conditions come from the derivatives or the smoothness of the derivatives so other conditions arise so this is the terminal point so these are all initial points of each polynomial then we have these conditions which are continuity conditions p i z i plus 1 is equal to p i plus 1 z i plus 1 so essentially we are saying that the initial value of this polynomial is same as the final value of the previous polynomial okay then you have condition dpi z i plus 1 by dz is equal to dpi i plus 1 and the third condition is so to ensure smoothness we make sure that the value of these two neighboring polynomials at the common point is same the first derivative of the two neighboring polynomials at the common point is same and the second derivative of the two neighboring polynomials at the common point is identical so these are the additional constraints which we impose okay these are additional constraints which we impose if you actually put all these constraints then you need two boundary constraints so there are different ways of specifying the boundary constraints one of them is called as you just look at the nodes here so free boundary conditions are we need two more conditions at the two boundary points at this point and this point so those are free boundary conditions would be d2p if you include these two boundary conditions that at the two boundaries the second derivative is 0 okay if you include these two boundary conditions then you get number of equations equal to number of unknowns what are the unknowns here alpha values okay how many alpha values 4 into n okay you can do some algebraic manipulations and finally reduce this problem to ax equal to b where a is a sparse matrix the manipulations are given here you can just have a look at it but interpolation when it is not possible to do high order interpolation this cubic polynomial interpolation will be used okay and that will that will lead to orthogonal collocations on finite element so we will get a method finite element method if you do you know polynomial interpolations on a smaller domain so the different ways of doing interpolation and it is a very very rich area I am just touching the tip of the iceberg by talking about these two basic ideas which I need to develop my further approach so the next lecture I will start using interpolation idea not the cubic polynomial not the piecewise polynomial interpolation simple interpolation okay nth order polynomial passing through all n plus 1 points this simple idea I am going to use to develop method of orthogonal collocations okay so by the method of orthogonal collocations we will discretize again boundary value problems and partial differential equations okay so you will again revisit the same problems but through a different approach okay so orthogonal collocations is a very popular method in fact it is a method which requires less number of grid points to get the same accuracy as that of the finite element finite difference method finite difference method would require large number of grid points this method gives better solutions with less number of grid points so this is computationally less expensive as compared to finite difference method okay that is why we want to go for polynomial interpolation