# Physics simulation

1. May 20, 2010

### Kynnath

I'm trying to write a program that will simulate basic gravity on a set of point particles. The problem I'm having is the program accumulates a lot of error very quickly, the particles start gaining speed and soon they're off to weird places. I can also track how the center of mass of the system starts moving all over the palce. Which shouldn't happen.

Not sure why error grows so quickly, though. I'm using the same equations I've seen in other simulators (actually a bit more refined, which should have led to less error I thought). But still they fly off in every direction.

The equations I'm using are:

Fg = G * ( mi * mj ) / ( d(i,j)^2 )

for the force of gravity, from which I get acceleration by dividing Fg by mi for the particle i and by mj for particle j. I run this for every pair of particles.

I add up the acceleration from each interaction, and use that to modify the velocity of the particle. I also store the acceleration from the previous iteration.

v(t) = v(t-1) + ( a(t) + a(t-1) ) / 2

and from this I get the new position for the particle

p(t) = p(t-1) + ( v(t) + v(t-1) ) / 2

The code itself is

Code (Text):
void nSistema::updateSistema1 ( )
{
long double G = 0.0000000001 ;
nLista<nVector<long double> > ainicial ; //Lista para almacenar las aceleraciones de la vuelta previa para cada objeto

centrodemasa_.setAceleracion( nVector<long double> () ) ;

for ( unsigned int i = 0 ; i < cuerpos_ ; ++i )
{
nCuerpo * cuerpoi ( sistema_.getObjeto( i+1 ) ) ;
if ( i == 0 ) // en la primer pasada tiene que guardarse la aceleracion previa y borrarse para calcular la nueva
{
ainicial.addObjetoEnd( cuerpoi->getAceleracion() ) ;  //aceleracion previa
cuerpoi->setAceleracion( nVector<long double> () ) ;  //setear a a 0 para el nuevo ciclo
}

for ( unsigned int j = i+1 ; j < cuerpos_ ; ++j )
{
nCuerpo * cuerpoj ( sistema_.getObjeto( j+1 ) ) ;
if ( i == 0 )// en la primer pasada tiene que guardarse la aceleracion previa y borrarse para calcular la nueva
{
ainicial.addObjetoEnd( cuerpoj->getAceleracion() ) ;
cuerpoj->setAceleracion( nVector<long double> () ) ;
}
nVector<long double> ij ( cuerpoj->getPosicion() - cuerpoi->getPosicion() ) ;
long double d2 ( ij.getNorma2() ) ;
long double d ( sqrt( d2 ) ) ; std::cout<< "d ( " << i+1 << " a " << j+1 << " ) = " << d << std::endl ;

nVector<long double> Fg ;
if ( d > 0.1 ) Fg = ij * G * ( cuerpoj->getMasa() * cuerpoi->getMasa() ) / ( d2 * d ) ;
else Fg = nVector<long double> ( ) ;

cuerpoi->setAceleracion( cuerpoi->getAceleracion() + ( Fg / cuerpoi->getMasa() ) ) ;
cuerpoj->setAceleracion( cuerpoj->getAceleracion() - ( Fg / cuerpoj->getMasa() ) ) ;

}
nVector<long double> a0 ( * ainicial.getObjeto( i+1 ) ) ;
//updatear v de cuerpoi ;
//nVector<long double> test ( cuerpoi->getAceleracion() ) ;
//std::cout << "test[" << i+1 << "] : a = ( " << test.getX() << " , " << test.getY() << " , " << test.getZ() << " )" << std::endl ;
nVector<long double> v0 ( cuerpoi->getVelocidad() ) ;
cuerpoi->setVelocidad( v0 + ( cuerpoi->getAceleracion() + a0 ) / 2 ) ;

//updatear p de cuerpoi ;
cuerpoi->setPosicion( cuerpoi->getPosicion() + ( cuerpoi->getVelocidad() + v0 ) / 2 ) ;

centrodemasa_.setPosicion( cuerpoi->getPosicion() ) ;
}
centrodemasa_.setPosicion( centrodemasa_.getPosicion() / cuerpos_ ) ;

}
I use a few classes I made myself, since the basis of the project is to get a bit of practice. In case they help, I have the definitions right here. I'd appreciate anyone who could offer some insight into why I'm having this problem, or how to fix it.

