[C++] Errore di calcolo impossibile...

[Jonny32]

Dopo 3 giorni di prove mi arrendo...
Il programma calcola le posizioni di 4 stelle che sono disposte in modo simmetrico in un immaginario piano cartesiano 2D. Il fatto è che disposte in questo mode le 4 stelle dovrebbero conservare, teoricamente la loro simmetria sia nelle posizioni, sia nelle velocità sia nelle accelerazioni...

Invece che succede? Il programma sbaglia il calcolo, piano piano sclera...
E nn è che sbaglia di tanto (come se mi fossi scordato di considerare una stella o di aggiornare una variabile...), sbaglia all'ultima cifra decimale!!

Vi posto il codice... Io mi vado uccidere

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
//#include <iostream>

#define G		6.6725985E-11				// Costante di gravitazione universale
#define TI		20	
using namespace std;

class Vect
{
public:
	double x, y, module;
};

class Bodies
{
public:
	double mass;
	Vect acc;
	Vect pos;
	Vect vel;
};

class Universe
{
public:
	int n;
	Bodies * body;
};

Universe * u;


void loop ( )
{
	int i, j;

	// Calcolo delle accelerazioni
	Vect dist;
	if ( u->n != 1 ) {
		for ( i = 0; i < u->n; i++ ) {
			u->body[i].acc.x = u->body[i].acc.y = 0;
			dist.x = dist.y = 0;
			for ( j = 0; j < u->n; j++ ) {
				if ( j != i ) {
					dist.x = u->body[i].pos.x - u->body[j].pos.x;
					dist.y = u->body[i].pos.y - u->body[j].pos.y;
					dist.module = sqrt ( dist.x * dist.x + dist.y * dist.y );
					if ( dist.module != 0 ) {
						u->body[i].acc.x -= u->body[j].mass * G * dist.x / 
								( dist.module * dist.module * dist.module );
						u->body[i].acc.y -= u->body[j].mass * G * dist.y / 
								( dist.module * dist.module * dist.module );
					}
				}
			}
		}
	}

	// Controlli anti-sclero
	printf ( "[0].acc.x = %.18lf\n[1].acc.x = %.18lf\n[2].acc.y = %.18lf\n[3].acc.y = %.18lf\n",
			 u->body[0].acc.x, -u->body[1].acc.x, u->body[2].acc.y, -u->body[3].acc.y );
	printf ( "[0].acc.y = %.18lf\n[1].acc.y = %.18lf\n[2].acc.x = %.18lf\n[3].acc.x = %.18lf\n",
			 u->body[0].acc.y, -u->body[1].acc.y, -u->body[2].acc.x, u->body[3].acc.x );
	system ("Pause");
	// Fine Controlli
	
	for ( i = 0; i < u->n; i++ ) {
		u->body[i].pos.x += u->body[i].vel.x * TI + u->body[i].acc.x * TI * TI / 2;
		u->body[i].pos.y += u->body[i].vel.y * TI + u->body[i].acc.y * TI * TI / 2;
		u->body[i].vel.x += u->body[i].acc.x * TI;
		u->body[i].vel.y += u->body[i].acc.y * TI;
	}
}

int main ( ) 
{
	u = new Universe;

	u->n = 4;
	u->body = new Bodies [u->n];

	u->body[0].mass = 2E+30;
	u->body[0].pos.x = - 1E+10;
	u->body[0].pos.y = 0;
	u->body[0].vel.x = 0;
	u->body[0].vel.y = -100000;
	
	u->body[1].mass = 2E+30;
	u->body[1].pos.x = 1E+10;
	u->body[1].pos.y = 0;
	u->body[1].vel.x = 0;
	u->body[1].vel.y = 100000;
	
	u->body[2].mass = 2E+30;
	u->body[2].pos.x = 0;
	u->body[2].pos.y = - 1E+10;
	u->body[2].vel.x = 100000;
	u->body[2].vel.y = 0;
	
	u->body[3].mass = 2E+30;
	u->body[3].pos.x = 0;
	u->body[3].pos.y = 1E+10;
	u->body[3].vel.x = -100000;
	u->body[3].vel.y = 0;
	
	for ( int i = 0; i < 1; ) 
	{
		loop ( );
	}
	
	system ("Pause");
	return 0;
}

Nn spaventatevi per la lunghezza del prog... nel main praticamente richiamo solo la funzione loop (oltre a dare i valori alle variabili... )

[Qu@ker]

Hai provato ad usare long double?

[Jonny32]

Originariamente inviato da Qu@ker

Hai provato ad usare long double?

Yes, stesso problema

[fek]

Originariamente inviato da Jonny32

Yes, stesso problema

