My planets went crazy - a Java related issue

Question

This morning I felt like writing a useless program and I ended up with this extremely minimal astronomic simulator in processing. MCVE version of code attached at the end of the post.

Then I added some planets, sit back and got ready to enjoy watching some planets orbiting around each other. However there turned out to be a problem.

The problem is, two originally static planets that getting too close to each other tend to fling each other to super ultra high speed and they will both disappear out of the screen forever. This is obviously against principle of momentum. Planets don't do this. I start to doubt there is something wrong with my Newton's Law implementation. But some other times the planets work just fine, provided that they keep a 'safe' distance.

If you wanna look into this problem, I have carefully included two setups for you. The first one, adds two planets and is beautifully stable. (You can do this by commenting out every line that has 'Three' in it.) You can test the second one by just input the default code. This one is where planets go crazy.

This issue seems like a mystery that origins from some computer science domain I have never explored(I am a noob though). For the sake of my hair, I would greatly appreciate it if someone could explain.

Start of code:

PVector planetOneLocation = new PVector(300, 200);
PVector planetOneSpeed = new PVector(0, -.1);
float planetOneMass = 1;

PVector planetTwoLocation = new PVector(100, 200);
PVector planetTwoSpeed = new PVector(0, .1);
float planetTwoMass = 1;

PVector planetThreeLocation = new PVector(200, 200);
PVector planetThreeSpeed = new PVector(0, 0);
float planetThreeMass = 10;

float g = 5;

void setup() {
  size(500, 500);
}

void draw() {

  updatePlanetOne();
  updatePlanetTwo();
  updatePlanetThree();

  planetOneLocation.add(planetOneSpeed);
  planetTwoLocation.add(planetTwoSpeed);
  planetThreeLocation.add(planetThreeSpeed);

  background(0);

  ellipse(planetOneLocation.x, planetOneLocation.y, 10*planetOneMass, 10*planetOneMass);
  ellipse(planetTwoLocation.x, planetTwoLocation.y, 10*planetTwoMass, 10*planetTwoMass);
  ellipse(planetThreeLocation.x, planetThreeLocation.y, 10*planetThreeMass, 10*planetThreeMass);
}

void updatePlanetOne() {
  PVector accDir = PVector.sub(planetTwoLocation, planetOneLocation);
  float a = g * planetTwoMass / (accDir.mag() * accDir.mag());
  accDir.normalize();
  PVector acceleration = accDir.mult(a);
  planetOneSpeed.add(acceleration);

  accDir = PVector.sub(planetThreeLocation, planetOneLocation);
  a = g * planetThreeMass / (accDir.mag() * accDir.mag());
  accDir.normalize();
  acceleration = accDir.mult(a);
  planetOneSpeed.add(acceleration);
}

void updatePlanetTwo() {
  PVector accDir = PVector.sub(planetOneLocation, planetTwoLocation);
  float a = g * planetOneMass / (accDir.mag() * accDir.mag());
  accDir.normalize();
  PVector acceleration = accDir.mult(a);
  planetTwoSpeed.add(acceleration);

  accDir = PVector.sub(planetThreeLocation, planetTwoLocation);
  a = g * planetThreeMass / (accDir.mag() * accDir.mag());
  accDir.normalize();
  acceleration = accDir.mult(a);
  planetTwoSpeed.add(acceleration);
}

void updatePlanetThree() {
  PVector accDir = PVector.sub(planetOneLocation, planetThreeLocation);
  float a = g * planetOneMass / (accDir.mag() * accDir.mag());
  accDir.normalize();
  PVector acceleration = accDir.mult(a);
  planetThreeSpeed.add(acceleration);

  accDir = PVector.sub(planetTwoLocation, planetThreeLocation);
  a = g * planetTwoMass / (accDir.mag() * accDir.mag());
  accDir.normalize();
  acceleration = accDir.mult(a);
  planetThreeSpeed.add(acceleration);
}

Update: After some effort I change float variables to double. But my planets are stilling bouncing out of the screen. I think besides the double/float problem, actually there is some issues with the resolution. I didn't define any time step, and that also contribute to inaccurate behaviour especially when speed is high.

