Solving N-Bodies Problem with Leap Frog & Runga-Kutta Methods

[Zaheer Ahmed]

Hello,

My Name is Zaheer Ahmed, I am Doing Masters In computer Engineering, I am currently working on Comparison of Methods to solve N-Bodies Problem, Firstly i Start with Sun-Earth problem, using Leap Frog Method, the behavior of Earth was not correct, as it was not making complete rotation, then i tried Runga-Kutta 4th order Method, againg the behavior of Earth was not Correct, it was making semi-circle and coming back to the original position, i select DT (time step) as 0.1,0.01,0.001 and run the program but the result was same, then i moved to Earth-Moon problem, the result is again same, either i am doing some calculation mistake or i am not using the correct procedure to perform the simulation. When i am taking DT (time step) as 1 and running the program for 100000 steps, the result is correct but when i am entering the value (604800) time for one week it is giving error. Can Anyone Advice Me, where i am going wrong, as it is really important for me as i have to submit this project by the end of this month. I can provide you the results for my program, if anyone willing to help me plese reply me on zaheer.mujtaba at gmail.com. Waiting for your reply. Zaheer Ahmed

 

[SpaceTiger]

It's hard to say what's wrong without seeing your code, but the leapfrog method is the exact equation of motion, so there's no reason it shouldn't work. As for the Runge-Kutta problems:

it was making semi-circle and coming back to the original position

That sounds more like a sign error than a problem with the integration method.

 

[Zaheer Ahmed]

Hello,

I can provide you the source codes of my program if you give me your email address. Waiting for the reply. Zaheer Ahmed

 

[tony873004]

Your time step is not the problem. I can simulate in Earth around the Sun with timesteps as large as several hours and have it complete millions of closed orbits. My program (you can get it here: www.gravitysimulator.com ), uses a simple routine, but I'm not sure what it's called. Maybe SpaceTiger will know . Euler first order maybe? But it basically works like this:

for every object in the simulator's universe compared to all other objects

compute distance using pothag
acceleration = GM / d^2
velocity = velocity + acceleration
position = position + velocity

repeat indefinately.

I'd have to agree with Space Tiger that you have a sign error if you're doing a half orbit, and then jumping straight back. Does it do this indefinately?

Something you've probably already figured out, but sure to screw you up if you havent: when you multiply your velocities by your timestep, you have to multiply your masses in the acceleration formula by the timestep ^2.