Code (Text):
#ifndef _nCuerpo_h_
#define _nCuerpo_h_

#include "nVector.h"

class nCuerpo
{
public :
nCuerpo ( ) ;
nCuerpo ( const long double & , const nVector<long double> & , const nVector<long double> & , const nVector<long double> & ) ;
nCuerpo ( const nCuerpo & ) ;
~nCuerpo ( ) ;
const long double getMasa ( ) const ;
const nVector<long double> getVelocidad ( ) const ;
const nVector<long double> getPosicion ( ) const ;
const nVector<long double> getAceleracion ( ) const ;
void setMasa ( const long double & ) ;
void setVelocidad ( const nVector<long double> & ) ;
void setPosicion ( const nVector<long double> & ) ;
void setAceleracion ( const nVector<long double> & ) ;

const bool operator == ( const nCuerpo & ) const ;
const bool operator != ( const nCuerpo & ) const ;

protected :

private :
long double masa_ ;
nVector<long double> velocidad_ ;
nVector<long double> posicion_ ;
nVector<long double> aceleracion_ ;
};

#endif // _nCuerpo_h_

Code (Text):
#ifndef _nLista_h_
#define _nLista_h_

#include "nNodoLista.h"

template < typename T >
class nLista
{
public :
nLista ( ) ;
nLista ( const nLista<T> & ) ;
~nLista ( ) ;

void addObjeto ( const T & ) ;
void addObjetoEnd ( const T & ) ;
void removeObjeto ( ) ;
// void removeObjetoAtras ( ) ;
const bool findObjeto ( const T & ) const ;
// const bool removeObjeto ( const T & ) ;
T * getObjeto ( const unsigned int = 1 ) const ;
const unsigned int getCantidad ( ) const ;

protected:

private:
nNodoLista<T> * primero_ ;
unsigned int cantidad_ ;

void deleteNodo ( nNodoLista<T> * ) ;

};

#endif // _nLista_h_

Code (Text):
#ifndef _nNodoLista_h_
#define _nNodoLista_h_

template < typename T >
class nNodoLista
{
public :
nNodoLista ( ) ;
nNodoLista ( const T & ) ;
nNodoLista ( const nNodoLista<T> & ) ;
~nNodoLista ( ) ;
T * getObjeto ( ) ;
void setSiguiente ( nNodoLista<T> * ) ;
nNodoLista<T> * getSiguiente ( ) ;

protected :

private :
T dato_ ;
nNodoLista<T> * siguiente_ ;
};

#endif // _nNodoLista_h_

Code (Text):
#ifndef _nSistema_h_
#define _nSistema_h_

#include "nCuerpo.h"
#include "nLista.h"

class nSistema
{
public :
nSistema ( ) ;
nSistema ( const nSistema & ) ;
~nSistema ( ) ;

void addCuerpo ( const nCuerpo & ) ;
void updateSistema1 ( ) ;
//unsigned int cuerpos ( ) { return cuerpos_ ; } ;
void drawSistema1 ( ) const ;

protected :

private :
nLista<nCuerpo> sistema_ ;
nCuerpo centrodemasa_ ;
unsigned int cuerpos_ ;
};

#endif // _nSistema_h_

Code (Text):
#ifndef _nVector_h_
#define _nVector_h_

