Velocity Verlet applied to solar system with C++

In summary, this conversation involved a programmer seeking help with their implementation of the velocity Verlet algorithm in C++ for simulating the solar system. They were specifically trying to plot the trajectory of an object in relation to the Sun and were encountering unexpected results. The conversation included discussion of the program's code, including a custom CelestialObject class and functions for calculating distance, force, and angle between objects. It was also mentioned that they had successfully used the algorithm for simulating a falling object, but were having trouble with a comet. Possible causes for the issue were suggested, such as incorrect time step or incorrect velocity components.
  • #1
Rebecca1990
1
0

Homework Statement


Hello, I am attempting to use the velocity Verlet algorithm of integration applied to the solar system in c++. My goal is be able to use the outputted position components in a plot to see if the trajectory of my object is elliptical/parabolic/hyperbolic resulting from the gravitational interaction with the Sun. I have coded it using c++. I then load the output file into MATLAB and plot. I did an example that mimicked a falling object that interacted with Earth, with a dt = 0.1. It worked fine when plotting. I attempted using x and y velocity components I found for Halley's comet using dt = 3.154e+7s (one year). When I load the data, and plot the x and y components of the position, I thought I would get an elliptical shaped object.Instead, I am getting a straight line that is decreasing as x increases. Can someone help me?

Homework Equations

The Attempt at a Solution


C:
[/B]
#include <string>
#include <iostream>
#include <vector>
#include <cmath>
#include <fstream>
#include <cassert>

   const double G = 6.67e-11;  // N*(m/kg)^2
   //double dt = 0.1; // seconds for falling
   double dt = 3.154e+7; // comet
class CelestialObject
{
private:
   std::string m_name;
   double m_mass;
   double m_positionx;
   double m_positiony;
   double m_velocityx;
   double m_velocityy;
public:
   CelestialObject(std::string name, double mass, double* position, double* velocity) // constructor
   {
      m_name = name;
      m_mass = mass;
      m_positionx = position[0];
      m_positiony = position[1];
      m_velocityx = velocity[0];
      m_velocityy = velocity[1];
   }

   // setters and getters

   void setPosition(double* position)
   {
      m_positionx = position[0];
      m_positiony = position[1];
   }
   void setVelocity(double* velocity)
   {
      m_velocityx = velocity[0];
      m_velocityy = velocity[1];
   }
   double getMass() const
   {
      return m_mass;
   }
   double getPositionX() const
   {
      return m_positionx;
   }
   double getPositionY() const
   {
      return m_positiony;
   }
  double getVelocityX() const
   {
      return m_velocityx;
   }
   double getVelocityY() const
   {
      return m_velocityy;
   }
  

