 Hi, we are learning numerical methods for initial value problems of a first order ordinary differential equations. In this, we have learnt some explicit methods as well as implicit methods. Explicit methods are easy to implement and they are also efficient in computation because their formulas are written explicitly for the unknowns in terms of all the known quantities. Whereas, implicit methods in general gives us a relation which may be a non-linear equation if our given ODE is a non-linear ODE in which case you may have to go for non-linear iterative method in order to get the approximate solution or you may have to go for a predictor character approach. This is what we have seen in the last few lectures. You may be wondering when we have explicit methods, why are we struggling with the implicit methods because we need to put more efforts in order to get the approximate solution through implicit methods? Well, there is a purpose for implicit methods. This purpose will be very clear if we understand the stability analysis behind these methods. In this lecture, we will try to understand the stability of these methods. In order to keep our discussion very simple, we will do the stability analysis for only forward Euler method and the trapezoidal method. Recall, forward Euler method is an explicit method and trapezoidal method is an implicit method. In order to see the importance of the implicit method, we will just do the stability analysis for these two methods. One can also extend the stability analysis for other methods like rangikutta methods and other multi-step methods that we introduced in our previous lectures, but because of the time constraints and also to keep our discussion simple, we will only cover the stability analysis for these two methods. For stability analysis, we will only consider the linear initial value problem of the form y dash is equal to lambda y post on the interval 0 to infinity with this particular initial condition y of 0 equal to 1. We say that the origin 0 is asymptotically stable if y of x that is the exact solution of this initial value problem tends to 0 as x tends to infinity. Remember, in our initial value problem, we will always take lambda to be non-zero. It may be a complex number or just a real number. You can see that the origin is stable if and only if the real part of the lambda that we have on the right hand side of our equation is strictly less than 0. A good numerical method should capture this decay property of the exact solution. Now, let us go to the forward Euler method and we will try to apply the forward Euler method to our initial value problem y dash equal to lambda y with initial condition as y of 0 equal to 1 and see whether forward Euler method captures the decay property of this initial value problem when lambda is less than 0. Let us keep in mind that the exact solution of this initial value problem is given by y of x equal to e power lambda x. Since we want to study the stability, we will always take lambda to be negative. To be more precise, let us take lambda equal to minus 20 in which case the solution y of x is given by e power minus 20 into x. The graph of the solution is shown in the blue solid line. The star here indicates the point at which we have specified the initial condition. You can see that as x increases, the solution drops down rapidly and tends to 0 as x tends to infinity. That shows that for lambda equal to minus 20, you have the asymptotic stability for the origin. Of course, we have seen that for lambda less than 0 for any value, you have the asymptotic stability because the solution has the exponential decay nature. Here, you can see that the solution drops down from 1 to 0 more rapidly within the interval 0 to almost a 0.3. Therefore, our numerical method should capture this rapid change in the solution. This is one of the aspects of the so called stiff problems. In stiff problems, the solution undergoes a rapid change in a small interval of the independent variable. Therefore, this particular initial value problem exhibits one aspect of the stiff problem. In general, explicit methods struggle to capture the solution when the problem is stiff. Let us see how forward Euler method captures this solution. Recall the forward Euler method for this initial value problem is given by this formula. If you recall, if you have y dash equal to f of x y, then the forward Euler method is given by y j plus 1 equal to y j plus h into f of x j comma y j. Here, we have f of x comma y equal to lambda into y. Therefore, the forward Euler formula is given by this expression. Let us take h is equal to 0.11 and see whether the Euler solution captures this decay property well or not. The approximate solution computed from the Euler method is shown in the red solid line. Here, blue solid line is the exact solution. Recall, the exact solution starts from 1 and it falls to 0 quite rapidly. However, since the y coordinate is shown with the scaling ranging from minus 3000 to 3000, therefore, in this y scaling, the exact solution is seen almost like a straight line just parallel to the x axis. Why it is happening? Because, the Euler solution is furiously oscillating around y equal to 0. In fact, as it oscillates, the amplitude also blows up to infinity. As x tends to infinity, you can see that this is no way near to the exact solution because exact solution starts from 1 at x equal to 0 and decays quite rapidly to 0. Whereas, the Euler solution started correctly from y equal to 1 at x equal to 0, but then it started oscillating and the amplitude of the oscillation keeps on increasing and as x tends to infinity, amplitude tends to infinity as it oscillates. It is quite surprising why Euler method gives such a drastic behavior and such a behavior is what we refer to as instability in the numerical method. When the exact solution is bounded, we at least expect the numerical solution to be also bounded. Whereas, in this case, the exact solution is nicely decaying to 0, but the approximate solution computed from the forward Euler method is blowing up. Let us try to understand why such an instability is observed in the forward Euler method. Let us go back to the Euler method formula. You have y j plus 1 equal to 1 plus h lambda into y j and in this particular initial value problem, it can be written as y j plus 1 is equal to 1 plus h into lambda to the power of j plus 1. How is this happening? It is very clear. You take y 1, it is given by 1 plus h into lambda into y naught, but y naught is equal to 1. Therefore, it is 1 plus h lambda. Now, y 2 is equal to 1 plus h lambda into y 1, which is equal to 1 plus h lambda into y 1 is equal to 1 plus h lambda. Therefore, if you multiply y 1 with 1 plus h lambda, you get y 2 equal to 1 plus h lambda square. Similarly, you have y 3 equal to 1 plus h lambda cube and so on. In general, y j plus 1 is equal to 1 plus h lambda to the power of j plus 1. We have already chosen lambda to be negative. Now, since we are with forward Euler method and we want the solution on the right side of the initial point x naught which is equal to 0, we have to take h to be positive. You have lambda negative and h positive. Now, with this the formula for Euler method which is given like this will show up the decay property only when mod 1 plus h lambda is less than 1 because you have 1 plus h lambda to the power of j plus 1 as j tends to infinity which is equivalent to saying that x tends to infinity y j plus 1 will tend to 0 if and only if 1 plus h lambda is less than 1. So, that is very clear. Now, you want this condition in order to capture the decay property which is there in your exact solution. Let us see what choice of your h will give this condition of course, from this inequality you can immediately see that you cannot achieve this inequality for any choice of h greater than 0, but you have to choose h to be less than minus 2 by lambda. Remember lambda is negative therefore, minus 2 by lambda will be some positive number. So, you have to choose h which is of course, greater than 0, but it should be less than this number. Now, when you take lambda equal to minus 20 you can see that this is equal to 0.1. Therefore, this small stability analysis says that your forward Euler method will give you the decay property which is there in your exact solution if and only if if you choose your h to be some number which is less than 0.1 of course, it should be greater than 0 also whereas, in the previous numerical experiment we have chosen h to be greater than 0.1. Thereby we had mod 1 plus h lambda to be greater than 1 that is why when you raised that to the power of j and as you go on increasing j your y tend to infinity that is what we can clearly see in the numerical solution here. Therefore, in order to have stability in the forward Euler method for this particular initial value problem you have to choose your h to be less than minus 2 by lambda. Let us choose h to be less than minus 2 by lambda which is in our particular choice of lambda it happens to be 0.1. We will now choose h to be less than 0.1 and let us see how the numerical solution looks like. Well, one good improvement is that the numerical solution is not oscillating widely and tending to infinity it is nicely decaying now. So, that one improvement is there in our numerical solution however, you can see that for small values of x it shows some spurious oscillation which is not there in your exact solution. This is purely a numerical noise which you do not want to have in the numerical solution, but Euler method is exhibiting it well if you decrease the value of h this oscillation will tend to decrease. Let us take h equal to 0.025 you can see at least at the points where we plotted the graph of the approximate solution we do not observe any oscillation for this value of h. However, you can see that at some points the accuracy is not so satisfactory you still have a big error that may be improved by taking smaller values of h. Now, I am taking h to be 0.01 you can see that the approximate solution is now better coinciding with the exact solution when compared to the approximate solution obtained for h is equal to 0.025. Let us further decrease the value of h and see at least for h is equal to 0.01 you can see that the red solid line is coinciding more well with the blue solid line. You can see at least visually that the accuracy is now better so from here you can see that to get a satisfactory accuracy like this you have to really choose your h very very small in the forward Euler method remember forward Euler method is an explicit method and it is a first order method. You can speed up the convergence by going for higher order methods like Rangikutta method and all. In fact, you can go for multi step methods in that way you not only gain higher order you also include more points into your formula in that way when you are working with stiff problems multi step methods work very well because it takes the information from more points and build the approximate solutions at the unknown points in that way multi step methods are more suitable for such stiff problems and for stability purpose explicit methods always comes with some extra conditions on h in order to achieve the stability this is what generally people call as conditionally stable methods like Euler method many other explicit methods are only conditionally stable and also they demand very small value of h in order to achieve the stability. Now let us go on with the implicit methods and see what is the stability situation in the implicit methods just for the illustration let us take the trapezoidal method recall that trapezoidal method is an implicit method because it gives an implicit relation for the unknown y j plus 1 and this implicit relation will be a non-linear equation if you are working with a non-linear ODE. However, in the present case we are working with a linear ODE therefore trapezoidal method in fact can be written explicitly and therefore we do not need any non-linear iterative method or we do not need any predictor character approach you can see that this expression can also be written like this and it is well defined for any lambda and h such that mod lambda into h is not equal to 2 right of course we are always taking lambda to be negative and h to be positive with this you can in fact see that the expression within the bracket will always lie between minus 1 and 1 that shows that whatever may be the choice of h as long as h is greater than 0 when you are working with lambda less than 0 you always have the expression within the bracket to be something lying between minus 1 and 1 therefore when you rise that to the power of j and as j tends to infinity you will always see that y j tends to 0 irrespective to any choice of h positive right. So, that shows that trapezoidal rule does not demand any condition in order to achieve the stability that is the good part of trapezoidal method not only trapezoidal method it also holds for many implicit methods that is why implicit methods are preferred especially for stiff problems. Let us see numerically how trapezoidal method works let us take h is equal to 0.11 recall that for this value of h the approximate solution computed from forward Euler method actually blows up as x tends to infinity right. Here you can see that for h greater than 0.1 the trapezoidal solution still decays and tends to 0 as x tends to infinity. However the accuracy is not that satisfactory especially for small values of x but that does not matter when you decrease h little bit you can see that trapezoidal method captures the exact solution more accurately if you recall Euler method captured the exact solution with this accuracy almost for h is equal to 0.001 which is very small when compared to the h that is given for the trapezoidal method right. Of course, you can prove that trapezoidal method is of order 2 and also it is an implicit method. So, trapezoidal method has two advantages one is it is an implicit method it is unconditionally stable second thing is it is a second order method. Therefore, even for a relatively bigger value of h you may get more accurate result when compared to the first order methods. Let us formally define absolute stable methods for the initial value problems of the form that we considered here. A numerical method for the above initial value problem possesses the property that y j tends to 0 as j tends to infinity for all h greater than 0. So, there should not be any restriction on h in order to achieve this property for any real part of lambda is less than 0. Remember if you have this property then the exact solution of this problem will have exponential d k and a good numerical method should capture that property in its approximate solution. And we say that a method is absolutely stable if this property is achieved for any choice of h positive. So, that is what is highly desirable for numerical methods if you want this property. In fact, the message from our previous example is that you have to go for an implicit method. If you go for an explicit method like Euler method they all come up with some condition under which you have this stability. Therefore, they are not absolutely stable whereas, implicit methods are often absolutely stable particularly we have seen that trapezoidal method is absolutely stable. In fact, there is one more important method which is absolutely stable is the backward Euler's method. If you recall in the very first lecture of this chapter we have introduced backward Euler method, but there we have used backward Euler method as an explicit method because we used it to compute the solution at the points on the left side of the initial condition. In that case the backward Euler method becomes an explicit method, but now what we are doing is we are using the backward Euler method in a different way. We are using it to compute the solution of the initial value problem at the points which are lying on the right side of the initial point x naught which is taken in our problem as 0. And we are interested in obtaining the solution for x greater than 0. Therefore, if you try to apply the backward Euler method for this problem you will see that the backward Euler method will also be an implicit method. More precisely for this initial value problem the backward Euler method is given by this expression and from here you can reduce this formula to y j plus 1 equal to y j plus divided by 1 minus h into lambda. You see lambda is negative and h is positive therefore, 1 minus h lambda will always be greater than 1 and therefore, 1 by 1 minus h lambda will be always less than 1. Again you see that backward Euler method is absolutely stable because you see that backward Euler method for this particular initial value problem is written as 1 by 1 minus h lambda to the power of j plus 1. And that tends to 0 as j tends to infinity irrespective to whatever may be the choice of h positive. Remember we always discuss the stability under the condition that lambda is less than 0. Therefore, backward Euler method is also absolutely stable. With this our discussion on stability is analysis for initial value problem is complete. From here you can clearly see the importance of using implicit methods when compared to the explicit methods. Before ending this lecture we will also spend some time on implementing the Euler's method for an initial value problem. Let us first recall the forward Euler method. The forward Euler method is applied to an initial value problem of the form y dash equal to f of x comma y post on the interval a comma b with the initial condition as y of a equal to y naught. The forward Euler formula is given by y j plus 1 equal to y j plus h into f of x j comma y j. And it runs for j equal to 0, 1, 2 and so on. We have to also give the value of h as the input. Apart from h we have to also give the expression for the right hand side function f as a function of x and y. Let us see how to code this method using python. As usual the first two lines are to import the libraries for plotting graph and also the numpy. The first line of the input is about taking the expression for the right hand side function f. Here I want to take the right hand side function as simply cos x. Therefore, I have np cos x because cos is not a function inbuilt in python, but it is inbuilt in the module numpy. That is why we have np dot cos x. Of course, I have given two parameters here x and y, but in this particular function f I have only x dependency and y dependency is not there, but still I have y because if you wish to change this function to something like x square plus y square, you can do that without disturbing the input parameter. That is why I just kept here. Otherwise y has no role for this particular definition of f. Next is the value of h. We are taking the value of h as 0.1 and then we are discretizing the interval a b. Here I am taking a as 0 and b as 6. Therefore, our interest in solving the problem is in the interval 0 to 6 and we are discretizing the interval 0 to 6 uniformly with step size as 0.1. Therefore, your array x will be 0, 0.1, 0.2 like that it goes on till 5.9, 6. The next input is the initial condition at x equal to 0. We are taking y naught to be minus 1. With this we have taken all the necessary inputs for our problem. Now, we are ready to apply the forward Euler method. The first line is to get the information about the length of the array x that will tell us how many grid points that is node points are there in our problem. It is x naught, x 1, x 2 up to x n and then we are initializing the approximate solution array which has n components y naught, y 1, y 2 up to y n. Its length is equal to the length of the array x and just to start with we are initializing the array y with 0 value and then we are putting y of 0 which is the first component of the array y with the initial condition. So, that is how it goes y naught is the initial condition that is what we are assigning here. From y 1 onwards we will compute using the Euler forward formula and that is given by y j plus 1 is equal to y j plus h into f of x j comma y j. That is what precisely we are writing here in the Python notation and that will go for each j is equal to 0, 1, 2 up to n minus 1 in fact. So, that is what is happening here. In fact, I should write it as n then only it will go up to n minus 1 that is the way range works therefore, it should be 0 to n. Then it will go from 0 to n minus 1 in Python and that will make the vector y 1, y 2 and so on up to y n right. So, once we have calculated the array y then we will go to plot the graph of the approximate solution y. You can see that the exact solution of the initial value problem that we have considered is y equal to sin x minus 1. Our interest is also plot the exact solution and compare it with the approximate solution. Remember approximate solution is stored in the array y and we are printing that in the first part of the plot command x comma y and in the same plot command we are also plotting x comma sin x minus 1. Therefore, in this single command we are actually plotting two graphs one is the approximate solution and another one is the exact solution and this is just the legend and with this the Euler code is complete. Let us run the code and see the output. The code is executed and now we can see the output of the code here you can see that the Euler solution is plotted with blue solid line and the red solid line shows the exact solution. You can see that the accuracy is not that good. However, if you reduce the step size here say you go and reduce it to 0.01 you can see that the approximation will be now better as you can see here. Let us also change the function say for instance let us take the function as minus 20 into y if you recall this is what the approximate we have discussed in our stability analysis part and let us take h is equal to 0.11 from our discussion you can see that the Euler forward formula should give unstable solution which is clearly seen here of course, it is not f of x equal to sin x minus 1 just because we asked it to print like this it is printing, but now it is minus 20 into y and you can clearly see that the Euler method shows furious oscillation and also it blows up as you go on increasing the value of x. You can see that the oscillation up to x equal to 6 is between minus 20,000 to 20,000 which is really big whereas, the exact solution is nicely decaying from 1 to 0 as x tends to infinity. From this lecture we can clearly see the importance of the implicit method and we also learn to implement the forward Euler method. Similarly you should also think how to implement the trapezoidal method and rangikutta methods with this note let us end this lecture. Thank you for your attention.