Solving Error Accumulation in Physics Simulation Program

In summary, the individual is trying to create a program that simulates gravity on a set of point particles. However, the program accumulates a lot of error quickly, causing the particles to move in unexpected ways. The individual is using a set of equations to calculate the force of gravity and acceleration, but is unsure why the error is increasing so rapidly. They are using a few custom classes for the project and are seeking insight on how to fix the problem.
  • #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

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_
 
Technology news on Phys.org
  • #2
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:
  • #3
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
http://en.wikipedia.org/wiki/Gravity_assist"
 
Last edited by a moderator:
  • #5
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
Kynnath said:
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.
 

Related to Solving Error Accumulation in Physics Simulation Program

1. What is error accumulation in physics simulation programs?

Error accumulation in physics simulation programs refers to the gradual accumulation of small errors in the simulation results over time. This can happen due to various factors such as rounding errors, numerical approximations, and limitations of the simulation model.

2. How does error accumulation affect the accuracy of the simulation?

Error accumulation can significantly impact the accuracy of the simulation results. As the errors accumulate, they can deviate the simulation results from the actual physical behavior, leading to inaccurate predictions and conclusions.

3. What techniques are used to reduce error accumulation in physics simulation programs?

To reduce error accumulation, various techniques can be employed such as using higher precision data types, implementing advanced algorithms, and refining the simulation model. Additionally, performing frequent error analysis and validation can also help in identifying and reducing errors.

4. Can error accumulation be completely eliminated in physics simulation programs?

No, error accumulation cannot be completely eliminated. However, it can be minimized by implementing techniques such as error control and correction, using more accurate data and algorithms, and regularly updating and refining the simulation model.

5. What are the consequences of ignoring error accumulation in physics simulation programs?

Ignoring error accumulation in physics simulation programs can lead to inaccurate and unreliable results, which can have serious consequences. It can affect the validity of scientific studies, engineering designs, and other applications that rely on accurate simulations, potentially causing financial losses, safety hazards, or other significant issues.

Similar threads

Why use polymorphism for grading activities?
    • Programming and Computer Science
2
Replies
36
Views
2K
Polymorphism using pointer to array
    • Programming and Computer Science
Replies
23
Views
1K
My destructor is causing my program to crash
    • Programming and Computer Science
Replies
3
Views
2K
Monte-Carlo integration does not work properly after implementing pdfs
    • Programming and Computer Science
Replies
1
Views
659
Why is my program crashing when using multiple classes and pointers?
    • Programming and Computer Science
Replies
5
Views
833
Understanding Inheritance and Polymorphism in C++
    • Programming and Computer Science
Replies
17
Views
1K
C/C++ Question about New syntax of C++
    • Programming and Computer Science
Replies
31
Views
2K
C/C++ Question in throwing a class object as exception in C++
    • Programming and Computer Science
2
Replies
35
Views
2K
Question on this overload program
    • Programming and Computer Science
3
Replies
89
Views
4K
Creating a 2D Vector Using C++
    • Programming and Computer Science
Replies
13
Views
1K

Hot Threads

Recent Insights

Back
Top