Celestial mechanics: Why is my planet orbiting faster and faster?

[S Holtom]

TL;DR Summary

Part comp sci part astronomy...not sure where to put this

I've tried to make a "naive" implementation of a planet orbiting a gravitational point source, in Processing (basically Java).
Gravity is a constant, and adds to the planet's velocity inversely proportional to the square of distance. I start the planet off with a tangential velocity.
I get an elliptical orbit, as desired, but the planet orbits faster and faster until it is ejected. What am I missing? Or is it really the case that very precise values are needed for a stable orbit?

(I'll post the source code shortly, it's on another machine)

 

[Hornbein]

I suppose it is roundoff error. You could try the continuity correction, which means averaging the acceleration in two places. Or maybe three places. That will reduce the peak value of gravity.

 

[kuruman]

It could be a matter of precision and cumulative errors. Do you check to see that the speed is the same when the planet returns to the starting point value after an integer number of revolutions? You could also verify that angular momentum is conserved when you move the planet to a new point and make the necessary corrections if it isn't.

 

[S Holtom]

Yes I was starting to suspect a similar thing. I have my planet orbiting at a distance of +1 unit, but then that means gravity needs to be something like 0.001, then dividing by 60 because 60fps...the numbers do end up pretty small.

 

[S Holtom]

Thanks for the very fast replies :)

 

[Baluncore]

I get an elliptical orbit, as desired, but the planet orbits faster and faster until it is ejected. What am I missing? Or is it really the case that very precise values are needed for a stable orbit?

You need to select a more accurate numerical integration technique.
Look for a symmetrical technique that, when run backwards, returns the body to the starting point.

Start here; https://www.ias.edu/ids/~piet/act/comp
See Time-Symmetric Codes … https://www.ias.edu/ids/~piet/act/comp/algorithms

 

[Vanadium 50]

.the numbers do end up pretty small.

Small numbers are OK.
Big numbers are OK.
Mixing big and small numbers is not OK.

 

[S Holtom]

Here's the (Processing) source code:

/* This world is just a simulation of an elliptical orbit */
// World coords are -1 to +1, both dimensions

boolean g_debug = true;

final float kGravity = 0.001f;
final float kTangent = 0.05f;
final float kDistance = 0.5f;
final float kPlanetSize = 0.03f;

float g_aspectRatio;

private PVector _planetPos;
private PVector _planetVel;

public void setup()
{
  size(800, 600);  
  background(0);
  
  g_aspectRatio = float(height) / float(width);
  
  _planetPos = new PVector(0.0f, -kDistance);
  _planetVel = new PVector(kTangent, 0.0f);
}public void draw()
{
  clear();
  
  int time = millis();
  float delta = float(time) * 0.001f;
  
  Update(delta);
  Render();
}

public void Update (float secsDelta)
  {
    if (secsDelta > 0.0f)
    {
      float sqDist;
      PVector diff = new PVector(-_planetPos.x, -_planetPos.y);
      sqDist = ((diff.x * diff.x) + (diff.y * diff.y));
      
      diff.normalize();
      diff = Vector_ScalarMult(diff, (kGravity / sqDist) * secsDelta);
      
      _planetVel.x += diff.x;
      _planetVel.y += diff.y;
      
      _planetPos.x += _planetVel.x * secsDelta;
      _planetPos.y += _planetVel.y * secsDelta;
    }
  }
  
  public void Render()
  {
    float halfX = float(width) / 2.0f;
    float halfY = (float(height) / 2.0f) * g_aspectRatio;
    float planetSize = kPlanetSize * float(width);
    
    stroke(200);
    fill(50, 50, 255);
    circle(_planetPos.x * halfX + halfX, _planetPos.y * halfY + halfY, planetSize);
    stroke(255, 0, 0);
    
    int middleX = width/2;
    int middleY = height/2;
    int lineWidth = width / 30;
    int lineLength = height / int(30.0f * g_aspectRatio);
    line(middleX, middleY - lineLength, middleX, middleY + lineLength);
    line(middleX - lineWidth, middleY, middleX + lineWidth, middleY);
  }

 

[S Holtom]

Oh and if anyone is compiling this, you'll need a function from my vector helpers class:

public PVector Vector_ScalarMult (PVector one, float mult)
{
  PVector vec = new PVector();
  vec.x = one.x * mult;
  vec.y = one.y * mult;
  vec.z = one.z * mult;
  return vec;
}

 

[S Holtom]

NB: In Processing, "width" and "height" are reserved terms for the width and height of the app window.

 

[bob012345]

Summary: Part comp sci part astronomy...not sure where to put thisGravity is a constant, and adds to the planet's velocity inversely proportional to the square of distance. I start the planet off with a tangential velocity.
I get an elliptical orbit, as desired, but the planet orbits faster and faster until it is ejected. What am I missing? Or is it really the case that very precise values are needed for a stable orbit?

