 In the previous class, I have introduced you to numerical methods of solving differential equations. In this particular lecture, we will continue that string and we will do a few examples on psi lab and to illustrate the basic properties of most of the commonly used numerical methods. The examples which I will choose today are examples which we have done before. We have tried to analyze them using some intuition as well as linear analysis in some cases. Today, let us just first revise what we have done in the previous class. In the previous class, we learnt a few numerical methods like Euler method. So, let us just quickly go through this and then begin with our examples. Euler method or is also known as forward Euler method is the first method which is taught to any student studying numerical methods for integration of differential equations. It is a very simple method. In fact, if you have got a differential equation x dot is equal to x comma t this x comma looks a bit big. So, I will just write it down again x of t. In that case, we evaluate the function using this algorithm at discrete time steps spaced by h. So, we have got x k plus 1 is equal to x k. This is x x k plus h into f into x k and that is the function f evaluated at the point x k and t k. Forward Euler method or simply Euler method is having order 1 that is if your numerical solution can be expressed simply as a first order polynomial in t then Euler method will give you exact answers. Typically of course, you know that our responses are not a first order polynomials in t. We have in fact, come across sinusoids and exponential functions which are not first order polynomials. So, Euler method is likely to have some error whenever you integrate most common differential equations. It is also an explicit method in the sense to evaluate x k plus 1. You really do not have to solve any differential equations. You have to just write the if you know the previous value of x k. In one shot, you can get x k plus 1. So, you get essentially an explicit algebraic equation wherein x k plus 1 is written explicitly in terms of x k. Forward Euler is also a single step method. So, you just require one previous the value at the previous time step in order to evaluate the value at the next time step. If you look at a special class of systems that are linear systems, then one thing we can attempt to do since we know the solution of linear equations, we can compare what a numerical method gives us and what is the actual solution of a linear system. Now, we cannot do this in general for non-linear system because you cannot write down the response of a non-linear system in terms of simple functions. In contrast, linear systems for example, x dot is equal to A x. We know it is stable when A is less than 0. If A is a real number, then if A is less than 0, then it is a stable. Euler method, if you apply Euler method, if you apply this algorithm where f of x is nothing but A of A x, in that case Euler method tells you that the system is stable if mod of the absolute value of 1 plus A h is less than 1. So, obviously the Euler method may give these two conditions in fact are not equivalent. If you look at for example, f h is greater than 2 by A, in that case even though A is less than 0, you may still have Euler method showing the system to be unstable. So, Euler method when you integrate the system may give you wrong answers. Even not only wrong answers, it may give a qualitatively wrong result. It may show a stable system to be an unstable one. So, this is basically one problem with Euler method that you will not get accurate answers and you may even get qualitatively wrong answers. Look at backward Euler method. Backward Euler method, the only difference is that the function evaluation is done at the new value of x. It is in fact done at the point which you want to evaluate. So, although backward Euler is also order 1, this is something you can prove. I am not prove this in this class. It is implicit. It is an implicit method in the sense to evaluate x k plus 1. You have to solve this equation x k plus 1 appears within this function. Therefore, you need a numerical method at every time step in order to get x k plus 1. So, this is if for example, for if x is f is a non-linear function, you will have to use some numerical method to solve algebraic equations like Newton-Raphson or Gauss-Seedle. Backward Euler also is a single step method. For a linear system, backward Euler is stable if this condition is satisfied. That is 1 upon the absolute value of 1 minus a h less than 1. This condition also is not equivalent to this condition. In fact, an interesting point which I will show you sometime later in one of the examples is that you can choose values of h. So, that an unstable system is shown as a stable system if you use backward Euler method. So, backward Euler method also can and sometimes depending on the choice of the value of h give you qualitatively wrong answers as far as stability is concerned of a linear system of this kind. Trapezoidal method on the other hand is an order 2 method. The basic algorithm is this. It takes the average of the function evaluation here at the previous point and the new point. It is an implicit method again and it is a single step method. It requires only the previous value of x in order to get the new value of x. For linear systems, trapezoidal rule says that the system is stable if this is less than 1, whereas the actual system is stable when a less than 0. But interestingly, this is a very interesting fact which I leave to you to work out. These two conditions are equivalent in the sense that if a less than 0, it also means that absolute value of 1 plus h by 2 divided by absolute value of 1 minus h by 2 is less than 1. So, that is an interesting thing. Trapezoidal rule whenever you are simulating linear systems will never give you a wrong result as far as stability is concerned. So, that is something interesting and useful. But of course, one important point is implicit methods like backward Euler and forward Euler are difficult to solve because especially when you come to non-linear systems because you will need to use a numerical method to compute x k plus 1 from x k. It is an iterative method at every time step itself you will have to iterate and solve a numerical or solve a non-linear algebraic equation. When you are talking of linear systems, of course, it would mean typically a matrix inversion or a solution of a linear system of equations. We also discussed in the previous class not trapezoidal. This method which I am talking of is rangekutta method. Rangekutta method of fourth order has the following evaluations. You have to evaluate k 1 at these points. I have done this in the previous lecture. You evaluate k 1 first, that is the function evaluation. k 2 is a function evaluation at a slightly different point depending on which depends on the first evaluation. So, f of x k plus h by 2 into k 1, k 1 is evaluated here. Similarly, you have got the evaluation of k 3 which requires the value of k 2 and k 4 which requires the value of k 3. So, what you see here is this is small error here. What you see here is the final value of x k is a combination of these. Of course, rangekutta requires a lot many more evaluations and it can be proven that this is a fourth order method. The one I have shown you is a fourth order method. It is of course, single step. The reason is that we require only the value of x k in order to compute all these functions k 1, k 2, k 3, k 4 and it is an implicit method. At least this one which I have shown you is an explicit sorry it is an explicit method. I am sorry rangekutta fourth order which I have shown you here is an explicit method. Explicit in the sense that all these functions evaluations are single shot. They are explicit evaluations of functions. You do not have to solve any non-linear algebraic equations to evaluate any value here. Of course, when you are trying to numerically integrate linear higher order systems, again let me point out here that the reason why I am studying the behavior of numerical methods using linear systems is not because linear systems cannot be solved. We saw in the lectures, previous lectures that linear systems in fact are amenable to exact solutions. You can write down the solution in terms of very well known functions. So, the only reason why I am using linear systems, a numerical integration of linear systems, the study of that is because I would like to know how the numerical methods behave. Now, if you have got a linear higher order system like x dot is equal to a x a is a matrix x is a vector. It is stable, you know it is stable if real of the any of the all eigen values should have the real part less than 0. Then we say it is stable, but if I numerically integrate this higher order linear system, in that case what you have is x k plus 1 is a vector is 1 plus a h times x k, which is again a vector this is a matrix. Now, this system is stable this is a discrete time system, it is stable if any the eigen values of 1 plus a h is less than 1. This is not something I am explicitly proving, but it is not difficult to prove you can use the same model analysis or linear systems analysis eigen values and eigen vectors to analyze discrete time systems as well. Now, if the absolute value of the eigen value of 1 plus a h is less than 1, then it is a stable system that is what Euler method says, the original system is stable if all eigen values have the real part less than 0. So, for a higher order linear system Euler method will give result that if 1 plus lambda h the modulus of lambda h is less than 1, then it is stable. So, where lambda is the eigen value of a. So, let me just repeat if the eigen values of 1 plus a h have the modulus or absolute value less than 1, then Euler method will show that the system is stable. It may be actually a stable system and Euler method may show to be unstable system if this becomes greater than 1. So, it really depends on the value of h here. You can show that mu is nothing but 1 plus lambda h, where lambda are the eigen values of a. So, the condition this condition can be rewritten as the condition for this discrete time system to be stable is 1 plus lambda h the absolute value should be less than 1. So, again let me just clarify what I am trying to say I am trying to say that even if the original system is stable if this condition is not satisfied then your numerical method will give you a qualitatively wrong answer about stability. So, you have to be careful when you do the simulation that you do not get wrong answers. Now, moving on we will do a few examples in psi lab it would be wonderful if you have got psi lab or mat lab or even if you program in C or Fortran and verify what I am trying to say here. What I will do is I will kind of make an algorithm for doing the numerical integration first of a system x dot is equal to minus 5 x. This is a stable system. So, I will first numerically integrate it using Euler method actually if Euler method is used remember that it will you know it is likely to give correct results one of the intuitive ways of choosing your h is that look at this eigen value this eigen value will probably will give in fact a solution of this kind. So, we will just try to verify whether we get that solution or not. In fact, you have to choose a h carefully if you look at this particular you know function it will decay almost to 0 in 5 times the time constant of the system. The time constant of the system is in this particular case 1 upon 5 that is 0.2 the time constant of the system. So, we expect that it will decay this particular thing will decay if I give a non equilibrium initial condition you will go to 0 which is the equilibrium in this case in about 1 about 1 second. So, I can just intuitively think that if I want to get a nice at least a nice picture of what is going to happen I should choose h at least you know something like 0.01 or so we will just try out with various values of h and see what we get. I done this in the previous class also, but it was somewhat in a hurried fashion let us do it a bit calmly in this particular lecture. So, this is the basic psi lab window I have already invoked psi lab what we will do is first of all let us look at the solution of. So, we will do a numerical integration using Euler method first of all we will clear this variable x we will choose a time constant let us say of time step of 0.01. Let us simulate or numerically integrate this system till three seconds because we know that it is probably to decay within that time. So, we have got a number of discrete time steps at which x has to be evaluated. So, let us assume the initial condition on x is 1. In fact, I should have written this x of 0, but psi lab does not allow 0 as an index. So, that is why I have called it 1. So, this is the value in fact at time t is equal to 0. So, I have given initial condition of 1 remember x dot is equal to minus 5 x has got an equilibrium value of 0. So, we are giving an initial condition which is not the equilibrium value. So, you will see a transient the value of a is 5 as I discussed some time back. So, what I will do is I will evaluate the behavior of the system using Euler method. So, x of i is equal to 1 plus a into h into x of i minus 1 this is what Euler method will do. So, if you recall we have got x k plus 1 minus x k upon h is equal to minus 5 of x k. So, this will effectively give us x k plus 1 is equal to 1 minus 5 h into x k. In fact, minus 5 h is nothing but plus a h this is nothing but this is 1 plus a h a itself is negative that is why you are getting 1 minus 5 h. So, what I will do is I will evaluate this and plot it. So, the way to do it is. So, we evaluate this it is giving a growing function. Now, why is this happening we will just look back perhaps I have done I have made some error here h is equal to 0.01 that is I am not saved it. So, I need to save it first I have saved it. So, I have given the value of h is equal to 0.01 which I changed from the earlier value and now I will run it again. So, let me just run this again. So, I am doing this of course, simulation for 3 seconds. So, now what we are getting is what we kind of expect exponentially decaying e raise to minus 5 t kind of response the initial value is 1 in case you are not able to see it clearly and this time period is from 0 to 3 seconds. So, what I will do next is I will change the time step to 0.1 and just to make it a bit easier to make see we save it plot it as a red. So, if you look at this you see that by changing the time step to 0.1 we have actually slightly the behavior is slightly change you see this red if you look at this red response it is slightly different. So, if I go ahead and let me let me make this 0.5 look at the response well you look what has happened your response is absolutely different from what is expected. So, if you choose a time step like 0.5 in which is in fact greater than 2 times 1 upon 5 I had mentioned that if h is greater than 2 by a the absolute value of 2 by a in that case your response is altogether different you really have a situation where you are getting an unstable response is not of course, the correct response we say the numerical integration is blowing up. So, though the actual system is stable the numerical method gives you a wrong answer. So, obviously, you should choose your time step carefully especially if you are using Euler method. Let us now go on and try to simulate the same system using backward Euler method. So, let us just have the same thing this is backward Euler method h is equal to 0.01 t final is 3 number of steps of course, is the depends on the duration as well as h the initial value of x is 1 a is minus 5. So, I just save this. So, as per backward Euler method you will have x of k or x of i is equal to x i minus 1 divided by 1 minus a h you can easily derive it from the basic formula for discretizing or you know using backward Euler method. Now, this particular sentence has been commented. So, I am not actually going to implement this particular you know this particular command what I am going to implement is this algorithm x of k is equal to x of k minus 1 divided by 1 plus a h. So, let us just numerically integrate and see what we get this is with a time step I believe of. So, we are getting the correct answer of course, the time step is 0.01 and a is minus 0.5. Now, I can make this also 5 we have done that before let us see what happens well the solution which you get is this the interesting thing is that of course, there is an a kind of error you know the earlier with a lower time step we got this response. Now, we are getting this slightly erroneous, but of course, it is still stable. So, what one small point which we should remember is that backward Euler a stable system will in fact, be shown as a stable system. So, that is no worry, but a very interesting thing is what if I got an unstable system a is equal to 5 and I choose my time step more than 2 times 1 upon 5 and interesting thing which you will see here is that this system which ought to be unstable you know after all a is equal to plus 5 this is the real system is unstable e raise to 5 t into x of 0 that is the solution what you will find is backward Euler when you use backward Euler it is actually showing the system to be a stable one that is very interesting. So, it really this something you should keep in mind that backward Euler in fact, is gives you can also give you wrong qualitative answers if your time steps are large this of course, not true if you choose a time step to be small. So, for example, if I choose 0.1 backward Euler of course, should give you the correct answer correct in the sense that it is an unstable system. So, you see that unstable system. So, it just grows up and if it is not visible on your screen this is almost 1.2 into e raise to 9. So, 1.2 into 10 raise to 9. So, that really means that it is really blown out in 3 seconds itself that is not very surprising. So, this is about backward Euler what if I do the same system simulation using trapezoidal rule trapezoidal rule has an interesting property that whether you are stable if you have a stable system it will show be shown as stable. So, if for example, this is for 0.01 just to say I will just check whether I have done the right thing I have to comment this statement and enable this statement which is essentially implementing trapezoidal rule this is x of k is equal to x of this x of i is equal to x of i minus 1 divided by 1 minus h by 2 into 1 plus a h by 2. So, you can actually work out this that you get this particular algorithm if you apply trapezoidal rule to x dot is equal to a x. So, if you have got a is equal to minus 5 and h is equal to 0.1 as we saw just now I save this run it remember although I have called this file backward Euler still it is doing trapezoidal it is using trapezoidal method of integration now. So, what if I increase the time step. So, I am using trapezoidal rule just remember that. So, for example, 5 this was unstable for example, Euler method when we try to integrate it with h is equal to 5 it was unstable what does trapezoidal rule is also stable of course, the errors in the response. So, if I actually go on increasing h the errors increase, but as we shall see I will just do something more sillier and make h is equal to 1 second even if I make h is equal to 1 second although I am going to get an error in the response it is not unstable this is my new response. Similarly, if a is equal to plus 5 if I choose h is equal to 1 this is an unstable system we will just do it on a fresh plot this is an unstable system and trapezoidal rule is showing it to be an unstable system. So, really trapezoidal rule you know at least in a linear system it kind of preserves the stability information. So, if I of course, if I make this time step smaller the answer is of course, closer to the right answer. So, you see this increasing here like this. So, to summarize whatever I have shown so far is that when we are trying to numerically integrate a system x dot is equal to minus 5 x. If you choose Euler method if you choose a h to be large in fact, if it is greater than the absolute value of 2 by a is equal to minus 5 in that case you in fact, get a wrong qualitative picture about the stability of course, as you go on increasing the time step a numerical methods accuracy versions. So, it is not unexpected that you will get better results if you choose h smaller, but the important point which I again emphasize is that Euler method may give you a qualitatively wrong result also it may show a stable system in fact, to be unstable if h is larger this does not happen with backward Euler or with trapezoidal rule. In fact, a stable system is always shown to be a stable system whether one uses backward Euler or trapezoidal rule, but interestingly with backward Euler method and essentially unstable system can be shown as a stable system if your h is very large. So, if h is greater than point 2 upon a the modulus of it and unstable system for example, a is equal to plus 5 is shown to be stable one whereas trapezoidal rule preserves the stability information of a linear system. So, that is one thing a very interesting point which we can we must note now the basic you know the three method which I have shown you to get accurate solutions in fact, for example, x dot is equal to minus 5 x it makes sense in fact, to use methods with time steps of 0.01 or so you will get a reasonably accurate picture you know it is another thing about you know if you choose h is equal to 0.5 etcetera not only you are going to get an inaccurate picture, but in case of Euler method you will get even a qualitatively wrong picture. So, you will use a it is a good idea if you have got a system with a time constant of 0.2 or so as in this case to use a time step of 0.01 and so on, but you know one you know for example, if your situations where a system is marginally stable you know you have got a system in which if the real part of the Eigen value of the continuous time system is 0 or near 0 Euler method is absolutely hopeless you will find that it does not give you a good you know it gives you a qualitatively wrong picture. What I will do now is try to simulate one of the systems which we have simulated before try to understand qualitatively before and the swing equations of a power system. So, the swing equations of a single machine infinite bus system are delta dot is equal to omega minus omega naught and the derivative of omega minus omega naught which is nothing but omega dot itself, because omega naught is a constant which is the frequency of the infinite bus is equal to omega b by 2 h p m minus k sin delta k is in fact we assume it to be a constant this is the internal voltage of the synchronous machine the infinite bus voltage divided by the reactance. So, if I got a single machine connected to an infinite bus e 2 angle 0 and this is e 1 well this is a synchronous machine which is modeled as an internal. So, this is our synchronous machine. So, this is e 1 angle delta this is the internal voltage this is the rotor angle x dash is a transient reactance. So, this these are the equations swing equations this is a non-linear system of equations. So, if I am trying to solve this system of equations by Euler method what is the kind of response I will get. We saw for small disturbances we have done the analysis of the system and we saw that in fact it is for the equilibrium corresponding to which is less than 90 degrees this one for small disturbances around the equilibrium we will get essentially an oscillatory response. So, we know that this if you are at this equilibrium you will get an oscillatory response. If I give of course, a very large disturbance from this then the system may even go unstable the synchronous machine may fall out of step or lose synchronism. So, we know the kind of responses which can be expected at least for small disturbances and qualitatively of course, we know that for large disturbances the system may go unstable. So, with this this is a non-linear system by the way this is not a linear system because your delta appears within the sine function. So, what I will do is try to integrate the swing equations by Euler's method. So, you just take let us assume P m is 1 h is 3 at this capital H is of course, the inertia constant k let us assume is e 1 e 2 divided by 0.4 this is omega and omega base we assume is the same as this frequency of the infinite bus. Let us take a time step of 0.01 the final time is 3 the number of steps of course, can be calculated from t final n h. Let us assume that delta is equal to sine inverse sine inverse P m by k in fact, this is the equilibrium value of delta omega is nothing omega let us assume is omega b plus 1. Now, remember one important point is that although delta is at the equilibrium value of the value of delta is actually its equilibrium value the value of omega is not its equilibrium value. So, as a result this we are not actually at the equilibrium because both states have to be at the equilibrium value if we are not to have any transient. So, by adding this plus 1 we have in fact, given an initial condition which is not the equilibrium value of the states. So, we are going to get a transient. So, I have given you a initial condition which is not equal to the equilibrium value if I had given you the equilibrium value the system would simply not move it would just remain where it is. So, if I apply Euler method you have got delta is equal to delta of i is equal to delta of i minus 1 plus h this is in fact, simply omega minus. So, we are just numerically integrating the system by Euler method omega i is equal to omega i minus 1 plus h times omega base divided by 2 times the inertia constant into P m minus k sine delta this easy to evaluate if I know what the initial values of delta and omega are I can find out the new values. Since, we are not at equilibrium since I have given an initial condition which is not the equilibrium condition of course, if I had made the 0 this 1 is 0 we would be at equilibrium. So, what I will do is I will just simulate the system. So, what I will do is again. So, I am evaluating this system using Euler method. So, if you look at the behavior of a system it looks as if it is an growing oscillation. In fact, if you look at we have done the small disturbance stability of this particular system and we know that the response is in fact, not a growing oscillation, but a steady oscillation if this disturbance is small. So, obviously, there is a problem in the method let us just try reducing the time step. What is the problem in reducing the time step? There is no problem you will have to wait longer for the simulation to get over and sometimes that becomes a bit that may become a limiting factor for doing any kind of study if you make the time step too small. So, now you are getting you know a different answer that shows that you have to actually reduce the time step. So, I will just redo this in fact, by using the smaller time constant time step we are in fact, getting a better answer, but you see that this is still growing it is not it is not you know you this value is slightly less than the value you have here. So, what we see is that if you have got you know marginally damped systems Euler method in fact, does not seem to be very nice because whatever I reduce it if I reduce the time step I can do that even further. Now, it is going to be painful because you will have to wait for a longer time now it is still it will take some time for this to come up. So, you see that there is slight difference in the what I got with the previous time step and the time step you are having now. One thing we can simulate for longer and longer time with this time step there is point triple 0 1, but the point you can actually try it out whatever time step you take you will find this eventually starts growing if you try to simulate beyond 3 seconds say 20 seconds even with this time step you will find that it is growing. So, Euler forward Euler method not very nice to simulate systems with which are marginally damped or have poor damping. So, let us just rerun the system using trapezoidal rule. I am using exactly the same system with the time step of point 0 1 and one problem comes here. Now, if I am going to try to do integrate this method using trapezoidal rule what is going to be the problem the problem which you see here of course, is that the statements in the program seem to have become much larger. One of the things you will notice is that delta i will not be easily obtained from the previous value that is delta i minus 1. So, let us just pay our attention on what we get if we try to discretize swing equations using trapezoidal rule. So, you have got delta dot is equal to omega minus omega naught. So, we will have delta k plus 1 is equal to delta k plus h by 2 times omega k plus 1 minus omega naught plus omega k minus omega naught. So, of course, we require the value of k plus 1 here omega k plus 1 will be equal to omega k plus h by 2 into omega b by 2 h is the into p m minus k sin k minus k sin k minus k sin delta k plus 1 plus p m we assuming of course, p m is a constant sin delta k. So, this is the discretization which you will have to use and that is what is really been programmed there. Now, the problem here of course, is that you do not know the value of delta k plus 1 you do not know what omega k plus 1 is, but you require that in order to evaluate delta k plus 1. So, if I if I given you delta k and omega k how do I get omega k plus 1 and delta k plus 1 you will have to use an iterative method. So, that is what exactly I have done in the program what I have done is try to use within every time step within every time step solve this using Gauss-Seidel method. So, what we do is get plug in this value get plug in this value of delta k and omega k you will plug it in here here here and here you plug it in omega k plus 1 is something you do not know. So, let us guess that it is equal to omega k in the first instance and this is delta k. So, you run this you you evaluate this assuming this is omega k and this is delta k you will get a new value of delta k plus 1 and omega k plus 1 plug in this here and you plug you plug in this here you plug in this here and re evaluate delta k plus 1 and omega k plus 1. Basically, we are following Gauss-Seidel method of solving this algebraic equation set of algebraic equation you go on doing this till there is no difference between a previously evaluated delta delta k plus 1 and the newest value of delta k plus 1 which you are getting. So, that is exactly what I am doing this here the program is a bit complicated and probably you will find it difficult to follow the steps, but I will just indicate what I am doing. So, for what I do here is I essentially in these steps I am in fact using Gauss-Seidel method using Gauss-Seidel method I have taken a maximum number of iterations per time step as 10 with the hope that of course, that will converge. So, in fact I am checking out the convergence criterion as well. So, I will just let me read go through the whole program slowly this is the check on the convergence. So, I will just run this of course, one small point will give the same disturbances before in the sense that the this was one in the previous. So, I have given a small disturbance we are away from the equilibrium save this. So, this is basically what you get we actually evaluated this we redo it this is not the what we get with trapezoidal rule. Now, we will get what we see. So, this is the response and actually this sinusoid does not grow with time. So, this is a this is a true in fact really reasonably accurate picture of how delta varies if the system is given a small disturbance. So, even if I increase the time step remember here this is this will not really be affected. So, I just change my time step to 0.1 10 times more save this rerun it because of course I have chosen a larger time step there is slight you know kind of pokey pokey solution here this particular sinusoid has got rough edges. But on the whole it does not give you a wrong picture of stability it is a undamped sinusoid delta is an undamped sinusoid. So, this is basically what you get with trapezoidal rule. If I increase the disturbance I increase the disturbance in that case you may even go unstable. So, I am now I am going to give a large disturbance relatively wrong you see that because of this larger disturbance the delta oscillations have become larger and one interesting thing that this is no longer sinusoidal looking you see this top part looks much fatter than the bottom part. So, obviously with the if you make the disturbance larger your behavior starts moving into the non-linear zone and you do not get the sinusoidal response as predicted by linear analysis. So, if I actually make this disturbance even larger you are going to get the system loses synchronism you see that this angle just goes on increasing with time. So, if I go on increasing the you know deviation from the equilibrium we go into non-linear behavior and this is of course, captured quite well by this numerical integration method. So, this is basically what I wish to tell you about using trapezoidal rule for small disturbances it gives you the correct response if I use backward Euler method for the same problem backward Euler method is also an implicit method remember. So, I have to use it an iterative solution to get this iterative method to get the values of the states at every time step. So, let us just do this again with backward Euler. So, so an interesting thing that I have just to observe here is that when I use backward Euler on a marginally stable system which is just having zero damping Euler backward Euler shows it as if it is a stable system as a damped oscillation. So, that is an interesting point about backward Euler it really changes the picture means especially if your system is marginally damped it is got very little damping. So, it is showing you the system to be a stable system whereas, actually it is an undamped oscillation. So, the system is a showing as if it is dying down. So, that is one interesting observation about backward Euler method of course, if I use backward Euler method with a very small time step you are going to get reasonably accurate solution. So, I just numerically integrated with a small smaller time step well we are coming the answer is slightly different and we are coming closer to the solution we have in mind. So, of course, reducing time step is one of the ways you can get accurate solution or a more accurate solution. So, with this I end my discussion about the numerical integration of swing equations. My discussion about numerical integration method would be quite incomplete if I do not consider what are known as stiff systems. We have learnt about stiff systems. Stiff systems are systems in which the Eigen values of the of linear systems especially are far apart in the sense that they are small as there is large Eigen values which translates to the fact that the system has got both fast and slow transients. And if you are actually interested only in the slow transients then you really need need not model it in great the fast transients in detail. In fact, I showed to you in one of the previous lectures that when you are trying to analyze a system with slow and fast transients the fast transient differential equations corresponding to the fast transients or the states associated with the fast transients can be converted to algebraic equations and we really do not lose outmatch as far as the information on the slow transients is concerned. A similar issue will kind of crops up when you have got even a numerically integrating stiff system. So, if you have got a stiff system you have got large and small Eigen values and the question comes what is going to be the choice of time step. So, if you have got a stiff system and you have got fast transients you will have to choose a time step which is compatible with the fastest transients of their system in order to accurately get the response. So, if you have got a system with slow as well as fast transients inherent then you have to choose a time step to be compatible with the fastest transients. So, that your fast transients are accurately obtained, but consider a situation where you are not really interested in the fast transients. So, how will you choose a time step? In Euler method it is very clear that if 1 plus lambda h the absolute value is greater than 1 then the numerical integration method is going to blow up. So, even if your fast transients are stable if I choose my time step which is not compatible with this particular condition then your numerical integration as I said will blow up, but what if I am interested only in the slow transients. Unfortunately with Euler method will be forced to take a time step which is small very very small which is compatible with the fast transients. So, that is one of the problems. So, even if you are interested in the slow transients you have to choose very small time steps. So, that is one of the problems which you will encounter when you try to use Euler method to integrate a stiff system. So, let us again look at a stiff system we have done this example before this is I L 2, this is I L 1, this is V c we know that this particular system has both fast and slow transients. So, this is 10 milli Henry, this is 1 Henry. In fact, the state associated with the slow transient is this and the states associated with the fast transient are these from the participation factors which we have evaluated in a couple of lectures before. So, if you have got this particular system this contains slower and faster transients. In fact, how do you know that it contains slow and fast transients. In fact, if you know the A matrix you have got x dot is equal to A x where x is equal to I L 1, I L 2 and V c then this A is equal to minus 10, 0, minus 100, 0, 0, 1. We have done this in a previous class. So, I am not reworking out this this is minus 10000 and 0. So, this is your A the Eigen values of the system are approximately minus 5 plus or minus j lambda 1 and lambda 2. So, they are two actually complex conjugate pair. So, this is a stable system, but one of the Eigen values is very small and the other one is very large. This is a example of a stiff system. If I try to numerically integrate this system, stiff system using Euler method I will have to choose my time step. So, that 1 plus lambda 1 h is less than 1. So, h should be such such that this is satisfied. In fact, it should be satisfied for all the three 1, 2 and 3 all the three Eigen values. So, one thing you will notice that if you want to satisfy it for all the three Eigen values then h will have to be very small. In fact, it becomes completely unstable in case you choose your h say as 0.001 or so. So, I will just do this numerical integration using Euler method and then we will conclude this particular lecture. So, let us just do this one integration. So, this is your system this is the A matrix I chose a very small time step. I will just simulate it for a short while because I know it is really not going to work out with Euler method even with this small time step. Even with this small time step we have a problem with Euler method. So, the first thing I will do is integrate it using Euler method. This is a remember higher order system. So, instead of 1 plus small case a h you have got identity matrix plus a into h into x plus h into u. So, I am implementing Euler method for this system of course, there is an input here. So, I have used Euler method I have generalized it to a situation where you have got an input as well. So, if I numerically integrate this system oops there is a error here. So, I will just check out what is the problem and redo it. So, the problem probably is there is some error here. So, we will what we will do it do is redo this example again in the next lecture. We will just debug the error which is there in the program and rerun this example in the next class. So, in the next class we will continue with this particular example see how it is in fact a very interesting and important thing we are going to learn in the next lecture that for stiff systems of this kind you have to use specific numerical methods even you know explore if you want to capture. In fact, both the fast and slow transition you may even have to explore things like variable time steps and so on. So, let us conclude this lecture at this point we will redo this example which did not work out in this particular lecture again in the next lecture.