Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Velocity Verlet C++ implementation

  1. Jul 29, 2017 #1
    I have been working on implementing a solar system system simulator in C++ - but am getting incorrect results on the order of 10^10 km, and it seems that the planets are moving directly away from the sun. I suspect that there is a mistake in the integrator (Velocity Verlet) that I have posted below. It would be appreciated if I could get someone to check that, assuming that the rest of the program is correct, the integrator is working. The code should be self-explanatory.

    Thanks!

    Code (C):
    void verlet(System& system, Trajectory& trajectory, double delta){

        long n = trajectory.get_number_of_rows();
        long n_bodies = system.get_number_of_bodies();

        double delta2 = pow(delta, 2);

        std::cout << "Starting integrator." << std::endl;

        for (int i = 0; i < n; ++i) {
            if(i == 0){
                std::vector<Vector3 > x0 = system.get_positions();
                std::vector<Vector3> v0 = system.get_velocities();

                trajectory.set_position(x0, v0);
                system.set_positions(x0);
                system.set_velocities(v0);
            }
            else {

                std::vector<Vector3> x0 = trajectory.get_positions_at_index(i-1);
                std::vector<Vector3> v0 = trajectory.get_velocities_at_index(i-1);

                std::vector<Vector3> a0 = system.get_accelerations();

                std::vector<Vector3> x1 {};

                for (long j = 0; j < n_bodies; ++j) {
                    x1.emplace_back(x0[j] + v0[j] * delta + delta2 * 0.5 * a0[j]);
                }

                system.set_positions(x1);

                std::vector<Vector3 > a1 = system.get_accelerations();

                std::vector<Vector3 > v1;

                for (long k = 0; k < n_bodies; ++k) {
                    v1.emplace_back(v0[k] + 0.5 *( a0[k] + a1[k]) * delta );
                }

                system.set_velocities(v1);

                trajectory.set_position(x1, v1);
            }
        }
        std::cout << "Integration finished." << std::endl;
    }
     
     
  2. jcsd
  3. Jul 29, 2017 #2

    Mark44

    Staff: Mentor

    Since your results are so far off, I suggest you use a debugger and put breakpoints in, and compare the computed results with results done on paper to see where the calculations are going wrong. I would put the first breakpoint on the last statement of the first if block in the outer for loop. If the calculations are OK when i == 0, I would put another breakpoint just after the calculation for a0 in the else block. Again, compare the results of your program with results obtained with paper and pencil. Possibly you'll find the error within just a couple of iterations.

    Since you are new to C++ (which you said in another thread), do yourself a favor and learn at least the rudimentary use of a debugger.
     
  4. Jul 31, 2017 #3
    Thank you, Mark44, for your reply. It was helpful in finding the bug.

    .. Mads
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?
Draft saved Draft deleted