1

To describe the problem, I am trying to use objects in my code to stream line a three body problem. I have the following code for the object:

#include <stdlib.h>
#include <cstdio>
#include <iostream>
#include <cmath>
#include <vector>
#include "star.h"

using namespace std;

Star::Star( double m, double x_p, double y_p, double x_v, double y_v )
{
    init( m, x_p, y_p, x_v, y_v);
}

void Star::init( double m, double x_p, double y_p, double x_v, double y_v )
{
    Mass = m;
    X_Position = x_p;
    Y_Position = y_p;
    X_Velocity = x_v;
    Y_Velocity = y_v;
    R_Position[0] = X_Position;
    R_Position[1] = Y_Position;
    R_Velocity[0] = X_Velocity;
    R_Velocity[1] = Y_Velocity;
}

double Star::potential( Star star2, double dx, double dy )
{
    double G = 3.0548e34;
    double Potential;

    double x_component = X_Position - star2.X_Position + dx;
    double y_component = Y_Position - star2.Y_Position + dy;

    double R = sqrt(x_component*x_component + y_component*y_component);

    Potential = G* Mass* star2.Mass / R;
    return Potential;
}

double * Star::compute_forces( Star star2 )
{
    double h_x = ( X_Position - star2.X_Position )/1000;
    double h_y = ( Y_Position - star2.Y_Position )/1000;

    double *F = new double[2];

    F[0] = ( potential( star2, h_x, 0.0 ) - potential( star2, -h_x, 0.0 ) )/2*h_x;
    F[1] = ( potential( star2, 0.0, h_y ) - potential( star2, 0.0, -h_y ) )/2*h_y;

    return F;
}

void Star::verlet( Star star2, double h )
{
    double *Force = compute_forces( star2 );

    X_Position += h*X_Velocity + 0.5*h*h*Force[ 0 ];
    Y_Position += h*Y_Velocity + 0.5*h*h*Force[ 1 ];

    double *Force_new = compute_forces( star2 );

    X_Velocity += 0.5*h*(Force[ 0 ] + Force_new[ 0 ] );
    Y_Velocity += 0.5*h*(Force[ 1 ] + Force_new[ 1 ] );
}

Now I believe that the velocity verlet algorithm is correct, but when I run the code using this main file:

#include <iostream>
#include <fstream>
#include <cmath>
#include <cstdio>
#include "star.h"

using namespace std;

int main()
{
    Star star1( 50, 0.0, 0.0, 0.0, 0.0 );
    Star star2( 1.00, 0.0, 1.0, -1.0, 1.0 );
    Star star3( 1.00, 0.0, -1.0, 1.0, 1.0 );

    Star arr[3] = { star1, star2, star3 };

    double h = 10/1000;

    //for ( double time = 0.0; time <= 10.0; )
    //{
        for ( int inst = 0 ; inst< 3; ++inst )
        {
            for ( int jnst = 0; jnst < 3; ++jnst )
            {
                if ( inst != jnst )
                {
                    arr[ inst ].verlet( arr[ jnst ], h );
                    double *pos = arr[ inst ].get_positions();
                    cout << " " << pos[ 0 ] << " " << pos[ 1 ] << endl;

                }
            }

                    }
        //time += h;
    //}

    return 0;
}

The values of members of the Star object are not updating :/. Is there something I am missing? the out put of the cout is this:

 0 0
 0 0
 0 1
 0 1
 0 -1
 0 -1

Thank you in advance!

Edit:

I tried implementing a std::vector<double> for my forces, but I ended up with a segmentation fault.

Edit 2: After checking my get_positions() method I noticed it was returning only the initialized values. So I tried implementing this:

std::vector<double> get_positions(){  std::vector<double> temp = { X_Position , Y_Position }; return temp; }

And it worked so i implemented the following into my main code.

        std::vector<double> p1 = star1.get_positions();
        std::vector<double> p2 = star2.get_positions();
        std::vector<double> p3 = star3.get_positions();

        cout << p1[ 0 ] << " " << p1[ 1 ] << " "  << p2[ 0 ] << " " << p2[ 1 ] << " " << p3[ 0 ] << " " << p3[ 1 ] << endl;