I numeri in floating point sono approssimazioni e come tali si portano dietro degli errori che aumentano ad ogni operazione. Se usi i risultati del calcolo come input per l'iterazione successiva, gli errori si propagano e si amplifcano perche' hai un feedback loop.

Hai due strade:

1) cerchi di calcolare gli errori e ne tieni conto nel risultato, di modo da cancellarli; non e' sempre possibile

2) fai i calcoli in virgola fissa usando gli interi, cosi' hai la sicurezza di non avere approssimazioni, ma perdi in range dinamico (non credo ti serva) e precisione vicino allo zero

Nel tuo caso sceglierei la seconda soluzione.

C'e' ancora una terza strada, ma dipende dall'applicazione: introduci "atrito" nel calcolo della tua fisica. L'atrito tende a dissipare energia e portare il sistema ad uno stato di equilibrio, quindi, col tempo, l'energia introdotta dagli errori di approssimazione viene anch'essa attenuata e il sistema si stabilizza senza esplodere.

[Jonny32]

Originariamente inviato da fek

I numeri in floating point sono approssimazioni e come tali si portano dietro degli errori che aumentano ad ogni operazione. Se usi i risultati del calcolo come input per l'iterazione successiva, gli errori si propagano e si amplifcano perche' hai un feedback loop.

Hai due strade:

1) cerchi di calcolare gli errori e ne tieni conto nel risultato, di modo da cancellarli; non e' sempre possibile

2) fai i calcoli in virgola fissa usando gli interi, cosi' hai la sicurezza di non avere approssimazioni, ma perdi in range dinamico (non credo ti serva) e precisione vicino allo zero

Nel tuo caso sceglierei la seconda soluzione.

C'e' ancora una terza strada, ma dipende dall'applicazione: introduci "atrito" nel calcolo della tua fisica. L'atrito tende a dissipare energia e portare il sistema ad uno stato di equilibrio, quindi, col tempo, l'energia introdotta dagli errori di approssimazione viene anch'essa attenuata e il sistema si stabilizza senza esplodere.

Quello di cui parti tu è un errore di approssimazione, che bene o male è genuino nel programma (nel senso che nn è colpa mia se c'è)...
Questo errore qui è qualcos'altro...

Guardate questi due eseguibili:

http://edt.altervista.org/con_printf.exe

http://edt.altervista.org/senza_printf.exe

Nel primo ho inserito, mentre facevo delle prove, dei printf nel ciclo delle accelerazioni...
Il secondo invece è normale.

Come potete notare i valori delle accelerazioni stampati ( quelli infondo ) sono diversi nei due programmi... Come mai?

I codice modificato (in cui ho aggiunto i printf) è questo:

...
if ( j != i ) {
	dist.x = u->body[i].pos.x - u->body[j].pos.x;
	printf ("dist.x = %.18lf\n", dist.x);
	dist.y = u->body[i].pos.y - u->body[j].pos.y;
	printf ("dist.y = %.18lf\n", dist.y);
	dist.module = sqrt ( dist.x * dist.x + dist.y * dist.y );
	printf ("dist.m = %.18lf\n", dist.module);
	if ( dist.module != 0 ) {
		u->body[i].acc.x += - u->body[j].mass * G * dist.x / 
					( dist.module * dist.module * dist.module );
		printf ("acc.x = %.18lf\n", u->body[i].acc.x);
		u->body[i].acc.y += - u->body[j].mass * G * dist.y / 
					( dist.module * dist.module * dist.module );
		printf ("acc.y = %.18lf\n", u->body[i].acc.y);
	}
}
...

[Jonny32]

2) fai i calcoli in virgola fissa usando gli interi, cosi' hai la sicurezza di non avere approssimazioni, ma perdi in range dinamico (non credo ti serva) e precisione vicino allo zero

Cosa intendi con range dinamico e precisione vicino allo zero?

Potrebbe essere una buona soluzione per l'errore di approssimazione...

[Jonny32]

Originariamente inviato da Jonny32

Cosa intendi con range dinamico e precisione vicino allo zero?

Potrebbe essere una buona soluzione per l'errore di approssimazione...

up!

[Jonny32]

up!

[Jonny32]

up

[andreaM]

Ciao, scusa ma ad ogni loop devono essere sempre visualizzati sempre gli stessi valori?

Questo sarebbe il primo loop, ti sembrano corretti i valori ?

[0].acc.x = 1.277277854497037124
[1].acc.x = 1.277277854497037124
[2].acc.y = 1.277277854497037124
[3].acc.y = 1.277277854497037124
[0].acc.y = -0.000000000000000026
[1].acc.y = 0.000000000000000026
[2].acc.x = 0.000000000000000026
[3].acc.x = -0.000000000000000026

I valori vicino a zero ma non zero mi lasciano perplesso.

Ciao.