[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... :muro: :muro: :muro: :muro:
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!! :mc: :mc:

Vi posto il codice... Io mi vado uccidere :incazzed: :incazzed: :incazzed:




#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... :stordita: )

[Qu@ker]

Hai provato ad usare long double?

[Jonny32]

Hai provato ad usare long double?

Yes, stesso problema :mc:

[fek]

Yes, stesso problema :mc:

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]

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]

Cosa intendi con range dinamico e precisione vicino allo zero?

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


up!

[Jonny32]

up! :mc: :mc: :mc: :mc:

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