 Welcome back, in the previous segment we discussed the basic cosmological simulation problem and we discussed the first order Euler method for it. In this segment we will talk about a second order method which is going to be a lot more accurate. So the idea is simple, instead of using two terms of the Taylor series we are going to use three terms. So f of t plus delta f is the variable we want to estimate at time t plus delta we are given the value at time t. So f of t plus delta is f of t plus f dash of t times delta and in the previous method we had stopped at this point but here we are going to add one more term. So plus f double dash of t times delta squared upon 2. This is still approximately equal to because as we know the Taylor series just keeps going and we are not taking all the terms of the series, we are just taking powers, second powers of delta and therefore by the way it is called a second order method. Now this can be written in a different manner. So let us look at this part, just this part. So this is what I have written down over here and I am going to rewrite it by taking delta common because delta is there in both. So this I brought out this delta so here I just get f dash t and here I get f double dash t times delta by 2. Now if you look at this you may observe that this is actually a quantity which can be written differently. So in what way? Well I can write it this way, why is that? So this is just a first order approximation of this. Now well if I want to know the value of a variable at t plus something, if I want to know the value of a variable at t plus something, how do I get it? Well I write it as the value of the variable at t then plus the derivative of this plus this something in the multiplication. So indeed in these parentheses what I have is the first order approximation of this. We are going in the other way of course that we are saying so we are going to substitute this for this over here. So instead of this we are going to get this term. So here is the alternate statement of the second order method. So we have the same thing f of t plus delta is equal to the first term does not change but this second term looks like this because this taken a first order approximation is exactly what this is. This alternate statement is nice in many ways. Why? First of all there are only two terms so that is good. Second the derivative, the second derivative has gone away. So it is a simpler thing and it almost looks like the first order method. The only difference is that in the second term we had t over here f dash of t times delta. So here we have t plus delta by 2 so that is the only difference. But this is still essentially the same expression and therefore in terms of error this is going to be as good. But because it is simpler we are going to use this in our algorithm. So we will use this form. But before we get to that I want to look at what this is geometrically just so that you get a better intuition of what is going on over here. So again let me draw a picture. So this is time and say let us say this is t and let us say this is t plus delta and on this axis we have f. So let us say so I am plotting the curve let us say something like this. So this is a plot of f. Now what does the first order approximation say? The first order approximation says that the increment that you are going to get is going to be the distance. So this distance is delta so the time difference which is delta times the slope and what is the slope? The slope of the tangent. So the slope of the tangent is going to be something like this. So this whole thing is going to be f dash at t times delta because this f dash is nothing but the slope of the tangent. So what is the second order term? This is the second order method saying. The second order method is saying that look do not take the slope at time t sorry do not take the slope at time t but take it at time t plus delta dash by 2. So take it somewhere over here so take it somewhere over here and now the slope looks something like this. So I want to multiply this slope by this delta and then add that up. So to do that I am just going to make it appear somewhere over here. So this slope I have just tracked it down and now I am going to multiply it by the delta and what has now happened is that the product is going to get me this. So this is this slope multiplied by this delta. This whole thing is this slope, see that this slope is going really sharp. So now it turns out that taking the slope in the middle is a good idea. In this picture also you can see that what I am getting to over here is much closer to the actual curve than this is and there is a different way of thinking about it and that is that to go from here to here if I take this so when I multiply this distance by the slope I am sort of assuming that the slope is going to be the same inside this entire interval. So if the slope was really the same in the entire interval then it would be just this straight line. In that case multiplying it would give me the exact answer but the slope is not the same. Should I take the slope at the beginning? Should I take the slope at the end or should I take the slope in the middle? Well the slope in the middle is likely to be closer to the average slope and therefore we should be taking the slope in the middle. So that is really the intuition behind this formula. It comes out from the Taylor series but the intuition is this that I want the slope, I want to multiply delta by the slope but I should pick a value which is sort of close to all the values in this interval. So this is t, this is t plus delta, I should pick a value of the slope which is close to all the values. So what are the values? Well the slope starts like this then it goes like this and it finally comes like this. So the slope in the middle which is this is kind of like this, it is closer to both and therefore we get less error. So we just did this, we got the geometric interpretation and this was just to informally see why this method should be doing better. So let us now try to apply this in the context of cosmological simulation. So we will do exactly like how we proceeded before. So let us say f is ri, so this f is ri and ri as we know is the position of star i. So what is f dash? Well f dash is the derivative of position which is the velocity. Again we are going to let t be 0, so what does this give us? So this tells us that the position at time delta is the initial position plus the velocity at the middle of the interval times delta. Now this may be another way to see why this new method is better because what we are saying over here is that we are multiplying delta by the velocity in the middle. So we are saying yes the velocity might be changing and we are hoping that the velocity in the middle is better representative of the velocities in the entire interval. The velocity gets lots of time to change from our selected velocity if we select it at the middle but if we select it, we select it at the beginning. If we select it in the middle then the initial velocity or the final velocity is not too far in time from the middle and therefore the middle is a better representative. So that is as far as the position is concerned. So if we know the position at time 0 and we know the velocity at time delta by 2 then we can calculate the position at time delta. Next let us see what happens if f is the velocity, so then the f dash is the acceleration and again let us take time equal to 0, sorry this time we are going to take time equal to delta by 2, you will see why. So what becomes of this, so t is delta by 2, so we are going to get delta by 2 plus delta. This is going to be v sub i of delta by 2, remember we are substituting v sub i for f and delta by 2 for t and f dash is acceleration and t is delta by 2, so this whole thing becomes delta times delta. So this tells us how we could calculate the velocity at time delta by 2 plus delta and for that we need to know the velocity at time delta by 2 and then the acceleration at time delta. Do you know the acceleration at time delta? Well as it happens we do. So this is the velocity term, this is for the acceleration what am I going to get, well for the acceleration we already know that it depends only on the position, oh I think I need a delta over here, so there is a delta over here which I have which I need to add. So this acceleration term is this term and this is the acceleration at time delta, do we know this term? Well we do if we know r sub i of delta and by the way when I say we can calculate r sub i I mean we can do it for all stars, so we can know it for rj and everything as well. So there is a delta term over here, so basically what has this led us to? So this says that if we know r i of 0 and v i of delta then we can calculate r i of delta and v i of delta plus delta by 2, let us just check again. So can we calculate r i of delta, well here is r i of delta, well do we have r i plus 0, yes we said if we know r i of 0 do we know v i of delta by 2, yes we said if we know v i of delta by 2, so yes so if we know these two things we will be able to calculate r i of delta, will we be able to calculate v i of delta plus delta by 2, well v i of delta plus delta by 2 we are going to update as v i of delta by 2, do we know this, well yes we said if we know v i of delta by 2 do we know these things, well all r i of delta we have calculated in this step over here, so therefore we know these terms as well and therefore we can calculate this. So what does this mean, what have we exactly accomplished, so I am going to draw a timeline over here, so this is 0, then this is delta, then this is 2 delta, then this is 3 delta and so on. So at this point we know r sub i and at this middle point we know v sub i and what that tells us that knowing these two things, knowing this and this we can calculate this and this. But now we can apply the same argument or let me write it down properly, we can calculate r i for this point and v i for this point. Now I can apply the same argument to these two terms, the same calculations, so I will be able to calculate the r i over here and the v i over here and I can keep going. So it is like the previous thing but the r i and v i values are not at the same time, at the same time instant. So they are at these staggered time instants. In fact what is going on is that the r i values are calculated at integral multiples of delta or multiples of delta by 2. So these are the r i calculations and these are the v i calculations. So you can fancifully think of the r i values sort of jumping across the v i values and in fact this method is also called the leapfrog method or it is also called the midpoint method. But anyway, so this is how the updates work, we calculate these two first, from these two we can calculate these two, from these two we can calculate these two and so on. So that is exactly our algorithm, what I am going to write down now is just this in a little bit more like code first of all. So we can repeat the steps to get things for 2 delta, 2 delta plus delta by 2 and so on. Well there is a little bit of a problem though, so we said if we know these two things, but do we know these two things? Well we do not, what was given to us were the positions and the velocities at time 0. So we know v i at this time and r i at this time. So our input specification said that we are given the positions and the velocity at time 0. We were not given this, so can we get this? Well by now you should be able to guess this, can we get this? Well given the velocity over here we can always get this if we know the acceleration, because that is and we can do this by a first order Euler expression, update. So to go from here to here we will use a first order update, but note that happens only once. So before we get the whole process started we are going to do this first order update to get the velocity at this point. So what is the expression needed over here? So we are going to have v i equals the original v i plus acceleration times delta. So that is the first order update, so that is what we are going to use here. So we were given these two things, from that we got the velocity over here and r i we already have and now we are going to use the method on this slide to get these, these, these and how are we getting these? So this is a second order update. This is also a second order update, so all of these are second order updates. So we do one first order update at the very beginning and subsequently we only do second order updates. So that is what the code or the algorithm is going to have. So as before we will read the simulation duration step size number of stars into t delta n. Then we will read the initial position velocity mass into those arrays and then this is going to be the first order update that we talked about. So calculating v i at this position. So at time delta by 2 and you can see that that is why we are multiplying by time delta by 2 and the values over here are at time 0. So these are values at time 0 and from that we are now getting v i at time delta by 2. And now we are just going to do this t by delta times and at the beginning of the loop r i will hold position at time s minus 1 delta, v i will hold velocities at time s minus half delta. So we are effectively at this point in our execution. So at this point we are we have advanced to this point. So r i is going to be over here we have not get got to the red values. So at this point r i is holding the positions at 0 and v i is going to be holding the positions at half but s is 1 over here. So this is 1 minus half so that is exactly the right thing. So exactly as we said we will calculated the updated values of the position and that will get us the position at this time. So from this we got to this time and then we will use the acceleration calculated earlier to go from this point to this point. So that will take us to acceleration oh sorry this should be no, so yeah so this is really the acceleration at s delta not earlier we are going to calculate the acceleration at this point. So we calculated the acceleration at this point and this acceleration we are going to add. So that will bring us to this point. So as a result our velocity will be available at this point which is s plus half times delta. That is it that is the entire algorithm. So as we said earlier the first order update is used only in this step and over here and over here we are using second order updates. So as a result the overall accuracy is high. This does turn out to be a good algorithm. So what have we discussed in this segment? So we have discussed the second order method and we have discussed sort of the details of the entire algorithm based on the second order method. In the next segment I will show you the program but we will take a short break.