Please express this relationship mathematically so we can see what you are doing. I can't make sense of the first sentence above.

 

[S Holtom]

Well I'm not absolutely sure how to do that...
Let's say that since the last time I updated the planet's velocity, then deltaTime seconds have elapsed.

Then my program says that:
new_velocity = current_velocity + (vector_to_gravitational_source * gravitational_constant * deltaTime / distance²)

And then of course the planet's position is just:
new_position = current_position + (new_velocity * deltaTime)

 

[S Holtom]

Oh and by the "tangential" comment, what I meant was:

1. The initial position of the planet is at some offset in the y axis, so coordinates are (0.0, 1.0)
2. The initial velocity of the planet is purely in the x axis, so e.g. (0.3, 0.0)

 

[bob012345]

Well I'm not absolutely sure how to do that...
Let's say that since the last time I updated the planet's velocity, then deltaTime seconds have elapsed.

Then my program says that:
new_velocity = current_velocity + (vector_to_gravitational_source * gravitational_constant * deltaTime / distance²)

And then of course the planet's position is just:
new_position = current_position + (new_velocity * deltaTime)

Are you trying to program a game of some kind? What is your goal? Did you have any physics courses?

 

[S Holtom]

It's involved in a kind of visualization of a galaxy that I am working on. I know what you're thinking: if I can't simulate one planet orbiting a gravitational point source, then simulating a whole galaxy is going to be quite a stretch.

But the visualization that I am working on is actually static. There are no orbits. The only reason that I am trying to model an orbit, is actually to do with how I am seeding my perlin flow fields. For which, in fact, a couple orbits of the "planet" is sufficient, so it doesn't actually matter that the orbits get faster.
But I would like to understand why my naive orbit maths fails.

 

[Vanadium 50]

Please express this relationship mathematically

Well I'm not absolutely sure how to do that...

That's kind ofa problem. If you don't know where you're going,, any road will take you there.

I won't say "The problem is..." but I will say "A problem is...". You are assuming that the steps are sufficiently small that the acceleration and position at the beginning of the step are the same as at the end. Or "close enough". This is known to accumulate errors.

You might want to look up Runge-Kutta.

 

[bob012345]

It's involved in a kind of visualization of a galaxy that I am working on. I know what you're thinking: if I can't simulate one planet orbiting a gravitational point source, then simulating a whole galaxy is going to be quite a stretch.

But the visualization that I am working on is actually static. There are no orbits. The only reason that I am trying to model an orbit, is actually to do with how I am seeding my perlin flow fields. For which, in fact, a couple orbits of the "planet" is sufficient, so it doesn't actually matter that the orbits get faster.
But I would like to understand why my naive orbit maths fails.

It fails because it lacks the correct physics behind why orbits are the way they are. The force of gravity is not constant in an elliptical orbit but total energy is conserved as is the angular momentum. You may not want a full primer on orbital mechanics to do your project but then what you could do is just implement the correct orbital equation for a hypothetical planet orbiting a hypothetical star. Or you could just model the Earth going around the sun as a circle and leave it at that.

I think specifically, you are adding energy to the orbit at every time interval so it will naturally spiral out of its orbit.

 

[S Holtom]

what you could do is just implement the correct orbital equation for a hypothetical planet orbiting a hypothetical star

Great. Can you recommend a good description?

 

[bob012345]

Great. Can you recommend a good description?

Make a model for the orbit of our good earth. You can get the data such as maximum and minimum distances, the period you know and you will have an equation of position vs. time. The Earth is an ellipse but nearly circular so you could approximate it as a circle for convenience then the equation is very simple and the velocity would be constant (in magnitude but the direction of course changes).

 

[vanhees71]

If I understand it right, you try to solve the equation of motion using a simple Euler integration. This is not very accurate, and it's likely that you accumulate a lot of numerical errors for simulations over a longer time. An easy check of the accuracy is to plot the energy and/or angular momentum as a function of time. Both should be conserved. If this is not the case, you can be sure that the numerics is not accurate anymore.

 

[bob012345]

If I understand it right, you try to solve the equation of motion using a simple Euler integration. This is not very accurate, and it's likely that you accumulate a lot of numerical errors for simulations over a longer time. An easy check of the accuracy is to plot the energy and/or angular momentum as a function of time. Both should be conserved. If this is not the case, you can be sure that the numerics is not accurate anymore.

It wasn't the correct equation of motion at all.

 

[S Holtom]

Well, this is embarrassing. It turns out that it was purely a programming error.
Basically, the function I was using for my timer does not return time elapsed since the last frame, but in fact time elapsed since the program started.
In fairness to me, I was thrown off by the fact that the planet got "ejected", when in fact the ejection was caused by an overflow, not the dynamics of the model itself.

With that bug fixed, the naive implementation works fine. I can have stable circular orbits, elliptical orbits and precessing orbits.

 

[S Holtom]

Thanks everyone for the replies, sorry if I wasted any time but I thought there were issues with the model.