 Good morning. In the last lecture, we studied numerical method for solution of ordinary differential equation in particular the initial value problem. In the initial value problem at one value of the independent variable, all the components of the dependent variable vector y is given, and that signifies an initial value problem. We will have quite a few further issues to discuss in the numerical solution of initial value problems, but today's lecture let us start with a brief introduction to the boundary value problem, after which we will return to the advanced issues of initial value problem. So, let us start with boundary value problem. From the initial value problem to boundary value problem in differential equations, there is a complete paradigm shift, which can be visualized with the help of these two questions. Suppose we ask this question, that is a ball is thrown with a particular velocity. Velocity means speed and direction combined. So, a ball is thrown with a particular velocity at an instant of time. What trajectory does the ball follow? That is in future, what will be the trajectory of the ball? That means at the starting time, all the conditions are given. The complete state vector is given position as well as velocity, and we want to know the position and velocity in instant of time in future. This is a typical initial value problem. Correspondingly, and in contrast, a typical boundary value problem will look like this. How to throw a ball such that it hits a particular window at a neighboring house after 15 seconds? That means, at time 0, we know the position. We have to throw the ball from here, and at time equal to 15 seconds, we want a particular position. That is, the conditions are given half at this time instant and half at a future time instant. So, at the initial time, we are giving only the position and not the velocities. That is, how to throw a ball? On the other hand, at the end point, we are giving the position. That is, position at this point and position at that point. So, this signifies a boundary value problem. Now, this is a difficult problem physically as well as computationally. Physically also, you can consider this issue that if you want to target a particular window where the ball goes exactly after 15 seconds, for that you need to make a very detailed preparation. And possibly, you have to try one or two times before perfecting the speed and direction. On the other hand, for the first problem physically, all that you need is a sophisticated measurement mechanism by which at every time instant, you will go on taking the snapshots of the ball. Now, physically also, this is a more tractable problem. Comparatively, this is more tractable. Computationally as well, we will find that solving boundary value problems will require more iterations and will require greater computational resource. A particular case, we will discuss here, which is called a two point boundary value problem. The kind of problems that is illustrated by this question. So, in this kind of boundary value problems, the conditions are given at two point, two values of the independent variable as here. So, boundary conditions at two values of the independent variable will signify a two point boundary value problem. The number of conditions could be dispersed even at more than two points. For example, we could have given some position at one time instant, some other positions at another time instant, some velocity values at another time instant that is also possible. So, a specific kind of problems involves two points typically at the two end points of the domain. So, that is that constitutes boundary of the domain and that is why it is a boundary value problem. But in general, boundary value problems could have conditions given at various points, but we will consider here this class of problems for two point boundary value problems. So, the methods of solution, one is shooting method which is conceptually the easiest and then you have finite difference method, finite element method and so on. So, we will consider two methods here. One is shooting and the other is relaxation or finite difference method. In the case of shooting method, let us take it, let us understand it with the help of an example. Suppose at x equal to a, we have got some values given and at x equal to b, we have some other values given. Now, let us consider a problem in which the differential equation is this, where y is 5 dimensional. That means, let us take y 1, y 2, y 3, y 4 and y 5 be the state variable. Now, at x equal to a, these are the values. If all these values, all these five values at x equal to a are given, then it would be an initial value problem. Now, in this particular case, suppose at x equal to a, all values are not given, but then at x equal to b, some values are given. Total how many conditions we need, how many values out of ten values, how many values we need to consistently solve the differential equations, we need five values out of these ten values. So, for the five values, suppose at x equal to a, y 1 is given, at x equal to b, y 2 is given. That means, the value of y 1 is given at x equal to a, the value of y 2 is given at x equal to b, then value of y 3 is given again at x equal to a and value of y 4 is given at both ends and value of y 5 is given neither at this end nor at that end, that is possible. So, a total of five conditions are given. Now, from these five conditions and from this differential equation, we want to solve the differential equation, solve the boundary problem over this domain. Now, we make this observation that if we knew these two values also, then it would be so nice, we would solve it as an initial value problem. So, suppose we make a guess of these two values. So, say we take p 1 as this value, which is not given. Similarly, another guess if we make, which is this value, which is also not given by five. The two values at a, which are not given for that, suppose we make some guess. It is like trying to hit the target, we try to make a guess of the direction and then shoot it. Similarly, we try to make a guess of y 2 and y 5 at x equal to a and solve the initial value problem. Solution of an initial value problem is equivalent to pulling the trigger and letting the bullet go. So, as the bullet goes, typically first round it does not hit the target, it goes somewhere else. What is meant here? By saying that the bullet does not hit the target, it means that the positions, which are prescribed at the end point, they are not correctly achieved. So, that means with these guesses of y 2 and y 5 at x equal to a, the point in the state space where we reach, from that point we can work out, we can figure out y 1, y 2, y 3, y 4, y 5 at x equal to b and we find that at these two points the value does not match. So, at these points whatever is the value we do not know, the true value we do not know. So, there is no way to match. These two values what it should be we know and where we have reached that we get from the solution of the initial value problem. So, between these two we can work out the discrepancies. First error is error in this. So, with these guesses for the two unknown state variable values at x equal to a, with these guesses for them the initial value problem that we have solved gives the value of y 2 here as y 2 tilde and the true value is this at b. Similarly, the other value that is known for that also we can find out the error. Because of our specific choice of these two values we have got these two specific values of errors. Now, the question boils down to finding out the suitable values here such that these two errors are 0. In a typical target hitting problem what we do after we find the errors that is compared to the bullseye whatever shift we have made in our first attempt we make an estimate of the direction change that we should give to the gun. So, that next round it hits the target correctly accordingly we change the direction of the gun and hit again. So, that means the initial condition whatever was assumed was slightly erroneous we try to make a correction in that right and corrections should be made in such a manner that eventually these errors turn out to be 0. Now, note that with small changes in p 1 p 2 there will be small changes in e 1 e 2 right. So, that way we can consider e 1 e 2 as functions of p 1 p 2 and then try to solve this system of equations. In this particular case this system of equations is actually a system of two equations e 1 equal to 0 and e 2 equal to 0 in two unknowns p 1 and p 2 right. So, we can try our standard Newton-Aubson method for solving these two equations in these two unknowns and for that we will need the Jacobian. Jacobian means derivative of these two errors with respect to these two variables that means a 2 by 2 Jacobian matrix we will get we will need and for that we need to give small change to p 1 and small change to p 2 and register whatever is the differential change in e 1 and e 2 and accordingly build the Jacobian matrix. Now, that means that for every iteration after solving one initial value problem in order to find the Jacobian we will need to solve two more initial value problems one with a little shift in p 1 another with a little shift of p 2 at p 2 with at p 2. So, then with so many initial value problem solutions we will get the error value at the end at these two values and its Jacobian its rate of change based on that we can apply one Newton-Aubson step one correction and this can go on iteratively eventually hopefully converging to e 1 equal to 0 and e 2 equal to 0 which will mean that now we have converged to the correct values of p 1 p 2 that is correct values of these two state state variables at the point x equal to a and the corresponding solution of the initial value problem then will be the correct solution of the boundary problem as well because initial values are matched and final values have zero errors. So, this is the theme of shooting method that is you make trials to shoot the target these two and depending upon error at every iteration you try to change the initial values and accordingly try to perfect step by step your attempt. Now, here we took a simplified case in which three values were given it could be somewhat more general rather than giving three values at this point x equal to a there could be three relationships between all the five values similarly here rather than giving two values it could be two relationships among five values then also conceptually the problem remains the same because in that case if there are three relationships between these five variables then we can find out three parameters or rather two parameters in terms of which these five can be expressed those two parameters will be p 1 p 2 because five variables three equations will mean two choices in our hands. So, in terms of two parameters p 1 and p 2 we could express all the five unknowns or five variables. So, in general the shooting method will look like this it will follow the strategy to adjust trials to hit a target. Now, in a two point one value problem where this is the differential equation and there are a few relationships given among the state variables at x equal to a and another a few relationships at x equal to b. Now, the number of this let us say it is n 1 and the number of this dimension of g 2 suppose that is n 2 then n 1 plus n 2 must be equal to n because if the dimension of the problem is n if the state space is n dimensional then we will need n conditions less than n will not be enough more than n may be inconsistent. So, in that case what we do with the help of these n 1 equations we parameterize y of a as a function of n 2 parameters that is n minus n 1 this is n 2 parameters. So, n 2 values p 1 p 2 p 3 p 4 up to p n 2 so n 2 becomes the dimension of p. So, in terms of those n 2 variables we parameterize y at this end. Now, if we make choices for these n 2 values then we get a complete state vector at x equal to a with such guesses with n 2 values of p we make guesses and then we define the initial value problem same differential equation and y a completely given with these n 2 guesses of p made p 1 p 2 p 3 up to p n 2. Now, this initial value problem can be solved over the interval a to b and we can evaluate y of b and then we try to put this y of b into this equation and that will give us n 2 equation which we try to solve. So, error vector is this and we try to solve e p equal to 0. So, that means that for the solution of this system from the current vector p n 2 perturbations having perturbation components in the e 1 e 2 e 3 e 4 direction that is coordinate shift we can evaluate the Jacobian matrix. That means that each Newton action step will require solution of n 2 plus 1 initial value problem 1 for getting e of p and another n 2 here to evaluate the Jacobian. After that we can execute 1 Newton action iteration which is next p is equal to current p minus that that Jacobian inverse into current error. So, this is the typical Newton step for this problem. So, shooting method this method has the advantage that it requires very few parameters to start with only n 2 variables we have to select and then we can start shooting. On the other hand the disadvantage is that there is no guarantee of convergence and the problem is very sensitive to the working of the method is very sensitive to the initial guess. That means initial guess becomes very important if the initial guess is far away then the iterations may go anywhere and that is the typical drawback of Newton's method. On the other hand it is a good method to start with that is if you have a problem for which you need the solution it is a good idea to first try shooting method quickly if it works you have your solution in hand. On the other hand if it does not work then a more sophisticated method can be tried which is finite difference method. Finite difference method adopts a global perspective. So, for this same problem let us see how a finite difference method will work. In a finite difference method first we will discretize this domain let us say further time being that we are discretizing this domain a to b into 8 steps. So, 8 interval this is x 0 and this is x 8. This is x 4 this is x 6 this is x 2 this is x 1 x 3 x 5 x 7. So, there are 9 values of x n is 8 and there are 9 values of x starting from x 0 up to x 8 and 8 intervals x 0 x 1 x 1 to x 2 and so on 8 intervals. Now, if we do that then what we can do that for every interval we can write this differential equation in discrete form. So, then from the differential equation we get what is called finite difference equation. So, for example suppose at x 3 if this is the value of y now here I am showing y directly along this axis, but actually that is not a one dimensional entity there is a 5 dimensional entity. So, suppose the vector y vector at x 3 is something and at x 4 it is this. Now, the differential equation written for this interval at what point we write the differential equation suppose we try to write it in the at the midpoint then it will look like difference in y that is y i minus y i minus 1 difference in y is equal to d y by d x which is this into difference in x h here x 8 minus x 0 by n or b minus a by n and this is same as x i minus x i minus 1. Now, this f evaluated at which point we are planning to evaluate it at the midpoint. So, that will mean x i plus x i minus 1 by 2 y i plus y i minus 1 by 2. So, now consider this this is an equation in x i x i minus 1 y i y i minus 1 and that is it right you can also say you bring the entire thing on this side and how many equations are here y has a dimension of 5 in this particular problem. So, actually you have got 5 equations here 5 equations here and like that you will get 5 equations over every interval. So, that means into 8 8 intervals are there. So, you will get 40 equations 40 equations you will get in these variables. The variables the unknowns involve the values of y values of x are known right. So, unknowns involve values of y. So, you have got 40 equations and how many unknowns are there there are 9 points at each point there are 5 unknowns that means 45 unknowns 45 unknowns 40 equations 5 equations less those 5 equations you get from here 1 2 3 4 5 5 conditions are given. So, 45 equations in 45 unknowns you can solve now there are 2 issues. If the differential equation is linear then this will be a linear function of y right and that means that this system of equation that you get from here will be a linear system of equation. On the other hand if the differential equation is non-linear then you will get non-linear system of equation and for solving them iterations will be necessary this is one issue. Second issue is of course in the first issue even if the differential equation is linear and the resulting algebraic equations that you get out of this are also linear still iterations could be used and at times iterations are found useful and relaxation method typically uses such iterative methods like Gauss Seidel method or Jacobi method to solve the linear system. Now, if the system is non-linear then the differential equation system is non-linear then this algebraic equation system is non-linear and the iterations will be necessary you cannot solve it in one step. So, you have to use iterative method and for that point you need Jacobian. Now, one good issue here good point here is that differential that the finite difference equation that you see here will involve variables only from here and here not from here not from here that means every bunch of equations like this every bunch of five equations that you get from here will involve only a few variables which are close by in around this interval and that means that the Jacobian that you will get out of it will be part and here you will take advantage of the percity of the matrix and you will use those kinds of method which work well with parse coefficient matrix for solution of linear system. Now, let us see the outline in a generalized version generalized form first of all finite difference relaxation method adopts a global perspective that is it tries to get the solution all over the domain at once and from this point to this point everywhere it tries to adjust the solution together not that it not the not like shooting method in which we try to build the solution from this end towards that end. Here in the global perspective all the parts of the domain are given due importance and they are adjusted all together. Now, one more issue when we try to frame this equation using this derivative function here then for developing these derivatives or the Jacobian there of we will need the values of x i x i minus 1 y i y i minus 1 and so on and therefore to begin with even at the very beginning of solution we need to assume one solution from this end to that end that may be wrong in general that will be wrong, but then one value one set of values for y at every point of the grid is necessary to start the solution process. So, what we do we discretize the domain a b with a grid of points x 0 to x n n was taken as 8 in this example then function values at x i will mean n into capital N plus 1 unknown capital N plus 1 because capital N plus 1 points are there into n which is the capital dimension of y. Now, we will replace the ordinary differential equations over intervals by finite difference equations which will look like this and this will give us n into n small n into capital N scalar equations that means capital N vector equations at so many interval and each vector equation will stand for small n scalar equations component wise. Now, we will assemble the additional n equations from the boundary conditions and then this total number of equations will be same as total number of unknowns and then starting from a guess solution over the grid we solve this system and at the end of the process we will get the complete solution. So, this method is called the relaxation method which works based on the finite difference equations. There is also a finite element method which is much larger than this course will allow us to discuss in that we do not consider finite differences like this, but we try to split the entire domain into small elements finite elements and we try to propose the type of functions with appropriate continuity requirements across boundaries and then the details of those functions we try to work out through iterations through the numerical method. So, finite element method constitute a full course in itself and so we do not go into the details of that towards the end of this course when we discuss variation of calculus there we will make a small introduction to the basic ideas of finite element method and now we go back to the initial value problems and consider a few important and interesting issues in the solution methodologies of initial value problems. In the previous lecture we discussed method for solving initial value problems of ordinary differential equations. Now, we try to address a few important questions regarding the way the method handles the differential equations. We try to make a stability analysis of the method. Adaptive Runge-Kutta fourth order method or fourth order fifth order method is an extremely successful method, but its scope is limited. Now, the focus of explicit method such as Runge-Kutta method is accuracy and efficiency of the method. For the purpose of efficiency we use adaptation of step sizes. Now, there is another important issue in the solution process of ordinary differential equation and that is stability. Stability of the method not the stability of the dynamic system for which the differential equation is framed that is another issue which we will discuss much later in the course. Now, stability of the numerical method means that if there is an error at one point then as the computation progresses does the error grow or does the error decay. If the error grows then we say that the process is not stable. Now, the issue of stability is handled in explicit method quite indirectly through the root of accuracy. Now, if we try to make an analysis of the stability then let us consider this differential equation and rather than starting our analysis from R K 4 method that is fourth order method let us start from Poisson method itself first order method the discussion will be simpler. The Poisson method gives us this from y n we get y n plus 1 to this relationship which we saw in the previous lecture. Now, Taylor series of the actual solution around y equal to x x equal to x n is this y at x n that is the true value of x n true value of y at x n plus f evaluated at this correct x and correct y into h plus higher order terms. You see what are the differences here? Here it is y n plus 1 the value of y at x n plus 1 given by Euler's method this is similarly the value at x equal to x n given by Euler's method and so on. Here this is the true value this is the true value and the function is evaluated at the true point the true solution of the ODE system. Now, if we try to see the discrepancy between these two points then this difference that is the value obtained minus the true value that is the error or discrepancy at x equal to x n plus 1 and this is the difference y n minus y at x n that is the value obtained from the method minus the true value that is this plus similarly this minus this. So, f evaluated here minus f evaluated here into h and then the rest are higher order terms. Now, this is error at nth step and what is this function at x y n and function at x y of x n this is the evaluated y at x n this is the true y at x n. Now, from mean value theorem we know that between y n and y of x n in the y space there will be a point in the straight line join in the straight line segment joining this erroneous y and true y in the straight line segment joining these two points there will be some point at which the derivative of f will be such that this derivative into difference in y values will give this difference of f values that we get from the mean value theorem. So, between y n and y of x n in the straight line join of these two points suppose that point the existence of which is guaranteed by mean value theorem is y n bar. So, at that point y n we are not making any such discussion for x because x is same at these two points. So, at y n bar then this turns out to be the Jacobian of f with respect to y the derivative matrix into delta n. Now if we combine these two and call this matrix as j the Jacobian matrix then we will get identity into delta n plus j delta n h that means i plus h j into delta n plus of course higher order terms. So, you see we get a relationship between the discrepancy at n plus 1 th step and the discrepancy at n th step. So, that shows that Euler's step magnifies the error by a factor of this. Now delta n is a vector delta n plus 1 is also a vector error in y at x n error in y at x n plus 1. Now this vector is mapped to this vector with this matrix i plus h j and that will give you the magnification produced in the error as a result of this one step. Now the question is is this magnification of size. Now we are talking of magnification of vectors. So, if the initial vector is in this direction delta n and the resulting vector is in this direction then we will be talking about magnification in terms of sizes. Now is this size magnification size that is size of magnification greater than 1 or less than 1. If it is greater than 1 that will mean that the size of the error will keep on increasing. On the other hand if it is less than 1 then that will mean the size of the error will keep on decreasing. So, using this j loosely as the representative Jacobian why we say loosely as the representative Jacobian because as we move from one point to another x will change y will change f will change and therefore, j will also change in general. But as a representative Jacobian using this same notation j we can say that from the error at step one as we have taken n steps and come to the step n plus 1 then. So, many times this matrix has been multiplied of course, in actual practice the matrix in between will keep on changing because j is changing. But using the rough Jacobian this will be the situation. Now for stability we need that this error at n plus 1 step goes to 0 that is it goes on reducing. So, for reducing the error step by step we will need this matrix to have all its eigenvalues of magnitude less than 1. So, we need eigenvalues of this matrix to fall within the unit circle around the origin. The eigenvalues may be complex that is why we are talking about the eigenvalues falling within the unit circle around the origin in the complex plane that is this. This is the unit circle and the eigenvalues of this matrix we need their magnitude we need less than 1. So, by shift theorem we know that eigenvalues of this matrix are related to the eigenvalues of h j because addition of the identity matrix is basically increasing the eigenvalues by 1. So, if the eigenvalues must fall within the unit circle with origin at with the center at origin for i plus h j then the eigenvalues of h j should fall within this circle which is shifted by 1 unit this is origin this is minus 1 0. So, within this circle all the eigenvalues of h j should fall and that gives us this condition 1 plus h lambda should be the magnitude of that should be less than 1 which will give this limit on h that is now the real lambda real part of lambda if it is negative then the this will be positive and then it will make sense that is h is a new a positive. So, then it will make sense and then the method will be stable that is the error will not go on growing, but it will keep on reducing over step that is if lambda is lambda has negative real part and if h is less than this this is what you will get from this figure also see this unit circle. If lambda is lambda has positive real part then anyway it the eigenvalue cannot fall within this circle if it has negative real part and h is so small that eigenvalues of h j that is h lambda comes within this then there is a hope that the eigenvalues of h j fall within this and the method is stable. On the other hand if you take an h step size which is larger so that the eigenvalues of h j fall somewhere here here here then the method will be unstable because in that case eigenvalues of i plus h j will have sizes which are larger than 1. So, the same result this you will get if rather than the matrix vector equation system if you consider a single o d which is w prime equal to lambda w allowing lambda to be complex. Now, the way we got this result from the first order method that is Euler's method for the second order Runge-Kutta method will get a similar relationship in which this matrix will appear this expression will appear for lambda which earlier for the Euler's method was 1 plus h lambda here it will be a quadratic factor. So, there the region of stability will be that this magnitude must be less than 1 which is this elongated figure. So, for R K 2 that is Runge-Kutta second order method the eigenvalues of h j that is h lambda should fall within this for fourth order Runge-Kutta method. Similarly, you can work out the expression of h and lambda that also will be in the form of an inequality and when you plot it you will find that the values of h lambda for stability should lie within this and so on. For higher order method there will be similar slight expansion still the problem remains that for large zones in the complex plane the instability problem remains that is if the eigenvalues of h lambda fall on this side then the method is unstable even if the values of h lambda fall on this side, but depending upon the value of h if they fall outside this figure here, here, here then also the method is unstable that means for unstable differential equations for which the matrix j has eigenvalues with positive real part the method is anyway unstable and if the differential equation is stable that is eigenvalues have negative real part, but if the value of h is large so that the eigenvalues fall here the method is still unstable. So, then how does the step size adaptation of Runge-Kutta fourth order method operate on this kind of a system as we discussed earlier the issue of stability in R K 4 or such explicit scheme is handled indirectly that is as the error goes grows up as the error grows up step by step then after a few steps it is identified that the error growth has been too much and that will be identified through the comparison of two estimates as we discussed earlier. Through the comparison of two estimates it will be found that at the seventh step say the error is too much that means over this interval we cannot use a single step we will subdivided into two steps and then any eigenvalue which was found here that will be drawn inside like this and it will be found that now it is fine the error is out of two estimates the error estimate will show that error is not too much and error is fine, but slowly because of the step wise magnification of the error after a point again the error will grow too much and then due to the error growth the step size adaptation will be performed and this goes on and finally after quite a few iterations quite a few steps the adaptation will be taking place so frequently that you will have extremely small steps, but that then the errors will be quite large this is the problem of explicit steps, explicit method. So, step size adaptation tackles instability by its symptoms it does not handle the actual melody what will happen if we do not use an explicit method, but use an implicit method what it means it means that rather than writing the differential equation in terms of the derivative d y by d x at x and y n we could write it in terms of x n plus 1 y n plus 1 at the other end of the small interval that is rather than taking the slope from here if we say that what about taking the slope at the point where we would finally read at that point the question is that where we will finally read that is y n plus 1 we do not know now. So, that means that here we do not get an expression for the unknown y n plus 1, but we get a system of equation a system of equations in y n plus 1. So, question is we still have n equations in n unknowns. So, we can solve it by Newton's iteration or whatever, but first we ask the question is it worth solving it that is the problem of instability will that be addressed if we solve this and accordingly find y n plus 1. So, let us first try to see whether it is worth solving this in order to evaluate y n plus 1. So, for that we try to make a similar analysis what is the discrepancy at x n plus 1 if we use this kind of a formula that will be y n plus 1 minus that true value that will be y n minus y x n plus 8 into the function value. The same equation which we used earlier except that now this f is evaluated at x n plus 1 y n plus 1 rather than x n y n similarly here and that will give us this is delta n and this will be 8 j this j is slightly different from that earlier j, but still it is the Jacobian matrix at a representative point into change in y difference in y. Now, notice that earlier in the explicit scheme we got delta n plus 1 is equal to delta n plus 8 j delta n. Now, here we are getting delta n plus 8 j delta n plus 1. So, in order to express delta n plus 1 in terms of delta n what we do we take this on the other side and then we get i minus h j here identity into delta n plus 1 minus h j delta n plus 1. So, i minus h j inverted and pre multiplied here will give us this. Now, this is the difference here. Now, for stability of the method we need delta to decrease over steps that is compared to delta n we want delta n plus 1 to be less. That means the i can values of this matrix the inverse of i i minus h j we want to be to fall within the unit circle. We want the i can values of i minus h j inverse to be within this. So, for that suppose the i can value is here. So, for i minus h j inverse i can value of i minus h j inverse. So, for the i can value to fall here for i minus h j inverse where should be the i can value of i minus h j that will be reflection. Reflection first by this circle and then by the real axis. So, if it is a plus i b then it is inverse will be 1 by a plus i b that means a minus i b divided by a square plus b square. This conjugacy is equivalent to reflection by the real axis and this division is equivalent to reflection by this circle. That means the i can values of i minus h j should fall outside the unit circle outside the unit circle i can values of i minus h j and that means that i can values of h j should be again accordingly shifted. So, if we analyze that that is i can values of i minus h j should fall outside the unit circle and that will give us the limitation on h as greater than this. Now, if lambda has positive real part then h larger than this quantity will be required for stability and if lambda has negative real part then this is always true because h is positive and then that means it will be always stable. So, that means we get absolute stability for a stable ODE system that is this i minus h j having i can values outside this unit circle will mean that h j must have i can values shifted i mean outside the shifted circle like this that is this. So, if h falls h j if i can values of h j fall within this circle then the method will be unstable for all the rest of the domain on the complex plane the method is absolutely stable. So, this is the advantage of the backward Euler's method. Now, we come back to this question should we solve this and if so how. Now, since we know that the stability issue will be guaranteed then we know that it is what solving this system to find out y n plus 1 because we are using the slope function at the other unknown end of the interval rather than the known end. So, we get this system of equations in the unknown y n plus 1. So, we need to solve that so for solving that what we need we need y n minus y n plus 1 minus y n that is this that is the difference equation we write it and this is a system of equations in y n plus 1. So, we need to solve g of y n plus 1 that is this expression equal to 0 for y n plus 1 Newton step Newton iteration will be this that is value at step k plus 1 is equal to value at step k into i minus h j inverse because it is Jacobian Jacobian of this function will be minus i plus h j that means h j minus 1. So, in the Newton step you have y k minus Jacobian inverse. So, that slip in sign has been taken care of here that minus has become plus. So, minus i plus h j has become i minus h j. So, i minus h j inverse into the difference of the function value. So, you have got this the error value at kth iteration. So, that is this this whole thing evaluated at the kth iteration. Now, this is an iterative method iterative step now in a semi implicit Euler's method what we say is that anyway we will use the guess value for this Newton Euler's step at the as the value of y n the previous value that is known at x equal to x n. So, if we use that as the guess value here and here then the first step in the Newton iteration that you will get will have y n here. So, that is same as this. So, this goes up and y n here. So, the function will be evaluated at x n plus 1, but y n which is known. So, the function is evaluated here and this is the first iteration of the Newton's step Newton's method. Now, if you utilize only the first iteration then what you get is called semi implicit Euler's method otherwise the complete implicit Euler's method will require you to have this iteration till convergence. Now, similar iterative schemes similar implicit schemes can be developed for second order Nguyen-Gekuta method also. Now, these implicit schemes will have the advantage that stability will be guaranteed, but the accuracy is another matter which has to be handled separately. Now, in the next lecture in the beginning we will consider what is called stiff differential equations and see how the application of implicit method helps in handling stiff differential equations. And further in that in the next lecture after discussing stiff differential equations we will raise more fundamental questions regarding the validity of solution of initial value problems in a physical and practical sense. And from that we will develop some deep mathematical results which will give us confidence to solve the problem use the solution of initial value problems as we solve in the typical problems. So, these two issues we will discuss in the next lecture. Thank you.