Newsgroup: it.scienza.fisica
Titolo: simulazione del sistema solare

Salve a tutti,
chiedo un aiuto circa una simulazione di un sistema solare al computer, non
di carattere informatico però :-)
Poiché non posso inserire i valori veri delle masse e dei raggi delle orbite
dei pianeti (per via della capacità del calcolatore di trattare numeri come
miliardi di miliardi), ho scalato insieme questi valori e la costante G
secondo la terza legge di Keplero, per mantere costante il periodo di
rotazione:
T²= r³/(G*(M+m))
Perciò ho diviso le distanze per 1 milione di Km, questo mi ha ridotto il
periodo di (10^6)^3 , perciò ho ridotto le masse di questo fattore essendo
sicuro di mantenere il periodo di rotazione costante. Quindi ho aumentato G
di 10^5 e ridotto ancora le masse di  10^5.
E' corretto procedere così per far rimanere costante il periodo di rotazione
dei pianeti attorno al sole?
Penso proprio di si, ma soprattutto serve questo ragionamento di far
rimanere costante T per simulare il sistema solare con piccoli numeri? O
dovrei semplicemente ridurre le distanze, e le masse al quadrato, lasciano
intatta G?
Grazie del vostro aiuto,
Daniele

"(¯`·.¸  D²  ¸.·´¯)" ha scritto:
 > Poiché non posso inserire i valori veri delle masse e dei raggi delle
 > orbite dei pianeti (per via della capacità del calcolatore di trattare
 > numeri come miliardi di miliardi),
Ci sono parecchie cose che non capisco nel tuo problema, ma la prima e
piu' fondamentale e' questa affermazione.
Che linguaggio stai usando?
Tutti i linguaggi, ormai da decenni, permettono l'uso del cosiddetto
"floating point", il che vuol dire che puoi benissimo usare numeri non
solo di miliardi di miliardi, ma amche molto piu' grandi.
Temo quindi che il vero problema stia nelle tue conoscenze, forse un
po' rudimentali... Sbaglio?
 > ...
 > Perciò ho diviso le distanze per 1 milione di Km, questo mi ha ridotto
 > il periodo di (10^6)^3 , perciò ho ridotto le masse di questo fattore
 > essendo sicuro di mantenere il periodo di rotazione costante. Quindi
 > ho aumentato G di 10^5 e ridotto ancora le masse di 10^5.
Io avrei fatto qualcosa di piu' usuale, come esprimere le distanze in
UA, le massa in masse solari (a proposito, ma sei sicuro che sia
importante tener conto delle masse dei pianeti? vedi dopo) e i tempi in 
anni.
(Tra parentesi, nella formula che haiscritto manca un 4pi^2...)
Il mio secondo problema e' pero' che non capisco che cosa intendi per
"simulazione del sistema solare".
Che cosa ne vuoi fare? Che grado di precisione vuoi ottenere?
Vuoi fare un presentazione grafica, vuoi dei risultati numerici, o che
cosa?
Perche' se ho capito bene stai trattando le orbite come circolari, il
che non e', e stai trascurando le interazioni tra i pianeti.
Quindi otterrai una simulazione molto molto grossolana...
Puo' darsi che per il tuo scopo basti, ma appunto: qual e' il tuo
scopo?
------------------------------
Elio Fabri
Dip. di Fisica - Univ. di Pisa
------------------------------

 T²= r³/(G*(M+m))
[sforbiciata]
Sono solo di passaggio e faccio qualche rapida osservazione.
Nella formula che rappresenta la terza legge di Keplero manca a sinistra un fattore
(2 Pi)^-2.
Che linguaggio di programmazione usi, non disponi di variabili floating point
in precisione multipla?
Se proprio devi "macinare" numeri piccoli, mi sembra che la miglior cosa sia
cambiare unità di misura, ad es. misurando le distanze in unita' astronomiche (UA),
il tempo in anni e le masse in masse solari, le grandezze T, r, G, (M + m) saranno
rappresentate da numeri all'incirca dell'ordine dell'unita'.
Se, come suppongo, risolverai le equazioni del moto con un procedimento iterativo,
con questa scelta delle unita' di misura gli incrementi delle coordinate spaziali e
temporale corrisponderanno a numeri comunque piccoli rispetto all'unita'.
Considera pero' che se usi nei calcoli variabili in un formato con poche cifre
significative, rischi dopo diverse iterazioni di perdere precisione nei risultati.
Ciao
--
Giorgio Bibbiani
Per rispondere inserire una sottolineatura
tra nome e cognome

