Mod - 01 lec - 05 foundation of scientific computing - 05

In today's lecture, we begin by describing the trapezoidal method for step differential equation. This has been used for both linear as well as non-linear equation, while this is straightforward for linear equation. We will see that for non-linear equation, we would require to have the solution strategy in an iterative fashion and that makes it difficult, and we will also find out that sometimes it would even not work in that context. In contrast to trapezoidal method, we will basically talk about a new method which is not very old; it has come about over the last couple of decades. It is called the compound matrix method, and in this compound matrix method, we will see its versatility by talking about instability problems, and in the context of flow instability, we will be introducing what is known as the Orr-Sommerfeld equation. This is basically a fourth order ordinary differential equation which has been known for its stiffness, and this viscous instability problems would require a special handling which was a very tough challenge and yielded solutions very efficiently by compound matrix method. We will talk about the Eigen values and Eigen function evaluation by compound matrix method. We were talking about a special class of ordinary differential equations or we can call it as stiff differential equations. Stiffness arises because of the feature of the differential equation which admits fundamental solutions, which have disparate length or time scales associated the variations, and for example, say, if we are trying to solve this vector equation, the unknowns being from u 1 to u m, and you have prescribed, let say some kind of an initial condition u 0, and what you look at in the analysis of the system, you work out this Jacobian matrix that is the partial of f with respect to these elements of u. This being a nth order system, so, you would get m Eigen values, and if you start them up in terms of the maximum and minimum of the real part, so, truly speaking, all these methods that you are talking about you can use it equally well for real equations or complex equations. So, that is what you need to look at the real part of it. The real part will tell you how this particular component is growing or decaying with time, if say time is the independent variable, and you look at the maximum of such values and the minimum of it; the ratio is what we call as the stiffness ratio, and as an illustration we talked about Stuart-Landau-Eckhaus equation, where this a 1s are basically complex quantities that you are looking at. So, are these coefficients alphas and beta i js; alpha i and beta i js, they are all complex quantities. So, you can see that if we remove the coupling term, the last term here, or this middle term is here, then you would note that the problem becomes somewhat easier. The stiffness does not come about because stiffness is associated with coupled ODE that is what we talk about here, that you have to have a multi-model system interacting with each other or you can have higher order differential equations. You can write them down as a set of first order equations, and then, you will perhaps get equations of this kind, and those ones would tell you whether you have a stiff system or not. So, if we remove these coupling terms, so, then it becomes somewhat simple and easy to solve. So, we talked about this example as a model which shows that it has two fundamental modes: one decays rather rapidly as e to the power minus 2000 t, and the second one gently decays with an exponent of minus 2 t, and if you try to solve this simple equation as it appears by the very best method that you can talk off, that admits a very small time step because you have to resolves this variation that restricts the time limit that you can advance your solutions, and this appears to be rather a small number implying that you have to do lot of work; even to solve a simple equation like this, for which you have the exact solutions. So, I may actually encourage you to get your hands dirty and try to really solve this equation yourself and find out what I am talking about, whether it makes sense or not. Now, of course, in this case the stiffness ratio is 2000 by 2. So, the stiffness is rather high, and one of the method to circumvent the problems of stiffness was traditionally to adopt what is called as the trapezoidal method and that is where we actually stopped. The trapezoidal method is really simple in concept that if you have this, let us say a 2 ODEs coupled here, so, the right hand side terms - f 1 and f 2 - continuously depends on time and this dependent variable y 1 and y 2, and we know those Jacobian matrix will give a very disparate time scale. So, what you would do is basically try to take an average of this right hand side at nth time step with the current time step. What it essentially does, it actually smears those Eigen values; the largest and the lowest averages out here and that gives you a little freedom. So, this was traditionally done, and if you look at the model problem that we talked about in the previous slide, this is essentially a linear equation. So, adopting trapezoid method for this is straightforward; you just simply take the right hand side, average it out, and it is easy to implement, if, because explicitly then you can get what they are, because this will appear linearly here, and so, this you can transpose on the left hand side and explicitly solve, and what happens by adopting this simple trick, will allow you to increase your time step by factor of 100; so, it is quite a big advantage that you can see. Now, what happens is when you have non-linear equations, the problem is not so straightforward, because then you end up with implicit equations and that would imply that you have to adopt some kind of an iterative procedure. For example, if I have the non-linear right hand side as symbolically written here, then what I would do is I will take the average of the previous time step where we have got the exact solution. So, I can explicitly evaluate this first term, whereas, this term what I would do, I will be evaluating it in a iterative sequence. So, iteration counter is basically m. So, what I do is I take some initial guess, work out some next step value, then use that back here and then keep doing it, and convince ourselves that we have reached some kind of a convergence criteria, and once we have done that, so, we are through for that time step, and this procedure has to be repetitively done for every time step, and you can understand that this will slow down the computation and you also, of course, have to do a little bit of a involved coding that would be a some additional chore for you to undertake. Now, while this was a situation for quite some time, a breath of fresh air actually came towards the late seventies - a new elegant method was evolved to tackle this stiff differential equation and that is what I am going to talk to you about, and this came out with this work of Ng and Reid, who were actually looking at solving problems of instability thus solving eigen value problems. So, essentially you may have some kind of a linear equation and fourth order equation happens to be a favorite for many physical systems. So, let us try to understand it with respect to this linear fourth order equation. You might think that linear equations are piece of cake, it isn't; in fact, let me tell you when it comes to fluid mechanics, where we spent quite bit of our time. We noticed that if we look at the generate non-linear governing equation, that is your Navier-Stokes equation, that is solved even by this calculated at the turn of last century. So, it wasn't that difficult while the corresponding eigenvalue problem involves having to solve this linear equation, and tell you what it was only 1951 the first sort of solutions were reported. So, you can understand that even though you may think that this would be easy to work on, it isn't so. Consider the fact that these coefficients a 1, a 2, a 3, a 4, they are functions of y. So, these are variable coefficient ODEs and they are quite difficult at times, and let us say your working domain ranges from y 1 to y 2, so, you would require auxiliary conditions; fourth order equation would suggest, we will require four of them. Suppose, we have to solve an eigenvalue problem, then we would surely have homogeneous boundary condition, right? Because this is a homogeneous differential equation and Eigen value relations require that you satisfy homogeneous boundary condition. Now, this is one of my favorite question, I keep asking whenever I have the opportunity - is this some kind of a black magic that you have the homogenous differential equation and you have homogeneous boundary condition and you get a non-trivial solution, how does it happen? You have been taught this from your, if not in the high school probably, surely from first year that you try to solve some kind of a Eigen value problem. So, you have the differential equation, there is no forcing, it is homogeneous. The boundary condition is also benign, there is nothing forcing the system still you get non-zero solution, how does it happen? I think somebody got really confused even comprehending the question, right? Shall I repeat? Yes sir The question is simply is this? I have a homogeneous differential equation, homogeneous boundary condition, how do i get a homogeneous solution? That is the basis for all your Eigen value problem, right? Because periodically it is repeated as you look at it What is repeated? Say, if you expand the vectors its function space, it is something like, it is a two layer shade. So, we can expect with a linear combination of infinitely many solutions. So, what happens? Still I am asking even a more fundamental question, it comes from your physics; talk about conservation principles. I have not done anything, I am, still I am getting non zero dividend, how does it happen? So, how does periodicity help there? And it is not necessary that I have to always use Fourier basis, I could use any other functions, any other orthogonal polynomial would do, right? Well, think about it. It is not very trivial; this is something that mathematicians and physicists teach wrongly in school. You give it to engineers, they will do better job, they will explain to you how this is done. If you look at it this way, suppose we are looking at some kind of a differential equation like what we have written there, say phi n is a function of all kinds of derivatives, all the way up to like this, and you are probably prescribing set of conditions, say phi, phi prime, phi double prime, etcetera, you would have n such conditions, that is what you do is define as zero for your Eigen value problem. Now, if I pick this problem and instead make this some kind of non-trivial boundary condition, that is what actually happens in a really life. Although you talk about homogeneous boundary condition, you would not find any system that works in beckon, right? You would have all kinds of environmental input. It is only that we cannot measure them, if you think back even tossing a coin ideally speaking, if you know all the inputs, and if you know the actual governing equation, you should be able to solve it as a deterministic system, you know, how to toss the coin and hope for heads or tails? It is only that lack of our knowledge of the dynamical system plus the lack of knowledge of the boundary conditions or auxiliary condition forces us to take a statistical route. So, what we are saying that just because we cannot quantify the background disturbances does not mean that we have this always equal to zero. So, in most of the cases, in real life system, you will always have something forcing there. Even in your differential equation, not all the processes have been modeled, and you know, in a dynamical system language, this is what is called as the process noise that our knowledge of the system itself is imprecise. Suppose I am tossing a coin, I do not know what happens if there is a temperature gradient in the air, or may be some kind of a gentle wind blowing, or some kind of a puff of gust is going on, they can affect the outcome. So, it is just that what we say borrowing some central limit theorem which you say that, if there are infinite such inputs to the system which are independent of each other, the outcome is what you get as a normal distribution, and then, you apply your statistics, and thereby, hangs the tail. So, now, if I look at it this way, what I can do here as we would shortly do. We will find that the solution that we are going to get. That would be something like the input times, what we call as the transfer function times the input. So, if I call this as the input, so, this is what we are going to get. So, phi will be some form of this kind, and now, what happens? Here, this is what we are writing, this happens to be your characteristic equation or Eigen value factor. This is what we will also call as the dispersion relation shortly. So, what you find that you try to find out those Eigen value which corresponds to the zero of this denominator, right? That is how you define your Eigen value problem. Now, what you can see is, what it says that, such a system is a kind of an amplifier, because even if you give B to zero, then you still have a 0 by 0 form and you can get a finite outcome, right? Finite solution, that is a 0 by 0form indeterminate form that will give you some kind of a finite value. So, what it actually tells you that even if B goes to zero, this ratio gives you a finite quantity, that's your mathematical answer. That is how you define your Eigen value problem. Physically it mean that you are not talking about the system in isolation, it would be that disturbance is drawing its energy from the lower order system. So, for example, if I am talking about flow instability, then I have a mean flow that is blowing, and then, in that mean flow, I am trying to see how, the, the disturbances are picking up. So, when I am talking about this growth of the disturbances, they have actually virtually access to a pool of energy coming from a mean flow. So, it is not violating any of your conservation principle, because the disturbance energy taps its energy from the mean flow, and mathematically, this is what you see, that is, what is the reason that you get that homogeneous differential equation. Homogeneous boundary condition give you inhomogeneous solution and that happens only for certain combinations of the parameters, that is, what you call is it your Eigen values, right? So, Ng and Reid were really guided in their approach in trying to solve such problems. So, two of the conditions, let us say are given here at the left extreme which are homogeneous, the other two conditions could be homogeneous or could be something like the variable decay into zero. So, phi and phi prime not exactly equal to zero but going to zero. So, that kind of a condition would be given, and what we would do is, what we have been telling all along that - even though, it is a higher order system, we would write it as a set of first order equations. So, it is a fourth order equation; so, we can define auxiliary variables like phi prime, phi double prime, and phi triple prime. So, this makes you a, gives a fourth order system, and what you can do is you can convert the 73 into a linear algebra equation. Phi prime is some A matrix times phi, and this is what would you get. So, of course, phi prime, that is the first component. If you differentiate it, you get the second, the, component, and second component, you differentiate it, you get the third component. The first three equations are straightforward from the definition itself. The fourth equation comes from the governing equation that we have written down. So, that is what you get. Now, as i say that, it is a fourth order system. So, will have four fundamental solutions, and let us say that two of them are physically admissible, that is what we say like we could, let us say create we have a domain here, and I am creating some kind of a disturbance here, and if we are talking about this as y, what I am saying that as i go away, far away from this the disturbances would decay. So, if I am following a constant, let us say x line and I have got that differential equation like this, like this. So, this is variation with respect to y for a fixed x. Then, as I go away from the source, of course, those two disturbances would decay and that is what we are, say the other two conditions would be prescribed as the limiting condition, and once we convert this differential equation in this form, then we will have two of the solutions which will be decayed with y that we will take it as admissible solution, that is what we are saying. Let them be phi 1 and phi 2. So, they actually satisfy the boundary condition, you may find that the other two would not satisfy your physical conditions that what we said that, it should decay with y. So, those two are of no concern to us. So, although, this is the problem of parasitic error growth. We talked about in the last class that two of them are admissible and two of them are non-admissible, and if you do not pay attention and care the numerical error, they could ride on those non-physical solutions and that could explode and give you a totally meaningless result. So, what we are trying to do is out of the four fundamental solutions. The ones which are admissible, we try to construct a solution matrix, which again we call it as phi consists of those two solutions phi 1, it is other three derivatives. So, basically, we have this four by two matrix. From there, what we could do is, we could construct various combinations 2 by 2 combinations. So, basically, what we are talking about? So, this method that we are talking about is called compound matrix method, and in the compound matrix method, we form the compound. For example, we take care of these two rows that will give you this variable Z 1 phi 1 phi 2 prime minus phi 1 prime phi 2. Then this is your Z 1, and now, if I couple the first and third row, I will get Z 2, and the first and fourth row will give me Z 3. Then, of course, there are other possibilities will be these two and these two, and finally, this one, well, you can very clearly see that it is going to be 4 c 2, because you have 4 y 2 matrix. So, you'll get six such combinations. And once you have obtained, let us say Z 1 as phi 1 phi 2 prime minus phi 2 phi 1 prime, I could differentiate it. I will get Z 1 prime as a phi 1 prime phi 2 prime plus phi 1 phi 2 double prime minus phi 2 phi 1, both of them like this, and phi 2 phi 1 double prime, of course, these two, we will cancel, and this one is, of course, nothing but your Z 2 itself. So, you get an equation from the definition itself that Z 1 prime equal to Z 2. You could keep doing that, keeping in mind that your governing equation is the fourth order equation. So, if i take Z 3 prime, of course, you can see that there would be phi 2 4, and here, we will get phi 1 4, and since they satisfy your governing equation, you can plug them in and work it out. You would see that you will get that Z 3 prime equal to a 3 Z 1 plus a 2 Z 2 plus a 1 Z 3. So, the same way we can work out on z 4 prime. Z 4 prime will not involve the differential equation. It will just have phi 1 prime and phi 2 triple prime and phi 1 triple and phi 2 prime, that would give you nothing but equal to z five, and if you look at z 5 prime, that will involve the differential equation and same way z 6 prime also will have this. So, you can see what we have done here, we have converted a fourth order system into a sixth order system. So, we have already decided to do little extra work. We have increased the order of the system, and what do we get we will tell you shortly through an example, that is the way I like to communicate that this compound matrix method is generally developed to study problem of instabilities with basically the modes that actually grow in space and time. So, what happens is for such a system, instead of working on the original problem, here phi prime equal to A of y into phi. If we solve instead this altered system compound matrix or variable system, that is much more easier. So, let me just demonstrate this by an example from fluid mechanics flow instability problem. This equation was derived by Orr and Sommerfield, and this is, as you can see a fourth order, well, d is the operated d d d y. So, you get a fourth derivative here, on right hand side, and those of you who are familiar with little bit of fluid mechanics would recognize this coming from the viscous dissipation, whereas, these are the terms which comes from convection, but let us not worry about fluid mechanics, let me state it as simply a problem of stiff differential equation, and then, let us see what we can do. So, basic problem is of this kind, let us say we are talking about flow where a simple geometry like even a flat plate and the flow is coming like this, and we have already done a few lectures ago that it creates a kind of a boundary layer and we have talked about that solution as u of y. And now, what we are seeing that giving that mean flow, we are trying to figure out what the disturbance is going to do and this disturbance quantity is defined by the stream function, and for this problem, if i am talking about a simple 2 D problem, then this is written down in terms of Fourier-Laplace transform with the basis given like this. See, that is what I just now told few minutes ago. In actual problem, if you really want to excite some kind of instability, you would have to sort of give some kind of an excitation and we will talk about it. Let us say I had put in a kind of a vibrator there to create a disturbance, and this vibrator that we have, that is oscillating, let us say at a frequency of omega naught, that is what this omega naught comes in here and we are just trying to find out. If that time dependent excitation creates some disturbance which grows in space, so, that is the alpha we are looking for. So, we are trying to look for complex alpha which will be governed by this Orr-Sommerfield equation, that is what you can now see that phi is that, well, that I think this is the lower case phi. So, let us keep it a lower case phi and then what we are going to do is obtain the solution given by that Orr-Sommerfield equation for all possible values of alpha and then integrate this. Well, those of you are from mathematics department. You would recognize or if you do not, this has to be evaluated along what is called as the Bromwich contour. So, not all contours in the alpha plane are admissible, you will have to find out which way the wave is going on. You have to construct a contour that does not violate physical principles. So, that is what we call as the Bromwich contour. But for the rest of you who do not bother, I am just stating the problem. Physically, this is the issue that there is a mean flow coming and that gives sets of this boundary layer, that is your mean flow, and then, we are creating a small disturbance at the wall and we want to study how this disturbance is going to manifest itself in the flow field, and that is what this equation tells us that the mean flow here u and its second derivative are involved. So, you can see, that is what I was telling you that in all disturbance evaluations equation, you would have input coming from the mean flow. It does not just evolve by itself, there has to be source of energy from where it will tap its energy. So, that is how this mean flow terms come into this equation, and these are various mechanism. That is what it does. and So, this is the statement of problem, and as you can see, we have just shown how this has to be converted into a sixth order system, and the boundary conditions, well, I have written it down phi and phi prime equal to zero. But if you give a excitation here, you do not need to do that you have a specific finite disturbance being given there, so, that will change to inhomogeneous path. Whereas, if you go far away from the source you would require those disturbance quantities to decay, and that would imply the phi and its derivatives to approach zero. It is not equal to zero; it is a different type of boundary condition asymptotic condition. Now, to get a glimpse of what happens in a problem, like this we have seen that this creates a kind of a boundary layer. So, outside the boundary layer, what happens? The flow would have reached its constant value, right? It would reach its constant value that is what we are saying that if non dimensional value would be equal to one and the derivatives are all zero. If the derivatives are zero, then although your governing equation is a variable coefficient ODE in the free stream, you will get this as a constant coefficient ODE, and once you get that, it is easy for you to start off with a given mode like this e to the power lambda i y plug it in this equation, and then, you would get a polynomial for lambda and two of it, you will see, they are given by minus and plus alpha and two of it is given by minus and plus of q where q is this. So, everything is obtained from this equation 77 and what you notice that what I told you that couples of the solutions are admissible couple of them or not that is why I particularly demonstrated this example to you. Now what happens is suppose alpha is a complex quantity, and if its real part is positive, then what happens one of the modes is growing with y, and that violates the principle which said that if we go far away, it should not decay to zero. So, what you immediately notice that if I take the modulus of alpha and if the real part of the alpha is positive, then, of course, lambda 2 is not admissible, lambda 1 is admissible, because that is the solution which will decay as e to the power minus alpha y. Same way you can see I am extracting a square root from this quantity, so, one of it will be positive and one of it will be negative, and the one with positive real part would again grow with y and that will be inadmissible. So, what you are seeing that, although it is a fourth order system, only two modes are admissible, that is one thing. The other thing is most of this problems of instability arises, when you are this parameter. Here, the Reynolds number happens to take large values, you know, for flows at lower Reynolds number, you do not see instability is a characterized by certain range of parameters, and in this case, it so happens this r e has to be large. Well, I have just shown you sort of estimate, it need not necessarily be have to be infinity but it should be fairly large, and if it is fairly large, you can see modulus of q is much larger compared to modulus of alpha, and this is the stiffness problem, because you have one solution going as e to the power minus alpha y. Another part is going as e to the power minus q y, and now, you can see the stiffness ratio would be modulus of q by modulus of alpha. So, that is why I chose this particular example to show you all the certain issues involved in solving stiff equations. So, what you do is, of course, you write down the fundamental solutions in terms of all those four quantities those are possible from mathematically, but out of those four, we have seen in the free stream this phi 2 and phi 4 grows with y and they have to be removed. So, they are not physically admissible, then what we are talking about, we should have the composite solution in terms of phi 1 and phi 3 which are admissible. Now, what we say that at the wall, this phi and phi prime are zero. So, that gives us this condition, right? So, these are two equations for these two coefficients c 1 and c 3, and as I told you that if it was not homogeneous boundary condition problem, then this right hand side would be non-zero. See, this will be non-zero, and then, you can straight ahead solve for c 1 and c 3. So, that is one thing, but how do we get the characteristic equation? We are looking for those combinations of complex alpha for which the determinant of the matrix formed by phi 1 phi 3 phi 1 prime, and phi 3 prime gives you a zero. So, that c 1 and c 3 would be right hand side vector will be zero, and if I take an inverse of that matrix, that will give you a determinant in the denominator. So that determinant in the denominator is this right If I switch back, you can see that is nothing but this phi 1 phi 3 and phi 1 prime and phi 3 prime. If I take that, that is what we are going get and that is your characteristic relation. We will also call that as a dispersion relation, and I will define what dispersion relation is, but looking at this expression, you can see what we are implying by dispersion relation. We are fixing omega naught and we are trying to find out that complex alpha which will satisfy 79. So, it is not all alphas are admissible, you will pick and choose a particular set of values of alphas which we will call as the Eigen values, and those complex alphas are the Eigen values for this fixed omega naught and r e right. So, this is something that we have to take home from here that, what is this quantity? Now, if I look at the definition of Z 1 that is that quantity. So, the dispersion relation that we wrote down there comes about the first compound Z 1 evaluated at the wall that has to be set equal to zero. So, that is your dispersion relation. So, this basically is a statement of Eigen value problem and we can study it. So, if we now look at what we have here, we can plug these particular values a 1 equal to zero a 2 will be this a 3 is zero and a 4 is this, and what we need to do is develop a methodology of solving this problem. The methodology for solving this problem is the following that your dispersion relation is, a, that Z 1 equal to 0 here at y equal to 0. So, what we tried to do is, then we tried to start the solution from faraway in this free stream, and then, we march downwards. So, what you do is - you assume some value of alpha, because you have already fixed omega naught and r e. You start from there and your governing equation is this Z dot or Z prime whatever we have written times this B matrix times Z. So, that is, that is the equation that we are solving. So, this Z prime will have, those Z 1 prime Z 2 prime that we have already evaluated. Now, what you have noticed? This is I am taking a little time to explain to you was Z 1 is this, and what we found that phi 1, in this case, goes as e to the power minus alpha y, and phi 2, in this case, goes as e to the power minus q y, we have worked it out. Now, if I look at this definition of Z 1, so, we all looking at Z 1 for y going to infinity. There, the, we have the fundamental solutions. You see, this phi 1 and phi 2 is not the same thing, everywhere it is at that part of the flow where you have the constant coefficient ODE, that is why we could get that exponent inside where u is not equal to one, and u double prime is not equal to zero. We do not know what these exponents are right. So, this is the issue that you have to understand. So, in the free stream, we know the variations like this and we can plug it in there. What we are going to get? We are going to get, we are going to get phi 1 would be something like this e to the power minus alpha y, and phi 2 prime will give us minus q or e to the power minus q y and minus e to the power minus q y, and if I differentiate it, I will get a plus and I will also get alpha e to the power minus alpha phi. So, this looks like this alpha minus q e to the power minus alpha plus q y. So, that is the variations of u in the free stream, and this is where you see the benefit of this compound matrix method comes out, because in the original problem, we had a admissible part which was going at this rate, but the solution was problematic, because mod q by mod alpha is very larger here. In defining these new sets of variable, we change the variation. Now, you see the variation is given by this quantity e to the power alpha plus q y with a minus sign embedded there. You would see that, that would be the situation for all the six variables with identical variations in the free stream. So, what it tells us that by adopting this procedure, we have been able to remove the stiffness of the problem. Once we have done that, then what we could do? We could adopt any of the methods which work for normal non-stiff equation. We can use it here, right? So, that is the whole idea in the compound matrix method and what you would find that what I have shown here for Z 1 y going to infinity, it goes like this, and this will be nothing but minus, and so on, so forth. So, what we could do is we have an estimate for the qualitative variations of these functions with y. So, what I could do is I can see a common factor here appearing. So, I could factor it out from all the quantities right. So, that is what I did here, we did Z 1 equal to 1 here. If I do that immediately Z 2 minus of alpha plus q and you check it out for yourself that Z 3 will be alpha square plus alpha q plus square and so on and so for, forth. So, that is your starting solution that is what I said that you start off with the solution. Those are given here Z 1 to Z 6, and then, you keep marching and come to the wall, and check what is the value of Z 1. If Z 1 happens to be zero, you have an Eigen value. If Z 1 is non-zero, then that alpha is not your Eigen value, so, you would reassign a value to alpha and do this process. And how do you reassign such a thing? We know one of the standard technique is to use Newton-Raphson method. I think you would have heard of even if you haven't pick up any elementary book on Numerical Analysis, and you will see, they talk about Newton-Raphason search and that will tell you how to obtain these Eigen values. Well, there are various other ways. We do not actually use Newton-Raphson method. We have a better way of finding it out but I am not going to spend any more time. But what you notice here that it is a clear example where stiffness of the problem has completely been eliminated, and you have the set of first order ODEs given like this, and you can use your RK 4 method that we just now talked about in the last class. So, you can see life can be made very simple by doing little bit of work. On the desk, evaluating all these compound matrix variables working out this initial conditions and write out a simple solver, and I have a people waited for seventy eighty years to get here. So, it was quite a bit of achievement by this group of people. So, in compound matrix method as I told you that, all the variables have equal growth or decay rate. That is why we can factor out the stiffness of the problem, and once the stiffness of the problem disappears, you could use any standard initial value problem routines to integrate, and as I also told you that by some other method, people have obtained this solution for the first time in 1951. So, I suppose that leaves us with the task of figuring out what the Eigen functions are. For example, if I have the Eigen value, alpha i could obtain the Eigen functions, which are nothing but the phi itself. So, that could be written in terms of c 1 phi 1 and c 3 phi 3. What you do is you could differentiate this few times to get phi prime, phi double prime, and phi triple prime, and then, you can eliminate c 1 and c 3, and if you eliminate those equations or those constants c 1 and c 3, you end up getting equations of this kind. I have just noted two out of four. There are four possible solutions. This is a fourth order system, so, we get four equations. Unfortunately, the problem becomes difficult though if the order increases. When we tried to do it for a fluid flow problem coupled with heat transfer, that makes it a sixth order ODE - the original equation, and then, the compound matrix variable becomes 62. So, you can see that it is somewhat more involved and you also get a very large number of candidate Eigen function equations. There are four possible solutions. So, out of those four possible Eigen function equations, only these two are the ones which work. The other two has some problem. I do not wish to getting to that part of the story. It is something which we have worked upon and we have suggested why those two equations are not to be trusted, because you see, what happens that our governing equation was a fourth order equation, but there are only two admissible solutions, right? So, if I look at the first equation, it is a second order equation. So, it has two modes and it happens that those two modes are the ones that we have retained. This one is a third order equation, but you can immediately see that there is a solution which does not vary with y, because it starts off with phi prime, right? If i write it has a sort of a polynomial for the operator d, so, it will be d cube d square and d. So, you can factor d out though corresponding to that, you will get a solution which does not vary with y. So, you can immediately see one mode is neutral. It does not vary with y where as the other two modes that you retain; they are the ones that you wanted them to be there. What happens, those other two equations actually are quite dangerous, because they involve some spurious numerical modes those which will be something like alpha plus q and all kinds of or plus alpha root. So, what happens is with this statement, I will stop talking about compound matrix method and stiff differential equation. You will also notice one curious fact is that if your system is defined by a second order ODE, your compound matrix method is very easy to work out. We will find out that it converts into a single equation. So, for a second order system, you have a benefit in a fourth order. You actually go from fourth order to sixth order system. For a sixth order system, you go to 62 equations and it complicates with a multi-physics problem, because suppose, you are talking about plasma control of flows, then you would be talking about the Maxwell's equation, and the order of the system keeps increasing. So, you probably go from sixth to eighth order, and then, you can see these problems of number of equations increasing rapidly can hurt. But at the same time, you do not have any alternative, because this is such a nice and clean method. The other possible way of solving stiff differential equations are fraught with danger, it does not come as easily as what you can get here. I think with this, we should stop our discussion on ODE. From the next class, we will be starting to talk about PDE, that is where the major interest lies, because the original governing equation in most of the cases are governed by partial differential equations. So, we will begin talking about classification of PDEs. The PDEs, of course, are not as you can get in a single basket that what we could, with for ODEs, although we have seen, there are variations, we have classified in terms of initial boundary value problem. We talked about stiff and non- stiff problems, but PDEs much more colorful, gives you lot of room to maneuver, and there are very interesting aspects of solutions of PDEs. Having said that, I would also add that we would not follow the classical route of solving, that is a promise which I would like to make today, and you will see, as we go along that, we will not be classifying them. As you may have done in your other math course in terms of parabolic elliptic and hyperbolic, those are old hat. We have to do something better and we can do. We will talk about that. So, I will stop here.