Solving Error Accumulation in Physics Simulation Program

[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

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.

#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_

#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_

#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_

#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_

#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_

 

[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, a_n, v_n, and p_n, are successive "guesses" that should converge quickly. F(...) calculates the acceleration based on p_n(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 a₁(t) is set = F(p₀(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).

v₀(t) = v(t-1)
p₀(t) = p(t-1)

a₁(t) = F(p₀(t)) (= a(t-1))
v₁(t) = v(t-1) + 1/2 (a(t-1) + a₁(t)) Δt
p₁(t) = p(t-1) + 1/2 (v(t-1) + v₁(t)) Δt

a₂ = F(p₁(t))
v₂(t) = v(t-1) + 1/2 (a(t-1) + a₂(t)) Δt
p₂(t) = p(t-1) + 1/2 (v(t-1) + v₂(t)) Δt

a₃ = F(p₂(t))
v₃(t) = v(t-1) + 1/2 (a(t-1) + a₃(t)) Δt
p₃(t) = p(t-1) + 1/2 (v(t-1) + v₃(t)) Δt

a₄ = F(p₃(t))
v₄(t) = v(t-1) + 1/2 (a(t-1) + a₄(t)) Δt
p₄(t) = p(t-1) + 1/2 (v(t-1) + v₄(t)) Δt

...

v(t) = v_n(t)
p(t) = p_n(t)
a(t) = F(p_n(t))

time += Δt
t += 1

This is a predictor-corrector type algorithm:

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

 

[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.

 

[Xitami]

http://en.wikipedia.org/wiki/Gravity_assist"

 

[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.

 

[rcgldr]

ways to correct the error

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.