 how to solve these problems. So, my lectures are going to be of course, centered on algorithms, but the derivation of algorithms, I am more about worried about how these algorithms are derived rather than what are these algorithms. So, I am going to derive Runge-Kutta methods, one class of algorithms which are based on Taylor series, the other class of algorithms based on polynomial or Weierstrass theorem, multi-step methods or predictor-corrector methods that is another class of algorithms. I will be if you touch upon orthogonal collocations, then we will move on to talk a little bit about convergence or numerical stability of these algorithms, why one algorithm has better properties than the other. So, comparing the algorithms, what is the basis for comparing the algorithms? After that, we will spend time on some special applications like there are methods of converting a boundary value problem into initial value problem, these are called as shooting methods. We look at shooting methods and then if time permits, we will have a peak at what are called as differential algebraic systems. So, differential algebraic systems are one in which you have to solve simultaneously differential equations and algebraic equations. In fact, a problem which is very, very often encountered in chemical engineering. Though it is very, very important, we can probably get time only to touch upon some algorithms, some ideas of DAE systems. So, let us begin our journey into ODE IVP. So, last class I talked about a general set of ODE IVP which I am, which is dx by dt is equal to f of x t. Initial condition corresponds to x naught, x belongs to Rn, n dimensional vector. And then I want to find out the time trajectories or I want to find out spatial trajectories. So, actually to generalize, I had made a small modification. I had just said here, some independent variable neta at time 0 or at spatial coordinate 0, whatever it is. At the initial condition is known to you, you want to integrate this differential equation over some finite interval. I want to integrate this differential equation over a finite interval in variable independent variable neta. That is my problem. Before we begin talking about the methods of solving, like we discussed a little bit about existence of solution, uniqueness of solution and so on. When we talked about solving linear algebraic equations or towards the end, I gave a lecture on existence of solutions for nonlinear algebraic equations. I am going to very, very briefly mention about this aspect, not really get into deep into this. Now, it is very, very important. So, there are three issues that are primary concern in the beginning. One is existence of solution, that is solution or a solution. First, essential thing that a solution for the given differential should exist. See, typically except for some system of differential equations which have been conceived by mathematicians, typically most of the differential equations that you want to solve are some models of a physical system. Now, a physical system exists and just like a physical system has different behavioral patterns for different starting conditions, a model also has different behavioral patterns for different starting conditions. For a specified initial condition, a solution should exist for the given differential equation that is very, very important. It may happen that you have done a wrong modeling, solution does not exist, physical system of course exists. You have developed a model for the physical system and there are some wrong assumptions or there is something wrong with your understanding of the system and no solution exist for the initial condition that you have specified. Now, if the solution exist, we are worried about one more thing. We are worried about the uniqueness of the solution. So, once I have a existence of a solution, I am worried about uniqueness. So, here concern is each set of initial condition. Uniqueness of solution is another important aspect of the solutions of differential equation. Now, why am I worried about existence and uniqueness which appeared to be very, very mathematical concepts. Actually, they are just abstractions of something that we know from the reality. First of all, a behavioral pattern exists for a real system. So, if the model represents a real system, a solution should exist for a given initial condition. Second thing that we know is that this I would call something like a principle of determinism that if you for most engineering systems, if you start with same initial condition, you will get the same identical behavior. Very, very important, you will not get different behavior under identical conditions. That is mathematical way of saying it is that solution should be unique. I should not get different solution for you know same initial condition does not make sense. I should get identical behavior for identical initial conditions. The third one is the third important aspect is continuity of solution with respect to initial condition. The third aspect is continuity of the solution with respect to initial condition. We want that the solution of the differential equation should depend continuously on the initial condition. Now, when I talk of continuity, you think of that epsilon delta definition and then you might think, my God, this is complex. It is not. If you try to understand what is it trying to quantify or abstract, then it is not so difficult. What we know from, what we know from, what do you know about continuity? About continuity, what we know is that will a small perturbation in the input will change, result in a small perturbation in the output. A bounded change in the input will it lead to a bounded change in the output is what we are worried about in the continuity. A function when it is continuous, in a very crude way of saying it is a small change should always, any amount of small change should lead to a small change in the function value. Now, we know from working with real systems that even though we demand uniqueness of the solutions, when you actually perform experiments, you can never repeat identical conditions. Even if you decide that, let us say you take a simple experiment that you have a pendulum and then you want to understand motion of the pendulum, you start from some theta. Let us say you want to start from some theta 0 and then you have developed equation for governing the pendulum. You want to understand those equations. You can, even if you want to repeat the experiment from same theta 0 every time, it is impossible to conduct an experiment which will be every time starting from same theta 0. If you do 10 experiments, it will be some theta plus delta theta every time. Now, what is important, what we know from physics, what we know from observations is that minor change in the initial condition leads to minor change in the solutions. It does not lead to significant drastic change in the solution. If a minor change in the initial condition leads to a drastic change in the, that does not happen. In the real systems, if you change the initial condition of this system slightly, the solution that you obtain for, see if I start from, first I start from theta 0. Next time I start from theta 0 plus some delta theta. Then next time I start with theta 0 minus delta theta. The solutions that I get for each one of these cases of the motion of the pendulum will be close. This you know. The way of quantifying this mathematically is saying that the solution is continuous with respect to the initial condition. If I perturb the initial condition a little, the solution will get perturbed a little. Very, very important. Perturbation in the initial condition can occur when you are solving a differential equation. Somebody might, I have given initial condition which has some variable is one-third. Somebody might decide to do it 0.33. Somebody else might do it 0.333. Somebody else might do it 0.3333. And then the solutions that you get should not drastically change because he takes 0.33 and she takes 0.333. They should be close by, similar solutions. So the third aspect which is important is that is the solution continuous with respect to initial condition. I am not going to go too deep into this. I am just going to state one theorem which talks about the continuity aspect, existence, uniqueness and continuity. So we will define something called a region. I am defining a region in n plus 1 dimensional space. n plus 1 dimensions because x is a n dimensional vector, eta is a scalar parameter time, space, whatever, dou by dou x, dou by dou t, whatever you are looking at. So this is, it could be time, it could be space depending upon what you are looking at. I am looking at a region. The region is defined using norms. This is using only one norm or this is using only absolute values because eta 0 is the initial point where you are starting. Eta 0 is the initial point. The value of x at eta 0, so this should be actually x eta 0, x minus x eta 0. Then let us see what this theorem says of existence of solutions. Now, this theorem gives a condition under which a unique solution exists. It says that if f function f and its first partial derivative that is dou f by dou x k, dou f by dou f is a function vector. It is partial derivative with ith component of x that is also a vector. This is also a vector. This derivative, partial derivative of the function vector with ith component of x is also a vector. So we are saying that if this function vector and its partial derivative and its partial derivative is continuous function on domain or on the set D that we have defined. If it is continuous then the x tilde here is the initial condition. I am starting with the initial condition and initial value of the independent variable. It also says that the solution is continuous function of the triplet. Detailed statement you can see here. What I want to point here is that how do I judge? See the problem is why what is the importance of this theorem? Importance of this theorem is that I am able to judge about existence and uniqueness of the solution and continuity of the solution without actually having to solve it, solve the problem. I do not want to solve the problem. I do not want to solve the differential equation problem. I just want to look at a differential equation and make a judgment whether solution exists, is it unique and if it is unique then the third thing is it continuous with respect to initial condition, time and spatial coordinate. Initial that is x tilde is my initial condition. Is the solution continuous with respect to x tilde? Is it continuous with respect to t naught? Is it continuous with respect to t? Now of course there is one more aspect which I am not talking about right now whether this solution corresponds to the real system behavior. I am not talking about that right now. I am just saying that in real system what I know that a solution exists. It is continuous with respect to perturbations in the initial condition. If I slightly start a little later with the same initial condition I will get a similar solution. If you think of two people doing experiments on two different pendulum which are identical, they start at different points in time from the same initial condition the solution will not be different if the systems are identical. So I can judge about the existence of the solution just looking at continuity of f and continuity of the first derivative with respect to each of the elements of x. I will just give one example then we will move on to solving the problems. Bifurcations we are talking about the behavior in the neighborhood of some steady state conditions. Even then uniqueness if these conditions are met solution will be unique. Bifurcations you are talking about the steady state behavior under some parametric variation. Do not confuse bifurcation with existence and uniqueness. In bifurcations we are looking at the local behavior of all the solutions in the neighborhood of some steady state point. Each one of them could be a unique solution starting from a unique initial condition. If you start from a initial condition it will get a unique from the same initial condition you will get the same solution. Correct but if you start from that same initial perturbation from the bifurcation point you will get the same solution that is important. See if it is unstable reactor if you perturb from the unstable point in the same way you will get the same solution you cannot get a different solution. See for example this is an unstable system I am trying to balance it on my right. Now this is an unstable operating point this is an inverted pendulum and if I put it like this it is a stable system. So this is an inverted pendulum. If I perturb it slightly the way it will behave if I do the same perturbation every time it will behave in the same manner every time that is what I am worried about. Let us see a brief application of this theorem on a specific example. See and what are these I am writing here in a very mathematical way that there is a region in which the solution region of interest region of interest is there in every physical system. For example if I am working with concentrations, concentrations cannot be negative. So you have to work with in the positive quadrant when you work with concentrations. There is no meaning to the solution or meaning to the differential equation when it goes to the negative quadrant. Same way about temperature pressure all physical variables all these physical variables unless you start talking about perturbation variables perturbation variables can be positive or negative but absolute variables have to be positive. So we have to worry about existence of solution in certain regions of state space where okay let us look at this one specific example okay. So this corresponds to f of x t I am giving you some differential equation right now. I am not worried about right now whether this corresponds to physical real system it is a mathematical example okay. Now I want to apply our theorem without having to solve I want to judge whether the solution will exist for every x and t okay for every initial condition say x naught t naught if I start from that will the solution exist will it be unique will it be continuous with respect to x naught t naught and t okay. Now this theorem asserts I just wanted to understand application of theorem if you are interested in knowing the proof of the theorem you can refer to books on differential equations where the proofs are given but for us as engineers we are just worried about you knowing whether a solution exist for the given problem okay. So first of all we have to see whether each one of them is it a continuous function is this a continuous function of x 1, x 2, x 3 and t if you examine all of them these are continuous functions of x 1, x 2, x 3 and t. So first condition is that f of x theta should be continuous function of the independent variable the second condition is I have to worry about three vectors what are the three vectors this is one is dou f by dou x 1, dou f by dou x 1 this is the first vector okay dou f by dou x 1 turns out to be 0 cos t and 1 dou f by dou x 2 there are three different vectors which I should look at this is my f of x t dou f by dou x 1 dou f by dou x 2 dou f by dou x 3 if this partial derivatives are continuous functions or continuous function vectors you can see for examining these function vectors you can see that these are continuous functions since these are continuous functions we are guaranteed that a unique solution exist for this particular set of differential equations starting from any initial condition any time t naught okay moreover if you perturb initial conditions by a small amount the solution will be perturbed by a small amount this we know from analysis of this particular f of x t okay an important aspect I am going to leave it at this point just the idea of doing this or sensitizing you about this fundamental issue you might be solving a problem which does not have a solution okay or in some region solution may not exist because in some region continuity of this might be lost okay so you should at least be aware of that one can examine without actually having to solve whether solution exist is it unique and will small changes in the initial condition will it lead to small changes in the solutions that you can just decide looking at the continuity of these vectors okay now let us move on to solving this let us say that we have a certain for the given problem solution exist it is continuous and all that so the so now we move on to the algorithms so before I start with the algorithms I want to talk about some basic concepts which are which will be using throat so first thing is first thing I would call is marching okay now what I am going to do is for the sake of you know a notational simplicity I am going to just just go back here and not worry about this the integration with respect to the spatial coordinates I am going to call this as T and I am going to say at time equal to 0 but the same things can be worked out when you are starting integration with spatial coordinate with you know space part equal to 0 so the spatial coordinate is at 0 so that is that is everything that I am going to say about time in the context of integration over a spatial coordinate same things can be applied okay now marching in time if you want to read it every time as marching in the independent variable is fine okay now I want to solve this problem typically I want to solve this problem from over a interval t belonging to some time 0 to some tf see I want to see how a reactor concentration profile is when I give a step change in the feed flow rate or feed concentration whatever variable of interest that is to me over a period of time from time 0 to next 30 minutes I have this model differential equation model I am playing with it it's like a toy with me I give a change in the input I record I want to integrate and find out concentration profile temperature profile as a function of time that's what is the assignment which I have given to you as computing assignment okay that's supposed to find it as a function of time so actually it's a entire function that I want to find out when you do this numerically okay it is not possible often to solve this problem over the entire interval that you intend to solve I might be wanting to find out the trajectory for next one hour okay but when I integrate okay I don't solve one initial value problem which is starting from initial condition at 0 and then entire trajectory over the entire interval what I do is I divide sub divide say I want to reach here so this is my tf okay I sub divide this into say t1 t2 t3 what I do is I sub divide this interval from 0 to tf into smaller intervals okay and I go marching in time so which means I integrate from time 0 to time t1 okay I solve one initial value problem with this as my initial condition at time 0 and I reach only up to time t1 this t1 could be for example I want to integrate the differential equation over one hour but I make a small step from 0 to one minute okay then the condition at the end of one minute is the initial condition for the next problem so I march from here to here one minute to two minute two minute to three minute and likewise I go on hopping in time okay now it is not necessary that these steps should be uniform they can be non-uniform depends upon depends upon the problem at hand so I go marching in time okay and instead of solving one initial value problem I solve multiple initial value problems one after another sequenced in time okay the end condition the final condition of the first initial value problem will be the initial condition for the next initial value problem okay and so on this is how even the problem which you are supposed to solve in the part of the assignment computing assignment is to be solved okay so you go marching in marching in time nine as I said the way you do this is t0 is equal to 0 is less than t1 is less than t2 is okay the difference between this that is t i minus t i minus 1 let us call this h i this we call as integration interval this we call as integration interval so we one of the major problems in solving ordinary differential equation initial value problem is how do I choose integration interval appropriately okay we will talk about this in some detail at a later point what is the basis for choose in integration intervals but remember that you are never going to solve the problem in one shot you are going to solve the problem by marching in time by subdividing into sequence of multiple initial value problems okay and the final value of the last problem is the initial value for the next problem okay that's how you are going to solve it so it's like saying that if you want to go from here to the gate you are never going to jump in one shot from here to the gate you are going to go in steps steps could be variable okay depending upon whether it's a you know up the slope or down the slope you know the pace could be variable you might be running you might be sometimes walking but you are never going to go in one shot okay we will go in steps and the difference between any two step will call it as integration step size the first important concept a notation that I am going to use throughout before I move on so this let's say function f evaluated at t i this function vector f evaluated at x at t i and time t i I am going to denote this as f i okay this is neither superscript nor subscript this is a vector it can have a superscript it can have a subscript I am talking about a time index okay here when I write f i which means function vector f evaluated at time t i using x t i and time t i okay I am going to use this shorthand notation f i for this okay similarly I am going to use a shorthand notation x t i as x i that means value of the vector x at time t i instead of writing every time t i t i makes my notations very very complex I am going to use this okay similarly whenever I get a Jacobian matrix I am going to say doh f by doh x evaluated at x t i t i when this matrix appears if it is evaluated at time t i I am going to just call this as doh f by doh x subscript i okay subscript i will indicate that it has been evaluated at time t i okay this is to simplify the notation as we go along because notation can become very very complex when you are trying to solve only initial value problem okay this is a notation please keep this in mind this is what you are always going to look at the second important concept that I want to mention here the second important concept that I want to mention here is explicit algorithm and implicit algorithm there are two classes of algorithms which we are going to look at some of them are called as explicit some of them are classified as implicit right now I am not going to derive after sometime we will derive algorithms some of them will turn out to be explicit some of them will turn out to be implicit right now I am not going to derive the algorithms while giving the idea I am going to take one simple algorithm which all of you probably know from your undergraduate mathematics is Euler's method I am going to illustrate what is an implicit Euler method and explicit Euler method this is just introducing the terminology one is the marching in time second is explicit and implicit this I need when I go along so second basic concept that I wanted to know is okay now I have this differential equation this look the differential equation I have this differential equation let us say I have done some integration and I have reached a point a time tn now tn need not be the last point tn is some intermediate point okay I want to go to some some some other final condition tn is some point where I have reached I want to go from I want to solve a initial problem which is over the interval t belongs to tn tn plus 1 okay see I have started from I have started from time 0 t1 t2 I have come to point tn I have come to point tn okay and then my problem is to go to tn plus 1 I have reached here I want to go to tn plus 1 so what I want to do is I have reached up to point tn and then I want to go to point tn plus 1 so I have broken down my initial value problem into sequence of initial value problems I have somehow solved up to point tn I want to go from tn to tn plus 1 okay the simplest algorithm you know from your undergraduate is Euler's method let us go to Euler's method I can I can approximate the left hand side okay using x tn plus 1 minus x tn this is an approximation of this is an approximation of the left hand side okay this is an approximation of the left hand side now let us take a simplified case where h is constant between the difference between let us take we are taking equal equal size steps okay so left hand side can be written as x n plus 1 minus xn divided by h with our new notation the left hand side can be approximated so I would say that dx by dt is approximated as dx by dt is approximated like this okay now comes the question of how do you first of all remember we are only solving this problem approximately in most of the cases when f of x is non-linear it is not possible to solve the problem exactly analytically the true solution is more than or often times is is not known so now question is how do I approximate the right hand side this is a function of x and t okay now the question arises is whether I should use value of x see this is my initial point I know what is xn here I know what is xn here okay at this point at this point I can evaluate the function vector f because xn is given to me it is known to me time tn is known to me okay but if you look at a differential equation carefully actually it requires the function derivative or the derivative vector that is f of x and t at each point in this isn't it when you integrate actually you should know at each point but we do not know what is the future value see you are currently at time tl you are advancing in time we do not know the future value so you make a simple approximation that the derivative over this interval can be approximated okay one one simple approximation is f x t over the entire interval is approximately equal to f x n tn you take the initial value find the derivative at the initial point local derivative and if you make this approximation okay now this is this we call as fn okay this we call as fn then with this approximation I can write explicit Euler algorithm this is xn plus 1 just see I am just rearranging the equations xn plus 1 is equal to xn plus h times you agree with me I have just combined these two approximations and I write this algorithm what it says is that vector x at time point tn plus 1 is function of xn and a derivative at computed with respect to xn okay this algorithm is called as explicit Euler algorithm because on the right hand side everything is known to me xn is initial value known to me I have arrived at xn by some some means and now I am getting new value xn plus 1 when you go from xn plus 1 to xn plus 2 what will happen xn plus 2 will be function of xn plus 1 xn plus 1 would be known to you so you start from n 0 see if you start from what will happen if you start from x 0 x1 will be known to you after you implement this step then when you go to x2 x1 is known so go to x2 and so on so you go on marching in time okay very very simple way of implementing the algorithm simple somebody might say well I do not really agree that you should take this approximation how can you say that the derivative in the time is going to remain constant okay so he would come up with an approximation x okay you would say no no I should take the derivative at the end point not at the beginning point okay well third person might say well no that is not the right thing you should take average of initial point and the end point okay now let us take that view I want to take average of this plus half let us say this is a better approximation okay if this is a better approximation than just taking the derivative at initial point what is my algorithm becomes xn plus 1 is equal to xn plus h by 2 f xn Tn plus xn plus 1 Tn plus 1 what is the trouble with this algorithm trouble with this algorithm is that xn plus 1 appears on the left hand side and on the right hand side how will you solve this problem in general f is a non-linear f is a non-linear operator you take a reactor or you take a distillation column f will be a non-linear operator how will you solve this problem I want to hear louder Newton Raphson or successes of solutions or optimization or it's a non-linear see now once you have developed this algorithm this equation is a non-linear algebraic equation with what is unknown xn plus 1 xn plus 1 is unknown xn is known to you so this is computable this part is computable what is not computable is this because xn plus 1 is not known to me I am marching in time I don't know what is the future value so how to guess my future value okay I have to guess my future value and do iterations to solve this problem I have to solve this problem iteratively and this particular algorithm is an implicit algorithm it is not an explicit algorithm here you have an explicit solution of xn plus 1 in terms of xn here you don't have an explicit solution so this is categorized as implicit while this particular algorithm is categorized as explicit algorithm okay this is implicit algorithm this is explicit algorithm well this is easier to solve we will see after some time does it not doesn't give you great solutions it's less accurate you have to use very very small h this one is difficult to solve every time you have to solve a non-linear algebraic equation gives great results numerical stability is excellent okay so what is difficult to solve actually gives you dividends it gives you more accurate better results easy to solve you have to be very very careful in choosing h I am not saying that you cannot use this unless you use h to be very very small this explicit order will not work here so well here the calculations are more because you have to choose small h your calculations are more because you have to do iterative calculations you get non-linear algebraic equations so inside solving so let's look at this I am solving an ordinary I am ordinary differential equation subject to initial conditions okay I am marching in time so one giant ODE IVP has been converted into sequence of ODE IVPs each one of them inside requires iterations a non-linear algebraic equation to be solved look at the complexity of calculations okay this is better than this so we tend to use this rather than this okay or you can use this no no doubt but you have to choose h to be very very small so you have to understand the differences and then use which one to use when and so on well of course get into those details when you use which one why and so on so this basic terminology will get on to the algorithm development from the next class there is one more concept called stiffness of algorithm of stiffness of differential equations but will come to stiffness of differential equations little later with these two are initially good enough to get started will develop these algorithms numerical methods of integrating differential equations and then we will move on to you know analysis part when is one better than the other and so