More precise orbital simulation program

grav · #1 (**permalink**) 15-April-2007, 04:36 AM

I thought this might be the best place for this thread. Actually, I don't know much about how simulation programs are made or how precise they are, but I make a few of my own, and I've noticed that adding a couple of simple lines to the calculations for the gravitation of a two-body system can make them much more precise, up to millions of times as much, so that they can be run faster as well. It works best for masses of similar quantity, though, such as binary star systems and such, while a small mass orbitting a much larger one will remain about the same.

The general program lines for gravity might normally look something like this:

d=((x(z)-x(3-z))^2+(y(z)-y(3-z))^2)^.5
vx(z)=vx(z)-G*M(3-z)*(x(z)-x(3-z))/(d^3)*ti
vy(z)=vy(z)-G*M(3-z)*(y(z)-y(3-z))/(d^3)*ti
x(z)=x(z)+vx(z)*ti
y(z)=y(z)+vy(z)*ti

which is repeated in a loop as the orbits proceed. Here, d is the total distance between the two bodies each time the loop is run through, x and y are the coordinates of each body, vx and vy the vectors of speed in the x and y directions, and z is the designated number for the bodies, 1 and 2, which alternate through these calculations to aquire their orbits. G is the gravitational constant and M is the mass of the body, which is given here as M(3-z) since we are considering the mass of the other body, so for z=1, 3-z=2 and for z=2, 3-z=1.

Okay, now here is what I have found to be a surprisingly precise program for at least a two body system of similar masses:

x(z)=x(z)+vx(z)*ti-G*M(3-z)*(x(z)-x(3-z))/(d^3)*(ti^2)/2*M(z)/(M(z)+M(3-z))
y(z)=y(z)+vy(z)*ti-G*M(3-z)*(y(z)-y(3-z))/(d^3)*(ti^2)/2*M(z)/(M(z)+M(3-z))
d=((x(z)-x(3-z))^2+(y(z)-y(3-z))^2)^.5
vx(z)=vx(z)-G*M(3-z)*(x(z)-x(3-z))/(d^3)*ti
vy(z)=vy(z)-G*M(3-z)*(y(z)-y(3-z))/(d^3)*ti

in exactly that order. I'm not sure how many programmers we have on this forum, but I think one will like the precision this program provides. If anyone else has any more suggestions for creating a faster or more precise simulation program, please let me know.

tony873004 · #2 (**permalink**) 15-April-2007, 07:43 PM

Why are you dividing by d^3? Your acceleration formula should be GM/d^2.

Also, if you want to make it faster, you computer probably does d * d much faster than d^2. Try them both and compare.

grav · #3 (**permalink**) 15-April-2007, 09:44 PM

Quote:

Originally Posted by tony873004

Why are you dividing by d^3? Your acceleration formula should be GM/d^2.

I am using vectors for the acceleration in the x and y directions, so a_x=(GM/d^2)*(x/d)=GMx/d^3 and a_y=(GM/d^2)*(y/d)=GMy/d^3. How do you do it?

Quote:

Also, if you want to make it faster, you computer probably does d * d much faster than d^2. Try them both and compare.

Yes, that does make it slightly faster also. I've found that out while experimenting with different variations of writing the same program in different ways, along with a few other "tricks", such as only printing text periodically, since that seems to take an incredible amount of time. But while these can increase it somewhat as far as the program itself is concerned, the lines for the program in this thread can make it hundreds to thousands of times faster (for two-body systems of similar mass, anyway) by allowing one to increase the time increments with the same precision, or up to millions of times more precise for the same running time.

tony873004 · #4 (**permalink**) 15-April-2007, 10:54 PM

Quote:

Originally Posted by grav

I am using vectors for the acceleration in the x and y directions, so a_x=(GM/d^2)*(x/d)=GMx/d^3 and a_y=(GM/d^2)*(y/d)=GMy/d^3. How do you do it?

The same way you did it, at least in the Euler integrator. I just broke it into 2 steps, that's why I didn't recognize your d^3.

Here's a link to a thread on my forum where the Euler integrator source code appears. I compute all the accelerations and velocities first, then in another loop, I update the positions from the newly-computed velocities.
http://www.orbitsimulator.com/cgi-bi...1160874415/6#6
Also, sqrt(x) is faster than (x)^0.5. These little things make a big difference since they're inside a loop that will be executed thousands to millions of times.

grav · #5 (**permalink**) 16-April-2007, 01:57 AM

