Sign in to follow this  
neoaikon

Solar System Simulation Issues

Recommended Posts

I'm creating a Simulation of our solar system using Newtonian Gravity, and Velocity Verlet Integration. I have the inner planets setup (Mercury, Venus, Earth, and Mars) and they orbit the sun as they should. However getting the moon to orbit the earth has been an issue, its orbit will slowly expand outward and become more elliptical in one direction. I'm going to refrain from posting any of my code here, as there is quite a bit, for those that are interested you can find a link to the source code at the bottom, its in C#. Allow me to expand upon what I am doing when I setup and integrate these bodies: - Setup Intial Positions, and Velocities Velocities are calculated using v = Sqrt(gc * M / r) For the moon I calculate a velocity to orbit the sun and add the velocity to orbit the earth to it. Where gc is the Gravitational Constant, M is the mass of the body we're orbting around, and r is the distance between their centers. To velocity verlet integrate I do the following: - X += (XVelocity * DeltaTime) - Y += (YVelocity * DeltaTime) - Get Direction Unit Vector to target - XAcceleration = ((gc * M) / (r * r)) * XDirection - YAcceleration = ((gc * M) / (r * r)) * YDirection - XVelocity += 0.5 * (XAcceleration + OldXAcceleration) * DeltaTime - YVelocity += 0.5 * (YAcceleration + OldYAcceleration) * DeltaTime - OldXAcceleration = XAcceleration - OldYAcceleration = YAcceleration This works fine for the planets, but the moon doesn't seem to work out.It seems that however I am setting the initial velocity doesn't seem to be correct. since the moon is the only body that I'm giving 2 initial velocities too. Is adding the velocity to orbit the earth and sun together for the moon the right way to do things? Can someone please point me on the right track to keep me from going crazy? Newton.txt Here is an example of the program for visual purposes: Newton.ex_ Change the extension to .exe to run it.

Share this post


Link to post
Share on other sites
Your problem is worse than setting the initial veocity correctly.

Your system is gaining energy because your integration isnt accurate enough. Likely, your planets are also gaining energy but much more slowly such that you simply havent noticed yet.

Use a smaller timestep, or use a better integration technique.

Share this post


Link to post
Share on other sites
You seem to be right, as increasing the time step causes the planetary orbits to slowly expand as well, although some seem to be alot more stable than others.I actually ahd to increase the timestep to about 5,000,000 to make this effect apparent. You can see in this screen shot(Screenshot), Mars is a lot further away than it should be, and the earth is orbiting at 187,000,000Km, 40,000,000Km farther than it should be. I understood the velocity verlet system was supposed to conserve energy, so it must be a mistake in my formulas.

I wonder if this could stem from a modification to the original formula, the original position integration I read from the examples was:

- X += (XVelocity * DeltaTime) + (0.5 * XAcceleration * (DeltaTime * DeltaTime))
- Y += (YVelocity * DeltaTime) + (0.5 * YAcceleration * (DeltaTime * DeltaTime))

But this caused my planets to just move toward the sun at a constant rate because the acceleration and thus the orbits were not circular...or at least that was the case when I was just using the earth, now when I'm experimenting uncommenting them in my code, and the orbits are still circular. I'll need to check into it to see whats going on since it wasn't behaving this way earlier.

EDIT. Adding back in the acceleration caused the system to be unstable compared to the other way around, so intuition tells me not to include it. Therefore I must be gaining that energy from somewhere else.

Someone said that starting the planets directly perpendicular to the sun (to the left side in my case) with velocities pointing straight down is not realistic, as this would give the body a slightly outward orbit that would over time lead to it being lost. Could this be where the extra energy is coming from?

[Edited by - neoaikon on October 12, 2008 3:36:43 AM]

Share this post


Link to post
Share on other sites
Your integration code looks funny to me. In particular, things like this:

//Merc Step
SolDir = NormalizeV2D(Sol, Merc);
SolAcc = new Vector2d((gc * (Merc_mass / (dist * dist))) * -SolDir.x, ((gc * Merc_mass) / (dist * dist)) * -SolDir.y);
SolVel.x += ((SolAcc.x + OldSolAcc.x) / 2) * dt;
SolVel.y += ((SolAcc.y + OldSolAcc.y) / 2) * dt;
OldSolAcc.x = SolAcc.x;
OldSolAcc.y = SolAcc.y;

SolDir = NormalizeV2D(Sol, Venu);
SolAcc = new Vector2d((gc * (Venu_mass / (dist * dist))) * -SolDir.x, ((gc * Venu_mass) / (dist * dist)) * -SolDir.y);
SolVel.x += ((SolAcc.x + OldSolAcc.x) / 2) * dt;
SolVel.y += ((SolAcc.y + OldSolAcc.y) / 2) * dt;
OldSolAcc.x = SolAcc.x;
OldSolAcc.y = SolAcc.y;