Rispondo a tutti:
Spiego meglio il funzionamento del programma, così chiarisco molte idee.
Ho fatto in C un risolutore di equazioni differenziali col metodo iterativo
di Runge-Kutta del secondo ordine, potendo scegliere un intervallo di
integrazione a mio piacimento, che per quanto riguarda la precisione sembra
soddisfacente, ma forse dovrei ragionarci meglio, per calcolare
effettivamente quanto errore mi da.
Immetto come input le posizioni x,y e le velocità x,y di ogni pianeta, e la
sua massa.
La forza gravitazionale agisce su ogni pianeta come nella realtà, tutto è
ben simulato, il sole è come un pianeta. Sostanzialmente è il problema degli
n corpi con un potenziale gravitazionale, dal quale ottengo come risultato
un file contenente le posizioni di ogni pianeta, da plottare.
Volendo simulare le condizioni reali del nostro sistema, dovrei immettere i
giusti rapporti delle masse e le giuste distanze, ma dopo vari tentativi
scopro che non si instaura un orbita chiusa. Questo solo inserendo un
pianeta come giove e il sole. Trovo strano che il giove si porti appresso il
sole nel suo viaggio :-)
In altre simulazioni che vedo infatti il sole si muove di molto poco,
verosimilmente.
Il problema però si evolve: anche se potessi inserire i valori veri (penso
che sia un inutile appesantimento del carico di lavoro del pc) io devo
inserire la velocità iniziale del pianeta (per comodità posto al perielio e
con componente della velocità solo lungo le y) che si esprime facendo uso
del tempo. Ora, visto che non ho un modo di controllare il tempo nel mio
programma, cioè non so come simulare il tempo che passa; penso quindi che la
velocità non abbia un significato reale, allora non potrei usare quella
reale per ottenere una simulazione veritiera: infatti non la ottengo.
Perciò mi chiedo: posso inserire il valore di circa 13,1 per la velocità di
giove? (Km /s) Ma se esprimessi tutto in anni cambirerebbe questo dato
numerico, variando così l'orbita, perché non ho modo di cambiare l'unità di
misura del tempo da nessun' altra parte. Ecco perché dico che l'informazione
della velocità non si applica bene al mio programma.
Spero di essermi spiegato bene, in caso contrario chiedetemi ulteriori
spiegazioni.
Grazie di nuovo dell'ascolto
Daniele

"(¯`·.¸  D²  ¸.·´¯)" ha scritto:
 > Spiego meglio il funzionamento del programma, così chiarisco molte
 > idee.
 > Ho fatto in C un risolutore di equazioni differenziali col metodo
 > iterativo di Runge-Kutta del secondo ordine, potendo scegliere un
 > intervallo di integrazione a mio piacimento, che per quanto riguarda
 > la precisione sembra soddisfacente, ma forse dovrei ragionarci meglio,
 > per calcolare effettivamente quanto errore mi da.
Si', ora e' un po' piu' chiaro...
Anche se io forse non mi sarei fermato alsecondo ordine: ho sempre
usato il quarto.
Per vedere che errore ti da', c'e' un modo molto semplice: prova per
es. a dimezzare il passo, e confronta i risultati.
Se sono suff. vicini, sei a posto.
BTW: Dunque il problema delle unita' e delle cifre l'hai risolto? Se
programmi in C, non vedoproprio che difficolta' ci fosse a usare
numeri floating, meglio in doppia precisione.
 > Immetto come input le posizioni x,y e le velocità x,y di ogni pianeta,
 > e la sua massa.
Beh, solo x,y? Le orbite non sono mica complanari...
 > La forza gravitazionale agisce su ogni pianeta come nella realtà,
 > tutto è ben simulato, il sole è come un pianeta. Sostanzialmente è il
 > problema degli n corpi con un potenziale gravitazionale, dal quale
 > ottengo come risultato un file contenente le posizioni di ogni
 > pianeta, da plottare.
OK, con una riserva che dico fra poco.
 > Volendo simulare le condizioni reali del nostro sistema, dovrei
 > immettere i giusti rapporti delle masse e le giuste distanze, ma dopo
 > vari tentativi scopro che non si instaura un orbita chiusa. Questo
 > solo inserendo un pianeta come giove e il sole. Trovo strano che il
 > giove si porti appresso il sole nel suo viaggio :-)
 > In altre simulazioni che vedo infatti il sole si muove di molto poco,
 > verosimilmente.