Well, I saw the program code you referenced. I like the way you saved time by running the loop for the coordinates separately so that it only needs to run through n times for the number of bodies rather than n^2. The thing is, though, I ran the configuration you are using, and it was way off compared to this one. Unless I put wrote yours incorrectly or something, I think you will be very pleased with the new configuration by itself, much less the improvement in the stability of the orbit. I had forgotten I originally experimented with the configuration (the order in which the calculations take place) to come up with the most stable orbit before I even tried to improve on that. The new lines will all have to be run together through the loop and so might appear that they would take more time, but you could probably run it up to hundreds of times faster and still gain a more stable orbit than that one seems to give, so I think you will like it.

I will rewrite the lines to be used and post it. Just to be sure I get it right, though, I have a couple of questions. How do you figure tM and the object mass. It appears that the increments of time are also figured in, but how exactly? And what is 398600440000000?

And just in case this turns out to be a million dollar deal or something (would be nice

), this thread constitutes a copyright.

grav · #6 (**permalink**) 16-April-2007, 02:48 AM

Quote:

Originally Posted by grav

Just to be sure I get it right, though, I have a couple of questions. How do you figure tM and the object mass. It appears that the increments of time are also figured in, but how exactly? And what is 398600440000000?

Okay, I saw on your forum that 398600440000000 is calculated as G times Earth's mass. But where is the time? I guess it must be figured in separately for the coordinates when setting the graphics, right?

grav · #7 (**permalink**) 16-April-2007, 05:50 AM

Well, Tony, I have to hand it to you.

It looks like you've got yourself a very good configuration after all. Perfect, even. I wrote it in wrong. You must have done some experimenting of your own or something. Mine has the same accuracy, but yours is apparently much more stable. Mine starts off just as precise, but then very slowly and steadily increases or decreases in size, so yours is actually better over the long run, and because of the way you've got your loops, faster even using the same increments for time. Apparently, by writing in the change in distance separate from the main loop, not only have you saved time in the program, but you've stopped any possible corruption of the orbits forthright. Guess I won't be making that million anytime soon

. So the only thing I've managed to improve upon is my own earlier program, but you've done it one better, so much thanks for that.

I'm going to keep trying to improve on it, though. Thanks for the tips.

As far as the timing in your program goes, then, about the only thing that might be improved upon to save more time is just the omission of a few lines, when figuring the acceleration. To run it just a little faster, you might try rewriting it as:

Code:

// CodeOptomizer.cpp : Defines the entry point for the DLL application.   
//   
   
#include "stdafx.h"   
   
BOOL APIENTRY DllMain( HANDLE hModule,    
   DWORD  ul_reason_for_call,    
   LPVOID lpReserved   
 )   
{   
    return TRUE;   
}   
   
int _stdcall RunLoop(int NumObjects, double * ObjMass, double * objx, double * objy, double * objz, double * objvx, double * objvy, double * objvz, double * objsize, long * CollisionObjectA, long * CollisionObjectB)    
{   
 if (NumObjects < 1) return 0; //Should be 1 or more   
   
 double tM;   
 double tM2;   
 double dx;   
 double dy;   
 double dz;   
 double D;   
 double f;   
 double fx;   
 double fy;   
 double fz;   
 int Collisions;   
 Collisions = 0;   
 for (int k=1;k<=NumObjects;k++) {   
  tM = ObjMass[k] * 398600440000000;   
  for (int j=k+1;j<=NumObjects;j++) {   
    dx = objx[j] - objx[k];   
    dy = objy[j] - objy[k];   
    dz = objz[j] - objz[k];   
    D = sqrt(dx*dx+dy*dy+dz*dz);   
    if  (((objsize[k] + objsize[j])/2) > D) {   
	Collisions = Collisions + 1;   
	CollisionObjectA[Collisions] = j;   
	CollisionObjectB[Collisions] = k;   
    }   
    if (tM > 0) { 
  f = tM/D/D/D;  
  objvx[j] = objvx[j] - dx*f;   
  objvy[j] = objvy[j] - dy*f;   
  objvz[j] = objvz[j] - dz*f;   
    }   
  tM2 = ObjMass[j] * 398600440000000;   
    if (tM2 > 0) {   
  f = tM2/D/D/D;
  objvx[k] = objvx[k] + dx*f;   
  objvy[k] = objvy[k] + dy*f;   
  objvz[k] = objvz[k] + dz*f;   
    }   
	
  }   
 }   
   
 for (int h=0;h<=NumObjects;h++) {   
  objx[h] = objx[h] + objvx[h];   
  objy[h] = objy[h] + objvy[h];   
  objz[h] = objz[h] + objvz[h];   
 }   
	
	
 //End of loop computations ------   
 return 1; //Finished successfully   
	
}

