Solar System Simulation Project (velocity verlet help)

Question

For my modelling and simulation class project, I want to simulate a solar system. I'm starting with just a star (sun) and a planet (earth), but I'm running into a few problems already. I've spent some time now just reviewing and learning about different formulas and way to simulate how the planet's orbits will be affected by the star and surrounding objects. I want to use velocity verlet and eventually look into the n-body problem. I'm having numerous issues with my velocity verlet function. Sometimes it acts as if it's making earth orbit normally and then it " warp drives" earth off into some random location. I've also noticed I never get a "negative" acceleration, so my x and y coord. are always increasing, so I don't see how the earth is suppose to wrap back around the sun. Any help is greatly appreciated. The AGK::Prints I have just so I can see how the different variables are changing.

double velocityVerlet(float positionCalc, double position2, 
                      float &velocity, double massCalc, double mass2)
//positionCalc is the position being updated, position 2 is position of 
// other object, same with mass
{
    float force = forceFunc(positionCalc, position2, massCalc, mass2);
    agk::PrintC("Force is: ");
    agk::Print(force);
    float acceleration = accelerationFunc(force,massCalc);
    agk::PrintC("Accel is: ");
    agk::Print(acceleration);`;

    double newAccel = 0;

    positionCalc = positionCalc + velocity*dt + 
                   (.5*acceleration)*pow(dt,2); //calculates new position
    agk::PrintC("New Position is: ");
    agk::Print(positionCalc);
    force = forceFunc(positionCalc,position2,massCalc,mass2);
    newAccel = accelerationFunc(force, massCalc);

    velocity = velocity + .5*(acceleration + newAccel)*dt; //new velocity
    agk::PrintC("Velocity is: ");
    agk::Print(velocity);

    return positionCalc;
}

Answer 1

The facts that your integrator accepts scalars and that your question is about 2-dimensional system makes me think that you are calling the integrator twice, once for each component. This simply won't work since your system will be taking unrealistic moves through the phase space. The integrator works with vector quantities:

X(t+dt) = X(t) + V(t) dt + (1/2) A(t) dt²

V(t+dt) = V(t) + (1/2)(A(t) + A(t+dt)) dt

Here X(t) is a column-vector that consists of the coordinates of all particles - this is the configuration subspace of the system's phase space. V(t) is a column-vector of the speeds of all particles, technically representing the momentum subspace. The same applies to A(t). Those have to updated simultaneously, not separately.

The whole velocity Verlet procedure translates as follows in code for force fields that do not depend on the speed (e.g. classical gravity):

Vector forces[num_particles];

// Compute initial forces
forces = computeForces(positions);

for (int ts = 0; ts < num_timesteps; ts++)
{
   // Update positions and half-update velocities
   for (int i = 0; i < num_particles; i++)
   {
      positions[i] += velocities[i]*dt + 0.5*(forces[i] / m[i]) * dt*dt;
      velocities[i] += 0.5*(forces[i] / m[i]) * dt;
   }

   // Compute new forces and half-update velocities
   forces = computeForces(positions);

   for (int i = 0; i < num_particles; i++)
   {
      velocities[i] += 0.5*(forces[i] / m[i]) * dt;
   }
}

Note that all positions are updated first before the next round of force evaluation. Also it is only necessary to evaluate the forces once per iteration since the positions do not change during the second update to the velocities. In the example code above Vector is a class that implements an n-dimensional vector and holds n components (e.g. 2 in your 2d-case). It also overloads the + and += operators to implement vector (component-wise) addition, as well as * and / to implement multiplication / division by a scalar. This is just to illustrate the case and can be replaced by inner loops over the components of each position/velocity vector.

The correct choice of time step is very important. Too small time step will slow down the simulation significantly. Too large time step will result in impossible physics, e.g. jumping planets.

Answer 2

There are some problems with the physics and some problems with the code.

First, the problem with the physics. Assuming that we're not modelling an alternative universe where the laws of physics are different, Newton's law of universal gravitation says that F=G*m1*m2/(r*r). However, force is a vector, not a scalar, so it has both magnitude and direction.

What the code calculates in forceFuncX is actually the magnitude rather than merely the component of the force parallel to the X axis. forceFuncY has the same flaw.

Next is the calculation for acceleration. Physics says that this is a=F/m. Mass is a scalar but both acceleration and force are vectors. So to calculate a_x, we could either use F_x/m or we could calculate F*cos(a)/m. Because cos(a) (with a being the angle from one CelesitalObject to another in 2D space) = dx/r, we can make this a_x = F*dx/(m*r) which is almost but not quite what you've got in your calculation (it's missing r in the divisor).

An alternative approach is to use std::complex but I'll not show that approach with the assumption that you may wish to extend this model to three dimensions.

That brings us to problems with the code. First, since you're using C++ and writing a simulation of a physical system of discrete objects, it makes sense that you've defined a CelestialObject class. What makes less sense is that your functions are called by picking out individual pieces of those objects and then calling C-style functions. We can improve the code by using those objects better. To start, since you haven't posted one, here's a CelestialObject class based on the interface I inferred from your code:

class CelestialObject 
{
public:
    CelestialObject(std::string name, float mass, float X, float Y, 
        float VX, float VY)
            : myname(name), m(mass), x(X), y(Y), vx(VX), vy(VY) {}
    void setPosition(float X, float Y) { x=X; y=Y; }
    void setVelocity(float VX, float VY) { vx=VX; vy=VY; }
    float getMass() const { return m; } 
    float getX() const { return x; } 
    float getY() const { return y; } 
    float getVX() const { return vx; } 
    float getVY() const { return vy; } 
    friend std::ostream& operator<<(std::ostream& out, 
                                    const CelestialObject& obj) {
        return out << obj.myname << '\t' << obj.x << '\t' << obj.y 
                   << '\t' << obj.vx << '\t' << obj.vy << std::endl;
    }
private:
    std::string myname;
    float m, x, y;
    float vx, vy;
};

Next, a few helper functions:

// returns square of distance between objects
float distance_sq(const CelestialObject &a, const CelestialObject &b)
{
    // distance squared is (dy^2 + dx^2)
    return pow(a.getY()-b.getY(),2) + pow(a.getX()-b.getX(),2);
}

// returns magnitude of the force between the objects
float force(const CelestialObject &a, const CelestialObject &b)
{
    //  F=(G * m1 * m1)/(r^2) in the direction a->b and b->a
    return G*a.getMass()*b.getMass()/distance_sq(a, b);
}

// returns the angle from a to b
float angle(const CelestialObject &a, const CelestialObject &b)
{
    return atan2f(b.getY()-a.getY(),b.getX()-a.getX());
}

Finally the actual verlet:

void updatePosition(CelestialObject &a, CelestialObject &b )
{ 
    float F = force(a,b);
    float theta = angle(a,b);
    float accela = F/a.getMass();
    float accelb = -F/b.getMass();

    // now that we have the acceleration of both objects, update positions
    // x = x +v *dt + a*dt*dt/2
    //   = x + dt * (v + a*dt/2)
    a.setPosition(
     a.getX() + dt * (a.getVX() + accela*cos(theta)*dt/2),
     a.getY() + dt * (a.getVY() + accela*sin(theta)*dt/2) 
    );
    b.setPosition(
     b.getX() + dt * (b.getVX() + accelb*cos(theta)*dt/2),
     b.getY() + dt * (b.getVY() + accelb*sin(theta)*dt/2)
    );
    // get new acceleration a'
    F = force(a,b);
    float thetap = angle(a,b);
    float accelap = F/a.getMass();
    float accelbp = -F/b.getMass();
    // and update velocities
    // v = v + (a + a')*dt/2
    a.setVelocity(
     a.getVX() + (accela*cos(theta) + accelap*cos(thetap))*dt/2,
     a.getVY() + (accela*sin(theta) + accelap*sin(thetap))*dt/2
    );
    b.setVelocity(
     b.getVX() + (accelb*cos(theta) + accelbp*cos(thetap))*dt/2,
     b.getVY() + (accelb*sin(theta) + accelbp*sin(thetap))*dt/2
    );
}

And finally some simple test code.

#include <string>
#include <iostream>
#include <vector>
#include <cmath>

const float G(6.67e-11);  // N*(m/kg)^2
const float dt(0.1);      // s
// all of the other code goes here...
int main()
{
    CelestialObject anvil("anvil", 70, 370, 0, 0, 0);
    CelestialObject earth("earth", 5.97e+24, -6.378e6, 0, 0, 0);
    std::cout << "Initial values:\n" << earth << anvil;
    std::cout << "Dropping an anvil from the top of a 370m building...\n"
              "It should hit the ground in about 8.7 seconds.\n";
    int t;
    for (t=0; anvil.getX() > 0; ++t) {
    std::cout << dt*t << '\t' << anvil; 
    updatePosition(anvil, earth);
    }
    std::cout << "Final values at t = " << dt*t << " seconds:\n" 
              << earth << anvil;
    return 0;
}

The test code uses a time step of 0.1s which is far too short for your solar system but fine for this quick test which is to see if we get a reasonable result for a known system. In this case, I've chosen a two-body system consisting of the planet Earth and an anvil. This code simulates dropping an anvil from the top of a 370m building which we can easily calculate will hit the ground in about 8.7s if we neglect air resistance. To keep the coordinates simple, I chose to place the origin (0,0) at the surface of the earth, and to consider the top of the building to be at (370,0). When the code is compiled and run, it produces the following:

Initial values:
earth   -6.378e+06  0   0   0
anvil   370 0   0   0
Dropping an anvil from the top of a 370m building...
It should hit the ground in about 8.7 seconds.
0   anvil   370 0   0   0
0.1 anvil   369.951 -4.27834e-09    -0.97877    -8.55668e-08
0.2 anvil   369.804 -1.71134e-08    -1.95754    -1.71134e-07
0.3 anvil   369.56  -3.85051e-08    -2.93631    -2.567e-07
   ...
8.3 anvil   32.8567 -2.9474e-05 -81.2408    -7.1023e-06
8.4 anvil   24.6837 -3.01885e-05    -82.2197    -7.18787e-06
8.5 anvil   16.4127 -3.09116e-05    -83.1985    -7.27345e-06
8.6 anvil   8.04394 -3.16432e-05    -84.1774    -7.35902e-06
Final values at t = 8.7 seconds:
earth   -6.378e+06  3.79705e-28 9.98483e-22 8.72901e-29
anvil   -0.422744   -3.23834e-05    -85.1563    -7.4446e-06

As you can see, this seems to work, but there are problems. The first problem is that since the objects should only move along the X axis, all Y components should be 0. They are not because this code is not very well-designed from a numerical analysis point of view. Doing additions and subtractions of floating point numbers when one number is large and another is small is one problem. Another is the use of the atan2f function which returns only a float and then using that result in cos() and sin(). Trigonometric functions are actually best avoided if possible.

Finally, this program currently only works with two objects. Adding a third would be painful with this kind of scheme, so a better design would be to process a std::vector<CelestialObject> by first calculating the net force on each object by considering the positions and masses of all of the others. I'll leave that to you, but this should at least give you a start in the right direction.