However now I am stuck on a completely new problem... Now I am getting the following numbers for the algorithm updates!

5.66002e-320 2.31834e-316
1.132e-316 4.63669e-313
1.698e-319 6.95503e-316
1.132e-316 4.63669e-313
5.66002e-320 2.31834e-316
1.132e-316 4.63669e-313
1.698e-319 6.95503e-316
1.132e-316 4.63669e-313
5.66002e-320 2.31834e-316
1.132e-316 4.63669e-313
1.698e-319 6.95503e-316
1.132e-316 4.63669e-313

Which means some where I am multiplying by zeros somewhere in my code. The problem is I cant for the life of me see where. Thanks if there is any help!

Adam Kalinowski
  • 61
  • 1
  • 6
  • 2
    Please don't take this as flippant advice: you need to learn how to use your debugger. You can use it to step through the code line by line and inspect your variables in real time. Using your debugger you'll be able to see where your calculations deviate from expected values immediately. I'm not familiar with the problem you're trying to solve, so it's hard for me to follow what's going on. – JohnFilleau Oct 04 '20 at 00:05
  • What do expect the output to be? – JohnFilleau Oct 04 '20 at 00:05
  • 1
    This line `double h = 10/1000;` should always set `h` to `0`, since the expression on the right is integer math. You can use `double h = 10.0/1000;` for a floating point math expression, or `double h = 10e-3;` for something a little more expressive. – JohnFilleau Oct 04 '20 at 00:08
  • 2
    Unrelated to the problem you're asking about, but EXTREMELY related to future `bad_alloc` errors you're going to encounter in the future, the function `compute_forces` generates a memory leak by allocated new memory on every call. `new[]` is a request to the OS for a block of memory. It MUST be accompanied (eventually) by a symmetric `delete[]` operation to return that memory to the OS. In modern C++ we try to avoid explicit dynamic memory allocation like this. We use one of the many C++ containers that handles memory allocation for us. – JohnFilleau Oct 04 '20 at 00:14
  • 1
    It's better if `compute_forces` returns a `std::vector forces`. `std::vector` handles dynamic memory allocation (and deallocation) for us. Much more compartmentalized. – JohnFilleau Oct 04 '20 at 00:14
  • 1
    It's obvious you're just using C++ as a tool to get some research done, which is totally fine, but you can decrease your coding time in the future if after you get this bare-minimum working (and functionally correct) you take the entire project over to [codereview.se] and ask for advice there on best practices. – JohnFilleau Oct 04 '20 at 00:16
  • I dont take it as flippant advice at all :) Here are some clarifications. I am using nvim as an editor on Arch linux. I expect that there is some change after I use the Velocity Verlet method; however, if i do something like `X_Position += 1.0` I dont get a change. So I am wondering If I am just not using the Class functions correctly. I will try to use `std::vector` and see how that helps and I will make an update. Thanks for the quick reply! – Adam Kalinowski Oct 04 '20 at 01:29
  • If `X_Position += 1.0` isn't doing what you expect, there may be an issue in the `get_positions` function. What's that look like? Please provide a [example]. – JohnFilleau Oct 04 '20 at 01:34

2 Answers2

2

Error

If you want to divide by 2*h_x, you need to write this as /(2*h_x), else you divide by 2 and multiply by h_x, giving miniscule values for forces and thus not moving the system by much.

To complement this, you defined the step size in the main program as

double h = 10/1000;

The value on the right is identified as result of an integer division, which is 0. With this step size nothing will change.

Style

Do not construct two data fields for the same value, you would have to ensure that these fields are always synchronized. Use getter methods to present data in a different format.

For science it would be better to use an established vector class that then also provides vector arithmetic, like the one of boost/Eigen.

Use initialization list syntax in the constructor, you do not need an init function to just assign the values.

Verlet

The Verlet method does not work this way. Even if everything goes right coding-wise, the result is a first order method that neither preserves energy nor momentum.

The short version is, the stages of the Verlet method are the outer frame. In each stage, all computations have to be carried out for all objects before changing to the next stage. That is, all velocities change, then all positions chance, then all forces are computed and accumulated, then all velocities change with the new forces/accelerations for all objects.