template < typename T >
class nVector
{
public :
const static int NVECTORSIZE = 3 ;
nVector ( ) ;
nVector ( const T & , const T & , const T & ) ;
nVector ( const nVector<T> & ) ;
~nVector ( ) ;

nVector<T> & operator = ( const nVector<T> & ) ;
nVector<T> & operator += ( const nVector<T> & ) ;
nVector<T> & operator -= ( const nVector<T> & ) ;
nVector<T> & operator *= ( const T & ) ;
nVector<T> & operator /= ( const T & ) ;
const nVector<T> operator + ( const nVector<T> & ) const ;
const nVector<T> operator - ( const nVector<T> & ) const ;
const nVector<T> operator * ( const T & ) const ;
const nVector<T> operator / ( const T & ) const ;
const T operator * ( const nVector<T> & ) const ;
const nVector<T> Xprod ( const nVector<T> & ) const ;

const bool operator != ( const nVector<T> & ) const ;
const bool operator == ( const nVector<T> & ) const ;

const T getNorma2 ( ) const ;
const T getX ( ) const ;
const T getY ( ) const ;
const T getZ ( ) const ;

protected :

private :
T data_ [ NVECTORSIZE ] ;
};

#endif // _nVector_h_

2. May 20, 2010

### rcgldr

The key here is the numerical integration used to calculate velocities and positions, which in turn are used to calculate acclerations for each time step. An iterative corrector method will improve the results. In the algorithm shown below, an, vn, and pn, are successive "guesses" that should converge quickly. F(...) calculates the acceleration based on pn(t), and Δt is the elapsed time per step. You may want to do 6 to 8 interations instead of the 4 shown in this example. The first step is essentially Euler (since a1(t) is set = F(p0(t)) (= a(t-1)), the remaining steps are trapezoidal. Even though each step of this algorithm will converge to a specific set of values, the algorithm is based on trapezoidal rule, a linear approximation, so you need to use small time steps (Δt).

v0(t) = v(t-1)
p0(t) = p(t-1)

a1(t) = F(p0(t)) (= a(t-1))
v1(t) = v(t-1) + 1/2 (a(t-1) + a1(t)) Δt
p1(t) = p(t-1) + 1/2 (v(t-1) + v1(t)) Δt

a2 = F(p1(t))
v2(t) = v(t-1) + 1/2 (a(t-1) + a2(t)) Δt
p2(t) = p(t-1) + 1/2 (v(t-1) + v2(t)) Δt

a3 = F(p2(t))
v3(t) = v(t-1) + 1/2 (a(t-1) + a3(t)) Δt
p3(t) = p(t-1) + 1/2 (v(t-1) + v3(t)) Δt

a4 = F(p3(t))
v4(t) = v(t-1) + 1/2 (a(t-1) + a4(t)) Δt
p4(t) = p(t-1) + 1/2 (v(t-1) + v4(t)) Δt

...

v(t) = vn(t)
p(t) = pn(t)
a(t) = F(pn(t))

time += Δt
t += 1

This is a predictor-corrector type algorithm:

http://en.wikipedia.org/wiki/Predictor-corrector_method#Euler_trapezoidal_example

Last edited: May 20, 2010
3. May 20, 2010

### Kynnath

Thanks! I'll give that a try, see how the program reacts and how many steps I'd need. I'd tried to figure decent ways to correct the error, such as trying to calculate the particle's energy at a given point and cut speed accordingly to keep energy constant, but hadn't gotten to figuring out how that would work.

4. May 20, 2010

### Xitami

http://en.wikipedia.org/wiki/Gravity_assist" [Broken]

Last edited by a moderator: May 4, 2017
5. May 20, 2010

### Kynnath

Right. Energy of the system doesn't change, not energy of each individual particle. As I said, was just an idea, didn't think it through.

6. May 21, 2010

### rcgldr

If you use the integration method I described in my previous post and keep the time step (Δt) small enough, you shouldn't have an issue with errors with a reasonable number of orbital cycles (100 to 1000). As mentioned you could calculate the total energy of the system every 10 to 100 cycles, and if it's increasing or decreasing, apply a small correction factor to the process.