tony873004 · #8 (**permalink**) 16-April-2007, 07:12 AM

You'll have to tell me where you made the change.

The reason you don't see time step in my code is that it is handled external to this code. When a user increases the time step by 2, the program increases the velocities of the objects by 2 and their mass by 2^2. This saves me from having to do the conversion in the loop, hence saving a lot of time.

tony873004 · #9 (**permalink**) 16-April-2007, 07:27 AM

Check that you're using double precision variables. I find it odd that you say your program gives the same results for short simulations, but not long ones.

grav · #10 (**permalink**) 16-April-2007, 08:01 AM

Quote:

Originally Posted by tony873004

You'll have to tell me where you made the change.

I just changed these

Code:

f = (1/D) * (1/D) * tM;   
  fx = (dx / D) * f;   
  fy = (dy / D) * f;   
  fz = (dz / D) * f;   
  objvx[j] = objvx[j] - fx;   
  objvy[j] = objvy[j] - fy;   
  objvz[j] = objvz[j] - fz;

and

 f = (1/D) * (1/D) * tM2;   
  fx = (-dx / D) * f;   
  fy = (-dy / D) * f;   
  fz = (-dz / D) * f;   
  objvx[k] = objvx[k] - fx;   
  objvy[k] = objvy[k] - fy;   
  objvz[k] = objvz[k] - fz;

to these.

Code:

 f = tM/D/D/D;
  objvx[j] = objvx[j] - dx*f;
  objvy[j] = objvy[j] - dy*f;
  objvz[j] = objvz[j] - dz*f;

and 

 f = tM2/D/D/D;
  objvx[k] = objvx[k] + dx*f;
  objvy[k] = objvy[k] + dy*f;
  objvz[k] = objvz[k] + dy*f;

Quote:

Check that you're using double precision variables. I find it odd that you say your program gives the same results for short simulations, but not long ones.

I'm using UBasic, which gives me up to 2600 digits, but the more I use, the slower it runs. I generally run it to about thirty or fourty decimal places. I think the reason my program runs off like that is because I was finding for the coordinates within the main loop, instead of doing it separately like you are, as well as running a loop for each body individually, so the slight change in coordinates of one body with each run affects the next.

HenrikOlsen · #11 (**permalink**) 16-April-2007, 08:14 AM

Ah, you've just explained what I couldn't understand.

I got confused when one of you wrote that x*x is faster that x^2, when all the compilers I know would result in the same calculation for both, (unless you force 2 to be real or use an explicit power function)

Try to learn a compiled language, you'll be amazed that the speed increase, and possibly confounded by the lack of effect many of the tweaks such as the ones you're doing has, since the compiler does most of them already.

grav · #12 (**permalink**) 16-April-2007, 09:03 AM

Quote:

Originally Posted by HenrikOlsen

Ah, you've just explained what I couldn't understand.

I got confused when one of you wrote that x*x is faster that x^2, when all the compilers I know would result in the same calculation for both, (unless you force 2 to be real or use an explicit power function)

Try to learn a compiled language, you'll be amazed that the speed increase, and possibly confounded by the lack of effect many of the tweaks such as the ones you're doing has, since the compiler does most of them already.

Thanks. I had been looking around for another program that might run faster and/or more accurately, but I was experimenting with different variations of formulas for gravity, and Basic is all I know, so I was hoping I could just change the few lines of calculations for gravity from another program to my needs, which wouldn't require much language skill, except for the layout I prefer.

grav · #13 (**permalink**) 16-April-2007, 09:13 AM

Well, it looks like changing the orientation of the calculations only exchanged one problem for another. Instead of the size of the orbits steadily increasing or decreasing with time, they now steadily precess counter to their lines of motion around each other, although the size of the orbits themselves don't change. It seems like there might be some simple way to help resolve one or the other of these.

tony873004 · #14 (**permalink**) 16-April-2007, 07:46 PM

Quote:

Originally Posted by grav

...Instead of the size of the orbits steadily increasing or decreasing with time, they now steadily precess counter to their lines of motion around each other, although the size of the orbits themselves don't change...