Mixing these steps destroys the order of the method and all conservation properties. (The first two stages can be interleaved, as there is no interaction between objects.)

I implemented some of the suggested changes, using the data of the Pleiades IVP test suite example, as the provided data lead to a rapid explosion of the system.

The main program solarsystem.c with the main Verlet loop

#include <iostream>
#include <cstdio>
#include <vector>
#include "star.h"

using namespace std;

int main()
{
    vector<Star> arr = {
        Star( 1, 3.0, 3.0, 0.0, 0.0 ),
        Star( 2, 3.0,-3.0, 0.0, 0.0 ),
        Star( 3,-1.0, 2.0, 0.0, 0.0 ),
        Star( 4,-3.0, 0.0, 0.0,-1.25 ),
        Star( 5, 2.0, 0.0, 0.0, 1.0 ),
        Star( 6,-2.0,-4.0, 1.75, 0.0 ),
        Star( 7, 2.0, 4.0,-1.5, 0.0 )
    };
    int N = arr.size();
    double dt = 0.001;
    int count = 10;
    for ( double time = 0.0; time <= 3.0; time += dt)
    {
        for ( int inst = 0 ; inst< N; ++inst ) {
            arr[inst].Verlet_stage1(dt);
        }
        for ( int inst = 0 ; inst< N; ++inst ) {
            for ( int jnst = inst+1; jnst < N; ++jnst ) {
                arr[inst].acceleration(arr[jnst]);
            }
        }
        for ( int inst = 0 ; inst< N; ++inst ) {
            arr[inst].Verlet_stage2(dt);
        }
        if( 10 == count) {
            count = 0;
            for ( int inst = 0 ; inst< N; ++inst ) {
                cout << " " << arr[inst].Position[1] << " " << arr[inst].Position[0];
            }
            cout << "\n";
        }
        count++;
    }

    return 0;
}

and the implementation of the Star class with header

#pragma once
#include <eigen3/Eigen/Dense>

typedef Eigen::Vector2d Vec2D;

const double G = 1;

class Star {
    public:
        Star( double m, double x_p, double y_p, double x_v, double y_v )
             :Mass(m),Position(x_p,y_p),Velocity(x_v,y_v) {};
        double Mass;
        Vec2D Position, Velocity, Acceleration;
        
        void Verlet_stage1(double dt);
        void Verlet_stage2(double dt);
        
        double potential(Star other);
        void acceleration(Star &other);
};

and corpus

#include "star.h"

double Star::potential( Star other )
{

    Vec2D diff = Position-other.Position;
    double R = diff.norm();
    return G * Mass * other.Mass / R;
}

void Star::acceleration( Star &other )
{
    Vec2D diff = Position-other.Position;
    double R = diff.norm();
    Vec2D acc = (-G / (R*R*R)) * diff;
    Acceleration += other.Mass * acc;
    other.Acceleration -= Mass * acc;
}

void Star::Verlet_stage1( double dt )
{
    Velocity += (0.5*dt) * Acceleration;
    Position += dt*Velocity;
    Acceleration *= 0;
}
void Star::Verlet_stage2( double dt )
{
    Velocity += (0.5*dt) * Acceleration;
}

This results in the trajectories below. The picture is very depending on the step size dt as near singularities of the potential function, that is, if bodies come very close together, the promise of symplectic methods of near conservation of energy and momentums breaks apart.

enter image description here

Lutz Lehmann
  • 25,219
  • 2
  • 22
  • 51
  • Thank you so much! But in order to make sure to understand you right, let me clarify some things. In order to fix my code should i get rid of the nested particle loops and calculate the each star's forces all at once? if so what do you suggest would be the best way. Thanks again. – Adam Kalinowski Oct 04 '20 at 13:23
0

I'm personally not against using raw pointers, but complications occur when they don't get managed properly. I have no idea what this code does, more so how it does! Nevertheless, I have tried ameliorating a few errors which I could observe, but evidently this code requires a serious overhaul. I suppose the shortcomings in this code are just due to inexperience, which can be understood.

https://gcc.godbolt.org/z/5zT5o9 Please do keep in mind that this code is still leaking due to usage(non-manage) of raw pointers in various function bodies.