 So, we have been having look at differential equations, ordinary differential equations with initial value specified and we have looked at different methods of different numerical schemes of solving ordinary differential equations subject to initial conditions. Then we also had some introduction to how does one do analysis of convergence, convergence of error between the true solution and the approximate solution. And finally, we saw one application of using a initial value problem solver for solving a boundary value problem, ordinary differential equation boundary value problem. So, this was shooting method and towards the end I also mentioned something about stiff differential equations. So, stiff differential equations are the ones in which different variables have different time scales. So, some variables are acting on a fast time scale, some variables are acting on a slow time scale and that gives difficulties in choosing integration interval. So, convergence and integration interval, choice of integration interval were tightly coupled and we had to look at Eigen values of the Jacobian matrix to get some clues into how to choose integration step size. We also looked at a fix that is variable step size method by which when we do not know how to choose the integration step size, we could get over the difficulties. So, today I just want to briefly touch upon one very very important topic. This is solving differential algebraic systems or these are popularly known as DAEs. In fact, in real chemical engineering applications, you are more likely to encounter differential algebraic systems than pure algebraic or pure differential systems because in many situations some phenomena have to be modeled as algebraic constraints. A classic example would be multi-component flash. An isothermal multi-component flash, I am just reducing the complexity by taking isothermal flash. When you are doing design of isothermal flash, you are typically concerned about the steady set operation. When you are looking only at the steady set operation, you only have to worry about algebraic equations. There are no derivatives coming into the picture. But if I want to simulate dynamics associated with this flash, more complex problem of course would be distillation column. Flash is only one tray, like a model of the flash might resemble just one tray. But just to give you an idea here of how differential algebraic systems arise, I am just going to write down differential equations and algebraic equations under the assumption that temperature is specified and constant. This is the component balance for the ith component. Overall mass balance, component balance. What is overall mass for ith component? Then total mass, liquid mass plus vapor mass. Then you have this yi is equal to ki xi, thermodynamic relationships. I am just writing this in an abstract way. There is some thermodynamic model for computing ki, ki values. And then the relationship between yi and xi is given by this. This is the thermodynamic model representation. Each ki will be found by some thermodynamic model. Well, then you have total volume in the flash. You will have models for the density in the liquid phase, density in the vapor phase. Density in the liquid phase will be function of composition, pressure, temperature. So it is the case for the vapor phase, vapor phase density. And then you will have two basic composition equations. Sigma xi equal to 1, sigma yi equal to 1. Or you can combine them and say that sigma xi minus sigma yi is equal to 0. Remember, I just want to do dynamic simulation of the isothermal flash. So I cannot ignore the two dynamic terms or two differential equations that I get. Not two differential equations. I am sorry. If I have n components, there will be n plus 1 differential equations. There are n components. There will be n plus 1 differential equation. One for each component balance. And so there are n components. So n differential equations for each component plus 1 overall mass balance. So there are n plus 1 differential equations and many more algebraic equations. Some of the algebraic equations you might be able to eliminate and substitute but not all. In such cases, when you cannot eliminate all the variables, you will get a differential algebraic system. And then I am just writing two more equations. You have to write the equation for the liquid flow out. Liquid flow out will be function of the mass here if it is gravity kind of flow. And vapor flow out will be function of pressure minus pressure in minus pressure out. Here also it will be function of pressure in minus pressure out. So I am just writing a general equation here. L minus these are the two flow equations. Now the variables that are specified in this particular case, these correspond to T, F, Z, V and P out. These are the variables which have been specified. Isothermal flash, so T is specified. V is specified because you are designing the flash for certain V or you want to keep the V constant. Then feed conditions are specified and the pressure out is specified. So these things are specified. And then you want to actually solve this coupled differential algebraic equations together. A typical problem in chemical engineering where if you have a reaction occurring there with two feeds coming in becomes more complex. You have conversion from one to another and then you have to maintain quite likely that your products are gas phase and the reactants are liquid phase. So a very common scenario. And then you are drawing the product. So differential equations and algebraic equations are coupled and this is what we have to solve together. So just algebraic equations, just differential equations. These situations arise only in certain cases but if you are doing dynamic simulations more often than not you get into DAE systems. Another example I had given previously of DAE systems was orthogonal collocation. It is not necessary that when we just get DAE systems, you get DAE systems when you get, when you are solving only dynamics. If you remember Tibler reactor with axial mixing, Tram problem. When we discretize this problem in space using orthogonal collocations, you get the boundary conditions or algebraic conditions. Of course in that particular problem was dynamics of Tram but if you were to do solve the Laplace equation by method of lines. So you have differential equations in one direction, in the other direction you have discretized and you might get algebraic conditions. So DAEs need not arise strictly when you have time. It can arise, we have only two special coordinates and you are discretizing in one special coordinate, you are solving as a OD in the other coordinate. So it can arise in some other context, not just time simulations but in time simulations quite often you get differential algebraic systems. So solving differential algebraic systems is much more complex than solving just the differential systems or algebraic systems and you should know at least something about them. So the idea of this one lecture is to give you a very, very brief introduction to some of the issues that are involved in solving DAEs. What I am going to do is, as I said just touch the tip of the iceberg but you should know that there exist separate methods for solving DAEs and also toolbox like the toolbox in MATLAB or some other thing like SILab would have solvers which can handle DAEs. MATLAB as a solver, OD15S stands for stiff. So OD15S solves differential algebraic equations simultaneously. So now let me abstract this, let me tell the most general form of differential algebraic equations. I can write as f y, y prime t equal to 0 where y is a vector and you have some coupled equations of derivatives and algebraics. So some of these could be algebraic states, some of these could be differential states. What I mean by algebraic states and differential states? We often write what is called as semi-explicit form. We often deal with what is called as semi-explicit form. The semi-explicit form is dx by dt is equal to, well let us use a little bit of a different notation here. Let us make this capital G and because we have been using f for derivative x, z, t and 1 minute or g. We will use here h. h is some non-linear function of y, y prime and t and then I can write a semi-explicit form which is 0 equal to g of x, z. The flash equation, you can rearrange into this semi-explicit form. dx by dt is equal to f x, z, t and 0 is equal to g of x, z where x we call as differential variables and z we call as algebraic variables. Differential variables as in you get dx by dt terms. You get dx by dt term. Algebraic variables because in the raw form you do not get derivatives of z. In the raw form you do not get derivatives of z. You may have to find out derivatives of z for some reason, we will come to that. But right now we are in this lecture we are going to look at this simplified form. Most of the cases it can be, most of the cases that we encounter can be written into this form. Then we want to solve them simultaneously. Now you can appreciate here that these are relatively difficult because initial condition for a ODE initial value problem, you can specify an arbitrary initial condition. Here you cannot specify an arbitrary initial condition. The initial condition must satisfy this constraint. The non-linear algebraic, linear or non-linear generally in chemical engineering problems these constraints will be non-linear. There are non-linear algebraic constraints. This must be met by the initial condition on x. So you have to give an initial condition on x and z which satisfies this constraint. That is very, very important. As you march in time, at all times this constraint has to be satisfied. So solving DA is a scale more difficult. In this classification, so this particular form is called as fully implicit. This particular form is called as fully implicit whereas this is called as semi-explicit form and in this x we call as differential variables and z are, and of course y is the differential is nothing but x, z put together. Actually the semi-implicit form can be written in the fully implicit form if you take this definition. If you take this definition then you can convert semi-implicit form into, so that is not really, so we will only look at semi-implicit form to simplify our life. Now how do you solve this? So solution approach is, one is called as nested approach and the other one is called as simultaneous approach. So the nested approach, so first thing we do here is, so basically we get g of, the first approach is that given x n you solve for this and the next step would be, now what it means is that, what it means when you are solving this differential equation you are going from n to n plus 1, internally when you are computing your function derivative you may have to solve every time inside you may have to call Newton Raphson method or Newton's method. Within your function evaluations, suppose you are doing Rangekutta, in Rangekutta you have to do evaluations of the function at intermediate points. At intermediate point when you do calculations you will have to solve for g x n. So you will have to solve for, you need an explicit way of or you need an implicit function. So this approach requires, this approach requires z to be obtained as a function of x, internally every time in let us say Rangekutta if you find an intermediate x for that x you have to compute corresponding z and then proceed you cannot do. So internally every time you have to compute or you have to solve the algebraic equations every time. So that is one approach. The second approach is simultaneous approach. The second approach is the simultaneous approach. In the simultaneous approach you use an implicit method, you use an implicit method to solve the differential algebraic equations together. So here I would say I am using an implicit multi step solver. For example simplest implicit multi step solver would be implicit Euler. In implicit Euler what I can do is I can solve for x n plus 1 minus x n minus f x n z n equal to 0 and sorry, x n minus f x n z n equal to 0. This is implicit method. So x n plus 1 z n plus 1 equal to 0 and g x the simplest implicit solver would be implicit Euler h times. h is my integration step size h times this. So this is my simplest method implicit method. See what is the advantage here? Advantage here is that both algebraic equations and differential equations get solved simultaneously and you get when you start from x n z n which are consistent you get x n plus 1 z n plus 1. Automatically x n when you are going to the new point x n plus 1 z n plus 1 they are consistent because you are making sure this happens. When you march again x n plus 2 z n plus 2 are consistent with the algebraic equations. So implicit Euler and how do you solve this if this g is f is a non-linear function you have to solve this using Newton's method or some non-linear algebraic equation solver. So this you have to treat them as a simultaneous this is this is these two equations have to be solved simultaneously x n is known to you x n plus 1 is not known z n plus 1 is not known you have to solve them simultaneously and then every time you solve you march in time you get consistent solutions. The other general class of multi step methods which are used here they are called as backward different solver BDF solver more popularly known as BDF. So these are in our multi step parlance these are methods of x n plus 1 is equal to h beta minus 1 f. See we had earlier seed Adams molten gears algorithm predictor corrector. So this is the class of this is general multi step algorithm here we do not use we only use one derivative that is f n plus 1 and all past x values and then here you solve this together with g of so we solve this equation and this equation simultaneously we solve this equation simultaneously. So the simplest example of this type is implicit Euler a BDF solver backward difference formally these are called as BDF formally. The nice thing is that you have x n plus 1 z n plus 1 appearing in derivative and here and then only use past x we are only using past x we are only using past x values. By the way some of you had this question about if you have multi step methods how do you initialize because in multi step methods you need suppose is a p step method you need first p values. So one approach suggested by Petzolt and Asher a very good book on ordinary differential equations by Petzolt and Asher. So they say that you start for a p step method you start one with Euler method then two step method then three step method four step method. So for first p steps you use one to p minus one step methods generate the initial points after that you continue using p step method. That is one way of initializing. So this backward difference formally so just to give you an example of this would be x n plus 1 is equal to two third f. So this is one of the backward difference formula and then you can write more you can get this from the textbook or you can derive you know how to derive this BDF formally you will be able to derive. So we are only going to use the derivative f n plus 1 and all x values from x n to x n minus p you can derive this coefficients. So likewise there are many methods the other way of solving this differential algebraic equations I am not going to do it on the board because we have done orthogonal collocations quite a bit on detail but other way of solving is orthogonal collocations. What is the advantage of orthogonal collocations? It converts a differential equation into algebraic equations. Then you can club the algebraic equation coming from orthogonal collocation discretization and algebraic equations which you have here and solve it together using Newton's method. So advantage of orthogonal collocations is to conversion to algebraic equations and then everything can be solved as one set of algebraic problem. So that I am just leaving it to you to work out yourself maybe in the exam or otherwise when it comes to. So these equations can be solved using these BDF approaches or using something like Lange Kutta if you have you know implicit function z as a implicit function of x then you can solve it using you know you can call internally inside your Lange Kutta you can call Newton's method and then solve it. So these are the way of solving this but solving these becomes much more complex than just solving the ODEs because in the function evaluation itself we have to call you know algebraic equation solver which could be non-linear so iteration it inside iterations and so on. So it can become quite complex okay. Now unfortunately just these equations you know it is not just unfortunately just using these methods for solving this because there are many other issues that come up when you have differential algebraic systems together. One of the issue that you must know is index of a differential algebraic system okay and what is difficult to solve is high index problem okay we will define what is index of a differential algebraic system index one systems are easy to solve high index systems are difficult to solve and you have to do some tricks. The main thing is that when you have a differential algebraic system you have to have the initial condition consistent okay that is X0 and Z0 should satisfy the algebraic constraint that is very very important okay. There is a textbook by Asher and Petzult on ordinary differential equations which a very good textbook which gives a very good mix of algorithms and theory talks in detail about differential algebraic solvers. There is a separate chapter dedicated to this systems and then Petzult has developed a software which is very very commonly used I think it is a freeware that is called Dasil okay. Differential algebraic equation solver called Dasil which is you can I think this is a freeware it is a full-time program you can download and then call convert into a MATLAB code and call it through MATLAB at least some small dimensional systems can be handled. Now what is index okay the index is number of times you should differentiate the algebraic equations okay to get a set of ordinary differential equations ODEs okay now you might solve it say well what is so looks like you can convert the system into set of ODEs and then you can solve it using ODEs solver. The problem is when you have algebraic equations okay which means some derivatives are 0 all the time okay this resulting system will be highly stiff system it is not at all easy to solve. It will be a stiff system though theoretically you can do this it will be a stiff system because typically when you arrive at differential algebraic system of equations you have chosen to neglect certain time derivatives because those time the time scale on which those derivative operate are very very small you neglect the time and you just look at algebraic constraints. So even though on paper you can do this by repeatedly differentiating the algebraic constraints this coupled set of equations will be difficult to solve. Anyway so number of times you have to do differentiation to arrive at a set of differential equations is called as index and in general index 1 problems are easy to solve okay if you differentiate the algebraic constraints only once and if you can convert them into set of differential equations great. If you cannot index 2 systems okay or index 3 systems or high index systems are difficult to solve very very difficult to solve and then you have to convert them typically the approach taken is to convert a high index system into a low index or index 1 system and then solve it as a index 1 system okay. Let us look at an example where we get a high index problem I will take a simple system just a mixing tank okay so v is the volume here v is the volume and one differential equation okay. Now for this particular problem for this particular problem I will just show you two phases of this problem okay in one case you can solve it as an OD in other case it turns out to be a differential algebraic system with index 2 okay. So one way is okay one way is you know C0 t inlet concentration is equal to some function beta t is specified okay and then initial value of C0 is specified of concentration inside the tank is specified simple OD problem okay you substitute this is a DAE but then you know you can substitute for beta t here and you are given initial condition it becomes one OD you can solve it that is very very easy problem okay. The second one actually the other problem okay so this one the first one is very very simple DAE that is dc by dt is equal to C0 t minus Ct by tau and C0 t minus beta t equal to 0. So this is a DAE and substitute for beta you can solve this using an OD solver not at all an issue simple index one system not at all difficult to solve the second one here is earlier I had specified C0 inlet condition and I wanted to do simulation okay now I want to specify C of t okay I want to specify C of t as a function of beta t as some function specified function of time is equal to beta t, beta t is generic way of writing a specified function I am just putting this as a beta t. See now I have to solve for dc by dt is equal to okay it is very very important how I choose the initial condition here this is an example in which I have specified Ct okay I have specified Ct this problem this differential algebraic system it turns out to be an index 2 system and this is much more difficult to solve. So well I will just quickly write down the solution and then we will also look at why it is index 2 system that maybe we have to do it in the next class but the solution for this problem is you can see that solution to this problem solution to this problem will not be Ct, Ct has been specified solution to this problem will be C0t the inlet profile okay so this is a concocted problem okay it actually tells you something that is intrinsic it difficult about da's well you might wonder where do you get this problem for example when I am doing control I am specifying how should the output behave I want to find out an input okay in the control problem I am specifying how the controlled variable should behave and then I want to reconstruct the input that will give me that particular control output okay so we will continue analyzing this particular system we will see that this particular system turns out to be an index 2 system and then we have to solve it by doing some more manipulations it is very very important that the initial conditions are consistent if you give a wrong initial condition you can get into trouble here.