Search code examples
c++oopobjectnumerical-methods

How can I update c++ class members with a function?

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!

Solution

  • 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