Di solito quando si trovano orbite nn chiusa, la colpa sta
nell'insuff. approssimazione dell'integrazione, ma nel tuo caso ho un
altro sospetto...
Se nel problema di due corpi parti col Sole fermo e dai una velocita'
al pianeta, di sicuro non troverai orbite chiuse, perche' il centro di
massa non sta fermo.
Non mi meraviglia che l'effetto si noti di piu' con Giove, visto che ha
massa piu' grande.
Infatti questi conti si fanno sempre in _coordinate relative_: si
assume fermo il centro di massa, e studia il moto del pianeta rispetto
al Sole.
 > Il problema però si evolve: anche se potessi inserire i valori veri
 > (penso che sia un inutile appesantimento del carico di lavoro del pc)
 > io devo inserire la velocità iniziale del pianeta (per comodità posto
 > al perielio e con componente della velocità solo lungo le y) che si
 > esprime facendo uso del tempo. Ora, visto che non ho un modo di
 > controllare il tempo nel mio programma, cioè non so come simulare il
 > tempo che passa; penso quindi che la velocità non abbia un significato
 > reale, allora non potrei usare quella reale per ottenere una
 > simulazione veritiera: infatti non la ottengo.
Sinceramente qua non ti capisco prorio: come sarebbe che nn sai come
simulare il tempo?
Hai scritto delle eq. diff. dove la var. indip. e' appunto il tempo;
ha fatto un'intgrazione discreta, in cui avrai per forza dovuto
scegliere un passo temporale.
Non so come ti sei comportato per le unita', ma come ho gia' detto
potevi liberamente usare unita' del SI oppure unita' astronomiche; non
fa nessuna differenza, purche' usi un sistema coerente.
Comunque sia, il tuo passo temporale avra' un certo valore, poniamo 10
giorni; basta che conti il numero dei passi, e vedi come "passa" il
tempo...


 numeri floating, meglio in doppia precisione.
...
Beh se non si sta attenti con l' ordine delle operazioni, coi float a 32 
bit e l' ormai standard rappresentazione IEEE754 il problema c'e': 
p.es
. 
se cerchi di moltiplicare la massa del sole per quella della terra 
(unita SI) trovi 2*10^30 * 6*10^24 che va in overflow in quanto supera 
il piu' grande float rappresentabile ( 3.4*10^38 ).
Naturalmente basta un po' di attenzione nella scritura del programma per 
evitare questi problemi.
Giorgio

 effettivamente quanto errore mi da.
...
Qualche commento in piu' oltre quelli che hai gia' ricevuto.
1. una misura diretta dell' errore di integrazione la puoi avere 
controllando la bonta' con cui le costanti del moto sonon conservate. 
Non aspettarti delle vere costanti ma, 
p.es
., se guanrdi l' energia 
totale dovresti avere fluttuazioni almeno un paio di ordini di grandezza 
minori delle fluttuazioni delle singole componenti dell' energia 
cinetica o potenziale del sistema (guardando ad almeno una rivoluzione).
2. Per quanto molto popolari tra i matematici numerici, i Runge Kutta 
non sono necessariamente i metodi numerici migliori.
Una possibilita' e' data da metodi adattativi alla Adam-Bashford (spero 
di non aver sbagliato troppo lo spelling) in cui il passo di 
integrazione non e' costante ma viene variato fino a raggiungimento di 
una precisione assegnata. Un 'altra (piu' semplice da programmare  e con 
notevoli vantaggi su altri fronti) e' data dai cosiddetti integratori 
"simplettici". Il piu' semplice dei quali e' dato dalle seguenti formule 
(per 1 grado di liberta', la generalizzazione a piu' gr di liberta' e' 
banale):
data l' accelerazione, velocita' e posizione al tempo t (a(t),v(t),x(t))
x(t+d) = x(t) + d*v(t) + d^2*a(t)/2
vp(t+d)= v(t) + d*a(t)/2
si calcola a(t+d) (funzione solo delle x(t+d) appena ottenute)
e infine
v(t+d) = vp(t+d) + d*a(t+d)/2
e quindi si itera per ottenete x,v,e a ai tempi t+2*d, t+3*d ....
Un algoritmo come il precedente ha due enormi vantaggi su molti altri 
(RK inclusi):
a. e' reversibile temporalmente;
b. e' in grado di garantire la conservazione della quantita' di moto 
totale e del  momento angolare totale  a precisione macchina (l' errore 
e' unicamente dovuto al roundoff e non alla discretizzazione).
3. naturalmente la simulazione su tempi lunghi e/o con oggetti con forte 
  variazione delle distanze sono intrinsecamente difficili.
Giorgio