Galaxy simulation with the velocity Verlet algorithm

[Adec]

Homework Statement

Simulate a galaxy in C++ using the velocity Verlet algorithm as the integration method

Relevant Equations

Newton's laws of motion, velocity Verlet

To simulate the trayectories of solar systems around a black hole (i.e. a galaxy) I have 3 classes in C++: cSystem, cBlackHole and cGalaxy. cSystem assigns initial values of position, velocity, etc to a solar system. cBlackHole does the same but just for the black hole. And cGalaxy mixes cBlackHole with an array of cSystem, it's where the algorithm works (as a method/class function)

I've implemented the algorithm as I was told:
(0) Give initial random positions and velocities
(1) Evaluate the initial acceleration $$\mathbf{a}_i=-G \sum_{j \neq i} \frac{m_j (\mathbf{r}_i-\mathbf{r}_j)}{|\mathbf{r}_i-\mathbf{r}_j|^3}$$
(2) Evaluate the change in position, and use an auxilary variable "w":
$$\mathbf{r}_i(t+h)=\mathbf{r}_i(t)+h\mathbf{v}_i(t)+\frac{h^2}{2}\mathbf{a}_i(t)$$
$$\mathbf{w}_i=\mathbf{v}_i(t)+\frac{h}{2}\mathbf{a}_i(t)$$
(3) Evaluate the change in acceleration using the new positions (see formula in (1))
(4) Evaluate the change in velocity
$$\mathbf{v}_i(t+h)=\mathbf{w}_i+\frac{h}{2}\mathbf{a}_i(t+h)$$
(5) t=t+h, go back to (2).

However, even when I use only 1 solar system and give it its orbital velocity, it ends up flying away, or collapsing into the black hole. So, my inclination is that the algorithm is the thing that is failing (but it could be other thing).

Here's the algorithm in my program (comments are in Spanish, ignore them if you don't understand it):

//Método de aplicación del algoritmo de verlet a la galaxia completa
void cGalaxia::verlet(cSistema (&Solar)[numsist], cAgujero Negro,
 double (&ax)[numsist+1], double (&ay)[numsist+1], double h)
{
    //Declaración de variables
    double wx[numsist+1],wy[numsist+1],xaux,yaux,denom, xauxBH, yauxBH, denomBH;
    int i,j;

    //Posiciones, y velocidades auxiliares
    for(i=0; i<numsist; i++)
    {
        //Cambio de r(t) a r(t+h)
        Solar[i].setPosX(Solar[i].getPosX()+(h*Solar[i].getVelX())+((h*h/2.)*ax[i]));
        Solar[i].setPosY(Solar[i].getPosY()+(h*Solar[i].getVelY())+((h*h/2.)*ay[i]));
        //Velocidades auxiliares
        wx[i]=Solar[i].getVelX()+((h/2.)*ax[i]);
        wy[i]=Solar[i].getVelY()+((h/2.)*ay[i]);
    }
    //Para el agujero negro
    Negro.setPosXBH(Negro.getPosXBH()+(h*Negro.getVelXBH())+((h*h/2.)*ax[numsist]));
    Negro.setPosYBH(Negro.getPosYBH()+(h*Negro.getVelYBH())+((h*h/2.)*ay[numsist]));
    wx[numsist]=Negro.getVelXBH()+((h/2.)*ax[numsist]);
    wy[numsist]=Negro.getVelYBH()+((h/2.)*ay[numsist]);
    

    //Aceleraciones
    for(i=0; i<numsist;i++) //Evalúo cada sistema
    {
        ax[i]=0;ay[i]=0; //Inicialización aceleraciones
        for(j=0;j<numsist;j++) //Evalúo la interacción con los demás sistemas
        if (j!=i && (abs(Solar[i].getPosX()-Solar[j].getPosX())<5)
        && (abs(Solar[i].getPosY()-Solar[j].getPosY())<5)) //Cuadrado de 5x5 para no hacer
                                                           //calculos inútiles, solo influencia cercana
        {
            //Variables auxiliares
            xaux=Solar[i].getPosX()-Solar[j].getPosX();
            yaux=Solar[i].getPosY()-Solar[j].getPosY();
            denom=pow(((xaux*xaux)+(yaux*yaux)),1.5);

            //Aceleraciones
            ax[i]-=G*Solar[j].masa()*xaux/denom;
            ay[i]-=G*Solar[j].masa()*yaux/denom;
        }   
        //Añado la contribución del agujero negro a cada sistema i
        xaux=Solar[i].getPosX()-Negro.getPosXBH();
        yaux=Solar[i].getPosY()-Negro.getPosYBH();
        denom=pow(((xaux*xaux)+(yaux*yaux)),1.5);
        ax[i]-=G*Negro.masaBH()*xaux/denom;
        ay[i]-=G*Negro.masaBH()*yaux/denom;
        //Calculo la aceleración del agujero negro
        xauxBH=Negro.getPosXBH()-Solar[i].getPosX();
        yauxBH=Negro.getPosYBH()-Solar[i].getPosY();
        denomBH=pow(((xauxBH*xauxBH)+(yauxBH*yauxBH)),1.5);
        ax[numsist]-=G*Solar[i].masa()*xauxBH/denomBH;
        ay[numsist]-=G*Solar[i].masa()*yauxBH/denomBH;
    }

    //Velocidades
    for(i=0;i<numsist;i++)
    {
        Solar[i].setVelX(wx[i]+((h/2.)*ax[i]));
        Solar[i].setVelY(wy[i]+((h/2.)*ay[i]));
    }
    //Para el BH
    Negro.setVelXBH(wx[numsist]+((h/2.)*ax[numsist]));
    Negro.setVelYBH(wy[numsist]+((h/2.)*ay[numsist]));
    
    return;
}

I've also changed the units of some variables, for the sake of productivity. I've used positions between 0 and 100 kilolightyears (the typical width of a galaxy), time is in mega years, and mass is in solar masses; for that I've also used $$G=6.351695379\times 10^{-10} \ \frac{kly^3}{M_{\odot} My^2}$$

You can see the full program by clicking here.
Here's a gif showing the behaviour of the system (made with gnuplot): gif

 

[DrClaude]

However, even when I use only 1 solar system and give it its orbital velocity, it ends up flying away, or collapsing into the black hole. So, my inclination is that the algorithm is the thing that is failing (but it could be other thing).

Since the initial positions and velocities are random, this does not surprise me. The conditions have to be just right to have closed orbits.

 

[Adec]

Since the initial positions and velocities are random, this does not surprise me. The conditions have to be just right to have closed orbits.

I'm supposed to give them random velocities as well. But, just for the sake of checking, I'm giving them the orbital speed (not random at all), and I do not get a circular orbit. Therefore, there has to be a problem.
If you read the code, line 87, there's the formula
$$v=\sqrt{\frac{GM}{r}}$$
Which is the speed that I'm using for every system, given its (random) position r.
That is why I said

However, even when I use only 1 solar system and give it its orbital velocity, it ends up flying away, or collapsing into the black hole.

 

[Adec]

I've already solved the problem. In line 351 I called "aceleracionx" twice instead of using the y component, now I get the circular orbits perfectly! :D