The correct way is to calculate the new acceleration from all sources and then update the velocity, like this:


SolAcc = new Vector2d(0,0);

//Merc Step
SolDir = NormalizeV2D(Sol, Merc);
SolAcc += new Vector2d((gc * (Merc_mass / (dist * dist))) * -SolDir.x, ((gc * Merc_mass) / (dist * dist)) * -SolDir.y);

SolDir = NormalizeV2D(Sol, Venu);
SolAcc += new Vector2d((gc * (Venu_mass / (dist * dist))) * -SolDir.x, ((gc * Venu_mass) / (dist * dist)) * -SolDir.y);

... and so on for the rest of the planets ...

SolVel.x += ((SolAcc.x + OldSolAcc.x) / 2) * dt;
SolVel.y += ((SolAcc.y + OldSolAcc.y) / 2) * dt;
OldSolAcc.x = SolAcc.x;
OldSolAcc.y = SolAcc.y;



Share this post


Link to post
Share on other sites
f?
I have problems with numerical accuracy with doubles.

http://www.youtube.com/watch?v=LoYFV_9En6I

Quote:
if (source.x == 0)
source.x = 1;
if (source.y == 0)
source.y = 1;
if (target.x == 0)
target.x = 1;
if (target.y == 0)
target.y = 1;


Are you sure about it?

Is it just me, or are you computing the acceleration for only one planet, and then after second step throwing the acceleration out? I forgot to hit reply and vorpy already pointed the problem.

Share this post


Link to post
Share on other sites
I want to thank you both! Especially for putting up with my code mostly uncommented! I always seem to get my questions answered when I can show people the gears of my programs.

I think you were correct Vorpy, I actually knew that was the correct way but the examples I was going by said different, and since I wasn't getting proper motion, I wasn't going to risk deviating from the examples.

if (source.x == 0)
source.x = 1;
if (source.y == 0)
source.y = 1;
if (target.x == 0)
target.x = 1;
if (target.y == 0)
target.y = 1;

As for this, I cannot explain it, I'm sure it was important at the time, but I've removed it now. But Raghar I think you hit the nail on the head about the floats. When I changed them to doubles, I got much better motion, additionaly this also solves an unrelated problem where too low of a time step would cause the moon to not move properly, just toward the earth and then fly off once the force reached high enough.

As for the system gaining energy, it still seems to be doing so, but nearly as much as before, much more stable, T=5,000,000 before the moon has gained enough energy to reach escape velocity at DT=200.

Isn't it true that ODE Integration Systems will tend to either gain or lose energy from the system? I remember reading this somewhere but I'm unsure of where, I believe gaining energy is called diverging from ideal and losing energy is converging to ideal.

I probably shouldn't expect to achieve perfect conservation (As much as I'd want it) right?, although I've heard good things about the energy conservation present in the Gear Predictor-Corrector method. I will be uploading a heavily commented version of my code at the same location above soon, once I get done doing that.

EDIT: Forgot to thank you as well Rockoon, your observation of my system gaining energy is something I wouldn't have discovered on my own.

[Edited by - neoaikon on October 12, 2008 2:30:36 PM]

Share this post


Link to post
Share on other sites
- X += (XVelocity * DeltaTime) + (0.5 * XAcceleration * (DeltaTime * DeltaTime))
- Y += (YVelocity * DeltaTime) + (0.5 * YAcceleration * (DeltaTime * DeltaTime))

Adding the bolded parts resulted in excellent conservation of energy, as you can see in this image
Newton2
the moon has remained orbiting the earth for once. This simulation is running at DT=50 and took 30min to run. The T value is inaccurate because the value actually maxed so it wasn't counting anymore. I've left trails on so I can see the path over time, the green spiral is the result of plotting a line with a magnitude of 10 Screen Units out in the direction the moon is in relationship to the earth, beautiful if I may say so :D. To calculate the accuracy of the simulation, I calculated how many orbits around our earth the moon should do in one year.

Orb = 365.26 / 27.3, Orb = 13.38. In the image the result is 14 orbits due to the fact the moons orbit isn't realistic. In reality it is more elliptical with Min/Max distance of 363,104/405696Km, Where in my simulation I'm starting it at 384,339, the semi-major axis. Adjusting the formula for elliptical, created accurate orbit plots. :D, Thank you all for your help! Once again these forums have come to my rescue.

Share this post


Link to post
Share on other sites

Create an account or sign in to comment

You need to be a member in order to leave a comment

Create an account

Sign up for a new account in our community. It's easy!

Register a new account

Sign in

Already have an account? Sign in here.

Sign In Now

Sign in to follow this