Update 2: Setting up time-step helps a lot. Some setups that went crazy without a time-step now works fine. But as long as there is the possibility that the center of two planets will be extremely close, there will be chances that the system go wild again. To solve this, a better integrator is needed.

Updata 3: In respond to @kevin-workman I ported his beautiful code here to replace the original project code. The same third planet in original post is added and the corresponding Newton math is updated. Contradictory to his tests, it looks even if mv.add(p.speed.mult(p.mass)); is commented out, the third planet still go crazy(Now since I have changed the demonstration code to a minimal version, there is no such line anymore). I think the error introduced by mult()do contribute, but concretely the unstable integrator also plays a major role.

Answer 1

This problem has nothing to do with float vs double accuracy. Float has more than enough accuracy for this, and in fact Processing uses float for everything by default, so any double values you try to use will just get lost anyway.

All of your problems are caused by this line:

for (Planet p: planets) {
    mv.add(p.speed.mult(p.mass));
}

Particularly this bit:

p.speed.mult(p.mass)

This line multiplies each planet's speed by its mass.

You might be thinking that p.speed.mult(p.mass) will leave the original p.speed PVector unchanged, but this isn't the case. PVector is not immutable, so calling the mult() function changes the underlying PVector instance.

Your first two planets have a mass of 1, so this line doesn't really affect them. But your third planet has a mass of 10, which means you're multiplying its speed by 10 every single frame.

You can test this by simply changing the mass of one of your original planets to 10 or even 2, or by changing the mass of the third planet to 1.

So, to fix your main problem, just get rid of this line, or change it so that it doesn't modify the p.speed PVector.

More info can be found in the Processing reference for PVector#mult():

Multiplies a vector by a scalar. The version of the method that uses a float acts directly on the vector upon which it is called (as in the first example above). The versions that receive both a PVector and a float as arguments are static methods, and each returns a new PVector that is the result of the multiplication operation.

Basically, you could change that line to this:

PVector.mult(p.speed, p.mass)

That will leave p.speed unchanged and return a copy with the multiplied value.

Taking a step back, you're going to have another issue because your planets can get arbitrarily close to each other. In other words, their distance can approach (or equal) zero. This doesn't happen in real life, and if it does, you can bet things are going to go crazy. So you're going to have crazy "gravity assists" if things get too close. You might consider limiting their distances.

If you have further questions, it'll be easier to help you if you work from this MCVE instead of posting your whole project:

PVector planetOneLocation = new PVector(300, 200);
PVector planetOneSpeed = new PVector(0, -.1);
float planetOneMass = 1;

PVector planetTwoLocation = new PVector(100, 200);
PVector planetTwoSpeed = new PVector(0, .1);
float planetTwoMass = 10;

float g = 5;

void setup() {
  size(500, 500);
}

void draw() {

  updatePlanetOne();
  updatePlanetTwo();

  planetOneLocation.add(planetOneSpeed);
  planetTwoLocation.add(planetTwoSpeed);

  background(0);

  ellipse(planetOneLocation.x, planetOneLocation.y, 10*planetOneMass, 10*planetOneMass);
  ellipse(planetTwoLocation.x, planetTwoLocation.y, 10*planetTwoMass, 10*planetTwoMass);
}

void updatePlanetOne() {
  PVector accDir = PVector.sub(planetTwoLocation, planetOneLocation);
  float a = g * planetOneMass / (accDir.mag() * accDir.mag());
  accDir.normalize();
  PVector acceleration = accDir.mult(a);
  planetOneSpeed.add(acceleration);
}

void updatePlanetTwo() {
  PVector accDir = PVector.sub(planetOneLocation, planetTwoLocation);
  float a = g * planetTwoMass / (accDir.mag() * accDir.mag());
  accDir.normalize();
  PVector acceleration = accDir.mult(a);
  planetTwoSpeed.add(acceleration);
}

Answer 2

The problem is in the integrating part. The one used is an Euler Integrator and it's not very accurate. Verlet integration is normally used for physics simulation.

Quoting Ilmari Karonen's Answer, verlet can be implemented as:

acceleration = force(time, position) / mass;
time += timestep;
position += timestep * (velocity + timestep * acceleration / 2);
newAcceleration = force(time, position) / mass;
velocity += timestep * (acceleration + newAcceleration) / 2;