 Hello and welcome to the NPTEL course on an introduction to programming through C++. I am Abhiram Ranade and in this lecture sequence I am going to talk about cosmological simulation and the reading for this is chapter 19 of the book. So let me begin by saying a few very very broad things. So you could say that the goal of science is to understand phenomena so that we can predict the future. Predicting the future is sort of a major test of science. And what this means? Well given the state of a system today we need to say what will happen to it tomorrow. In contrast the goal of engineering is to decide how should we interact with the system so that we get the state we want. And these questions are actually interesting in every field, pretty much every field. So in cosmology the questions would be how did the universe start and what happens tomorrow. In meteorology the question is what will the weather be like tomorrow and even in something like biology these questions are very important and for example these questions effectively are asking well how does an organism grow and even how can we make it grow well. The goal for today is to write a program that can predict the position of stars after some time given the position and velocities today. So this is the central question of cosmology and we are going to look into it. But the methods we use in principle are applicable to almost every other problem. So what is cosmological simulation? Well we are given positions velocities and masses of some n stars and we will use r sub i v sub i for the positions and velocities of the i th star and m sub i for the weight mass of the i th star. And since the positions and velocities are functions of time we may put 0 after r i and p i and I should point out that r i and p i are vectors. So they have 3 coordinates each. What do we want? Well we want the positions and velocities at some time t which we will specify. So we want to know how our universe has evolved and we also want a movie of this evolution exactly what happens to the stars. Well of course this requires us to know some physics and the main thing it requires is how one star exerts force on another star. So the force on star i due to star j has magnitude g times m i m j upon r j minus r i magnitude squared. Here g is a constant and let me explain what that r i minus r j means. So I said r i and r j are vectors so let me draw a picture of the vectors. So let this be any origin and let us say this is r i and let us say this is r j. So our mass m i is over here and our mass m j is over here and we really want this distance. So r j minus r i so this is r j going in this direction minus r i is really this vector. So this vector is r j minus r i and what is being done in the slide is that we are taking the magnitude of this. So magnitude of r j minus r i so we are dropping the direction and we are squaring it. So this is the famous inverse square law that the force varies as the product of the masses and inversely as the square of the distance between them. And what is the direction of the force? So the direction of the force is from star i to star j so it is an attractive force. So the force on star i due to star j is attractive so that the star i is going to be attracted towards m j. So we can write this as the unit vector in this direction. So the direction can be indicated by the unit vector from star i to star j and alternatively the unit vector may be written as this which is the vector divided by so it is r j minus r i which is the vector itself in the direction i to j divided by magnitude of r j minus r i. So this would be the unit vector because we are dividing by its magnitude. So now we can substitute this quantity into this so to say to get a vector. So what happens? So the force due to star j on star i as a vector is g times m i m j times r j minus r i upon r j minus r i cubed. So what has changed? So this was the magnitude to it we are multiplying by the unit vector. So now this becomes a vector quantity and its magnitude is still going to be this. So we now have the force due to star i due to star j on star i as this vector. Now in doing this cosmological simulation we want the total force on star i due to all other stars. So this is simply obtained by adding up this quantity for all j which are not equal to i. And in fact we are most interested in the acceleration so that will be this divided by the mass of i. So it is the same expression except that the mass of i has been divided out. So this is what we need to know in order to calculate the forces and the acceleration and in addition of course we know we need to know Newton's laws Newton's laws of motion. So or we need to know basic kinematics basic kinematics is that if there is a certain velocity then that velocity is the derivative or the velocity is the rate of change of the position and the acceleration is the derivative or the rate of change of the velocity. So here is an outline of this lecture sequence. So I am going to start by discussing the basic method or the first order method as it is also called. Then I will talk about an algorithm for cosmological simulation using the basic method. Then I will talk about a more accurate second order method and an algorithm based on that. After that I will show you the code for the second order algorithm and a demo and then I will conclude. So what is the basic method and this basic method is due to Euler. So the input to the basic method is the value of a state variable at time t. It could be position, it could be velocity, it could be temperature. The method really is applicable to pretty much anything not just cosmology. So and we want to know what the value is at time some t prime bigger than t. So the main idea is to approximate it, approximate the whole process using the derivative. So what does that mean? So it means that if f is the variable you want to estimate and you are given this value f of t. So we said that the value of at time t was given to you so f of t is given to you and you want to know the value at time t dash but let us say for a minute that we want to know the value at some t plus delta. Then the value at time t plus delta can be written as f of t plus f dash t the derivative times the time difference. So maybe I should draw a picture over here. So say time is going in this direction and I am plotting the variable in this direction. So suppose the curve of that f, the plot of f versus time looks something like this and so let us say this is t and let us say this is t plus delta. So we do not know this point but we want to know this, this is the point that we want to know, this is the point that we do know. So what the method is saying is that at this point you can approximate the curve by the tangent. So this is the tangent and then you are going to multiply the tangent by this. So this is delta and notice that the slope of the tangent. So f dash t is the slope of tangent at time t. So and the slope is nothing but this distance h divided by delta. So all that is saying is that h is equal to f dash of t plus sorry f dash of t times delta. So this h is f dash of t times delta. So what is this entire thing? This entire thing is this distance and that is f of t. So f of t plus f dash t times delta that is what we have written over here. So it just, this is just a simple thing which comes about by approximating the curve by a tangent at that same point and another way to think about it is that we have already talked about the Taylor series expansion. So this is just the Taylor series expansion to one term or in addition the increase is expressed in using one more term. So this is called the first order estimate and it is called so because delta appears with power 1. Now this expression turns out to be reasonably accurate for small delta and delta is often called the step size. So how do we, what do we do if we want to go all the way till time t dash, well we divide this interval into steps of small size delta that is step size and then take as many steps as needed to go from t to t plus t dash. So how many steps will be needed? Well we want to cover the distance t dash minus t so that many divided by delta that many steps are needed. So we repeat this t dash minus t upon delta times so that we will first get t of, we will get f of t plus delta then f of t plus 2 delta f of t plus 3 delta and so on substituting t plus delta in place of t over here and keep going in that manner. So this is the basic first order method and it requires us to know the derivatives at all points and I should also say that this applies to scalar or vector state variables. So f could be a scalar, f could be a vector. So what if f is a vector, well then we are taking the derivative of each of the three components separately. So how do we use this for our cosmology problem? Well let us say that f denotes r i, f f is r i or the position of star i then what is f dash? Well f dash is the derivative of the position so it is the velocity so it is v i. So of course these are functions of time but I have not put time here right now just to keep my notation simple. And let us say for simplicity that t is 0, so then what does this equation look like? Well so instead of f we are putting in r i the position of star i and t is 0 so this is just delta and then this f is just this r i and t is 0 and this is v sub i because we said that f dash is v sub i and the delta is there. So what have we said that if you want to know the position at time delta and given you know the position at time 0 then to the position at time 0 you just have to add the velocity times the whatever time interval that is going to elapse and this is going to be the velocity at the beginning of the interval but the point is that if the interval is small enough we can say that this is the velocity for the entire thing. So this is an approximation that is why we have put in this approximation approximately equal to sin over here but if the interval is small then this is not bad this will be essentially constant during that interval. So that is it that is what this expression that we had on the last slide means as far as the position is concerned. And it is reasonable for small delta as we just said because in a small time interval we can sort of assume that look the velocity will not change all that much. Now the velocity actually changes so we also have to apply to the velocity so we are now going to ask what if this f were to be representing the velocity well then what does f dash represent? So we said earlier that the derivative of the velocity is the acceleration so this is nothing but the acceleration. So what does that get us well it says that V dash of delta or the velocity at time delta is the original velocity so it is the original velocity plus whatever acceleration that might have accumulated in that time difference delta. So again it is kind of a reasonable thing. And now let us use what we know about acceleration well we mentioned earlier that the acceleration on of star i due to every other star is proportional to the mass of the star times the distance between them upon the square of the distance and of course the q of the distance well so this is a vector quantity and this is a scalar quantity so it really is mass upon square of distance in the direction in this direction. And in our earlier expressions these zeros were not there so time was not mentioned but here we are we want to apply this to time 0. So we actually know this let me point out that we actually know this quantity why because we are given the positions at times at time 0. So therefore we know the accelerations at time 0 we are also given the velocities at time 0 so therefore we can calculate this as well. So we know the velocity at time 0 so we can calculate this we know the velocity at time 0 and we know the acceleration at time 0 because the acceleration is this which depends only on the positions and we know the positions. So basically what has happened is that given the position and the velocity at time 0 we are able to determine the position and the velocity at time delta for all delta sorry for all i so for all the stars. So basically we can predict what happens to our system after delta time units and of course delta has to be a small number in order for this prediction to be reasonably accurate but once we can do it for a small number we can just keep on repeating it and then we will get we can get to whatever time we want. So here we got we went from 0 to delta but then we can start with this and apply these two updates and then we will get to t equals 2 delta and 3 delta and so on. So the first order algorithm is quite simple but I am going to state it nevertheless you probably many of you probably have understood what exactly we are doing here. So in the first step we are going to read the simulation duration step size and the number of stars into the variables t delta and n. Then for each star we are going to read its initial position and velocity and mass into these variables so these are arrays actually arrays or vectors but let us say arrays because we also have the physics vector right now so the ith element of the arrays will get the initial position of the velocity and the initial position and the velocity and of course the mass is not changing but the mass of the ith star will be put in the ith element of the vector m. And then we have to go from let us say the time 0 to some target time t so that is what this the simulation duration is so from 0 the current time to the target time t. So we have to take steps s going from 1 to t upon delta so we have to take t upon delta steps. So what happens in each step? So basically we are going to do what we said on the last slide. So here what I have done is I have calculated the acceleration of star i due to all the other stars and this is at the current time. I have put down instead of the subscript j I have put down index j because the mass of the jth star is now in m of j and similarly r of j minus r of i to get the difference between the vector difference the vector distance between the jth star and the ith star and this is of course the cube of the magnitude of the distance. So this we are going to calculate with whatever values we have in r i r j so these are at time 0. So this is just the acceleration at acceleration of star i at time 0 because in the first iteration this these are values at time 0 and in general you will see that they are going to be at time s minus 1 times delta. So s in the first round is going to be 1 and so this is indeed going to be in the first iteration it is going to be the acceleration at time 0. Then we are going to apply the position update rule. So the position update rule was that we take the current position add to it the current velocity times delta and we will get the new position at time delta ahead of whatever time we started. So this is going to give us position at s delta because at the very beginning over here we knew the quantity r for all stars at time 0 or in general at time s minus 1 s minus 1 times delta. So when we do this update we will know the quantity at time s delta and this is the velocity update. So the velocity update says take the acceleration at time s minus 1 delta or 0 in the first iteration and multiply it by delta and add it to the velocity at 0. So this is actually 0 or s minus 1 times delta. So this is exactly what I have written down in the previous slide but we have put it in a loop and we have sort of put it in terms of variables rather than formulae. So that is it, that is the entire algorithm. So at the end of this we will have r i and v i contain the velocities and the positions and velocities of the i star and at time t because we would have taken the right number of steps over here. Now the performance of this is actually not that great and let me just give you a clue to this. So even if I take two bodies let us say the sun and the earth. So what this ends up doing is that the earth will very quickly spiral out from the sun. So this is the sun, see this is supposed to be the orbit of the earth and if the earth is somewhere over here in the first step it will go a little bit outside then it will go further outside and by the time it comes over here it will have gone out substantially and within 2, 3 turns it would go out very far and this would not work at all. Or let me put it another way if I want this to work I would have to take that step size so small that it would take an excruciating amount of time to make the simulation go through any reasonable amount of time period. So that step size would have to be really small and if the step size delta is very small then that really means that the number of steps I have to take is going to be huge. So we are going to fix that though but let me just remind you what we have discussed. Well we discussed the basic cosmological simulation problem and we discussed the first order method due to Euler and next we are going to talk about a second order method but before that we will take a quick break.