# Galaxy simulation with the velocity Verlet algorithm

• Comp Sci
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 blackhole (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):
Algorithm:
//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;

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]));
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
{
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)
//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;
}
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;
}

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 megayears, 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
Mentor
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.

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.

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

DrClaude