Velocity Verlet C++ implementation

[madsmh]

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!

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;
}

 

[Mark44]

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.

 

[madsmh]

Thank you, Mark44, for your reply. It was helpful in finding the bug.

.. Mads