 We are considering stability of various methods. So, single step methods such as Euler's method or Runge-Kutta method will not have any problem about stability, because when we consider an approximate method, initial value problem is a differential equation of order 1. So, that differential equation you replace by difference equation. So, in case of single step methods, you replace differential equation of first order by difference equation of the first order. Whereas, for multi-step method such as say midpoint method, you were differential equation of order 1 is replaced by difference equation of order 2. So, then the difference equation will have solutions and this creates some problem. So, we are going, we are studying those things. So, as a representative case, we are looking for a differential equation looking at the initial value problem y dash is equal to lambda y with y of 0 is equal to 1. So, we know it is exact solution, the exact solution is going to be e raise to lambda x. For this particular differential equation, we want to see what happens when you apply say Euler's method, Runge-Kutta method or midpoint method. So, y dash is equal to lambda y, y at 0 is equal to 1. Let us look at solution over interval 0 to b. The interval 0 to b will be subdivided into n equal parts. So, h is going to be equal to b by n, x n point is going to be equal to n h because our x 0 is equal to 0. In the Euler's method, we have got y n plus 1 is equal to y n plus h times f of x n y. Our f x y is lambda times y. So, that is why you have y n plus h lambda y n. So, that means it is equal to 1 plus h lambda y n. Now, replacing y n by 1 plus h lambda y n is equal to 1 plus h lambda y n minus 1 and so on, you get y n plus 1 is equal to 1 plus h lambda raise to n plus 1. So, this is what I was talking that y dash is equal to lambda y. It is a differential equation of first order. Here you have got y n plus 1 is equal to 1 plus h lambda y n. So, again it is a difference equation of first order. Look at the Range-Kutta method. In the Range-Kutta method, y n plus 1 is given by y n plus k 1 plus k 2 by 2, where k 1 is nothing but h times f of x n y n, k 2 is h times f of x n plus h y n plus k 1. Function f x y is lambda y. So, that is why k 1 becomes h times lambda y n and k 2 will be h, then lambda times y n plus k 1, but y n k 1 is h lambda y n. So, that is why k 2 becomes h lambda into 1 plus h lambda y n. Substitute in this formula. So, we have got y n plus 1 is equal to y n plus h lambda y n plus h square lambda square by 2 y n. So, again it is a difference equation of first order. y n plus 1 is going to be equal to 1 plus h lambda plus h square lambda square by 2 into y n and the solution will be given by this. Now, let us look at the difference or let us look at the error. So, in order to find the error, we look at expansion of e raise to x. So, e raise to x function is equal to 1 plus x plus x square by 2 plus x cube by 3 factorial and so on. Hence, I can say that e raise to x is equal to 1 plus x plus terms of higher order that means terms of the order of x square. If I written three terms, then e raise to x will be 1 plus x plus x square by 2 plus terms of the order of x cube. So, when we want to compare the exact terms of the order solution with the approximation which we have obtained by using Euler's method, we are going to retain two terms e raise to x is equal to 1 plus x. When we want to compare the error in the Runge-Kutta method, we will keep three terms. So, we have e raise to x is equal to 1 plus x plus higher order terms. So, that gives you e raise to n x will be equal to 1 plus x raise to n and then again terms of the order of x square. So, this is what one uses and one gets e raise to n x is 1 plus x raise to n plus term of the order of x square e raise to put x is equal to n h or put x is equal to h lambda. So, you are going to have e raise to n h lambda is equal to 1 plus h lambda raise to n plus term of the order of h square lambda square. This 1 plus h lambda raise to n that is the approximation in the Euler's method. So, that is our y n our x n is n h. So, e raise to n h lambda because this exact solution is e raise to lambda x, this is nothing but y at x n. So, we have y at x n minus y n is equal to term of the order of x square lambda square. So, this is consistent with what we have been saying that the local discretization error in Euler's method is going to be h square. In Runge-Kutta method of order 2 the local discretization error will be of the order of h cube. So, now, let us look at the error in the Runge-Kutta method. So, in case of Runge-Kutta method you have got y n to be 1 plus h lambda plus h square lambda square by 2 whole thing raise to n that is what we obtained y n plus 1 is equal to 1 plus h lambda plus h square lambda square by 2 raise to n plus 1. So, if I consider y n then it is going to be nth power of this bracket. So, thus this is going to be your y n e raise to n h lambda as before it is going to be y at x n. So, y x n minus y n is of the order of h cube lambda cube. So, truncation error to be of the order of h cube as one expect. So, these are the single step methods and now let us look at the midpoint method. We are looking at disinitial value problem. The midpoint method is y n plus 1 is equal to y n minus 1 plus 2 h f of x n y n very much similar to the Euler's method where instead of y n we have got y n minus 1 and instead of h f of x n y n we are saying 2 h f of x n y n and last time we saw that the local discretization error is h cube. So, in the case of Euler's method it is h square for midpoint method it is h cube, but as we will see that there are some problems with this method. Now, look at this formula. So, you have got y n plus 1 is equal to y n minus 1 plus 2 h lambda y n. If I take it on the other side you have got y n plus 1 minus 2 h lambda y n minus y n minus 1 is equal to 0. So, this is the difference equation and in the difference equation y n plus 1 is in terms of y n and y n minus 1. So, this is differential equation of first order, this difference equation of the second order. In order to determine the solution of this uniquely you will need two conditions like in the case of differential equation. If you have got differential equation of first order then the unique solution is determined by one condition. If you have got differential equation of second order then you need two conditions. These two conditions they can be initial conditions or they can be boundary conditions. So, now, we have got a difference equation of second order. So, we will need two conditions. One condition is going to be y 0 is equal to 1 the same as in the case of initial value problem. Now, here we need to supply one more condition. Now, since in this initial value problem we know the exact solution. The exact solution is e raised to lambda x. What is y 1? y 1 is approximation to y at x 1. Our interval is 0 to b. We are subject to x 1. So, we are dividing it into n equal parts. So, our x 0 is 0, x 1 is going to be equal to h. So, then y 1 which is approximation to y at x 1 that is going to be where we will supply it to be e raised to lambda h. So, now, let us find solution of this difference equation with two conditions y 0 to be equal to 1 and y at 1 to be equal to e raised to lambda h. Now, the solution of the difference equation assume that it is of the form r raised to n where r is to be determined. Then, y n plus 1 will be r raised to n plus 1 will be minus 1 minus 2 h lambda y n is r raised to n minus y n minus 1 which is r raised to n minus 1. This will be equal to 0. This will imply that r square minus 2 h lambda here there should be r minus 1 is equal to 0. So, it is a quadratic equation r square minus 2 h lambda r minus 1 is equal to 0 and our y n is r raised to n. So, the solution the roots of this quadratic equation they are given by r 1 is equal to 2 h lambda plus square root of 4 h square lambda square plus 4 divided by 2 which is nothing but h lambda plus 1 plus square root of 1 plus h square lambda square. The second root will be given by h lambda minus square root of 1 plus h square lambda square. General solution y n is given by beta 1 r 1 raised to n plus beta 2 r 2 raised to n. Beta 1 and beta 2 they will be determined by the conditions y 0 is equal to 1, y 1 is equal to e raised to h lambda. So, thus in case of midpoint method the solution is given by y n is equal to beta 1 r 1 raised to n plus beta 2 r 2 raised to n out of which one part beta 1 r 1 raised to n plus beta 1 r 1 raised to n is going to correspond to our exact solution. Our y n is approximation to y at x n, y is the exact solution. The solution of our difference equation is given by beta 1 r 1 raised to n plus beta 2 r 2 raised to n. One part corresponds to the exact solution, other part is extraneous solution. It appears because you are replacing differential equation of first order by difference equation of the second order. So, this extraneous solution this is this is what causes the problem like depending on the value of lambda it corresponds to the exact solution. Now, how this happens let us look at it in more detail. So, we have y n is equal to beta 1 r 1 raised to n plus beta 2 r 2 raised to n, y 0 is equal to 1, y 1 is equal to e raised to h lambda. We have seen that r 1 the first root is given by h lambda plus square root of 1 plus h square lambda square and this is the second root, y 0 is equal to 1, put n is equal to 0. So, this will become beta 1 plus beta 2 is equal to 1, then put n is equal to 1. So, y 1 that means beta 1 r 1 plus beta 2 r 2 will be equal to y 1. So, the constants beta 1 and beta 2 they are determined by these two condition. So, the first condition tells us that beta 2 is going to be equal to 1 minus beta 1. The second equation tells us when I substitute for beta 2 I will get beta 1 into r 1 minus r 2 plus r 2 is equal to y 1 that gives us beta 1 is equal to y 1 take r 2 on the other side. So, y 1 minus r 2 and divide by r 1 minus r 2, but look at here r 1 is h lambda plus square root of 1 plus h square lambda square r 2 is equal to 2 is h lambda minus square root of 1 plus h square lambda square. So, when I consider r 1 minus r 2 h lambda will get cancelled and you are left with 2 into square root of 1 plus h square lambda square. Now, as h tends to 0, so when h tends to 0, what is going to happen is y 1 is e raise to h lambda. So, when h tends to 0 e raise to 0 will tend to 1 then r 2 this part will tend to 0 this part will tend to minus 1 and hence when I consider y 1 minus r 2 that is going to converge to 2 in the denominator you have got term 2 into square root of 1 plus h square lambda square. So, that will tend to 2, so beta 1 will tend to 1 as h tends to 0 and beta 2 will tend to 0. So, our by n which consists of 2 parts the part beta 1 r 1 raise to n that corresponds to our exact solution and this is something extra. Now, let us look at this in more detail carefully. What is our r 1? r 1 is h lambda plus square root of 1 plus h square lambda square as we did in case of e raise to x we can do for f x is equal to root of 1 plus x its derivative is given by 1 upon 2 square root of 1 plus x. So, f x is equal to 1 f of 0 plus x times f dash c. So, what I am interested is saying that square root of 1 plus x is 1 plus terms of the order of x then our r 1 r 1 will be h lambda plus square root of 1 plus h square lambda square. So, this I can say that this is the derivative of r 1 r 1 is equal to 1 plus terms of the order of h square lambda square. So, r 1 is equal to h lambda plus square root of 1 plus h square lambda square. So, that is h lambda plus 1 plus higher order terms for r 2 the second root we have h lambda minus square root of 1 plus h square lambda square. So, using similar argument you will get this to be equal to h lambda minus 1 plus higher order terms. Now, r 1 raise to n that is 1 plus h lambda raise to n plus higher order terms r 2 raise to n will be minus 1 raise to n minus 1 raise to n 1 minus h lambda raise to n plus the higher order terms. So, now our solution is beta 1 r 1 raise to n plus beta 2 r 2 raise to n. So, let us substitute in that we are trying to look at what our y n looks like when we considered Euler's method or Runge-Kutta method things very easier. We got a formula for y n plus 1 in terms of y n then the same formula tells us the relation between y n and y n minus 1. So, like that we could go up to y 0 y 0 is equal to 1 that is given to us. So, that is how in case of Euler's method we got y n to be equal to 1 plus h lambda raise to n. Similar was the case with Runge-Kutta method for the midpoint method we have to solve a quadratic equation get its roots and then our solution is y n is equal to beta 1 r 1 raise to n plus beta 2 r 2 raise to n. Now, this r 1 raise to n and r 2 raise to n we have shown that r 1 raise to n will be 1 plus h lambda raise to n plus higher order terms. Our exact solution is E raise to lambda x. So, when I consider y at x n. So, this is the exact solution x n is going to be n times h. So, y at x n is E raise to lambda n h. So, as we compared in case of Euler's method or Runge-Kutta method y at x n minus y n we want to do similar thing here. So, that is why r 1 raise to n we write as 1 plus h lambda raise to n. So, then one can compare E raise to lambda n h and 1 plus h lambda raise to n. So, that is the idea. So, we have y n to be beta 1 r 1 raise to n plus beta 2 r 2 raise to n r 1 raise to n is 1 plus h lambda raise to n plus terms of the order of h square lambda square. For r 2 raise to n we have got minus 1 raise to n and 1 minus h lambda raise to n. This 1 plus h lambda raise to n will be equal to E raise to n h lambda plus terms of the order of h square lambda square. Similarly, here you will get E raise to minus n h lambda and thus our y n is equal to beta 1 E raise to lambda x n plus beta 2 minus 1 raise to n E raise to minus lambda x n plus terms of the higher order. So, we have this is our initial value problem. This is the exact sum of the solution y x is equal to E raise to lambda x y n is approximation to y at x n. So, y at x n is going to be equal to E raise to lambda x n. So, this presence of this term is normal y n is approximation E raise to lambda x n is the exact solution you have got this error h square lambda square, but this term beta 2 minus 1 raise to n E raise to minus lambda x n. Now, if 1 is doing exact calculations then the initial value problem our initial values they take care like you have got 2 terms, but the coefficient of the term is extra which should not be there its coefficient tends to 0 as h tends to 0. So, when h is small enough the second term should not matter because we have got beta 1 E raise to lambda x n plus beta 2 E raise to minus lambda x n plus term which is of the order of h square lambda square. We have shown that beta 1 will tend to 1 as h tends to 0 the term order of h square lambda square as h tends to 0 that term also will tend to 0. The second term is beta 2 E raise to minus lambda x n of which the coefficient beta 2 is tending to 0. So, the exact computations will not pose any problem, but we do the computations in finite precision. So, that is why there are going to be round off error and when you look at the term beta 2 E raise to minus lambda x n. Suppose your lambda is less than 0 we are looking at y dash is equal to lambda y. So, if lambda is less than 0 E raise to minus lambda x n that becomes significant. So, we have y n is equal to beta 1 E raise to lambda x n plus beta 2 minus 1 raise to n E raise to minus lambda x n which is extraneous solution and the terms of the higher order. If lambda is less than 0 E raise to minus lambda x n will increase you are going away from 0 you are looking at the interval 0 to b. So, the term will increase exponentially our beta 2 because of the round off error it will not tend to 0. So, this term will start dominating this term our approximate solution consists of two parts. So, this is the part which corresponds to the exact solution, but if lambda is less than 0 then this solution becomes dominant which is corresponding to the error and you will get the error to be increasing. Last time we looked at a specific example it was y dash is equal to minus 2 y plus 1 and there I have shown you the numerical results where the Euler's method if you apply to that particular example the error goes on decreasing. Whereas in the case of midpoint method the error instead of decreasing it went on increasing we had looked at the interval to be 0 to 4. So, this is the problem with the midpoint method to the same problem whether lambda is less than 0 or lambda greater than 0 Euler's method Dranke-Kutta method there will not be any such problem. So, that is why if your lambda is bigger than 0 then you should use midpoint method because then in that case even though you have extraneous solution beta 2 e raise to minus lambda x n your beta 2 whether it tends to 0 or not e raise to lambda x n with lambda bigger than 0 that will become smaller and smaller and then it would not have it would not affect our true solution. Now, I want to consider or I want to rather state result for the Adams-Bashford method and then we will go to boundary value problem. So, in case of Adams-Bashford method again it is a multi step formula. So, let us look at the formula. So, we have y n plus 1 is equal to y n plus h by 24 55 g n minus 59 g n minus 1 plus 37 g n minus 2 minus 9 g n minus 3 where g n is value of g at x n g is our function f x y. So, it is f x n y n. So, here we have got a formula for y n plus 1 in terms of y n then y n minus 1 y n minus 2 y n minus 3. So, this is going to be or if we if you consider y dash is equal to lambda y you will get a difference equation of order 4. In case of midpoint method we had difference equation of order 2. So, we got a solution of the type beta 1 r 1 raise to n plus beta 2 r 2 raise to n. In this case it is equal to 0. So, it is going to be difference equation of order 4. So, that is why our solution will be of this type beta 1 r 1 raise to n plus beta 2 r 2 raise to n plus beta 3 r 3 raise to n plus beta 4 r 4 raise to n of which one corresponds to the true solution. All the three remaining terms they are extraneous they should not be there, but in this case what happens is what one can show is modulus of r 2 modulus of r 3 and modulus of r 4 they are less than 1 for h small enough. And then because they are less than 1 their effect will go on diminishing when you increase n and it they would not affect our true solution. So, the Adams-Bashforth method that is going to be stable for this y dash is equal to lambda y y at 0 is equal to 1. So, now this completes our study of initial value problem. Now, we want to consider boundary value problem. So, as the name suggests the boundary value problem will be the conditions will be given at the two end points of the intervals or at the boundary. So, we are going to look at second order differential equation with two boundary conditions given. Now, what I am going to do is just describe what is the finite difference method. So, it is a classical method. So, we will describe the finite difference method and then that will complete our solution of a progress. And then that will complete our solution of approximate solution of differential equation. There are other methods for solving the boundary value problems which we will not be considering. So, we have got this is our boundary value problem y double dash x plus f x y dash x plus g x y x is equal to r x for x belonging to interval a b. You need to have two conditions. So, the conditions are y a is equal to alpha y b is equal to beta y dash is d y by d x y double dash is second derivative d 2 y by d x square. So, our assumption is the coefficient functions f and g and the right hand side are these are going to be continuous functions on interval a to b. So, this is second order boundary value problem. We can have fourth order boundary value problem which can have the fourth derivative or the derivatives up to order 4 appearing and then the boundary conditions they will involve not only the function values, but may be also the derivative values or their combination. So, now, here is our boundary value problem where the derivatives they are appearing. Once again our aim will be to find approximation to the exact solution y at finite number of points. As we had done before we will divide our interval a b into n equal parts and our aim is to find approximation at these n plus 1 discrete points. So, we have this is the uniform partition of interval a b, h is the length of sub interval b minus a by n and interior mesh points they are given by x n is equal to x 0 plus n h, n is equal to 1 2 up to n minus 1. When x n is equal to x 0, it is x 0 is equal to a the end point. When n is equal to capital N then it is going to be the right end point b. The value of our function y is given at a and at b. So, those values they are known. So, our aim will be to find approximation to y at x n at these interior mesh points. Now, this is satisfied for every x belonging to a b. So, in particular it will be satisfied at x n. So, we have got y double dash x n plus f of x n y dash x n plus g x n y x n is equal to r x n, n is equal to 1 2 up to n minus 1. I do not write for the end points because we already know what is y at x 0 and y at x n. So, now the idea in the finite difference method is this y dash x n and y double dash x n these are the derivative values. So, you replace these derivative values by finite difference. We have considered numerical differentiation. So, if you have got a function f it is derivative at a. So, we looked at the forward difference formula that is f dash a is approximately equal to f of a plus h minus f a divided by f of a minus h divided by 2 h. So, you choose a point a plus h a minus h and then look at the values there, take their difference and divide by 2 h. Now, the advantage is in that case the error was of the order of h square whereas, you get a error to be of the better order if you consider central difference formula. Same thing is true for the second derivative f double dash a will be approximately equal to f of a plus h minus 2 of f of a plus f of a minus h divided by h square. So, what we are going to do is this y dash at x n and y double dash at x n we will replace the values by central difference formula and that is going to give us a tridiagonal system. So, our unknown is going to be y n which is approximation to y at x n. So, here is our equation. So, this equation is exact there is no approximation. So, y x n its approximate value we denote by y n this is the central difference formula. So, we have got y dash x n to be approximately equal to y at x n plus h minus y at x n minus h divided by 2 h. So, this will be nothing but y n plus 1 minus y n minus 1 divided by 2 h. So, this is for the first derivative for the second derivative y double dash x n will be approximately equal to y at x n plus h minus 2 times y x n minus y at x n minus h whole thing divided by h square. So, this will be equal to y n plus 1 minus 2 y n here it should be plus plus y n minus 1 divided by h square. So, these are the values we substitute in this formula. So, when you substitute you you are going to have a relation which involves y n y n minus 1 and y n plus 1. So, we have this is the substitution for y double dash x n y n minus 1 minus 2 y n plus y n plus 1 divided by h square plus f x n and then y dash x n. So, that is y n plus 1 minus y n minus 1 divided by 2 h plus g x n plus y n y at x n. So, that is y n is equal to r x n n is equal to 1 2 up to n minus 1. So, we get n minus 1 equations in n minus 1 unknowns our n minus 1 unknowns are going to be y 1 y 2 up to y capital N minus 1 y 1 will be approximation to y at x n. Let us simplify multiply throughout by h square collect the coefficients of y n minus 1 y n and y n plus 1. So, it will become y n minus 1 you are multiplying by h square. So, from here it will be 1 from here it will be minus h by 2 f n f n is nothing but f of x n into y n minus 1 plus coefficient of y n from here it will be minus 2 and from here it will be h square g n plus coefficient of y n plus 1 from here it is 1 and from here it is going to be h by 2 f n y n plus 1 is equal to h square r n n varying from h square minus 1 2 up to n minus 1. So, thus this for the sake of simplicity let me call this as a n the diagonal entries as d n and these values as c a. So, we have got a system of the type a n y n minus 1 plus d n y n plus c n y n plus 1 is equal to h square r n. If I look at n is equal to 1 when n is equal to 1 y 0 is given to us it is alpha. So, I can take it on the other side. So, I get for n is equal to 1 it becomes d 1 y 1 plus c 1 y 2 is equal to h square r 1 minus a 1 times alpha. Similarly, when I consider n is equal to capital N minus 1 y n plus 1 will be y capital N and y capital N is given to us y capital N is y at b that is equal to beta. So, let me take it on the other side. So, this is going to be the last equation a n minus 1 y n minus 2 plus d n minus 1 y n minus 1 is equal to h square r n minus 1 minus c n minus 1 beta. So, in each equation you have got diagonal entry and then two more entries and that is it. So, thus our system is going to be tridiagonal. So, you have got this is going to be your coefficient matrix along the diagonal d 1 d 2 d 3 up to d n minus 1. You will have one super diagonal and then you will have one sub diagonal and then we have looked at various methods of solving system of linear equations. So, one solves this system after solving the system you will get y 1 y 2 up to y n minus 1 that was our aim in solving that we said that y is our unknown solution and then we will try to find its approximation at some discrete points. So, y at x n we wanted to find. So, when you replace the derivatives by the difference formula or by numerical differentiation then you get a system of linear equation and then by solving that you will get the approximate values. Now, I just want to recall the result about numerical differentiation. For numerical differentiation what we had done was f x is going to be equal to f of a minus h plus divided difference based on a minus h a plus h into x minus a plus h. So, this is a linear polynomial which interpolates the given function at a minus h and a plus h and this is going to be the error term. You take the derivative. So, f dash x will be equal to f of a minus h a plus h. The derivative of this is given by f of a minus h a plus h x x multiplied by the product of these two terms plus the divided difference and then you have to take the derivative of this. So, that is x minus a minus h plus x minus a plus h when you put x is equal to a this term will go away. So, you get f dash a to be equal to f of a minus h a plus h minus h square f of this a minus h a plus h a a. So, thus f dash a is approximately equal to this with the error to be of the order of h square and then similar is the case for the second derivative f double dash a plus h a plus h a. So, if you consider the central difference formula then again the error is of the order of h square. Now, here is the example y double dash plus y is equal to 0, y at 0 is equal to 0, y at 1 is equal to 1. So, that means in our notation f x is identically 0 there is no y dash term g x which is coefficient of y that is 1 and q x is equal to 0. So, the diagonal entries they are going to be equal to 0 d n they are going to be equal to h square g n minus 2. So, that will be nothing but h square minus 2 whereas a n and c n these are equal to 1. So, that is our tridiagonal system. So, thus so we have looked at now solution of differential equation the boundary value problem by finite difference method finite difference method is a classical method. Now, there is a finite element method. So, which has been now used successfully it competes now with finite difference method. The finite difference method has advantage of simplicity that you just replace the derivatives by numerical differentiation. Here what we have done is we have replaced the derivatives by numerical differentiation with the error to be of the order of h square. So, whatever approximations y in which we obtain again the error is going to be of the order of h square. So, y at x n minus y n will be of the order of h square. So, if you want whatever accuracy you want you may have to choose h very small. If you choose h small it increases the size of your system of linear equations. So, then other way is use some better formulae. So, for the numerical differentiation use some better formulae but then that becomes complicated. Now, the finite element method it is one of the important method, but we will not be doing it in this course. So, now our next topic is going to be Eigen value problem. Eigen value problems like when you consider matrix of size 5 or bigger you cannot have a formula for Eigen values. So, that is why one tries to find as much information of Eigen values as possible indirectly that is without finding and then we are going to consider some approximate methods for finding Eigen values and Eigen vectors. So, that is the topic which we will start in our next lecture. So, thank you.