 methods for solving only initial value problems we have looked at the algorithms we have looked at the way algorithms are derived and now one of the most crucial aspects is how do I select the integration step size. I am at some point T n I want to move to T n plus h what should be integration step size h it is a crucial decision when you have to implement or when you have to use a ODIVP solver. Now I was describing one approach in which you are not really worried about the analysis but a practical approach to deal with this problem is variable step size solver or variable step size implementation this variable step size approach can be used only when you are going from T n to T n plus 1 without using the information from the past ok. So, Runge-Kutta can be used alternatively you could even tailor orthogonal colocations if you want using variable step size. So, variable step size is something that is possible when the method involves only going from only you need information at time T n you do not need past information about derivatives you do not need. So, multi step methods are difficult to tailor for variable step size, but the other Runge-Kutta methods can be easily tailored for variable step size implementation. So, here I do not have to worry about what should be the right steps I can just work with a changing step size I make a on the spot on the fly decision as to what should be my integration step size and next time it could be different next time it could be you know smaller or larger it depends upon how the system is behaving in a particular region. So, let me just write down the steps in the algorithm that way you will understand what is really happening yesterday I tried to explain the philosophy now let us look at. So, my problem is I am given x T n is equal to x n at T is equal to T n and epsilon is epsilon is the tolerance this could be 10 to the power minus 8 or it should be some small value. Now, my first step is to choose step size let me let us say we have chosen a step size h 1 then now using this step size h 1 what I do is I compute x n plus 1 I am choosing to index it as 1 because this is using h 1. So, using whatever method I have chosen it could be RK method it could be Euler method Rangekutta or Euler method whatever whatever Rangekutta class of method I have used as I said you could also tailor this for orthogonal collocation. So, whatever method you have chosen to march from T n to T n plus h h 1. So, you can use that and compute x 1 n plus 1. Now, step 3 now I am calling h 2 which is h 1 by 2 I am calling h 2 which is h 1 by 2 and now what I am going to do is. So, here just to put the notation correctly here. So, this is this one is computed at T n plus h 1 I am calling it as T 1 n plus 1 n plus 1 see n is the counter n is the counter right n is the index of time T is actual time and we have marched from in this first step 2 we have marched from T n to T n plus h 1. Now, after I define this h 2 my counter will see 2 steps of h 2 will be equal to 1 step of h 1 isn't it because I have divided this by 2. So, now I have to adjust and redefine my counter. So, now I will call T n plus 1 2 is equal to the counter we have to advance twice because you are going half step here. So, counter would have to advance twice. Now, what I do is now what I do is I compute x 2 n plus 1 x 2 now 2 here this index 2 is given because this is using h 2 this is using h 1 this is using h 2 that is why I given index 2 here. Now, as far as the time is concerned T 2 n plus 2 is same as T 1 n plus 1 right this this these 2 are equal this is equal to this is equal to this. So, what should happen is why you will be wondering why I given the tolerance the tolerance is to check whether the value that you get here and the value that you get here are the 2 different ok. If they are very close what do you mean by close is using tolerance if they are very close I accept h 1 as my step ok. If not I reset my h 1 to be equal to you know some fraction of initial h 1 you can halve it or you can make it 0.75 or whatever ok. So, now what I am going to do here is I am going to check whether this relative error this is remember this is x 2 n plus 2 minus x 1 n plus 1 because they are time matched ok you cannot use matching of n n and n here ok. So, if this relative error is very very small ok. So, whether you make whether you make 2 steps and go to that point or whether you make 1 single step it is not making too much difference. So, we accept we accept this as the next value else what we do is we shrink h ok. So, we shrink h 1 ok. So, what I do is if the tolerance condition is not met I shrink my h 1 ok. So, I shrink my h 1 to h 2 that is h 1 by 2 ok and then I redo the whole thing I keep repeating till this condition is satisfied. So, essentially I keep shrinking I start with some guess the initial guess you have to give I start with some guess ok and then you go on shrinking till you get tolerance satisfied. In MATLAB you have this Runge-Kutta solver there are 2 Runge-Kutta solvers which can be very easily used R K 2 3 and R K 4 5 ok. These solvers are actually reasonably very good solvers they can be used for large class of problems except some stiff equations what is the stiff business I will come to that. So, these are variable step size solvers they will ask you maximum integration interval ok and they will also ask you tolerance ok. The solution when it comes it will not be at one solution at the end point the solution will be profile because it will have shrunk the integration step size to the point where tolerances are met and then only it proceeds what is the advantage here ok. The advantage of this this is a practical way of solving the problem when you do not know how to choose a integration step size unfortunately this will not give you insights into what is happening but it it solves your problem does not. Now will of course will go to the point where we start getting insights or we start understanding why why this is happening before that I have given or I have described an algorithm which practically solves the problem ok. So, when you do not know what the step size to choose just use variable step size implementation it is safe even if you use Euler with variable step size it will work because even if you make a wrong guess of the step size it will keep shrinking the step till this condition is met ok. So, you do not accept till this condition is met ok. So, that is a practical way out now let us start getting into the insights earlier we were talking about linear difference equations arising out of iterative schemes right x k plus 1 or E k plus 1 was some matrix S inverse T times E k where E was the error and k was iteration index. Now I am going to talk about difference equations in time ok we are going to use difference equations in time. So, this concept of Eigen values is so fundamental that you know it pervades the analysis that we do and those of you who are doing the analytical methods for partial differential equations and ODEs that course also you will appreciate more and more that how this Eigen values play a role why they are so important of course there you get Eigen functions and Eigen values whereas finite dimensional case you get Eigen vectors and Eigen values. Now let us start getting into this business of convergence analysis. So, this convergence analysis is tightly coupled with selection of step size. So, two things cannot be really separated how do you choose integration step size will decide whether convergence to the true solution occur or not ok. So, these two things are not really separate entities ok though now we have a fix of how to go about if you do not know anything we would still like to get insights into what is what is really going on ok. It also allows you to compare different methods ok. Now what I am going to do is to do the analysis ok. First of all I need to know see when can you analyze whether a particular method is converging to the true solution only when you know the true solution right. If you do not know the true solution how will you analyze whether given method is going to the true solution or not ok it is very difficult to do that. So, the first criteria is that you should start looking at a system for which true solution is known ok. Then only we can start looking at the difference between the true and approximate solution that is the first first thing. So, which other which other systems for which we know the true solutions exactly whether it is scalar or whether it is vector case is linear differential equations initial value problem ODE initial value problem for linear ordinary differential equations could be scalar or could be vector. So, I am going to use linear ordinary differential equations initial value problem as a benchmark to understand to get insights into convergence analysis ok. How to extend it to a non-linear case will worry about it little later let us first get some insights as to what is happening and then let us let us see how to extend it to a non-linear case ok. So, I am going to consider two kinds of systems one is dx by dt is equal to A x, x belongs to R and x tn is equal to xn or I am given some you could consider this initial condition you could consider x0 is x0 is initial condition whichever way ok. The nice thing about this is that A is a constant A is a scalar x is scalar variable one variable differential equation initial value problem I know the solution I know the true solution ok. Now, I can use Euler method or whatever method and see whether you know whether truth goes to the or the approximation goes to the truth and so on. The second thing which I am going to look at the second benchmark which I am going to look at is dx by dt is equal to A x, x belongs to Rn and x tn is equal to xn and A here is a n cross n real valued matrix real valued in the sense or containing all real entries I am taking a shortcut and A is a real scalar ok. What is the true solution here for the first case if I am given x0 what is the true solution the true solution x star t let us call it let us call it stars as a true solution x star t is e to the power At x0 if it is x at tn it will be xn ok whether it is x0 or xn does not matter this is the true solution this is the true solution ok. What about this case if you are attending the other course this must have been done analytical methods what is the true solution here e to the power ok if you are not done it I will derive it in the class not an issue. So, the true solution here well what she is saying also correct but the true solution can be written here as x star t is e to the power At where A is the matrix I suppose you have been introduced to this e to the power exponential of a matrix no yes or no how many of you do not know this ok we will talk about it not an issue into x0 ok I derive this not an issue. You can write this as e to the power At x0 I derive it for a special case but the result holds for general case and is that I derive this for a case where you have Eigen values or Eigen vectors are linearly independent. So, we can derive this very easily other cases it requires little more work ok. So, this is my true solution and now I can use this as a benchmark to compare whether my approximate solution is going to the true solution or not ok. Let us not worry right now about the multivariable case let us start looking at the scalar case it is easier to get understanding of or to understand the convergence behavior. So, this will visit later ok. So, my true solution is this ok typically when I am using some numerical integration or some solver like Rka 45 or whatever I want to go from xn to I want to go from xn to xn plus 1 right. So, I will write this solution in that format ok. So, x star tn plus 1 let us take right now step size to be constant. Now, I am not worried about variable step size ok let us take constant step size step size h ok x tn plus 1 is e to the power a tn plus h x naught right and x star tn is e to the power a tn x naught is this ok. I can combine these two and say that x star tn plus 1 is equal to e to the power a h x star tn this check is everyone with me on this I am just expanding this substituting for I have just expanded and substituted this. So, these two equations I can combine to write this ok in our notation this is same as x star n plus 1 is equal to e to the power a h the notation that we have adopted with that this is x star n plus 1 is this e to the power a h x star n this is the true solution how it should actually evolve ok. Now let us look at explicit Euler method first ok how will you write explicit Euler method explicit Euler method is x n plus 1 now when I am writing just x no star this is approximate solution x star is the true solution x is the approximate solution ok. Now explicit Euler explicit Euler ok what is x n plus 1 x n plus 1 is at least today you should know because you are supposed to submit Euler's method and h times today or what is the submission did fourth ok. So, but it is too close right fourth is not too far from. So, I can expect that you should know what is Euler method today. So, h times f n right but what is f n in this case a x n. So, this I can write here x n plus h a x n right. So, I can write this as 1 plus a h x n everyone with me on this ok. Now let us define an error let us define error which is e n this is x star n minus x n ok x star n minus x n this is my error. So, what is this error this is error between the truth and the approximation ok. Now if I take this as my equation b and if I go here and call this equation as my equation a then I can subtract the two equations and then get the difference equation that governs the error dynamics ok. So, a minus b equation a minus b will give me ok. Now do not confuse between this e exponential and or let us do like this to reduce the confusion I will call this epsilon I will call this error epsilon because you should not confuse with e to the power something ok this is what I get. Now I have to convert this into a difference equation ok I am going to play a little bit of a trick I am going to substitute for x star n ok in terms of epsilon plus x n and going to do some read read adjustment ok. So, this equation this equation let us call it c and this equation let us call it d. So, I am going to combine c and d and get a difference equation combining these two you get e n plus 1 is equal to 1 plus a h plus e to the power a h minus 1 plus a h ok 1 plus a h sorry I forgot about e n I forgot about e n ok. What I have done is I have eliminated x n on the right hand side of that equation I have written x star n I have written x star but I have removed all that I have done is I have written x n is same as x star n minus epsilon n. If I just use this this is my definition this is my definition I am just using this definition substituting and rearranging that equation I get this equation ok I get this equation ok. So, this gives us some idea about how the error behaves how the error behaves it also gives us some insights we will look at the insights and then move on and do some more analysis. But this one simple insight which you will get here is that error dynamics ok will depend upon this difference you see something familiar here 1 plus a h how do you expand e to the power a h 1 plus a h plus a square h square by 2 factorial plus whatever. If h is very very small ok higher order terms can be neglected and this difference will be very very small ok. If this difference is very very small then if this difference is very very small what will govern the error dynamics e n plus 1 is equal to 1 plus a h right 1 plus a h how 1 plus a h behaves will decide ok whether the error goes to 0 whether error goes to infinity or whatever happens. If this is very very small if this can be neglected ok can you just look at this term and say what should happen this is a linear difference equation scalar linear difference equation what is the Eigen value 1 plus a h 1 variable there is no matrix is only a scalar how will the difference equation behave how will the difference equation behave if you if this if you can neglect e to the power a h minus 1 plus a h this is approximately equal to 0 this can be neglected then you can say that e n plus 1 is equal to 1 plus a h e n ok 1 plus a h e n let us make a simplifying assumption just to get our understanding correct ok let us assume that a is strictly less than 0. If a is strictly less than 0 how should the true solution behave what is the true solution true solution x star t is e to the power a t x naught what should happen to x star t as t goes to infinity if a is less than 0 exponential minus something will exponentially exponential decay it will go to 0 ok if a is less than 0 ok when will error go to 0 how will this equation behave see this equation this equation will behave as e n you can derive what I am writing here very very easily e n will be 1 plus a h rest to n e 0 rest to e 0. So what will happen to 1 plus a h when will it go to 0 when will this go to 0 this quantity should go to 0 error should go to 0 right that difference between the truth and approximation should go to 0 when will it happen louder I cannot hear no something more than that minus 5 is less than 1 plus less than 0 then minus 5 between more of 1 plus a h not 1 plus a h more of 1 plus a h should be less than 1 the way you should put more of 1 plus a h should be less than 1 ok. So this will go to 0 more of a h if more of a h 1 plus a h if this is strictly less than 1 what will happen e n will tend to 0 as n tends to infinity yeah. So you start with e 1 there will be a mistake committed one step you move right women to take one step there will be a mistake committed from x 0 is 0 is some initial point ok which need not be same as see because when if you start from the same point ok women to take one step ok the solution of approximate solution will be different from the true solution you start analysis from e 1 ok. So let us start from e 1 this will be n minus 1 yeah good question is 0 will always be 0 now let me start from even ok. So this will be n minus 1 right and this is even now even is different even if I start from same x 0 ok x 1 which is calculated using approximate solution is not same as x star 1 ok now let us do the analysis fine yeah good thinking I appreciate your e n minus 1. So this is is this fine so this has to go to 0 right otherwise we have we have trouble. So e n will go to 0 only when this is strictly less than 1 ok. So this being strictly less than 1 what does it mean in terms of choice of see this condition is giving you a way to choose h this condition is giving you a way to choose h for explicit Euler method. If you violate this condition the difference between the true and approximate will not vanish ok. If you violate this condition the difference between the true and but if you obey this condition what you know is that asymptotically the difference between the truth and approximate will vanish ok. So this is this is a circuit condition for explicit Euler can we derive an condition for implicit Euler what is implicit Euler just do it. So for implicit Euler now one thing one before you proceed I just want to point out one more thing here see look at this look at this expression here whether h is small or whether h is large it depends upon a what is whether whether see that this is very very important ok what is small h somebody might give you some number 10 to the power minus 6 is small h but that is not the thing whether h is small or h is large it depends upon what is value of a it is not. So it depends upon the characteristic equation of or characteristic constant of the equation a ok which appears in the equation when we go to the vector case this a here will get replaced by Eigen values ok will get replaced by Eigen values. So what is step size it depends upon the system it depends upon the system parameters a that is also an important take home message that it is not independent of what value of what value a has ok. Of course this you can write in terms of you can expand this and write this as inequality in terms of but now let us go more on to x implicit Euler method implicit Euler method is x n plus 1 will be f n plus 1 ok for the linear case you do not have to do iterations to compute the solution. So this is nothing but x n plus h a a x n plus 1 right ok. So I can write that x n plus 1 is equal to 1 upon 1 minus a h I just taking this one left hand side ok 1 minus a h x n right. If a is less than 0 if a is less than 0 1 minus a h this will always be positive your integration step size is positive ok 1 minus a h is always positive. So 1 plus a positive number this is always going to be a fraction ok this is this is very nice ok. Now if I am doing analysis for this method how will how will this equation change what will appear here 1 upon 1 minus a h 1 upon 1 minus a h ok there are things called paray approximants and this is a paray approximation a type of paray approximation of e to the power a h 1 upon 1 minus a h is a type of paray approximation this is a better approximation than Taylor series and truncating ok. So now again if this is very very small you could just look at this ok you could look at this and the equation for error will change to see this is explicit Euler this is my implicit Euler ok a is less than 0 when a is less than 0 ok the solution should the true solution asymptotically goes to 0 e to the power minus a t that asymptotically goes to 0 what can you say about the error in this case what is our criteria for convergence this coefficient should be less than 1 strictly less than 1 right for implicit Euler will this always be less than 1 for any choice of h you see why implicit Euler is better even if you make a wrong guess of h you are still you know you do not have to worry whereas here you have to be very very careful explicit Euler you better be careful about how you choose your integration step size implicit Euler when if you make a wrong choice eventually Euler will go to 0 ok that is why implicit methods are many times perform much better than the explicit methods yeah if a is positive no matter what will happen if a if a is positive it is a difficult problem because the solutions are going to infinity so matching two infinities both will go to infinity actually but the difference between the two infinities can grow below so difficult to get understood insights into what is happening if you make an assumption about if you say something about if you take a to be greater than 0 that is why I have taken the case where a is less than 0 we basically want to get insight into how Euler behaves this if the the system is truly unstable system ok then getting the true unstable solution is a difficult problem conversionally not an easy problem even a small difference can blow up in this case not here why so I have to I have to make sure in this case so I have to make sure that minus 1 is less than 1 plus a h is less than 1 right so I have to choose my h carefully so that this inequality is not violated can you do analysis for the trapezoidal rule just do it see what what expression that you get see trapezoidal rule is xn plus 1 is equal to xn plus h by 2 fn plus fn plus 1 this is the trapezoidal rule ok what will be equivalent difference equation for E n if I just if I do the algebra ok I will get this thing here is just check whether you are getting the same thing 1 plus h by 2 into a divided by 1 minus h by 2 a ok and again the same same thing we can say no we can argue that if this 1 plus h a by 2 upon 1 minus h a by 2 if this is close to 0 ok you could only look at how well I will do the I will write down the exact equation when it is not close to 0 but right now this is approximate analysis just so in that case in that case this is 1 plus h a by 2 this is 1 minus h a by 2 and the condition becomes more of this should be less than 1 so this is this is trapezoidal rule or Simpson's method or whatever ok this is trapezoidal rule or Simpson's method and in this case it will become right in fact this ratio is always going to be a fraction the ratio is always going to be a fraction why if a is of course if a is less than 0 then this ratio is going to be a fraction this is always going to be higher than this ok so we are assured that this is a fraction so trapezoidal rule will converge the error will you know error will go to 0 asymptotically will become smaller and smaller ok actually the way you should do analysis is I am doing a very hand-wavy kind of analysis the way I should look at this problem is I will write it down a little more I will write an exact expression the exact expression I should write is E n plus 1 E n plus 1 X star n plus 1 is equal to 1 plus a h by 2 1 minus a h by 2 e to the power a h minus 1 plus a h by 2 1 minus a h by 2 0 and e to the power a h well right now when I am writing this equation I am not doing any hand-waving business I am not saying this is close to 0 and all that ok this is a linear difference equation ok this is b sorry this is my let us call this z let us call this vector as z z n plus 1 is equal to if you call this matrix as b and if you call this as z n this is a linear difference equation z n plus 1 is equal to b times z n ok what is the criteria for convergence spectral radius of this should be equal to you can find out the spectral radius of this is a triangular matrix ok and then you can argue about the so basically we are using this fundamental concept of analyzing qualitative behavior of linear difference equations in this context to understand how the error goes to 0 same idea in fact those of you who will probably continue to work in area of digital control will find the same idea you come back comes back there when you talk about you know design of controllers digital controllers they will there again you will use spectral radius of the closed loop dynamic should be inside unit circle same idea linear difference equations that is the fundamental applications are different iterative processes or analyzing behavior of integration methods or you know solving some control problem the same idea in the context of iterative methods we had k plus 1 as b times k where k was iteration index here n is n is time ok and my condition still is spectral radius of b should be strictly less than 1 ok one thing you can you can see is that if a is less than 1 if a is less than 0 ok this part is going to go to 0 anyway asymptotically because x star is going to 0 ok everything will depend upon this and this of course this is going to 0 a is less than 0 so e to the power minus a h is going to 0 so everything finally boils down to this this term if this term is strictly less than 1 approximation error will go to 0 asymptotically ok this is what gives us insight into what happens now when I extend this to multivariate case eigenvalues of a matrix will appear and then we have to analyze them we will see that later.