 We are considering approximate solution of initial value problem. So, we have already considered Euler's method and the local discretization error in the Euler's method is h cube and the total error it is of the order of h square. So, now today we are going to consider in the Euler's method we had the local discretization error to be h square. Now, today we are going to consider Runge-Kutta method of order 2 and we will derive this method along with the error and then I will state Runge-Kutta method of order 2. After that we will consider multi step formulae for example, Adams-Bashforth method. So, we are looking at the initial value problem y dash is equal to f of x y x x belonging to a b with initial value y at a is equal to y 0 given, y is the unknown function to be determined and the notation is y dash is equal to dy by dx. What we do is we subdivide this interval a b into n equal part. So, a is equal to x 0 less than x 1 less than x n is equal to b each sub interval is of length h which is going to be equal to b minus a divided by capital N y at x n is the exact value of the unknown function and y n is going to be an approximation to this. So, our aim is to find approximation to y at x n at these n plus 1 points. The right hand side function it is f is a function of x and y and y is a function of x. So, we have got this function g x defined on interval a b when we want to calculate the derivative of g. So, that is going to be equal to g dash x is equal to by chain rule partial derivative of f with respect to x plus partial derivative of f with respect to y and now y is a function of x. So, it is y dash of x we know that y dash is equal to f. So, g dash x is going to be f x plus f y into f evaluated at x y. So, this is the first derivative. Now, let us look at the second derivative. So, the second derivative we have calculated the first derivative. In the first derivative we have got partial derivatives of f with respect to x partial derivative of f with respect to y and the function f. Function f partial derivative with respect to x and partial derivative with respect to y they are going to be function of x and y. So, when we want to calculate the second derivative we will again apply chain rule and then when you are considering partial derivative with respect to y into y dash we will be substituting for y dash is equal to f. So, taking these things into consideration the second derivative of g is going to be given by g double dash is f x x then f x y partial derivative of f x with respect to y multiplied by y dash. So, that is going to be equal to f. So, this is derivative of f x with respect to x these two terms. Now, we will be applying product rule. So, here f y x into f then f y y into y dash that is f and then there is one f here. So, it is f y y f square then f y into f x plus f y f y. So, that will be f y square into y dash that is f and all these functions they will be evaluated at x y x. So, we have calculated the derivatives of the right hand side with respect to x. Now, we are going to look at the Taylor series expansion for unknown function y. So, value of y at x n plus 1 will be equal to y x n plus h times y dash x n plus h square by 2 y double dash x n plus h cube by 6 y triple dash x n plus terms of the order of h raise to 4. So, this is the Taylor series expansion. I truncate this series by keeping the first three terms then we are going to have the local discretization error to be of the order of h cube. The formula which you obtain involves the function f and its partial derivatives with respect to x and with respect to y. What we want to do is get a formula which involves only function f and not its partial derivative. So, we try to determine y n plus 1 to be equal to y n plus a k 1 plus b k 2 where k 1 is going to be h times f of x n y n, k 2 is going to be h times f of x n plus alpha h y n plus beta k 1. So, the constants a, b, alpha and beta they are at our disposal. We want to determine these constants so that the first three terms they are going to match. Here, you have got y dash at x n y dash is equal to f. So, you have got a term f of x n y n. We also want to match the term which contains h square and then the formula which we will obtain will have a local discretization error to be of the order of h cube. So, the way we are going to do this is we have got f of x n plus alpha h y n plus beta k 1 use Taylor series expansion for function of two variables. So, this f of x n plus alpha h y n plus beta k 1 will be equal to f of x n y n plus alpha h and then partial derivative of f with respect to x evaluated at x n y n plus beta k 1 partial derivative of f with respect to y evaluated at x n y when we look look at the Taylor series expansion for y at x n plus 1 that also involves the terms f f x f y. So, then we will be able to compare the terms and determine a, b, alpha and beta so that the constant term the term which contains h and the term which contains h square they are going to match and that will give us a formula which involves only the function f, but the local discretization error is of the order of h cube. So, let us write down the Taylor series expansion for function of two variables f of x n plus alpha h y n plus beta k 1. The Taylor series expansion is f of x n y n plus alpha h into partial derivative of f with respect to x plus beta k 1 into partial derivative of f with respect to y all the functions to be evaluated at x n y n plus the next term will be alpha h square by 2 the second order partial derivative with respect to x plus alpha h beta k 1 f x y plus beta k 1 square by 2 f y y evaluated at x n y n plus term of the order of h cube. So, this is the Taylor series expansion for the function of two variables. Now, this is the Taylor series expansion for the unknown function y, y dash is f y double dash is nothing but derivative of this f. So, that is going to be equal to f x plus f y into f plus h cube by 6. Now, derivative of this and we have seen the derivative of this to be equal to f x x plus 2 times f x y f assuming that the mixed partial derivative they are equal that means f x y is equal to f y x plus f y y into f square plus f y into f x plus f y square into f x n and then y at x n. So, here this h square by 2 this function also will be evaluated at x n y x n and same thing here. So, our y n plus 1 the formula is going to be of this form y n plus a k 1 plus b k 2 where k 1 is h times f of x n y n k 2 is h times f of x n y n plus a k 1 plus x n plus alpha h y n plus beta k 1. Then we have seen this to be the Taylor series expansion. So, this Taylor series expansion we are going to substitute here and then we get y n plus 1 is equal to y n plus h times a plus b f of x n y n because our y n plus 1 is y n plus h times k 1 plus or a times k 1 plus b times k 2. Then so, we have y n plus 1 is equal to y n plus a k 1 plus b k 2 k 1 is h times f of x n y n k 2 is h times f of x n y n plus alpha h y n plus beta k 1. For this we write Taylor series expansion and then when we substitute in this formula. So, you are going to have h times a of f plus from here you will get b of f. So, that is why we have got y n plus 1 is equal to y n plus h times f of x n y n plus a plus b f of x n y n then h square times b alpha f x plus beta f y into f x n y n plus the coefficient of h cube is going to be this term plus term of the order of h raise to 4. So, it is just substitute the Taylor series expansion and then write down what you get then this is our Taylor series expansion for y at x n plus 1. We are going to compare these two we are going to neglect the error between y x n and y n we are going to look at only local discretization error. That means you have calculated approximation y n. Now, from y n I am going to y n plus 1. So, at this particular step for going from x n to x n plus 1 whatever is the error that is our local discretization error and that we are trying to determine. So, we will assume that y at x n is equal to y n. So, with this assumption when you compare these two y n is going to be same as y at x n. So, the constant terms match then here you have h times f of x n y x n which is same as f of x n y n and here you have f of x n y n. So, if we want these two terms to match our a plus b should be equal to 1. So, we get the first condition on a and b that a plus b should be equal to 1. Now, look at the next term you have got h square here h square here then here you have f x plus f y f divided by 2 and here you have got b alpha f x plus b beta f y f x n y. So, if I want this term and this term to match then what I should do is b into alpha should be equal to 1 by 2 b into beta should be equal to 1 by 2. So, we have a plus b is equal to 1 b alpha is equal to 1 by 2 b beta is equal to 1 by 2. So, thus if I choose my a b alpha and beta such that a plus b is equal to 1 b alpha is equal to 1 by 2 b beta is equal to 1 by 2 then constant term the term which contains h and the term which contains h square they are going to match with the assumption that y at x n is equal to y n. Now, we have got four constants we have got one more say degree of freedom at our disposal like I can impose one more condition. Now, whether by imposing that condition if I can also match the coefficient of h cube then that will be good because then our local discretization error will be of the order of h raise to 4. But no matter how you choose your a b alpha beta it will not be possible because when you look at y n plus 1 in the coefficient of h cube you have got f x x f x y f and an f y y f square. Whereas in y at x n plus 1 you have got these three terms and in addition f x f y and f y square f. So, no matter which way you choose your a b alpha beta you cannot hope to match the coefficients of h cube because of the presence of these. So, now we have got the condition a plus b is equal to 1 b alpha is equal to half b beta is equal to half four constants and only three conditions. So, there are infinitely many ways of choosing this. So, let us choose which is something symmetric. So, choose a is equal to b is equal to half alpha is equal to beta is equal to 1. Then the local discretization error is going to be of the order of h cube. So, this is the known as Runge-Kutta method of order 2. What you need to do is evaluate your function f at x n y n and then evaluate your function f at x n plus h because r phi is equal to 1 and y n plus k 1 you have calculated k 1. So, you calculate that and then you take the mean of these two a and b is equal to half. So, it is going to be h times f of x n y n plus h times f of x n plus h y n plus k 1 divided by 2. So, that is the Runge-Kutta method of order 2 and local discretization error is of the order of h cube. So, if I consider the error now y x n plus 1 minus y n plus 1. So, this is the term in y at x n plus 1 which contains h cube. The corresponding term in y n plus 1 which contains h cube is given by this I am subtracting and then you have got terms of the order of h raise to 4. So, this is equal to when you simplify it will be minus h cube by 12 times this factor plus the term of the order of h raise to 4. So, this is one of the undesirable property of Runge-Kutta method when you look at y at x n plus 1 minus y n plus 1. So, we look at the coefficient of h cube. So, that coefficient of h cube it is quite complicated and then you have higher order terms. I will talk about this little more in detail when we compare Runge-Kutta methods with other methods. Now, what I am going to do is state Runge-Kutta method of order 4. Whatever technique we have used we can try to do it for higher order methods like here we had tried to match the coefficients in the Taylor series expansion for the unknown function y which contain the constant terms the term which contain h and term which contain h square. Now, in the based on a similar principle one can derive a method which is known as Runge-Kutta method of order 4. The principle is the same the computations they become messy. So, let me just state what is Runge-Kutta method of order 4. So, that is y n plus 1 is equal to y n plus 1 by 6 k 1 plus 2 k 2 plus 2 k 3 plus k 4. So, what is k 1? k 1 as before is h times f of x n y n k 2 is going to be h times f of x n plus h by 2 and then y n plus the increment for y n is going to be this k 1 which we have found k 1 by 2 k 3 will be h times f of x n plus h by 2 as before and then y n plus this k 2 divided by 2 and k 4 is equal to h times f of x n plus h by 2. So, this is the principle which is known as h times f of x n plus h by n plus k 3 k 3 you have found here. It can be proved that the local error is of the order of h raise to 5. So, we have got now 3 methods one is Euler's method then we have got Runge-Kutta method of order 4 for the Euler's method we had y n plus 1 is equal to y n plus h times f of x n y n that means per step there was one function evaluation. When you look at Runge-Kutta method of order 2 we have got y n plus 1 is equal to y n plus k 1 plus k 2 by 2 where k 1 is h times f of x n y n k 2 is h times f of x n plus h y n plus k 1 that means for each step you are evaluating your function twice and then the local discretization error it improved from h square to h cube. So, your function evaluation is doubled, but the local discretization error improves from h square to h cube. In Runge-Kutta method of order 4 you have got y n plus 1 is equal to y n plus k 1 plus 2 k 2 plus 2 k 3 plus k 4 divided by 6 for k 1 k 2 k 3 k 4 you need to evaluate function once. So, you have got in all 4 function evaluations and the local discretization error is improved from h cube to h cube. Now, these 3 methods which I have described they are known as single step methods because in order to evaluate y n plus 1 what comes into picture is y n. Using y n you calculate k 1 and then that k 1 gets introduced in k 2 and so on, but they all depend on y n. So, these methods they are known as self starting method because y 0 is given to us that is the condition in the initial value problem the initial condition. Once you know y 0 using the formulae you can calculate y 1 then once you have y 1 you can go to y 2 and so on. So, these are known as single step methods and we have calculated or we have evaluated their local discretization error. Now, these methods these are classical methods. There are methods which are of relatively recent origin in which case instead of this Taylor series expansion we are going to use numerical quadrature. Now, what we are going to achieve is for example, in Runge-Kutta method we had local discretization error to be h raise to 5 and we had 4 function evaluation. So, whether I can reduce the function evaluations so, in the Adams-Bashforth method we will show that essentially doing only one function evaluation per step we will achieve the local discretization error to be h raise to 5. So, then that is very good that instead of 4 function evaluations we have got only one function evaluation. So, you are reducing your work by means it becomes one fourth by retaining still the local discretization error to be of the order of h raise to 5. Unfortunately, there is a trade off because if we had such a method available then Runge-Kutta method would have disappeared because the methods which we survive those are the ones which are better in some respect. So, here the method which we are going to derive which essentially needs one function evaluation per step it will not be self starting that means you will need not only y 0, but you will need y 0, y 1, y 2, y 3. So, these are the things that is not such a big disadvantage, but there are also stability considerations that single step methods they are always stable whereas, the multi step methods they can have stability problem. So, those things we will be studying in detail. So, now let us first see what is the principle behind this numerical quadrature. So, we have got our initial value problem y dash x is equal to f of x y x. Let me write this function f x y x as g x x belonging to a b. This is our uniform partition a to b integrate both the sides between x n to x n plus 1. So, you have integral x n to x n plus 1 y dash x dx is equal to integral x n to x n plus 1 f of x y x dx or integral x n to x n plus 1 g x dx. Integration on the left hand side is y at x n plus 1 minus y at x n plus 1 minus y at x n. So, if we can evaluate this exactly or if we can integrate this exactly, then value of y at x n plus 1 will be equal to y at x n plus this integral. Now, in general it will not be possible to integrate this exactly. So, we use numerical quadrature rule. So, the simplest numerical quadrature rule which we have considered is rectangle rule. In which case we approximate the given function by a constant polynomial and that constant polynomial is value of your function at the left hand point. So, then we considered also midpoint rule, then trapezoidal rule, Simpson's rule and so on. So, the other rules we will consider, there will be some problem. So, let us first look at what is the rectangle rule what it gives. So, we have g x that is our function f x y x is equal to g of x n and this is the error term divided difference of g based on x n x into x minus x n integral x n to x n plus 1 g x d x is approximately equal to h times g x n the error is given by integral x n to x n plus 1 divided difference x minus x n d x. The divided difference is continuous function provided your y is sufficiently differentiable. So, this is continuous x minus x n will be bigger than or equal to 0 on the interval x n to x n plus 1. So, you can use mean value theorem for integrals take this out and then replace it by g dash of c n integral x n to x n plus 1 x minus x n d x will be x n plus 1 minus x n square divided by 2 that means h square by 2. So, this is the rectangle rule and then y at x n plus 1 minus y x n this is the equation satisfied by the unknown solution y. When you apply numerical quadrature rectangle rule to the right hand side you are going to have y n plus 1 minus x n square divided by 2. So, this is the rectangular rectangle rule to the right hand side you are going to have y n plus 1 minus y n where y n plus 1 is approximation to y at x n plus 1. So, y n plus 1 minus y n is equal to h times g x n that means h times f of x n y n this is nothing but our Euler's method and the error is h square by 2 g dash c n which is same as h square by 2 y double dash c n. So, now the rectangle rule gives us back the Euler's method. So, when we derived Euler's method what we did was we looked at the Taylor series expansion for y truncated it keeping only one term or two terms we had kept and that gives us Euler's method. Now, the another way of looking at Euler's method is apply numerical quadrature to the right hand side use rectangle rule and then we get Euler's method. Now, next method is midpoint rule. So, in the midpoint rule what we will be doing will be looking at the midpoint of interval x n to x n plus 1 and then evaluating your function f at that midpoint. So, remember you have we had integral a to b f x d x to be approximately equal to b minus a by 2 f of or rather b minus a not a by 2 b minus a f of a plus b by 2. So, this is the midpoint rule. Now, if you consider integral x n to x n plus 1 f of x y x d x this will be approximately equal to h times f of x n plus x n plus 1 by 2 and then you will have y n plus x n plus 1 by 2 and then you will have y n or it will be y at x n plus x n plus 1 by 2. So, see we are looking at this uniform partition of interval a b. So, you have got a this is x n this is x n plus 1 the value here is y n value here is y n plus 1. If you apply midpoint rule then you are introducing one more value here. So, what we want is we want to determine approximation to y at the node point. We do not want to introduce some point in between and then those many values also coming into picture. So, as such the midpoint rule directly it is not applicable. Now, what one can do is instead of integrating between x n to x n plus 1 you integrate from x n minus 1 to x n plus 1. In that case the midpoint will be x n and that will give you a rule, but that part we are going to look at it little later. What will happen if I use trapezoidal rule? In case of trapezoidal rule you consider the value at the two end points. Now what we are saying is we have come up to y n and you are going to determine y n plus 1. So, if you apply trapezoidal rule your y n plus 1 will be determined only implicitly. So, again that comes under predictor character formally. So, that we will do little later. Now, what I am going to do is I want to consider interval x n to x n plus 1. I do not want to introduce new mode nodes like the midpoint of x n to x n plus 1 and so on. I want to consider a cubic polynomial approximation. So, when you consider cubic polynomial approximation when we had considered Newton quotes formally. What we did was we subdivided our interval a b into we want to consider four points that means three equal parts and then we fit a cubic polynomial. I do not want to subdivide interval x n to x n plus 1 into three equal parts because that will introduce new nodes. So, that I do not want to do, but what I can do is I have come up to x n. So, I know value at x n, but I also know value at x n minus 1 x n minus 2 x n minus 3. So, I can look at these four values fit a cubic polynomial and now integrate that cubic polynomial over the interval x n to x n plus 1. So, this is something different to numerical quadrature we have done so far. In the Newton quotes formally we would have subdivided we would have chosen four equidistant points in the in our interval of integration x n to x n plus 1. So, there and then we would have fitted a polynomial. Now, we are doing something else what we do is we go outside of our interval of integration. So, we have got we are integrating over x n to x n plus 1. So, you consider x n plus 1 and then x n minus 1, x n minus 2, x n minus 3. So, look at g of x n minus 3, g of x n minus 2, g of x n minus 1 and g of x n fit a cubic polynomial based on these four points. So, once you get a cubic polynomial it is defined over whole of r. So, use its value over x n to x n plus 1 and integrate because the cubic polynomial you know how to integrate. So, that is going to give us a formula which is known as Adams-Bashforth formula. So, let me derive that formula. The computations they are straight forward, they are bit more complicated or if you want to really work out the details. Now, what I am going to do is I am going to consider the cubic polynomial interpolation. The points are equidistant. So, we have got y at x n plus 1 minus y x n to be integral x n to x n plus 1 g x dx, g of x n plus 1. So, g is function f x y x, p 3 is a polynomial of degree less than or equal to 3 interpolating our function g at x n x n minus 1, x n minus 2, x n minus 3. So, the interpolation points except for this x n they are outside our interval x n to x n plus 1. Now, this is a formula you are familiar with g of x n. So, this is the cubic polynomial g x n plus divided difference based on x n x n minus 1 into x minus x n plus divided difference based on x n x n minus 1 x n minus 2 into x minus x n x minus x n minus 1 plus divided difference based on the 4 points x n x n minus 1 x n minus 2 x n minus 3 multiplied by x minus x n x minus x n minus 1 x n minus 2. So, this is the cubic polynomial and this is the error term. So, in the error term we have got x here and then w x and then w x and then x where w x is product of x minus x n x minus x n minus 1 x minus x n minus 2 x minus x n minus 3. Next, the divided difference based on x n x n minus 1 will be g x n minus g x n minus 1 divided by x n minus x n minus 1. This is h and this is the Newton's forward difference delta g n minus 1. So, delta g n minus 1 is precisely g of x n minus g of x n minus 1. The divided difference of g based on these 3 points you have got x n minus x n minus 1 is equal to h x n minus 1 minus x n minus 2 is also equal to h. So, this is equal to g of x n minus 2 times g of x n minus 1 plus g of x n minus 2 and 2 h square. The quantity in the numerator is delta square g n minus 2. So, it is the Newton's forward difference of second order upon 2 h square and now we will also write a similar formula for divided difference based on x n x n minus 1 x n minus 2 x n minus 3. That formula is g x n minus 3 times g x n minus 1 plus 3 times g x n minus 2 minus 2. This is equal to minus g of x n minus 3 equidistant points. So, upon 6 h cube. So, this will be the numerator is delta cube g n minus 3 divided by 6 h cube. So, thus our polynomial p 3 x is g n plus delta g n minus 1 into x minus x n by h. Because the term is x n minus 2. So, the next term here was divided difference based on x n x n minus 1 into x minus x n. So, we have the divided difference to be delta g n minus 1 divided by h. So, that h I am writing with x minus x n. The next term is delta square g n minus 2 divided by 2 h square into x minus x n x minus x n minus 1 plus delta cube g n minus 3 divided by 6 h cube. This is for the divided difference of g based on x n x n minus 1 x n minus 2 x n minus 3 multiplied by this. So, we have got y x n plus 1 minus y x n is equal to integral of x n minus 3 integral x n to x n plus 1 f of x y x d x. This we are going to replace by y n plus 1 minus y n is equal to integral x n to x n plus 1 p 3 x d x. If I make change of variable here as s is equal to x minus x n by h, then I get y n plus 1 minus y n to be h times x n plus 1 integral 0 to 1 p 3 s d s. So, we had our p 3 x to be given by this. So, our x minus x n by h is going to be equal to s x minus x n minus 1 by h that will be s plus 1 and so on. So, we have this h times integral 0 to 1 p 3 s d s. This is our p 3 x n minus x n by h is s x minus x n by h is s and x minus x n minus 1 by h is going to be equal to s plus 1 because we have x minus x n minus 1 is equal to say upon h is equal to x minus x n plus x n minus x n minus 1 divided by h. x minus x n by h is s plus h upon h. So, that is going to be equal to s plus 1. So, thus we have p 3 x n minus 1 by h is equal to be g n plus delta g n minus 1 s plus delta square g n minus 1 s into s plus 1 by 2 plus delta cube g n minus 1 g n minus 3 x minus x n by h will be s x minus x n minus 1 by h will be s plus 1 and x minus x n minus 2 by h will be s plus 2 by 6. So, distribute this h cube 1 h to each of the bracket and now we have got y n plus 1 minus y n is equal to h integral 0 to 1 p 3 s d s. So, you substitute. Now, you integrate. So, you will have first will be g n then delta g n minus 1 integration of s will be s square by 2 between 0 to 1. So, that becomes delta g n minus 1 by 2 plus delta square g n minus 1 by 2 integration of s square will be 1 by s cube by 3 between 0 to 1. So, it will be 1 by 3 integration of s again s square by 2. So, it is 1 by 2 plus delta cube g n minus 3 by 6 and now integration of this will be 1 by 4 plus 1 plus 1. When you simplify you will get it to be equal to h multiplied by g n plus delta square g n minus 1 by 2 plus delta square g n minus 2 multiplied by 5 by 12 plus delta cube g n minus 3 by g n minus 3 multiplied by 9 by 24. g n is g of x n delta g n minus 1 is g n minus 1 delta square g n minus 2 is g n minus 2 times g n minus 1 plus g n minus 2 and a similar expression for delta cube. So, we substitute that and then we get it to be equal to after simplification you get the formula h by 24 55 g n minus 59 g n minus 1 plus 37 g n minus 2 minus 9 g n minus 3. So, this is Adam's Bashforth method. Now, here when you look at the formula the coefficients they are not very nice numbers, they are not very good, but it does not matter you are going to write a program for calculating. So, whatever are the coefficient it is once for all one does the derivation, you write the program and then you know how to calculate the approximations to y n. Now, the formula which we obtain for y n plus 1 it involves y n plus 1 plus 1 plus y n minus 1 y n minus 2 y n minus 3. So, initially y 0 is given to us. So, you cannot calculate y 1. So, you will be able to use this formula say y 5 onwards. So, we have suppose y 0, y 1, y 2, y 3 they are given to us then using these formulae you can calculate y 5 or y 4 onwards. So, for y 4 you will need y 0, y 1, y 2, y 3 and then you can continue. So, that is why this method is not self starting you have to somehow get y 0, y 0 is given y 1, y 2, y 3 and then you can continue. In our next lecture we will see that here the discretization error is of the order of h raise to 5 same as in case of Runge-Kutta method of order 4, but I will explain what I mean by essentially one function evaluation per step. So, we will do this and then the in our next lecture we are also going to look at what are known as predictor-corrector formulae. So, one of the method which we are we are going to consider is known as Adams-Moulton method and then there are similar methods and we will compare these methods. So, thank you.