   friend std::eek:stream& operator<<(std::eek:stream& out, const CelestialObject& object)
   {
      return out << object.getPositionX() << '\t' << object.getPositionY() << '\t' << object.getVelocityX() << '\t' << object.getVelocityY() << std::endl;
   }
};

   // returns square of distance between objects
   double distanceSquared(const CelestialObject &a, const CelestialObject &b)
   {
      // distance squared is (dy^2 + dx^2)
      double rSquared = pow(a.getPositionY()-b.getPositionY(),2) + pow(a.getPositionX()-b.getPositionX(),2);
      return rSquared;
   }
   // returns distance between objects
   double distance(const CelestialObject &a, const CelestialObject &b)
   {
      // distance is (dy^2 + dx^2)^(1/2)
      double r = sqrt(pow(a.getPositionY()-b.getPositionY(),2) + pow(a.getPositionX()-b.getPositionX(),2));
      return r;
   }

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

   // returns the angle from a to b
   double angle(const CelestialObject &a, const CelestialObject &b)
   {
      double angle = atan2f(b.getPositionY()-a.getPositionY(),b.getPositionX()-a.getPositionX());
      return angle;
   }

   // Velocity Verlet algorithm

   void updatePosition(CelestialObject &a, CelestialObject &b )
   {
      double F = force(a,b);
      double theta = angle(a,b);
      double accelerationA = F/a.getMass();
      double accelerationB = -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)
   double newPositionA[2] = {a.getPositionX() + dt * (a.getVelocityX() + accelerationA*((b.getPositionX() - a.getPositionX())/distance(a,b))*dt/2), a.getPositionY() + dt * (a.getVelocityY() + accelerationA*((b.getPositionY() - a.getPositionY())/distance(a,b))*dt/2)};
   a.setPosition(newPositionA);
   double newPositionB[2] = {b.getPositionX() + dt * (b.getVelocityX() + accelerationB*((b.getPositionX() - a.getPositionX())/distance(a,b))*dt/2), b.getPositionY() + dt * (b.getVelocityY() + accelerationB*((b.getPositionY() - a.getPositionY())/distance(a,b))*dt/2)};
   b.setPosition(newPositionB);

   // get new acceleration a'
   F = force(a,b);
   double thetaNew = angle(a,b);
   double newAccelerationA = F/a.getMass();
   double newAccelerationB = -F/b.getMass();
  
   // update velocities
   // v = v + (a + a')*dt/2
   double newVelocityA[2] = {a.getVelocityX() + (accelerationA*((b.getPositionX() - a.getPositionX())/distance(a,b)) + newAccelerationA*((b.getPositionX() - a.getPositionX())/distance(a,b)))*dt/2, a.getVelocityY() + (accelerationA*((b.getPositionY() - a.getPositionY())/distance(a,b)) + newAccelerationA*((b.getPositionY() - a.getPositionY())/distance(a,b)))*dt/2};
   a.setVelocity(newVelocityA);
   double newVelocityB[2] = {b.getVelocityX() + (accelerationB*((b.getPositionX() - a.getPositionX())/distance(a,b)) + newAccelerationB*((b.getPositionX() - a.getPositionX())/distance(a,b)))*dt/2, b.getVelocityY() + (accelerationB*((b.getPositionY() - a.getPositionY())/distance(a,b)) + newAccelerationB*((b.getPositionY() - a.getPositionY())/distance(a,b)))*dt/2};
   b.setVelocity(newVelocityB);
   }int main()
{
   // example of halley's comet
   double positionHalley[2] = {48696201283, -6.8734238e+10};
   double positionSun[2] = {5.6378e7,5.6378e7}; // in metres, take origin at centre so equal component on both axes
   double velocityHalley[2] = {-9.096111, -6.916686};
   double velocitySun[2] = {0,0}; // may have to switch
   CelestialObject halley("halley", 2.2e+14, positionHalley, velocityHalley);
   CelestialObject sun("sun", 1.99e+30, positionSun, velocitySun);
   std::cout << "Initial values:\n" << sun << halley;
  
   std::eek:fstream write_output("halley2.txt");
   assert (write_output.is_open());

   int t;
   for (t=0; t < 80; t++)
   {
      write_output << dt*t << '\t' << halley;
      updatePosition(sun, halley);
   }
   std::cout << "Final values at t = " << dt*t << " seconds:\n";
   std::cout << "Sun: " << sun << "\n";
   std::cout << "Halley: "<< halley << "\n";
   write_output.close();
/*
   // example of mass falling off building file output
   double positionMass[2] = {0, 370};
   double positionEarth[2] = {0, -6.378e6};
   double velocityMass[2] = {0, 0};
   double velocityEarth[2] = {0,0};
   CelestialObject mass("mass", 70, positionMass, velocityMass);
   CelestialObject earth("earth", 5.97e+24, positionEarth, velocityEarth);

   std::eek:fstream write_output("falling.txt");
   assert (write_output.is_open());

   int t;
   for (t=0; mass.getPositionY() > 0; t++)
   {
      write_output << dt*t << '\t' << mass;
      updatePosition(mass, earth);
   }
   std::cout << "Final values at t = " << dt*t << " seconds:\n";
   std::cout << "Earth: " << Earth << "\n";
   std::cout << "Mass: "<< mass << "\n";
*/
   return 0;
}
[/B]
 
Last edited by a moderator:
Physics news on Phys.org
  • #2
Did you try smaller steps?
Did you check the size of all the relevant numbers? A straight line could mean you underestimate the strength of gravity.
 

1. What is Velocity Verlet and how is it applied to the solar system in C++?

Velocity Verlet is a numerical integration algorithm commonly used in physics simulations. It is used to calculate the position and velocity of objects in a system over time. In the context of a solar system simulation in C++, Velocity Verlet would be used to update the positions and velocities of planets and other celestial bodies based on their gravitational interactions with each other.

2. Why is Velocity Verlet a popular choice for simulating the solar system?

Velocity Verlet is a popular choice for simulating the solar system because it is a simple and efficient algorithm that accurately calculates the trajectories of objects in a system. It also conserves energy and momentum, making it a more accurate choice for long-term simulations.

3. How does Velocity Verlet handle the gravitational interactions between objects in a solar system?

Velocity Verlet uses Newton's law of gravitation to calculate the force between objects in a solar system. The algorithm then uses this force to update the positions and velocities of the objects in the system.

4. Are there any limitations or drawbacks to using Velocity Verlet for solar system simulations?

One limitation of Velocity Verlet is that it assumes a constant force between objects, which may not be accurate for highly eccentric orbits. Additionally, the algorithm can become unstable if the time step used is too large, leading to inaccurate results.

5. Are there any alternative methods for simulating the solar system besides Velocity Verlet?

Yes, there are several other numerical integration methods that can be used for simulating the solar system, such as the Runge-Kutta method and the Leapfrog method. Each method has its own advantages and limitations, and the choice of which one to use depends on the specific requirements of the simulation.

Similar threads

  • Engineering and Comp Sci Homework Help
Replies
2
Views
1K
  • Engineering and Comp Sci Homework Help
Replies
2
Views
2K
  • Engineering and Comp Sci Homework Help
Replies
2
Views
3K
  • Engineering and Comp Sci Homework Help
Replies
1
Views
2K
  • Engineering and Comp Sci Homework Help
Replies
1
Views
935
  • Engineering and Comp Sci Homework Help
Replies
15
Views
2K
  • Engineering and Comp Sci Homework Help
Replies
5
Views
3K
  • Engineering and Comp Sci Homework Help
Replies
3
Views
1K
  • Programming and Computer Science
Replies
8
Views
1K
  • Engineering and Comp Sci Homework Help
Replies
5
Views
2K
Back
Top