My code does that if you use too large of a time step. Take Mercury for example. If you use a time step like 64K seconds, Mercury's orbit noticabely precesses. And it's precessing backwards from the direction it precesses in real life. If you slow the time step down, the precession gets less and less. If you slow it below 2k seconds, you've slowed it past 0, and it is now free to precess in the natural direction. As you slow it further, the rate of precession asymptotically approaches the correct Newtonian-predicted value of about 530 arcseconds per century. Here's a graph I made of time step vs precession in arcseconds per century.

Try running it to 15 decimal places. I doubt you need much more than that. 15 is what my program uses (Double precision). With a slow enough time step, if you propogate Apophis' orbit from today's vectors, to its 2029 approach of Earth, it only misses JPL's predicted close approach distance by about 200 km. That's not too bad considering both Apophis and Earth travelled over 20 billion kilometers each, on completely separate trajectories during this 22 year period to get to their encounter. For comparison, if someone could hit a golf ball from San Francisco to Tijuana, Mexico, this would be like predicting within 1/10th of a millimeter where it would land. This shows you just how powerful a Newtonian-only model, using a simple Euler algorithim with double precision variables and a slow time step can be.

I suspect that the 200 km difference has little to do with my use of only 15 digits. Their simulator takes GR into account, and perhaps more gravitational perturbers from the asteroid belt, as well as non-spherical gravity.

Since you say that UBasic lets you choose how many digits you use, if you get your code working, I'd be interested in seeing a plot of how different number of digits effects the results. How many are required for good results? How many are overkill?

publius · #15 (**permalink**) 16-April-2007, 08:06 PM

Tony,

Is the source code for your simulator available? If it's "trade secrets", then no problem.

The reason I ask is I'd like to run it through Intel's C++ compiler and compare it. That sucker does vectorizing and parallelizing for multi-CPUs (and multi-core and hyperthreading) and I always like to see if it's as good as they like to brag about. It has impressed me with other stuff.

-Richard

tony873004 · #16 (**permalink**) 16-April-2007, 09:13 PM

Quote:

Originally Posted by publius

Tony,

Is the source code for your simulator available? If it's "trade secrets", then no problem.

The reason I ask is I'd like to run it through Intel's C++ compiler and compare it. That sucker does vectorizing and parallelizing for multi-CPUs (and multi-core and hyperthreading) and I always like to see if it's as good as they like to brag about. It has impressed me with other stuff.

-Richard

Gravity Simulator is written in Visual Basic 6.0, with the exception of the number-crunching routine. It's written in C++ since C++ is much faster than VB 6.0, and this is the only part of the program where speed is advantagous.

The C++ code is compiled into a .dll. The VB code calls this routine once per iteration. I provided a link to the C++ code above. No trade secrets here. It's just a basic Euler's method propogating formulas from your High School Physics book. I guess it would be as easy as copying it, pasting it, and compiling it into a .dll using Intel's compiler. I used Visual C++.Net.

Assuming that works, and assuming you have Gravity Simulator installed on your computer, you simply need to replace the file "codeoptomizer.dll" with your replacement .dll. Just rename the old one incase you want to restore it. Then to test speeds, you want to run Gravity Simulator with the graphics turned off. This forces the program to only crunch the numbers.

I'd be very interested to see if a different compiler would make a big difference. So if you attempt this, let me know.

BTW, if you want to see a comparison between the speed of VB 6.0 and .NET compiled C++, simply rename codeoptomizer.dll, so VB can't find it. This forces it to use VB to do the math. I believe the difference in speed is a factor of 6.

publius · #17 (**permalink**) 16-April-2007, 10:05 PM

Tony,

I see. I doubt the Intel compiler would do much better on that simple function in the .dll that MS's now. What it would shine at if the actual time-stepping, the big loop was done as well. I don't know see enough if the above function to do much.

If the whole thing were there, Intel's compiler has vectorizing and parallelizing, plus profile guided optimization. The latter is something else. You do an "instrumented build", where the compiler puts metric code in that keeps track of the data flow and stores the results in a data file. You turn that instrumented build loose on typical input data conditions. It's slow as heck, but it is carefully keeping track. Now, the compiler then takes that data set to better optimize the final build.

That can do wonders. But to do that here, the Intel compiler would need to see the whole picture, with the main time step loop in there as well. And it would need to know the better data picture of where the inputs to the loop are stored, etc, etc.

-Richard