- #1
Kynnath
- 17
- 0
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
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.
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:
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:
#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:
#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:
#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:
#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:
#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_