I wrote a Runge-Kutta 4 once for this program, but never debugged it enough to use it for multiple objects. It could only do a 2-body ), where a massless planet circled a stationary Sun (or would this be called 1 body since the planet's effect on the Sun is never computed). I wasn't that impressed with the results compared to the way I do it now. I could do a timestep twice as fast as the current method before the orbit fell apart, but it took it a lot longer to do that timestep.

 

[enigma]

Welcome to PF, Zaheer!

It's possible you have a cosine issue somewhere.

cos^-1 is the same in quadrants 1-2 as it is in 4-3, and in my orbital dynamics code I was constantly making mistakes ensuring I was operating in the correct quadrant to the point where I wrote my own arccos function which took quadrant into account.

 

[DuncanM]

Hello, Zaheer.

You mentioned several aspects, but said nothing about trying to simulate elliptical motion; after all, wouldn't a 2-body problem (Earth revolving around the Sun) simply reduce to plotting an elliptical orbit?
There are several programs available for solving such problems. For example, here is a numerical utility (it is just a numerical root-finder):

http://www.akiti.ca/KeplerSolver.html

The code is in Javascript, whose syntax is almost identical to "C", so translating it should be easy. In fact, you could easily write your own code using the bisection method (not elegant, but it works).
Adding more bodies should simply be a matter of tracking more elliptical orbits, with perturbations caused by the other bodies.

On the GravitySimulator.com "Links" page is listed the book, "Solar System Dynamics". I have this book and have found it to be quite good. There are others I could recommend too.

Also, the Satellite Toolkit is available for free to students, schools, and aerospace professionals. You might want to get a copy. It would probably have several powerful tools that would assist you in your work (if nothing else, to confirm the results of your own work).

 

[DuncanM]

Here's another good book:

"An Introduction to the Mathematics and Methods of Astrodynamics, Revised Edition"

by Richard H. Battin, Senior Lecturer at MIT

It is part of the AIAA Education Series
J. S. Przemieniecki, Series Editor-in_chief

 

[NEOclassic]

Wellcome Zaheer and Duncan to PF

Hello, Zaheer.
You mentioned several aspects, but said nothing about trying to simulate elliptical motion; after all, wouldn't a 2-body problem (Earth revolving around the Sun) simply reduce to plotting an elliptical orbit?

Isn't it so that classical mechanics tells us that in the situation of a single planet orbiting around a massive body (2-body problem) that the centers of both bodies circle in the same direction about an axis somewhere between those centers (and perhaps very much closer to the more massive body) but that the distance between those centers remains constant since the gravitational attraction between them is constant only if the distance between them is constant. We also know that tangential momentum must be conserved and further that the angular momentum can only be conserved if the distance between the bodies remains constant which has already been stipulated. In short the orbit of the single planet must be circular relative to the center of the massive body - clearly the addition of a second planet to the 2-body problem could perturb the circular orbit of the original planet etc, etc for the n-body situation that is evidenced by the "ovate" (not elliptic) behavior of the solar system's reality. Quod Erat Demonstrandum. Cheers, Jim

 

[SpaceTiger]

Isn't it so that classical mechanics tells us that in the situation of a single planet orbiting around a massive body (2-body problem) that the centers of both bodies circle in the same direction about an axis somewhere between those centers (and perhaps very much closer to the more massive body) but that the distance between those centers remains constant since the gravitational attraction between them is constant only if the distance between them is constant.

The force need not remain constant between the two bodies. Rather, the position of the total center of mass must remain constant in some inertial frame. Elliptical motion satisfies this condition and is part of the generic solution to the two-body problem.

 

[Zaheer Ahmed]

Hello,
Thanks for your precisious messages, I found the mistake, i was entering the mass of sun as 41.989e30 Insted of 1.989e30 kgs, Now after correcting the mistake my simulation results shows that Earth is making one complete circle and in the form of an ellipse, this results are same for both Runga-Kutta 4th order as well as Leap Frog Method, But the problem is " Earth is making one complete circle in 23 weeks, where as when i am running the simulation for 52 weeks, it is making more that two circles. Can Anyone Tell me what could be the problem. The distance between Earth and Sun varies from 1.49 e11 to 1.45 e11. Please Can anyone tell me where i am going wrong, waiting for your precious suggestions. Zaheer A.

 

[tony873004]

What time step are you using? What are you initial Position (x,y,z) and Velocity(x,y,z) for Earth. Can you post your code?

Also, use 1.98911e30 for the Sun. This won't be the source of you error as it is too small, but later, when you get the big error sorted out, it will keep the Earth from drifting over several decades.

 

[Zaheer Ahmed]

Hello Tony,

I am posting my source codes for Runga-Kutta 4th order method, Please Verify and let me know where i am going wrong. Thank you.
================================================================

#include<iostream>
#include<iomanip.h>
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<time.h>
#include<fstream.h>


// Total number of bodies
const unsigned int N = 2;

const double G = 6.67e-11;
const double PI = 3.14159265358979;

double mass[N];
double position[N][2];
double velocity[N][2];
double force[N][2];
double sim_time;

const int NWEEK = 53;
double DT = 1;

// Function Prototypes
void moveBodies();
void calculateForce();
void computeForces(int index, float);
void write2File(char* filename);

typedef struct _runga{
	double matrix[4][2];
}runga;

runga deltap[N];
runga deltav[N];
runga deltaf[N];


int main()
{
	unsigned long i,j,k;
	time_t			start,end;
   struct tm *tmstart, *tmend;

	// Initializing the mass for the planets
	mass[0] = 1.98911e30;
	mass[1] = 5.97e24;


	//Initializing the position of the planets according to given data
	//Sun
	position[0][0] = 0.0 ; position[0][1] = 0.0;

	//Earth
	position[1][0] = 150e9			      ;   position[1][1] = 0.0;



	//Initializing the velocity of the planets according to given data
	//Sun
	velocity[0][0] = 0.0 ; velocity[0][1] = 0.0;

	//Earth
	velocity[1][0] = 0.0 ; velocity[1][1] = 2.9786e+004;

	//Initializing the velocity and force

	// Loop for week
	for(j=0;j<N;j++)
	{
		force[j][0] = 0.0;
		force[j][1] = 0.0;
	}

	//Initialization
	for(i=0;i<N;i++)
	{
		for(j=0;j<4;j++)
		{
			for(k=0;k<2;k++)
			{
				deltap[i].matrix[j][k] = 0.0;
				deltav[i].matrix[j][k] = 0.0;
				deltaf[i].matrix[j][k] = 0.0;
			}
		}
	}


	ofstream outClientFile("solarRungaKutta4.txt",ios::out);
	if(!outClientFile)
	{
		cerr << "File could not be opened" << endl;
		exit(1);
	}

	//srand((unsigned) time(NULL));
	// Starting time for the program
	start = time(&start);
 	tmstart = localtime(&start);
   cout << "The  simulation is started at: ";
 	printf("%ld/%ld/%ld  ",tmstart->tm_mday, tmstart->tm_mon + 1,
            tmstart->tm_year + 1900);
 	printf("%ld:%ld:%ld\n\n",tmstart->tm_hour,tmstart->tm_min, tmstart->tm_sec);
   outClientFile << "Starting time for the program:\t";
   outClientFile << tmstart->tm_mday <<"/"<< (tmstart->tm_mon + 1) <<"/"<<
   							(tmstart->tm_year + 1900)<<"  ";
	outClientFile << tmstart->tm_hour <<":"<< tmstart->tm_min <<":"<<
   							tmstart->tm_sec<<endl;
   outClientFile << "\nMass of Body 0 is: "<< mass[0]<< endl;
   outClientFile << "Mass of Body 1 is: "<< mass[1]<< "\n" << endl;
   outClientFile << setw(7)  << setiosflags( ios::left ) << "Body";
   outClientFile << setw(20) << setiosflags( ios::left ) << "X position";
   outClientFile << setw(20) << setiosflags( ios::left ) << "Y position";
   outClientFile << setw(20) << setiosflags( ios::left ) << "X velocity";
   outClientFile << setw(20) << setiosflags( ios::left ) << "Y velocity";
   outClientFile << setw(20) << setiosflags( ios::left ) <<
   																"Distance from Origin\n\n";


	for(k=0; k<NWEEK; k++)
	{
		cout << "Week: " << (k+1) << endl;
		for(i=0;i<604800;i++)
		{
			calculateForce();
			moveBodies();
		}

		for(j=0;j<N;j++)
		{
			outClientFile << setw(7) << setiosflags( ios::left ) << j;
			//outClientFile << setw(20) << setiosflags( ios::left ) << mass[j];
			outClientFile << setw(20) << setiosflags( ios::left ) << position[j][0];
			outClientFile << setw(20) << setiosflags( ios::left ) << position[j][1];
			outClientFile << setw(20) << setiosflags( ios::left ) << velocity[j][0];
			outClientFile << setw(20) << setiosflags( ios::left ) << velocity[j][1];

			if(j==0)
				outClientFile << setw(20) << setiosflags( ios::left ) <<
                  (sqrt(pow(position[0][0],2) + pow(position[0][1],2))) << endl;
			else
				outClientFile << setw(20) << setiosflags( ios::left ) <<
         			(sqrt(pow(position[j][0] - position[0][0],2) +
                  pow(position[j][1] - position[0][1],2))) << endl;
		}

		outClientFile << endl;
	}


	// Ending time for the program
	end = time(&end);
   tmend = localtime(&end);
   cout << "The  simulation is ended at: ";
 	printf("%ld/%ld/%ld  ",tmend->tm_mday, tmend->tm_mon + 1,
            tmend->tm_year + 1900);
 	printf("%ld:%ld:%ld",tmend->tm_hour,tmend->tm_min, tmend->tm_sec);
   outClientFile << "Ending time for the program:\t";
   outClientFile << tmend->tm_mday <<"/"<< (tmend->tm_mon + 1) <<"/"<<
   							(tmend->tm_year + 1900)<<"  ";
	outClientFile << tmend->tm_hour <<":"<< tmend->tm_min <<":"<<
   							tmend->tm_sec<<endl;

	sim_time	= difftime(end,start);
   outClientFile << "\nThe time taken by the simulation is: " << sim_time <<
   						" seconds" << endl;
	cout << "The time taken by the simulation is: " << sim_time <<
   						" seconds" << endl;


	//Closing the file
	outClientFile.close();

	return 0;
}
/************************************************************************/
void calculateForce()
/************************************************************************/
{
	double distance, magnitude, temp;
	double direction[2];
	long i, j;

	for(i=0;i<(N-1);i++)
	{
		for(j=(i+1);j<N;j++)
		{
			temp = pow(position[i][0] - position[j][0], 2) + pow(position[i][1] - position[j][1], 2);
			distance = sqrt(temp);

			magnitude = ( G*mass[i]*mass[j] )/temp;

			direction[0] = position[j][0] - position[i][0];
			direction[1] = position[j][1] - position[i][1];

			force[i][0] = force[i][0] + magnitude*direction[0]/distance;
			force[j][0] = force[j][0] - magnitude*direction[0]/distance;
			force[i][1] = force[i][1] + magnitude*direction[1]/distance;
			force[j][1] = force[j][1] - magnitude*direction[1]/distance;
		}
	}
}
/************************************************************************/
void moveBodies()
/************************************************************************/
{
	int i,j,k;

	for(i=0;i<N;i++)
	{
		//Position
		deltap[i].matrix[0][0] = DT*velocity[i][0];
		deltap[i].matrix[0][1] = DT*velocity[i][1];

		//Velocity
		deltav[i].matrix[0][0] = DT*force[i][0]/mass[i];
		deltav[i].matrix[0][1] = DT*force[i][1]/mass[i];
	}

	for(i=0;i<N;i++)
	{
	
	 //Position
		deltap[i].matrix[1][0] = DT*( velocity[i][0] + 0.5*deltav[i].matrix[0][0] );
		deltap[i].matrix[1][1] = DT*( velocity[i][1] + 0.5*deltav[i].matrix[0][1] );

		//Velocity
		computeForces(0,0.5);
		deltav[i].matrix[1][0] = DT*deltaf[i].matrix[0][0]/mass[i];
		deltav[i].matrix[1][1] = DT*deltaf[i].matrix[0][1]/mass[i];
	}
	for(i=0;i<N;i++)
	{
	
		//Position
		deltap[i].matrix[2][0] = DT*( velocity[i][0] + 0.5*deltav[i].matrix[1][0] );
		deltap[i].matrix[2][1] = DT*( velocity[i][1] + 0.5*deltav[i].matrix[1][1] );

		//Velocity
		computeForces(1,0.5);
		deltav[i].matrix[2][0] = DT*deltaf[i].matrix[1][0]/mass[i];
		deltav[i].matrix[2][1] = DT*deltaf[i].matrix[1][1]/mass[i];
	}
	for(i=0;i<N;i++)
	{
	
	//Position
		deltap[i].matrix[3][0] = DT*( velocity[i][0] + deltav[i].matrix[2][0] );
		deltap[i].matrix[3][1] = DT*( velocity[i][1] + deltav[i].matrix[2][1] );

		//Velocity
		computeForces(2, 1.);
		deltav[i].matrix[3][0] = DT*deltaf[i].matrix[2][0]/mass[i];
		deltav[i].matrix[3][1] = DT*deltaf[i].matrix[2][1]/mass[i];
	}


	for(i=0;i<N;i++)
	{

		position[i][0] = position[i][0] + (deltap[i].matrix[0][0] + 2*deltap[i].matrix[1][0] + 2*deltap[i].matrix[2][0] + deltap[i].matrix[3][0])/6;
		position[i][1] = position[i][1] + (deltap[i].matrix[0][1] + 2*deltap[i].matrix[1][1] + 2*deltap[i].matrix[2][1] + deltap[i].matrix[3][1])/6;


		velocity[i][0] = velocity[i][0] +  (deltav[i].matrix[0][0] + 2*deltav[i].matrix[1][0] + 2*deltav[i].matrix[2][0] + deltav[i].matrix[3][0])/6;
		velocity[i][1] = velocity[i][1] + (deltav[i].matrix[0][1] + 2*deltav[i].matrix[1][1] + 2*deltav[i].matrix[2][1] + deltav[i].matrix[3][1])/6;

		force[i][0] = 0.0;
		force[i][1] = 0.0;
	}

	// Equating the force to zero
	for(i=0;i<N;i++)
	{
		for(j=0;j<4;j++)
		{
			for(k=0;k<2;k++)
			{
				deltaf[i].matrix[j][k] = 0.0;
				deltaf[i].matrix[j][k] = 0.0;
				deltaf[i].matrix[j][k] = 0.0;
			}
		}
	}
}

 

[tony873004]

I think if you put code tags around your code then it keeps your indenting. Without it, its kinda hard to tell which { goes with which }

It's still hard to decipher someone else's code, but here are a few comments:

// Loop for week
for(j=0;j<N;j++)

I think this is just a comment error. Should be Loop for number of bodies?

temp = pow(position[0] - position[j][0], 2) + pow(position[1] - position[j][1], 2);
distance = sqrt(temp);

magnitude = ( G*mass*mass[j] )/temp;

This looks a little suspicious. You used temp as a temporary value in your pothag formula so you could add the squares before you square rooted them to solve for distance. But now you're using temp in the magnitude formula instead of the distance variable you computed? Also, don't you want to square that since
[tex]F= \frac{GMm}{distance^2}[/tex]

magnitude =  ( G*mass[i]*mass[j] )/distance^2;

or for faster execution time:

magnitude =  ( G*mass[i]*mass[j] )/(distance*distance);

 

[Zaheer Ahmed]

Hello Tony,

First of all thanks for your reply, Let me explain why I wrote "Loop for Week" in my program, with this I am selecting the simulation period For Example 13 weeks or 26 weeks like that. Secondly, yes I am using temp as a temporary value in my formula so that I could add the squares before I square rooted them to solve for distance. But look when I am calculating the magnitude I am using magnitude = (G* mass * mass[j])/temp; Where Distance = sqrt(temp);
Anyway yes what you are telling is right, "It's still hard to decipher someone else's code" But I appreciate you for coming forward and helping me with your valuable advice. Any suggestions, According to the results generated by my program, Earth is revolving around the Sun in 23 weeks as against 52 weeks. What could be the reason for this kind of behavior; if you have any idea please let me know. Thank you. Zaheer A.

 

[Zaheer Ahmed]

I think the absence of other planets matters while simulating 2-Bodies problem

Hello,

I was waiting for the valuable suggestions from your sides, Anyways what I think is when I am simulating the Sun-Earth Problem, Earth is compl- eting one revolution in 23 weeks as against 53 weeks, this is happening be- cause the distance between Sun and Earth is gradually decreasing and one more thing I notice is when the Earth is entering the second quadrant and third quadrant it is covering the distance very fast as compared to the first and fourth quadrant. What I think is because of absence of other planets, the velocity is increasing something like that, Still confused and need your ideas, Comments from your side are please welcome. Zaheer A

 

[pervect]

Your program is obviously still buggy. I would suggest using your RK4 method to solve an easier problem first, such as linear first order differental equation with one variable, dx/dt + x = 0.

With proper programming, it is possible to separate out the equation(s) that you are solving from the procedure that does the integration. You can find one way of doing this in the book "Numerical Recipies in C". This book includes the RK4 method, it also includes something I didn't notice in your code, timestep control. It has complete source code for the problem you are working on. [Clarification - complete source code for the RK4 method with associated timestep control and error bounds control.] Check it out at your local library - and last I looked, you can even find portions of the book on the internet.

 

[tony873004]

...What I think is because of absence of other planets, the velocity is increasing something like that, Still confused...

The absence of the other planets will have no effect on your simulation. Without the other planets, this is just a solveable 2 body. Your Earth should orbit in exactly 1 year.

If you included the other planets, you'd have to run your simulation for thousands of years (sim time) to see very small changes such as fluxuating inclination and semi-major axis.

...it also includes something I didn't notice in your code, timestep control...

I think DT is his time step variable. It is hardcoded as DT=1.

Zaheer, is it possible that you are changing this to a faster DT when you execute the code? This could be causing a problem because I believe the masses in your force formula you would have to use M*DT^2. But if you're using DT=1, this won't make a difference. It's probably not your problem since that would give you bigger slower orbits instead of the smaller faster ones you're getting. Just my guess.

 

[Zaheer Ahmed]

My Simulation Results Shows that the presence of other planets matters

Hello,

My Simulation results shows that the presence of other planets matters, as I performed the 2-bodies experiment for 53 weeks and found that though Earth was rotating around the sun in an elliptical orbit, but it was making more that two circle’s in 53 weeks time, where as I performed the experiment considering nine planets and the results shows Earth is revolving around the sun in 53 weeks time, this says that the absence of other planets matters, Any questions or any suggestions are please welcome. Thank you. Zaheer A.

 

[tony873004]

Hi, Zaheer

There's no way that the other planets presence or absence is going to effect your result in any meaningful manner, such as shortening your period by half.

Kepler's 3rd law used to compute the period of a planet has no inputs for additional planets. All the other planets combined are only 1/1000 of the mass of the Sun. They are insignificant in the short term (< 1000s of years).

You can try it in my simulator. Start with the "9 planets simulation", verify that Earth takes 1 year to circle the Sun, then delete the other 8 planets and try it again. It will still take 1 year.

I think your error may lie in the way you are adding your numbers. You may be accidently advancing the Earth's position twice per time step.

Next week, after my finals are over, I'll look closer at your code. Unfortunately, I had no luck opening it in my VC++. I think you are very close and just the order of things are causing you problems.

*** edit ***

One way you may want to verify your findings is to set the masses of the additional planets to something very low, like 1 kg. Then see if they still affect Earth.

 

[Zaheer Ahmed]

It is showing correct behavior for 3 or more bodies But with 2-bodies it seems a prob

Hello,

My program is giving correct results for 3 or more bodies, for example i consider "Three bodies problem" as well as "solar system", and when i plot the graph, Earth is making one complete rotation around the Sun in 53 weeks, where as when i am considering 2-bodies, Earth is making more than two rotations in 53 weeks, it seems some where i am doing mistake, let me work on this and try to rectify my mistake. Thank you, Specially to you Tony . Zaheer A