 So, till now we have been looking at problem discretization we transform the problem from original problem which was let us say a boundary value problem or a partial differential equation into a computable form that computable form could be a set of linear algebraic equations linear different or non-linear algebraic equations the transform form could be O the initial value problem and now we need ways of actually solving constructing the solution. So I just want to draw the picture that we started with when so we had this original problem this original problem was non-linear algebraic equations boundary value problem is partial differential equations all kinds of things that arise in modeling of engineering systems they could be differential algebraic equations. So we have different kinds of equations that we need to solve so this is coming from your physics this is coming from your modeling transport phenomena heat and mass transfer and so on. So this will give you set of equations most of the times these equations are non-linear not at all amenable to constructing analytical solution these are in multiple dimensions and you have to solve them using some numerical techniques okay. So the next step was problem approximation so we constructed we used concept of approximation theory we got this transform problem this is transform computable problem. So in approximation theory we use three different tools one was Taylor series approximation basically we use polynomial approximations that was one of the fundamental tools but of course we use also function approximations later when we talk about least squares right. So we had we used concepts of approximation theory to come up with a transform problem. So here either interpolation Taylor series approximation or these squares these were two these were three basic ideas that were used to transform the problem. So I would just list them here is either Taylor series or interpolation. So Taylor series interpolation or these square approximations so these were three predominantly used tools to transform the problem to a computable form and now we want to attack this problem so we are and as you could see that this problem which is started with was a partial differential equation you ended up here sometimes with non-linear algebraic equations sometimes with initial value problem ordinary differential equation. So you the transform problem from the view point of the structure of equations could be completely different okay so it did not be it did not resemble the original problem original problem is a partial differential equation in space and time here you get ordinary differential equation only in time original problem is partial differential equation in space you just get algebraic equations either linear or non-linear so transform problem is completely different. We hope that if we solve the transform problem we get something close to the true solution not the exact solution. Now how do you attack the transform problem that is so the next thing is look at the tools so I am going to talk about three different tools ax equal to b I would put this as solving linear algebraic equations this is one of the fundamental tools that we are going to use so this tool will be used to attack this problem well the other one other tool that I am interested is actually solving f of x equal to 0 so this is another tool which is used to solve the transform problem okay the third tool which I am going to be looking at is od initial value problem initial value problems is again the fundamental block and Euler method all those methods will come under this here we will revisit Newton's method will look at its other nuances and see how you can enhance the convergence and well the fourth tool which is quite commonly used is stochastic methods but I am not going to deal with this okay stochastic methods this is a fourth tool what we expect after that once we attack this problem with these tools or their combinations it is not always that you will use only one tool you might use this and this because Newton Raphson will require this okay to solve this and so on so their combinations and finally we what comes out is the approximate solution so this is my end result this is my end result so I start with the original problem I use the approximation theory come up with a transform problem which is computable I have only three four tools with me which I use ingeniously to concoct a solution you should get this feel at the end of the course that actually it is like a belpuri shop you know you only have few things with you okay and with the same things you can make sev puri you can make they better the puri you can make you know belpuri so it is the same thing is the same ingredients okay and by just making mixtures you can create different dishes the same idea is here that we do not have too many tools we do not understand how to solve very very complex equations exactly we just know after all that we have done progress in mathematics we just know how to solve ax equal to b we just know how to solve non-linear algebraic equations and all these you will realize we know how to solve this approximately we know how to solve this approximately so ultimately the you know what you land up with is a third order approximation you approximate from here to here okay then you try to use these tools but they themselves are approximate okay and then there are errors in computation because of inherent limitations of a computer so all three combined together okay you get a third order approximation which is which you hope is which you hope is the close to the reality what is what the real solution is okay well there is one more approximation which I forgot if you start with a real system okay so this is the real process or a system when you write a mathematical model okay that itself is an approximation well when I say that reactor is like a CSTR it is not CSTR or when I say it is a PFR it is not the perfect PFR PFR and CSTR these are models these are idealized models which try to approximate the reality so from here to here itself there is an approximation what we formulate here we cannot solve exactly so there is one more level of approximation then when you try to attack them using different tools those of those tools are also approximate okay because all these ideas of Taylor series then polynomial approximation interpolation all these again we will revisit when we are doing this okay so they are not going to leave us because again many of these things are not solvable so you go back to again Taylor series you go back to again and then you solve an approximation of approximation and so on so you could say it is a fourth order approximation and in between you know you have computer which is having its own limitations so all things put together what you get is the approximate solution now this approximate solution should represent something that is happening in this problem and in the reality and this is where your inputs as engineer or a physicist or a chemistry will come into picture you should make you should know whether this makes sense often times you need to actually give good guess you must have realize this when you are solving some of the problems numerically which and if you do not give a good guess the solution can be absurd and so giving good guess is where your background in engineering or science is very very critical okay so now what we do is we have come up to this point now we look at a b and c we will not be able to get into stochastic methods hopefully we will be able to cover give up I will be able to do some justice to these three tools okay so let us begin with a x equal to b and then you might wonder well I have been doing this since my 12 standard so what more to it what more to a x equal to b I have been doing this for ages a lot more to it I will at least need 8 or 10 lectures that too we will just touch that tip of the iceberg when using or when you are solving problems in school takes books with some two equations in two unknowns or three equations in three unknowns you know simple methods work well the most the first thing that you learn is Kramer's rule right now what I show is that Kramer's rule beyond four equations becomes unwieldy you cannot use it for computer implementation the most practical method is of course Gaussian elimination again credited to Prince of mathematics Gauss this method is probably one of the most efficient methods of directly solving linear algebraic equations but then again these methods will be good if you have thousand equations in thousand unknowns once you start going into ten thousand twenty thousand one one like equations in one like unknowns you have problems okay you still get into roadblocks of time constraints you can use iterative methods and we will talk about iterative methods so iterative methods are you start with a guess solution and then try to converge to the close to the true solution as quickly as possible again iterative methods we have to look at when you are guessing and trying to go to a solution you should know whether you will converge or not okay so under what conditions you will converge is going to be one of the main themes there are also optimization based methods which are again iterative methods I will briefly touch up on optimization methods and then finally I will also talk about a fundamental thing what is ill conditioned systems and well conditioned system so how do you classify ax equal to b well it depends upon this a matrix so we will we will talk about properties of a matrix like a condition number in detail I suppose you may have been introduced to this idea of condition number but I will derive the the basic ideas that lead to a condition number and we will talk about well conditioned matrices ill conditioned matrices when the matrices are ill conditioned you cannot get reliable solutions through numerical computations and you should be aware of that because nowadays you have tools you know you just give a matrix it will pop out a solution you should know whether this solution makes sense or not okay one is from physics viewpoint other is from computation viewpoint if the matrix is ill conditioned okay you cannot get a good solution and you should know should be aware of that okay so let us begin by talking about existence of solutions when does the solutions exist b lies in the column space of a so there are two actually pictures there are two geometric pictures that are typically used I will just use this you should actually look at the book by Gilbert Strang he gives very very nice introduction to solution of linear algebraic equations this is example from Gilbert Strang I am just taking one simple example which is 2-1 1 1 x y is equal to 1 and 5 okay now there are two viewpoints by which you can geometrically visualize the solution well one viewpoint which is taught in our schools is you know two lines intersecting or when it goes to three variables in three unknowns we talk about three planes intersecting at one point and so on in general we can say that hyper planes intersecting so one viewpoint I would say is this that is 2-1 x y is equal to 1 this is equation number 1 and second equation is 1 1 x y is equal to 5 okay so this is we would say that well so we draw these two lines and you know we this is this is my x y plane this is very commonly given example so this line is 2 x minus y is equal to 1 and the other line is x plus y is equal to 5 and then this line and this line intersects at this point this is the solution okay commonly one example means not this particular example but intersecting lines is or intersecting planes so this is a solution that is what we know the other interpretation is of course the column space interpretation which is coming from vector addition so I could view this equation as 2-1 sorry 2-1 times x plus y times y times minus 1 1 is equal to 1 5 okay this is the second equation and here we are not drawing pictures in x y plane any longer x and y are coefficients by which these two vectors add to give you the third vector so this is your vector addition okay and this picture is I will call this as so this is row picture row picture becomes very very difficult to visualize in some sense beyond three dimensions okay but my claim is that the column picture is it will be easier to visualize after three dimensions four dimensions five dimensions we are just talking about linear combination of vectors so here what we are saying is we have these two vectors we have this vector 1 minus 1 and I have this another vector this is my vector 2 1 this is my vector 1 minus 1 okay and then my third vector is this 1 5 we are just trying to use the law of vector additions okay and then you can complete the geometric picture by drawing the parallelogram here so actually it's just linear combination of these two vectors gives me the third vector okay gives me the third vector and when will you when will you be able to get solution for any right hand side when if you want to specify any right hand side okay then the linear combination of the columns put should span the entire space that means that means these two these two should be linearly independent if these are linearly independent all possible linear combinations of these two vectors will span R2 and then any vector well there will be a problem if suppose if this was if this was 2 minus 2 suppose I change this equation to 2 minus 2 what is the trouble these two are these two are vectors in the same direction okay no linear combination of these two vectors can ever give me 1 5 okay so the column picture would say that equation doesn't have a solution that is because you have one equation you have two vectors which are aligned along the same direction there is no linear combination which is ever going to give me this vector okay no linear combination is ever going to give me this vector so now when column picture tells you that there is no solution same thing will happen for row picture if you take row picture here okay for this changed case you will see that there are no two different lines see they are parallel lines which will never meet so 2 x minus y is equal to 1 and minus 2 x plus y is equal to 5 so this will be the row picture and you will get two lines which never intersect okay so when this picture shows that there is a problem you cannot get a solution the row picture also will give you a problem there is no solution no hyper planes meet or no planes or no lines meet okay we can generalize this to any n dimensional case we can talk about linear combinations of column giving you a solution which lies in the column space okay and row picture will be hyper planes and then the solution when it exists is all the hyper planes meet in one point okay so actually line is a one dimensional hyper plane in two dimensional space plane is a two dimensional hyper plane in three dimensional space and in four dimensional space there will be a three dimensional hyper plane and so on and those three dimensional hyper planes should meet well difficult to visualize how three dimensional planes meet at one point but easier to look at the column picture and talk about linear combinations of columns in n dimensional spaces okay so we will not get beyond a certain point here what is important is because all of you know this I am just revising this very quickly when does the solution exist solution exists when the columns are linearly dependent independent okay and same thing is true about rows if the columns are linearly independent the rows are really independent and number of the fundamental theorem of linear algebra is number of linearly independent rows equal to number of linearly independent columns is equal to rank of the matrix okay so with this brief brush up we will just get into the solving problems using okay now Gaussian elimination is something which I am not going to touch upon because I assume that you already have a sufficient background of Gaussian elimination I have uploaded some nodes on Gaussian elimination you should have a look at those nodes I just want to compare two things here to begin with to give a motivation for methods which are more efficient one is number of computations that are number of multiplications and divisions that are involved so two methods that you know for solving a x equal to b the first method you learn is Kramer's rule okay the first method that you learn is Kramer's rule now I am going to use this psi to denote number of required to arrive at the solution if psi represents the number of multiplications and divisions that are required to arrive at the solution then for Kramer's rule psi estimated is approximately equal to I am talking about n cross n matrix so my concern here is a is in general an n cross n matrix where n could be large okay n could be thousand ten thousand or whatever a large number and we have seen that these kind of matrices arise when you do let's say finite difference method problem discretization of boundary value problem a two-dimensional boundary value problem will give rise to even if you take hundred hundred grid points okay along x and y you will get a huge matrix which has to be solved so large matrices are not uncommon when you do approximate solution so here it comes out to be n minus 1 into n plus 1 into n factorial plus n which is approximately or this is equal to this is exactly equal to which is approximately equal to n square into 1 this is n factorial plus n this is n square into n factorial okay estimated number of operations divisions and multiplications for Kramer's rule I am not going to derive this you just accept this number right now I am just going to give you estimate of what can happen if you try to use Kramer's rule for a matrix of moderate size hundred by hundred okay if you take n equal to hundred okay this data I have taken from process book so this is you can show that the size is close to 10 to the power 162 okay and he says that a deck computer deck 1090 would take 10 to the power 149 years to solve this problem okay so absurd number so to solve a problem of hundred cross hundred matrix involving hundred cross hundred matrix using Kramer's rule is impossible okay a computer would take the the deck computer when he probably wrote the book was the fastest one and then he has given this number that it would take 10 to the power of 149 years okay so Kramer's rule is ruled out it just good for you know 10 standard or 12 standard whatever you learn it and good for solving three equations in three unknowns maybe four equations if you have patience to do determinants then so this is not a way to go certainly not the way to go what about Gaussian elimination Gaussian elimination seems to be the hope here so in Gaussian elimination we start with ax equal to b then we transform it to an upper triangular matrix ux is equal to b prime or b cap let's call it we transform by elimination we transform it to an upper triangular matrix and then we do the backwards sweep okay so what is the number of operations that you need here the number of operations that you need here is n cube plus 3 n square minus n by 3 which is approximately for large n which is approximately equal to n cube by 3 so n cube by 3 is far far far less than you know factorial so this the reason why Gaussian elimination is so popular and is computer savvy method is because this is one of the most efficient ways by which you can solve ax equal to b okay so that's why it is introduced right at 10 standard or 12 standard and then you start building on it as you go into your higher grades okay so now so obviously we want to we want to work with Gaussian elimination okay and now what I want to do is before I move to iterative methods I want to see whether Gaussian elimination can be made more efficient okay by exploiting certain structures that appear in discretization problem discretization so at some time back when we did this discretization I told you that finite difference or orthogonal collocations you know or orthogonal collocations on finite element and so on we get some special matrices these are called as sparse matrices sparse matrices are ones which are filled with large number of zeros and they may have some nice structures okay if you exploit these structures okay you can come up with direct methods or so these these methods which I am talking about these are classified as direct methods okay when I say direct methods as against iterative methods iterative methods in which you start with a guess and from the old guess you construct a new guess and so on okay like Newton Raphson the iterative method direct methods in Gaussian elimination you do not start with a guess right you directly solve the problem okay other variants of Gaussian elimination are Gauss-Jordan method okay LU decomposition so all these are variants of Gaussian elimination but basically the foundation is Gaussian elimination okay also I have uploaded something about LU decomposition LU decomposition is Gaussian elimination represented as you know sequence of two triangular matrices my triangular matrix calculations one is upper triangular other is a lower triangular so you can have LU decomposition also you can you can have a look at LU decomposition but that is derived from Gaussian elimination so the order of computations are again of the order of n cube by 3 it is not not too different okay so in direct methods now what we are going to look at is direct methods for solving when I say direct method not the Kramer's rule of course the Gaussian elimination sparse linear systems are those set of equations ax equal to b where matrix a has large number of zeros and I want to exploit spatial structure of that matrix a to come up with a efficient solution which will significantly reduce the number of multiplications and divisions that you need for solving the problem okay that is the motivation for these sparse methods again in this course I cannot do a justice to sparse methods I am just going to introduce you to sparse methods okay everything is a iceberg and we just touch the tip of the iceberg if you look at sparse method in Matlab there is a sparse matrix toolbox and it will list so many you know algorithms I cannot cover all of them but what is important is to sensitize you that there exists something called sparse methods and you should look at you should try to see whether there is a nice structure in your problem which can be exploited to reduce the computations that's what is more important so I am going to touch upon three or four methods and we will move on to something else but at least this will give you a flavor of what it is what is this sparse matrix business okay so we are going to look at some things which are some matrices some sparse matrices this is one is tridiagonal and block tridiagonal matrices we will look at block diagonal matrices we will look at well in general one can talk about banded matrices so banded matrices with only m diagonals okay I am going to touch upon this I will talk about these two then I will also mention about triangular matrices block triangular matrices so these are the some these are some nice forms actually if you start looking at these forms just in a linear algebra textbook you might wonder where do I get all these forms but all these forms appear in problem discretization okay and once we write some of these forms will realize where they appear and how you so the motivation for these looking at these sparse matrices becomes clear if you look at the problem discretization methods okay those will give you these nicely patterned matrices which you can actually exploit to come up with efficient solutions let's look at first the the simplest one I would say is block diagonal matrix a block diagonal matrix would appear a block diagonal matrix would appear in orthogonal collocations on finite elements OCFE a block diagonal matrix would look something like this this original equation which is ax equal to b I am writing it in a block diagonal matrix form okay here a1 a2 up to am are matrices which are full they are not they are they don't have too many zeros these matrices are typically dense matrices not sparse matrices but all these square zeros are actually matrices of appropriate dimension which contain only zeros so this is a sparse matrix because nonzero nonzero elements appear along these diagonals appear along these diagonals and then here this is a dense matrix this is a dense matrix this is a dense matrix but all these are zero matrices null matrices so this is a huge matrix filled with large number of zeros only here on the diagonal you have some small sub matrices which are dense okay the small sub matrices which are dense what I have done here x x here I have partitioned here okay I'll talk about this partitioning very soon and even b vector I am partitioning into b1 b2 okay now what are these partitions what are these partitions well my matrix ai is actually a matrix which is ni cross ni okay so this could be see this could be a 3 cross 3 matrix this could be a 5 cross 5 matrix next one could be a 3 cross 3 matrix if you are doing orthogonal collocations on finite elements on the first element you take three internal collocation points next one you take five internal collocation points and so on you will get this matrices of different dimensions only these matrices on the diagonal are dense because if you remember when you do orthogonal collocations on finite elements only variables in that small segment appear appear in that equation okay so you have this funny structure that you will get well will you get this or you will get some overlap probably in OCFE you will not get exactly this structure you might get some overlapping and you will have to do some some more tricks to bring it to this form okay but there are some applications where you get this nice banded structure and then each one of them is a ni cross ni matrix this vector xi x i belongs to R ni so it's a ni dimensional vector okay and same is true about vector bi which belongs to R ni so these are subsystems these are subsystems which are arranged in a diagonal form okay now what do I do when I want to solve this well one way is the brute forces take this entire matrix and use Gaussian elimination okay I can use Gaussian elimination on this entire matrix okay so what is n here total n my total n is summation i goes from 1 to m ni right my total n all number of how many number of what is the dimension on this entire matrix n cross n what is n see these are all subsystems I am just adding them so the entire matrix this entire matrix is n cross n where n is equal to summation of all ni okay these are ni subsystems which I want to solve together so if I use Gaussian elimination here directly then the number of multiplications and divisions if I use one without without trying to exploit the structure of this matrix if I try to use blanket Gaussian elimination then number of divisions and multiplications we just listed this will be n cube plus 3 n square minus n divided by 3 the other option is I could solve subsystems so I could solve for ai xi is equal to bi i going from 1 to m I have this m subsystems okay I could solve each small sub problem okay see this this is 3 by 3 or 4 by 4 you know each one of them is a smaller smaller dimensional system I can solve this each one of them using Gaussian elimination okay now what is the what is psi i what is psi i here that will be ni cube plus 3 ni square minus ni by 3 right so this ni r 3 4 5 you know the smaller dimensional matrices whereas this one this this matrix is a large matrix this matrix is a large matrix which has large number of rows and columns so instead of using a Gaussian elimination on entire one if I do this trick of solving each subsystem separately then the total number of total number of this thing is sigma ni cube plus 3 ni square minus ni by 3 is actually far far less than this will be far far less than this because see what is n n is summation ni summation ni cube will have lot many terms than okay so this simple trick of identifying sub blocks which are dense and solving those sub blocks okay can actually reduce your computations multiple folds okay so within Gaussian elimination also if you are instead of instead of solving one giant problem in which there are lots of zeros and you are going to waste time in eliminating zeros because they are already zeros in what do you do in Gaussian elimination you first make it low upper triangular matrix so you have to bring zeros and then if you just apply Gaussian elimination without thinking you will waste time in bringing zeros from zeros doesn't make any sense and but the computer will just do computations without you know understanding that there is a zero there so you have to tell the computer well there are zeros don't worry about those zeros just do the subsystem calculations and then you can solve the problem much more efficiently than okay so this is very very important particularly when you have iterative calculations suppose you have to solve ax equal to b kind of equations inside a Newton-Raphson where a has this banded structure then you know just imagine see when you are doing this repeatedly say thousand times or ten thousand times well you better solve it efficiently okay so suppose this ax of ax equal to b this banded structure appears inside a loop okay iterative loop okay each time if you do these many calculations it is much more expensive than doing these many calculations okay that's where this banded matrix structure helps well the next one we are going to look at is the Thomas algorithm okay so the Thomas algorithm is for a special case of Thomas algorithm maybe some of you have done this in your undergraduate because it is often taught so Thomas algorithm is for block is first for a tridiagonal matrices okay I am going to write this notation I am going to change a little bit here the right hand side this elements I am going to call them as d1 d2 dn my matrix is now going to be written with b1 c1 0 then a2 b2 c2 0 then 0 a3 b3 c3 0 so this is a tridiagonal matrix you have zeros here and then finally you will get an bn so you have this b1 b2 b3 up to bn that is the main diagonal okay there is one more diagonal about this which is c1 c2 c3 and one more diagonal below this a2 a3 so this starts from a2 this starts from 1 but this ends here at cn-1 and this starts from a2 ends in n you have seen this kind of a matrix finite difference method you have seen this kind of a matrix okay so in finite difference method you get this tridiagonal matrix okay and then question is question is can I exploit the spatial structure this has so many zeros filled with so many zeros how many elements are non-zero 3n-2 3n-2 n elements of diagonal n-1 of first diagonal above and n-1 of first diagonal below so 3n-2 elements are non-zero rest all elements are zero okay so I am going to do simple Gaussian elimination here can you do it can you try this obviously do Gaussian elimination you should just eliminate I don't have to do Gaussian elimination for these elements here I just have to eliminate a2 okay when I come here I just have to eliminate a3 okay and so on so I just have to eliminate one element below the diagonal main diagonal okay that will give me a upper triangular matrix and then I just go on doing in the backwards sweep so if you just write the steps here it will be something like this so I am going to define this gamma1 which is c1 by b1 I am just writing it algorithmically what I am doing is basically Gaussian elimination okay and gamma k is equal to ck bk minus ak gamma k minus 1 this is for k equal to 2 3 up to n minus 1 then I have this beta1 is equal to d1 by b1 and beta k is equal to dk by dk minus ak beta k minus 1 bk so what I have written here it looks probably at the first site little complex but all that I am doing is eliminating eliminating one element below the diagonal okay and these gamma and beta are the elements which appear in the so what I am doing is after doing this operations I am going to get this new matrix which looks like this new matrix looks this new matrix after doing all this will look something like this 1 gamma 1 0 0 0 1 gamma 2 0 of course you have this matrix x here vector x here and the right hand side is right hand side is beta1 beta2 betan these equations written here are for computer implementation you can put it in a for loop and just do calculations to get this matrix this is this is just these are the steps of Gaussian elimination written in an algorithmic way so that you can put it in a computer program too so now I get this matrix this is upper triangular matrix when you do upper triangular matrix here you just have to worry about one row above the diagonal right and now what you do backwards sweep okay you do backwards sweep you get x1 directly from this x2 you get right all of you know back substitution so now I do backwards sweep and then solve the problem what is the number of multiplications and divisions that are required in Thomas algorithm this is called Thomas algorithm now I will do backwards sweep on this should I write those equations or you are clear about it I think all of you know about this right how do you solve from this very very simple okay what is the point to take home is that number of multiplications and divisions for this case reduced to 5n-8 this is just linear function of n as against cubic function of n if I were to use if I were to use Gaussian elimination without thinking about the structure I would be doing n cube by three calculations roughly here I just do 5n-8 significantly lower than what I need to do as you know by normal Gaussian elimination so Thomas algorithm is one example where you significantly reduce the computations required by exploiting the structure the next what we are going to see is of course a simple extension of this is if you have this I will do in my next class is what if see here we are looked at this algorithm where a bc are scalars what if these are matrices what if b is a mate b1 is a matrix c1 is a matrix a2 is a matrix b2 is a matrix those things are those matrices are called as block tridiagonal matrices okay block tridiagonal matrices so we just simply extend this idea to block tridiagonal matrices and still we can do much less computations than doing entire solving the entire thing okay so that is what is going to be the next extension so these are some of the well-known sparse matrix algorithms and you can significantly reduce the computations even in Gaussian elimination we are still doing Gaussian elimination we are not going to do something new okay so we will continue on more of this in my in the next lecture