Perturbazioni in Meccanica Celeste *Stefano Mandelli 14 settembre 2010 *Osservatorio Astronomico Città di Seveso e-mail: [email protected]Sommario In questo breve articolo si presenteranno alcune nozioni teoriche sull’ultimo teorema di Poincaré che saranno completate con delle applicazioni teoriche di particolare interesse come per esempio: Standard Map e Twist Map. Successivamente si porterà l’attenzione sulla risonanza 3:1 nella fascia di ateroidi principale. Vedremo che la risonanza con Giove porta gli asteroidi presenti in quella zona ad avere una evoluzione della loro orbita caotica. Non potendo usare una mappa di Poincaré e non potendo usare altri tipi mapping, per mettere in luce questa caratteristica di caoticità delle orbite si opererà in questo modo: Scrittura di un software del problema N-Corpi. Tramite l’uso del software simulare l’evoluzione dell’orbita di un asteroide posizionato nella fascia di risonanza 3:1. Plottare l’andamento di eccentrità e semiasse maggiore in funzione del tempo e dare le relative conclusione. Successivamente sfuttare anche le mappe di A.Morbidelli del 1996 per dare maggior senso a quanto ricavato dalla simulazione numerica. 1 L’ultimo teorema di Poincaré 1.1 Introduzione Verranno introdotti degli esempi molto particolari per capire la vera essenza del problema. Il percorso si svolgerà attraverso la presentazione dello standard map classico imperturbato. Suc- cessivamente si discuterà il caso perturbato cercando di trarre alcune conclusioni. Verrà presentato anche l’esempio della twist map che ci introdurrà all’ultimo teorema di Poincaré e al teorema della varietà stabile. 1.2 Lo standard map E’ una trasformazione puntuale estesa, prossima all’identità (quindi canonica) del toro bidimen- sionale in se stesso, dipendende da un parametro ∈ R x 0 = x + y + sin(2πx) (mod1) y 0 = y + sin(2πx) (mod1) (1) Riscrivibile come: x 0 = x + y 0 (mod1) y 0 = y + sin(2πx) (mod1) (2) Questo modo di riscriverla mette meglio in evidenza che lo standard map è una trasformazione canonica, infatti è facile trovare la generatrice della trasformazione che risulta essere: F (y 0 ,x)= y 0 x + 2π cos(2πx)+ (y 0 ) 2 2 (3) 1
SommarioIn questo breve articolo si presenteranno alcune nozioni teoriche sull’ultimo teorema di
Poincaré che saranno completate con delle applicazioni teoriche di particolare interesse comeper esempio: Standard Map e Twist Map. Successivamente si porterà l’attenzione sullarisonanza 3:1 nella fascia di ateroidi principale. Vedremo che la risonanza con Giove portagli asteroidi presenti in quella zona ad avere una evoluzione della loro orbita caotica. Nonpotendo usare una mappa di Poincaré e non potendo usare altri tipi mapping, per mettere inluce questa caratteristica di caoticità delle orbite si opererà in questo modo: Scrittura di unsoftware del problema N-Corpi. Tramite l’uso del software simulare l’evoluzione dell’orbitadi un asteroide posizionato nella fascia di risonanza 3:1. Plottare l’andamento di eccentritàe semiasse maggiore in funzione del tempo e dare le relative conclusione. Successivamentesfuttare anche le mappe di A.Morbidelli del 1996 per dare maggior senso a quanto ricavatodalla simulazione numerica.
1 L’ultimo teorema di Poincaré
1.1 Introduzione
Verranno introdotti degli esempi molto particolari per capire la vera essenza del problema. Ilpercorso si svolgerà attraverso la presentazione dello standard map classico imperturbato. Suc-cessivamente si discuterà il caso perturbato cercando di trarre alcune conclusioni. Verrà presentatoanche l’esempio della twist map che ci introdurrà all’ultimo teorema di Poincaré e al teoremadella varietà stabile.
1.2 Lo standard map
E’ una trasformazione puntuale estesa, prossima all’identità (quindi canonica) del toro bidimen-sionale in se stesso, dipendende da un parametro ε ∈ R
x′ = x+ y + ε sin(2πx) (mod1)y′ = y + ε sin(2πx) (mod1)
(1)
Riscrivibile come:
x′ = x+ y′ (mod1)y′ = y + ε sin(2πx) (mod1)
(2)
Questo modo di riscriverla mette meglio in evidenza che lo standard map è una trasformazionecanonica, infatti è facile trovare la generatrice della trasformazione che risulta essere:
F (y′, x) = y′x+ε
2πcos(2πx) +
(y′)2
2(3)
1
che soddisfi le equazioni di hamilton:
p =∂F
∂q, Q =
∂F
∂P(4)
dove la trasformazione Q = q, P = p è una trasformazione prossima all’identità. 1 Quello che sivuole fare di questo modello è lo studio delle orbite nello spazio delle fasi. La tecnica per costruirele orbite nello spazio delle fasi è quello (per lo standard map) del time mapping cioè considerareil punto iniziale x0 e vedere come esso evolve al tempo t atraverso un mapping di questo genere:
gt : Rn → Rn (6)
Nel caso dello standard map, si utilizza il time-one map cioè una trasformazione del tipo:
g1 : Rn → Rn (7)
che consiste nel mappare al tempo 1 tutti i vari punti iniziali x0. Poi potrei partire da un ben notopunto iniziale e a quel punto avrei una successione di punti {xn} che mi identificano il movimento.Questo è un tipo di mapping, oltre a questo ve ne sono altri due:
• Time-1 map;
• Dipendeze periodiche dal tempo;
• Sezione di Poincaré.
Nel caso imperturbato con ε = 0 prendendo vari dati iniziali sulle y, notiamo che la loro evoluzioneè lineare. Il caso con ε > 0 mostra invece delle sensibili diversità. In particolare nei punti (0,0) e(1/2,0). Le rette deformate che hanno preso il posto di quelle continue del caso integrabile sonosempre dei tori KAM, invece la natura delle orbite cambia sensibilmente in prossimità dei puntiprima citati. Il punto P1(0, 0) per ogni valore di ε si può mostrare che è un punto fisso iperbolicoe quindi il fenomeno di caoticità sarà legato al famoso punto omoclino di Poincaré.Analizziamo in modo profondo lo standard map imperturbato. Dalle osservazioni che verrannofatte cercheremo di trovare somiglianze e sostanziali differenze col caso perturbato. Per ε = 0 ilmio mapping diventa:
y′ = y , x′ = x+ y (mod 1) (8)
Il valore in y resta inalterato, infatti fissato il dato inziale y0, y = y0 è una costante del moto(funzione invariante del moto). Una costante del moto è detta ance funzione invariante se: datoil flusso legato all’equazione differenziale si ha che:
F (gt(x)) = F (x) (9)
Per tutti i punti iniziali x e per tutti i tempi t. In altri termini il valore di F è costante sull’orbita.E’ interessante considerare le superfici di livello di F (perchè indicano, nelle fasi piatte, che sicomporta come una costante del moto ?). Se definiamo la varietà:
Γc : F (x) = c (10)1Cos’è un toro: T 1 non è altro che la retta reale in cui si identificano i punti che differiscono di un intero quindi
il toro T 1 è il segmento [0, 1]. Il toro T 2 è il piano reale in cui si identificano i punti le cui cordinate differisconoper interi. In questo caso è il quadrato unitario: [0, 1]× [0, 1] con i lati opposti ”incollati”. Tipicamente si ha a chefare con un toro quando le variabili in gioco sono angoli. Dev’essere chiaro che se prendiamo α ∈ R : |α| < 1 ecalcoliamo la seguente trasformazione:
x′ = x+ α (5)
Vuol dire prendere un x ∈ [0, 1] calcolare x′ = x+ α e poi eventualmente riportarlo ancora in [0, 1] aggiungendo osottraendo 1.
2
allora questa è una varietà invariante ed è costituita da orbite. Se ogni punto dell’orbita appartienea Γc e si ha che:
Rn =⋃c
Γc (11)
allora dico che lo spazio delle fasi è stratificato o in lnguaggio tecnico fogliato. Il GRANDERISULTATO di Kolmogorov fu proprio quelli di mettere in luce il fatto che a differenza deisistemi integrabili le superfici invarianti ci sono ancora e sono delle piccole deformazioni di quelleprecedenti ma non si mostrano per tutti i valori di c si presentano soltanto per valori ben precisi.Questo produce nello spazio delle fasi delle isole di stabilità, in cui tori KAM vengono mappati intori KAM. Al di fuori delle isole di stabilità il moto è caotico.
1.2.1 Trasformazione del toro di Arnol’d-Avez
Sia x ∈ T 1. Sia Φ definita come:
Φ : x′ = x+ α (12)
Ogni punto inizale x0 da luogo ad orbite periodiche se e solo se α = pq con p, q interi : q > 0
se α è irrazionale non abbiamo più orbite periodiche ma abbiamo un’orbita densa sul toro T 1.Dimstriamo questa ultima proprosizione.
Teorema: sia data la trasformazione Φ definita precedentente con costante α irrazionale, al-lora si ha un’orbita nello spazio delle fasi densa sul toro T 1
Dimostrazione: Suddividiamo l’intervallo [0, 1) in q intervalli uguali. Prendiamo un punto ini-ziale x0, considero i primi q+1 puntidell’orbita. Poichè sono tutti distinti ne esistono almenodue (chiamiamoli il punto xn e xm relativi all’interazione n ed m del punto x0) ad un medesimosubintervallo e quindi distano tra loro di una quantità più piccola di 1/q:
|Φnx0 − Φmx0| = ε <1q
(13)
riscalo tutto ponendo p = m− n:
|x0 − Φpx0| = ε <1q
(14)
la trasformazione conserva le distanze tra i due punti, quindi la relazione vale per tutti i punti:
|x(k+1)p − xkp| < ε <1q
(15)
Per l’arbitrarietà di q è dimostrata la densità dell’orbita sul toro T 1
q.e.d.
Questa proprietà del sistema (cioè quella di dare luogo ad un’orbita densa su T 1) viene spessoracchiusa nel termine ergodicità topologica. Questa proprietà dei sistemi ergodici si osserveràpiù avanti avere radici ben più profonde. Se consideriamo un insieme aperto I e un sistema iteraton volte, denotiamo con N(I, n) il numero di volte che l’orbita attraversa I, si può dimostrare che:
N(I, n)n
→ m(I) per n→ +∞ (16)
In questo modo possiamo dire che l’orbita è asintoticamente distribuita sul toro T 1.
3
1.2.2 Standard map imperturbato - Conclusioni
Tutto qanto detto fin’ora è relativo al caso dello standard map imperturbato. Riassumendoabbiamo visto che:
i. La dinamica del toro T 2 viene descritta sullo spazio delle fasi da tori T 1 invarianti fogliati;
ii. Esiste una costante del moto ed è F (x, y) = y. Ogni toro inviariante corrisponde la costantedel modo y=c;
iii. Ogni orbita è una traslazione di una fattore α;
vi. Se α ∈ Q allora ho orbite periodiche;
v. se α ∈ R \Q allora ho orbite dense.
1.3 Standard Map Perturbato
La dinamica più facile da studiare e quella ovviamente in cui non c’è moto, quindi andare astudiare i punti fissi2 del mapping Φ. I punti fissi dello standard map imperturbato abbiamo vistoche era tutto il toro y = 0 infatti la trasformazione essendo definita come:
Φ0 =
{x′ = x+ y
y′ = y(17)
per y = 0 otteniamo x = x quindi l’identità 0 = 0 quindi abbiamo un continuo di punti fissi sulsegmento [0, 1] di y = 0. Nel caso perturbato invece le cose cambiano. Il nostro mapping chechiameremo Φε per quanto sia piccola la perturbazione il continuo di punti fissi si ”sgretola” dandoorigine ad un diagramma delle fasi costituito da due soli punti fissi. Il mapping è:
Φε =
{y + ε sin(2πx) = y +m
x+ y + ε sin(2πx) = x+ n(18)
Dalla prima equazione ho che:
sin(2πx) = 0 x ∈ [0, 1] =⇒ x = 0 (x = 1) e x =12
(19)
La seconda equazione ci da y = n che riducendosi al quadrato fondamentale3 equivale a considerarey = 0. In questo modo i due soli punti fissi sono:
P1(0, 0) = (1, 0) ; P2 =(
12, 0)
(20)
Si osserverà subito che la natura dei due punti fissi considerati è completamente diversa da quelliprecendenti. Prima avevamo punti fissi di tipo parabolico mentre ora si osserverà che P1 è ditipo iperbolico mentre P2 è di tipo ellittico. Alla base di questo sgretolamento dei punti fissi c’èil teorema di Poincaré che metterà in luce le vere cause profonde che ne sono alla base. Ora persemplificare il problema, lo linearizziamo , nel senso che vado a studiare una trasformazione Φl incui la Φε è sostituita dalla sua matrice Jacobiana valutata nel suo punto fisso:
Φl(x) = Φε(x+ h) = Φ0(x) +∂Φε∂x|x · h+R (21)
2Un punto fisso di una trasformazione (o un punto critico di una equazione differenziale) è un punto tale percui si annulla il campo vettoriale, quindi x è punto fisso se f(x) = 0
3he?
4
1.3.1 Trasformazioni lineari simplettiche piane
Sia data una trasformazione lineare Φ identificata dalla matrice A. Supponiamo che A ∈M2×2 eche detA = 1 Se gli autovalori di A sono λ1 e λ2 la relazione sul determinate unitario può essereespressa come:
λ1λ2 = 1 (22)
Essendo detA = 1 allora la trasformazione Φ conserva le aree nello spazio delle fasi4. In fun-zione alla natura degli autovalori (complessi o reali) abbiamo questa prima classificazione dellatrasformazione Φ:
i. Caso iperbolico λ1, λ2 ∈ R e sono entrambi positivi o negativi;
ii. Caso ellittico λ1, λ2 ∈ C : =λ1 6= 0,=λ2 6= 0,=λ1=λ2 < 0;
iii. Caso parabolico questo è un caso molto particolare e si verifica quando i due autovalorisono uguali. λ1 = λ2 Questa uguaglianza implica che gli unici due valori che possonoassumere λ1 e λ2 sono ±1
Il nome iperbolico, ellittico e parabolico non sono scelti a caso ma sono in funzione al fatto chein tutti e tre i casi gli iterati di un certo (arbitrario) punto iniziale giacciono su curve invariantiche nel primo caso sono iperboli, nel secondo ellissi e nel terzo parabole. Alcuni esempi che sipossono portare sono i seguenti:
Caso1: Un sistema tale che λ1 = 2 e λ2 = 12 lo spazio delle fasi è fogliato da curve invarianti che
sono delle iperboli.
Caso2; Un sistema tale che λ1 = eiα e λ2 = e−iα lo spazio delle fasi è fogliato da curve invariantiche sono circonferenze, aventi come centro l’origine.
1.3.2 Genericità ed Eccezionalità
Tenendo conto che gli autovalori sono funzioni continue degli elementi di matrice, allora possiamopostulare la seguente affermazione:
Teorema:Sia A una matrice di tipo ellittico, se vario di poco gli elementi di matrice i suoi auto-valori sono ancora sul cerchio unitario a distanza fissa dall’asse reale. In altri termini si ha cheper un punto corrispondente ad una matrice ellittica esiste un intorno U tale che ∀x ∈ U è ancoraun punto di una matrice ellittica.
Definizione: Se una matrice soddisfa la condizione sopra esposta si dice che il sistema che essarappresenta è strutturalmente stabile.
Lo stesso vale in generale anche per i punti iperbolici, ma non vale per quelli parabolici. Variandogli elementi di una matrice parabolica, si ricadrà in un sistema ellittico o iperbolico. Sono ecce-zionali quelle variazioni di parametri che tengono inalterato il carattere parabolico della matrice.Lo studio delle variazioni del tipo qualitativo in funzione ai parametri è l’argomento chiave dellamoderna teoria delle biforcazioni
1.3.3 Il caso di trasformazioni non lineari
Come detto nella sezione precedente, quando abbiamo a che fare con un sistema non lineare comelo standard map perturbato, tentiamo di linearizzarlo usando uno sviluppo in serie sui puntifissi. Il punto fisso su cui viene valutata la matrice Jacobiana prende poi la qualità di iperbolico,
4Teorema di Liouville (?)
5
parabolico o ellittico a seconda della natura della stessa matrice Jacobiana.5Studiamo ora i punti fissi del mapping perturbato. Precedentemete avevamo trovato che i puntifissi dello standard map perturbato era in P1(0, 0) e in P2( 1
2 , 0). Calcoliamoci la Jacobiana neidue punti: matrice:
JacAP1 =(
1 + 2πε 12πε 1
)(23)
Mentre in P2
JacAP2 =(
1− 2πε 1−2πε 1
)(24)
risolvendo l’equazione agli autovalori
λ2 − TrAλ+ detA = 0 (25)
ottengo che P1(0, 0) è sempre di tipo iperbolico. Mentre P2(1/2, 0) è ellittico solo per ε < 2π . Per
valori maggiori diventa parabolico e per valori grandi di ε diventa iperbolico. Quindi in (1/2, 0)si ha una biforcazione.
1.4 Twist Map - Ultimo teorema di Poincaré
Introduciamo un esempio molto interessante per verificare quanto detto fin’ora. L’esempio ciporterà automaticamente ad enunciare e dimostrare l’ultimo teorema di Poincaré. L’esempio chesi vuole portare è quello del Twist Map. Iniziamo a scriverne la trasformazione:
Ψ0 =
r′ = r
φ′ = φ+ 2πα(r) (mod2π) :dα(r)dr
6= 0(26)
Ogni circonferenza di raggio r ha la sua rotazione, che è in funzione del raggio. Per quandoriguarda la dinamica si ha una condizione del tutto simile allo standard map imperturbato, cioèuno spazio delle fasi fogliato da tori T 1 invarianti.Si hanno quindi gli stessi casi che si avevanoper lo standard map cioè che se la traslazione avviene con un coefficiente irrazionali allora leorbite andranno a riempire in modo denso tutte le superfici invarianti (cioè i tori T 1), mentre perrotazioni con coefficiente razionale si hanno delle orbite periodiche. Ora prendiamo un segmentoγ trasversale alla circonferenza questo viene trasformato in una linea curva a causa del fatto cheogni punto ha una rotazione diversa che è in funzione di r 6 però il punto di intersezione della lineacurva con il toro T 1 invariante è sempre lo stesso! Ora invece consideriamo il mapping perturbato.
Ψε =
r′ = r + εf(r, φ)
φ′ = φ+ 2πα(r) + εg(r, φ) (mod2π) :dα(r)dr
6= 0(27)
Il segmento viene perturbato di una quantità piccola a piacere avendo come risultato una suatralazione verso l’alto. Il punto di intersezione della linea curva con segmento γ c’è sempre perònon è più il punto P ! Questo nuovo punto lo chiamiamo Q′. Questo nuovo punto non è un puntofisso del mapping e l’unica cosa che possiamo dire è che appartiene a γ quindi il risultato ottenuto ènotevole, abbiamo visto che esiste sempre un punto Q che nel passaggio da mapping imperturbato
5Un’osservazione interessante riguarda il fatto che se abbiamo un continuo di punti fissi allora questi punti sononecessariamente parabolici. Questo è abbastanza intuito da osservare perchè se fossero iperbolici o ellittici una lorovariazione mi dovrebbe portare ancora su un punto fisso, invece proprio perchè i punti sono o su curve invariantiche sono ellissi o iperboli tendono a ruotare. L’unico modo per far si che una loro variazione porti ancora ad unpunto fisso è che siano di tipo parabolico.
6in particolare per x < r∗ i punti ritarderanno, mentre per x > r∗ i punti hanno maggior angolo di rotazionequindi anticiperanno sulla posizione.
6
a mapping perturbato viene trasformato radialmente nel punto Q′! Questo è un risultato moltoimportante infatti ci spinge subito a definire il luogo dei punti Q che hanno qusta particolareproprietà. Chiamiamo il luogo dei punti Σ (che sarà rappresentata da una curva chiusa) e oraconsideriamo le intersezioni tra le curve Σ e Ψε(Σ) La trasformazione è di tipo simplettico linearequindi ha determinate uno e mantiene le aree, quindi la Ψε(Σ) rispetto alla stessa Σ non potràessere tutta esterna o tutta interna, bensì ci dovranno essere dei punti di intersezione in qunto learee devono uguagliarsi. Quindi ecco che si riscontra ancora lo stesso fenomeno che avevamo vistoper lo standard map. Un logo denso di punti fissi, tramite una piccola perturbazione, si trasformain un luogo di punti fissi discreti7
1.5 Teorema 1 di Poincaré
Per il Twist Map possiamo dire: Nel passaggio dal T0 al Tε in generale avviene che tutti i toriinvarianti di T0 corrispondenti a traslazioni razionali α ∈ Q si ”sgretolano” rimangono solo unnumero finito di orbite (2k) periodiche di periodo q sotto Φε cui corrispondono punti fissi delmapping Ψε = T qε alternativamente ellittici ed iperbolici.
1.5.1 Numeri irrazionali di classe (K,σ)
Un numero irrazionale α viene detto di classe (K,σ) se è soddisfatta la seguente relazione:∣∣∣∣α− p
q
∣∣∣∣ > K2
|q|2+σK > 0, σ > 0 (28)
Teorema:Per ogni numero irrazionale α esiste un’approssimazione razionale esatto quanto sivuole, il cui errore è inferiore all’inverso del quadrato del denominatore∣∣∣∣α− p
q
∣∣∣∣ < 1q2
(29)
1.6 Teorema 2 (KAM)
Sempre per il Twist Map possiamo ora enunciare il teorema di Kolmogorov-Arnol’d-Moser:Consideriamo un toro invariante per il twist map imperturbato, con corrispondende costante dirotazione α irrazionale. Se α è di classe (K,σ) allora per ε abbastanza piccolo esiste un toroinvariante anche per il mapping perturbato che è una piccola deformazione di quello imperturbato.Infine la misura di Lebesgue dell’insieme dei tori invarianti del sistema perturbato tende allamisura piena dell’insisme nel limite che ε→ 0
1.7 Teorema della varietà stabile
Per spiegare i meccanismi della dinamica delle orbite intorno ai punti fissi iperbolici dobbiamofar riferimento al concetto di punto omoclino e a come Poincaré aveva descritto lo spazio linearestabile e instabile relativo ad un certo punto fisso (iperbolico).Consideriamo un mapping non lineare così costituito:
φ : M →M (30)
supponiamo che sia invertibile su tutto M e che x sia tale che: φ(x) = x. Sia A la matriceJacobiana della trasformazione φ allora se la matrice A è iperbolica (cioè ha autovalori di modulodiverso da uno) posso definire la varietà stabile e la varietà instabile legate al punto fisso x. Orasviluppiamo il discorso facendo riferimento ad un esempio semplice. Nel caso della trasformazionedel piano (lo standard map per esempio) la matrice era iperbolica di tipo simplettico, in più lacondizione λ1λ2 = 1 faceva in modo tale che un autovalore fosse maggiore di 1 mentre il secondo
7in particolare un numero pari di punti fissi, alternativamente iperbolici ed ellittici
7
fosse minore di 1. Quindi è ovvio che se prendo l’autospazio di A relativo al primo autovaloreavrò uno spazio ”contratto” 8 mentre l’autospazio di A relativo al secondo autovalore sarà dilatato.Dunque le nostre varità lineari stabili ed instabili sono così definite:
E+ = {x ∈ Rn : Anx→ 0 per n→ +∞} (31)
mentre la varità dilatante è definita come:
E− = {x ∈ Rn : Anx→ 0 per n→ −∞} (32)
Il punto geniale della trattazione viene ora. La linearizzazione del problema ci ha portato a definirelo spazio lineare stabile ed instabile associato al punto fisso iperbolico. Ora però questa proprietàdei punti di E+ ed E− bisognerebbe riuscire a trasportarla anche alla trasformazione non lineareφ. Questo è quello che fa Poincaré, prima in senso locale ad un intorno U(x) ⊂ M , poi vedremoche può essere esteso in forma globale a tutta la varietà M .
Teorema: Sia x ∈ M un punto fisso di φ iperbolico. Sia E+ la relativa varietà lineare sta-bile associata al punto fisso. Consideriamo U(x) ⊂ M sufficientemente piccolo e consideriamol’insieme W+ identificato da:
W+ = {x ∈ U : φn(x) ∈ U ∀n ≥ 0 , φn(x)→ x per n→ +∞} (33)
Allora si ha che W+ è una varietà che ha la stessa dimensione di E+ passa per il punto x ed ivitangente allo spazio lineare stabile E+. In più la varietà W+ non può ripiegarsi sopra E+ essendoquindi rappresentabile attaverso un grafico
Quindi abbiamo visto come una proprietà vista per il caso lineare è stata trasportata al casonon lineare. Ora il punto consiste nel rendere i concetti fin’ora esposti generali e quindi estenderele proprietà a tutta la varietà M .
Passaggio da locale a globale: Se noi definiamo:
W+ = {x ∈M : φn(x)→ x per n→ +∞} (34)
e per quando riguarda la varietà instabile:
W− = {x ∈M : φn(x)→ x per n→ −∞} (35)
La varietà stabile è definita come l’insieme di punti per cui per t → +∞ tendono a cadere sulpunto fisso x allo stesso modo per la varietà instabile con t→ −∞. Quindi sono insiemi invarinati,cioè variatà comporte da orbite. Se un punto appartiene alla varietà allora gli appartiene anchetutta la sua orbita. Questo quindi fornisce un metodo costruttivo della variatà stabile ed instabileda poter usare col calcolatore. Posso costruire la curva della varietà stabile semplicemente agendoin questo modo:
γ2 = φ−1γ1 = φ−2γ0 (36)
Aumentando il numero delle iterazioni vado ad aggiungere sempre più intervalli, ottendo infinela varità. La varietà stabile non può intersecare se stessa come quella instabile, ma nulla vieta ilfatto che possano avere dei punti in comune. Il punto di incontro tra varietà stabile ed instabileè proprio l’origine dei moti caotici, ed è il famoso punto omoclino di Poincaré che denotiamocon ν. Il caso più semplice è considerare in cui varietà instabile e stabile siano coincidenti. Peròil punto più interessante e quello che capita maggiormente è quando W+ 6= W− e si ha cheν ∈ W+eν ∈ W− quindi ν ∈ W+
⋂W−. Il punto è che quando varietà stabile ed instabile si
8Infatti lo spazio lineare stabile è detto anche ”contraente”, mentre quello relativo alla varità lineare instabile èdetto anche ”dilatante”
8
incontrano danno luogo nel punto di intersezione a contuinui intrecciamenti dovuti al fatto cheW+ e W− sono due varietà invarianti cioè se contengono un punto, allora ne contengono tuttala sua orbita che quello genera. Se ν ∈ W+
⋂W− allora ν1 = φ(ν) ∈ W+
⋂W− e così via per
tutti i suoi iterati. Questo spiega quindi lo spezzamento della sepratrice del pendolo. Si viene acreare un punto omoclino in cui la verietà stabile e instabile si intersecano trasversalmente, quindila derivata a seconda delle direzioni è diversa, quindi si ha lo splitting della separatrice.
1.8 Teorema di Arnol’d-Jost
Sia H(q, p) un’Hamiltoniana che ammette n integrali primi in involuzione {φi}ni=0. Sia c ∈ Rntale che la varietà:
M =n⋃i=0
{q, p : φi(q, p) = ci} (37)
contenga una componente connessa e compatta Mc allora in un intorno U(Mc) esiste un sistemadi variabili canoniche di angolo-azione (θ, I) tali che: Tn × G → U(Mc) (dove G è un aperto diRn) tale che l’Hamiltoniana trasformata dipenda solo dalle Ii e quindi le equazioni di Hamiltonhanno soluzione:
θi(t) = θi,0 + tωi(I1,0, . . . , In,0) (38)
dove abbiamo che Ii,0, θi,0 sono i nostri dati iniziali e ωi = ∂H∂Ii
.Il teorema di Arnol′d − Jost formalizzato qui sopra è un modo per dire in maniera più generalequello che nelle sezioni precedenti si era osservato nel caso dello standard map e del twist mapimperturbati. Infatti afferma che in presenza di n integrali primi le orbite nello spazio delle fasigiacciono tutte su tori invarinti parametrizzati dalle sole azioni Ii, come visto precedentemente ilparametro fondamentale che ci indica se avremo delle orbite dense o periodiche è il parametro α.Se è irrazionali avremo delle orbite dense sul toro, se è razionale avremo delle orbite periodiche.Per generalizzare questo fatto in modo topologico ed ampliare la definizione ad un caso generico didimensione n possiamo introdurre il concetto del modulo di risonanza. Il modulo di risonanzaè così definito:
Mω = {k ∈ Zn : (k, ω) = 0} (39)
La dim(Mω) è detta molteplicità di risonanza. I casi estremi sono dim(Mω) =0 (non risonante) edim(Mω) = n− 1 (completamente risonante). Quindi l’orbita è densa sul toro: Tn−dim(Mω). Nelcaso non risonante l’orbita è densa su tutto il toro Tn.
1.9 Teorema di Poincaré
Iniziamo ora a trattare quello che nelle scorse sezioni è stato solo enunciato, cioè ilTeorema di Poincaré.Vedremo bene la formalizzazione topologica e tutte le questioni annesse. Quello che Poincaré sipropose di fare era studiare e dare dei risultati generali sulla dinamica di un sistema descritto daun’Hamiltoniana di questo tipo:
• H(q, p, ε) è funzione analitica di tutti i suoi argomenti, in particolare sviluppabile in serie dipotenze di ε in un determinato intorno;
9
Ora potrei andare a verificare se il sitema analizzato ammette n integrali primi e quindi siaintegrabile. Per fare questo cerchiamo di costruire degli integrali nella forma:
per essere un integrale del modo bisogna che sia soddisfatta la relazione:
{φ,H} = 0 (42)
Sostuituendo lo sviluppo di H e dell’integrale primo otteniamo che:
{φ0, H0}{φ1, H0} = −{φ0, H1}. . .
{φr, H0} = −r∑j=1
{φr−j , Hj}
La condizione di non degenerazione implica direttamente che dalla prima equzione Φ0 non puòdipendere da variabili di tipo angolo, ma solo da quelle di azione infatti:
Φ0 =∑k∈Zn
φk(p)ei〈k|q〉 =⇒ {Φ0, H0} =∑k∈Zn
i
⟨k
∣∣∣∣∂H0(p)∂p
⟩φk(p)ei〈k|q〉 (43)
per soddisfare la condizione iniziale di annullamento si deve avere allora che:⟨k∣∣∣∂H0(p)
∂p
⟩= 0
oppure φk(p) = 0. Differenziando la prima ottengo che:
∑j
kj∂2H0(p)∂pj∂pl
(44)
ma la condizione di non degenerazione mi impedisce di annullare il determinante dell’Hessianoquindi ciò che rimane è che le φk(p) = 0 ∀k 6= 0Come sono fatte le φk(p) ? Sviluppiamo in serie di Fourier e ugualiamo i termini:
φ1 =∑k∈Zn
φkei〈k|q〉 φk(p) = φ∗−k
H0 =∑k∈Zn
hkei〈k|q〉 hk(p) = h∗−k
(45)
ponendo ωj = ∂H0∂pj
ed uguagliando i termini della serie otteniamo che:
ck(p) := i
⟨k
∣∣∣∣∂H0
∂p
⟩hk = i 〈k |ω(p)〉 φk (46)
da cui:
φk(p) = −i ck(p)〈k |ω(p)〉
(47)
Il procedimento può essere reiterato per tutte le equazioni ottenendo tutti i coefficienti. Si pongonotuttavia due problemi fondamentali:
• ConsistenzaLa validità della (39) nel caso in cui k = 0 impone che sia c0(p) = 0. Questacondizione è soddisfatta per Φ1 non è detto in generale;
10
• Piccoli denominatori : la soluzione trovata non vale sulla varietà F = {p ∈ G : 〈k |ω(p)〉 = 0}densa in Rn e quindi in G. In più anche se si riuscisse a risolvere l’equazione (38) ottenendodelle soluzioni non appartenenti ad F è possibile che i denominatori siano così piccoli da fardivergere ugualmente l’integrale.
Il teorema di Poincaré è quindi così enunciato:
Teorema: Sia H(1, p, ε) una Hamiltoniana che soddisfi le ipotesi:
• Non degenerazione: det dell’hessiano sia diverso da zero;
• Genericità: i coefficienti hk(p) dello sviluppo in serie di Fourier di H1(q, p) sono tutti NONnulli sulla varietà F
Allora non esistono integrali primi analitici Φ(q, p, ε) indipendenti da H(q, p, ε).
2 Un semplice algoritmo
2.1 Il problema di Keplero
Il sistema di equazioni differenziali usato è quello sviluppato da P.H.Cowell e A.C.D Crommelinnel 1911 di cui si riporta nota in un articolo di F.L.Whipple presente a pag 388 Vol 111 di Apjnell’anno 1950 intitolato : ”The Comet Model: The acceleration of Comet Enke”. Il sistema diequazioni differenziali proposto da Cowell e Crommelin è il seguente:
..−→ri= −µ−→ri||ri||3
+G ·N∑
j=0 ; j 6=i
mj
[ −→r0j
||r0j ||3−−→rj||rj ||3
](48)
Dove per esteso abbiamo che:
−→r 0j = −→r j −−→r 0 (49)
−→r 0 è la posizione dell’oggetto su cui si stanno calcolando le perturbazioni generate degli altri Ncorpi di massa mj situati nelle posizioni −→rj .Facendo riferimento al fatto che le scrittura r ed r0j senza la freccia di vettore simboleggiano ilmodulo dei vettori, possiamo separare in componenti l’equazione (1) ottenendo:
..xi= −µ
xi||ri||3
+G ·N∑
j=0 ; j 6=i
mj
[xj − xi||r0j ||3
− xj||rj ||3
](50)
..yi= −µ
yi||ri||3
+G ·N∑
j=0 ; j 6=i
mj
[yj − yi||r0j ||3
− yj||rj ||3
](51)
..zi= −µ
zi||ri||3
+G ·N∑
j=0 ; j 6=i
mj
[zj − zi||r0j ||3
− zj||rj ||3
](52)
L’integrazione numerica che viene proposta in questa relazione, avviene proprio applicando ilmetodo di (RKF45) alle equzioni (2) e (3).
11
2.2 I metodi di integrazione
Sono stati usati due metodi di integrazione differenti. Runge-Kutta (RK) di ordine 4 e Runge-Kutta-Fehlberg (RKF45) di ordine 5. Il metodo RKF45 è così enunciato:9
k1 = f(tj , yj)
k2 = f(tj +14h, yj +
14k1)
k3 = f(tj +38h, yj +
332k1 +
932k2)
k4 = f(tj +1213h, yj +
19322197
k1 −72002197
k2 +72962197
k3)
k5 = f(tj + h, yj +439216
k1− 8k2 +3680513
k3 −8454104
k4)
k6 = f(tj +h
2, yj −
827k1 + 2k2 −
35442565
k3 +18594104
k4 −1140k5) (53)
Il punto yj+1 è definito nel seguente modo:
yj+1 = yj +[
16135
k1 +665612825
k3 +2856156430
k4 −950k5 +
255k6
]· h (54)
2.2.1 Verifica dell’ordine del metodo
Per verificare che il metodo usato sia effettivamente di ordine 5, si è risolto prima il seguenteproblema di Couchy: {
f ′(x) = f(x)f(0) = 1 (55)
Questo problema di Couchy ha la banale soluzione:
f(x) = ex (56)
che è stata usata per verificare l’affidabilità del metodo valutandone l’errore in funzione del passodi integrazione e quindi ricavando in modo diretto l’ordine dell’algoritmo di integrazione.
9Questo metodo è stato trovato sul seguente testo di analisi numerica: NUMERICAL METHODS forMathematics, Science and Engineering, 2nd Ed, 1992 Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A
12
In verde è plottato il metodo standard RK e come si può notare ha ordine 4. Mentre in rossoè stato plottato il risultato ottenuto con il metodo RKF45. Come si può notare qeust’ultimo haordine di integrazione pari a 5.
2.3 L’algoritmo
Inizialmente è stato necessario riproporre in C++ l’equazione differenziale del nostro sistema per-turbato. L’equazione differenziale è stata scritta nel seguente modo:
Altro punto fondamentale del programma consiste nell’implementazione dell’algoritmo RKF45.Qui di seguito ne riportiamo il codice:
void rungekuttafehlberg (double& t, double h, int N, double *x,double *xt,double *xg,double*masse,int p, double* (funzione)(double *,double *,double *,double *,int,double)){double *k1 = (*funzione)(x,xt,xg,masse,p,t);double *temp = new double[N];for(int i = 0; i < N; i++) temp[i] = x[i] +0.25*h*k1[i];double *k2 = (*funzione)(temp,xt,xg,masse,p,t+0.25*h);for(int i = 0; i < N; i++) temp[i] = x[i]+(3.0*h/32.0)*k1[i]+(9*h/32)*k2[i];double *k3 = (*funzione)(temp,xt,xg,masse,p,t+(3.0*h)/8.0);for(int i = 0; i < N; i++) temp[i] = x[i]+ (1932.0*h/2197.0)*k1[i] -(7200.0*h/2197.0)*k2[i]+ (7296.0*h/2197.0)*k3[i];double *k4 = (*funzione)(temp,xt,xg,masse,p,t+(12*h)/13);for(int i = 0; i < N; i++) temp[i] = x[i] + (439.0*h/216.0)*k1[i] - (8.0*h*k2[i])+(3680.0*h/513.0)*k3[i] - (845.0*h/4104.0)*k4[i];
13
double *k5 = (*funzione)(temp,xt,xg,masse,p,t+(4*h)/5);for(int i = 0; i < N; i++) temp[i] = x[i]- (8.0*h/27.0)*k1[i] + (2.0*h*k2[i]) -(3544.0*h/2565.0)*k3[i]+ (1859.0*h/4104.0)*k4[i] - (11.0*h/40.0)*k5[i] ;double *k6 = (*funzione)(temp,xt,xg,masse,p,t+h);for(int i = 0; i < N; i++)x[i]+=(16.0*k1[i]/135.0 + 6656.0*k3[i]/12825.0 + 28561.0*k4[i]/56430.0 - 9*k5[i]/50.0+2.0*k6[i]/55)*h;t+=h;delete k1;delete k2;delete k3;delete k4;delete k5;delete k6;delete temp;}In questo modo si riconduce il problema non lineare di partenza a tanti problemi lineari. Il numeroN che viene dato in input permette di applicare i vari passaggi N volte, risolvendo a cascata anchesistemi complicati di equazioni differenziali. Nel caso degli N corpi abbiamo un sistema di 3 equa-zioni differenziali (una per ogni coordinata x,y,z) del secondo ordine non omogenee, che vengonorisolte tutte tramite l’ausilio di questa unica funzione. Ovviamente è necessario che arrivi ancheun vettore di dati iniziali di lunghezzaN∗g doveN è il numero di eq. differenziali e g è il loro grado.
2.4 Troncamento e Tolleranza
In questo tipo di algoritmi ad un passo (one-step) bisogna prestare molta attenzione a quantoincide l’accumulazione dell’errore di troncamento. Il buon senso ci porterebbe a dire che facendotendere h → 0 (o comunque a valori molto piccoli in quanto h 6= 0) la precisione con cui andia-mo a stimare la nostra soluzione del problema dovrebbe essere sempre migliore. la precisione ciattendiamo che segua l’andamento descritto nella figura precedente. Ora si mostrerà che questadeduzione è errata in quanto per piccoli valori di h aumentano moltissimo anche il no di iterazioni,ad ogni iterazione è associato un errore che si somma passaggio dopo passaggio. Il grafico prece-dente è stato costruito con opportuni valori di h scelti in modo tale che l’errore dipendesse solodal grado di accuratezza quindi questo effetto non si nota. Ora verrà presentato il comportamentodell’errore con h ”critici” in cui questo comportamento sarà ben evidente.
14
Dal grafico si può notare che diminuendo il valore di h oltre ad un certo limite critico si perdeprecisione a causa dell’accumulazione dell’errore di troncamento dovuto alle troppe iterazioni.In punto cruciale sta quindi scegliere la giusta h in modo tale che l’errore sia minimo e che cheaiuti a compensare l’accumulazione dell’errore di trocamento su grandi scale di integrazione
A supporto di questo problema intervengono i metodi nominati come ”Predictor-Corrector”.Si usa un metodo di ordine p come predizione e con un metodo di ordine p − 1 si effettuanole stime degli errori. Il tipico esempio è la coppia Predictor Corrector fatta con gli algoritmidi Adams-Bashforth e Adams-Multon. In questa breve relazione però si è scelto di usare comePredictor il metodo RKF45 mentre come corrector il semplice algoritmo di RK. Il sistema consistenel modificare il passo di integrazione ad ogni iterazione in modo tale che l’errore si mantengacostante di una certa tolleranza scelta inizialmente. La modifica del passo di integrazione h avvienetramite il calcolo della ”sezione di h” che chiameremo s. Ne risulta che il nuovo passo sarà:
hnew = s · hold (57)
Il calcolo esplicito di s, viene sempre presentato da John H. Mathews & Kurtis K. Fink nel lorolibro : Numerical Methods Using Matlab, 4th Edition, 2004 :
s =(
tol · hold2 |zk+1 − yk+1|
) 14
=⇒ hnew = s · hold (58)
dove tol è il massimo errore che noi tolleriamo in tutto il calcolo. Ovviamente questo metodo, senon adeguatamente controllato potrebbe rendere particolarmente instabile l’algoritmo. Proprioper questo è necessario definire un hmax e un hmin entro il quale non si puù uscire. Se calcolo uns molto minore di 1 avrò un h troppo piccolo e il calcolo oltre ad essere pesantemente influenzatodall’accumulazione dell’errore di troncamento, rischia di non convergere.
2.5 Il sistema solare
All’inizio per ogni pianeta viene definito il suo vettore di posizion x, posizione y, massa, velocitàx e velocità y nel seguente modo:
while (t <= tmax) {//Mercuriorungekuttafehlberg(t,h,4,posmer,pianetix,pianetiy,masse,0,asteroide);mercuriox[i]=posmer[0];mercurioy[i]=posmer[1];pianetix[0]=posmer[0];pianetiy[0]=posmer[1];t-=h;//Venererungekuttafehlberg(t,h,4,posven,pianetix,pianetiy,masse,1,asteroide);venerex[i]=posven[0];venerey[i]=posven[1];pianetix[1]=posven[0];pianetiy[1]=posven[1];
15
ecc...}
Dove i vettori pianetix e pianetiy contengono tutte le informazioni circa le posizioni aggior-nate degli 8 pianeti maggiori. Dato che il programma è pensato nel particolare per calcolare leperturbazioni dei pianeti maggiori sugli asteroidi, nell’algoritmo viene lasciato libero un vettore”asteroide” che puù essere inizializzato a piacere coi parametri iniziali del’asteroide di cui si vuolecalcolare la dinamica. L’agoritmo è pensato in modo tale che l’utente possa aggiungere anche Nasteroidi, ma questi non rientrino nel gioco delle perturbazioni con i pianeti maggiori del sistemasolare. Quindi i pianeti risentono delle mutue perturbazioni, gli asteroidi risentono delle pertur-bazioni degli 8 pianeti maggiori, ma loro non perturbano il moto dei pianeti. L’orbita viene poiscritta su un file .dat nel formato opportuno per essere letto da gnuplot e quindi essere in gradodi disegnare anche le orbite.
2.6 Effetti particolari sull’orbita di mercurio
Un effetto decisamente particolare che si è osservato rigurda l’orbita di mercurio. Qui di seguitoinseriamo due rappresentazioni della sua orbita calcolata la prima con il metodo di runge e kuttadi ordine 4 la seconda con il metodo RKF45.Runge-Kutta (ordine4):
Runge-Kutta-Fehlberg (ordine5):
16
Il fatto particolare consiste nel SOLO allargamento dell’orbita di mercurio. Sia nella prima im-magine che nella seconda le orbite di Venere e Terra sono praticamente comparabili. Passandodall’ordine 4 all’ordine 5 invece l’orbita di mercurio subisce un netto cambiamento. Nel caso diordine 4 si nota un evidende allargamento dell’orbita dato presumibilmente dall’accumulazionedell’errore di troncamento. Questo effetto di allargamento sparisce passando all’ordine 5.
2.7 Modifiche effettuate
Quello presentato che è stato presentato è un abbozzo abbastanza gergale di un problema ben piùcomplicato. Per semplificare i conti, per prima cosa, non si usano più km e s come unità di misuradi lunghezza e di tempo ma si usano unità di misura molto più comode che sono le UA e i day. Inquesto modo la costante di gravità universale va leggermente modificata e si adotta il parametrodi Gauss µ con relativi riscali ad Hoc. Nella presentazione precedente venivano fatti partire daposizioni ”belle” , invece il problema finale conivolge direttamente l’epoca iniziale da cui si vuolepartire. Un programma a parte estrae dal file DE406.dat del JPL le posizioni iniziali in formatocartesiano dei pianeti del sistema solare, vengono scritte su un file che viene successivamenteletto dal programma, che imposta le condizioni iniziali del problema agli N-Body. Noi non siamointeressato solo alla dinamica generale, ma anche a quello che noi osserviamo per esempio dallaterra e quindi dare in output delle cordinate equatoriali di posizione in funzione al nodo ascendentee all’equatore celeste! I due angoli di posizione vengo calcolati effettuando banali conti matricialicontenuti nel pacchetto lapack. Viene effettuata la correzione topocentrica e viene calcolata laParallasse sfruttando un pochettino di relatività ristretta. Il programma totale non è liberamenteaccessibile al pubblico.
3 Lacune di Kirkwood
Si è preso in considerazione l’asteroide Musa 600 si è integrato il moto considerando il problemaa 4 corpi Terra Sole Asteroide Giove su un periodo di 10000 anni. Il più è stato plottato anche unasteroide FINTO situato nella fascia di risonanza 3:1 con giove . I dati delle variazioni dei loroelementi orbitali kepleriani su identiche scale temporlai sono qui ora messi a confronto:
Variazione semiasse maggiore
17
Variazione eccentricità
18
Variazione Argomento del Perielio
Variazione Longitudine del Periastro
Come si può ben vedere dai grafici nella prima colonna di sinistra sono plottate le variazionidegli elementi orbitali dell’asteroide 600 Musa che ha un semiasse maggiore di circa 2.6UA. Glielementi orbitali variano in modo previsto dalla teoria perturbativa di Lagrange. Il semiasse mag-giore oscilla intorno ad un valore fisso, l’eccentricità pure lei oscilla molto lentamente e argomentidel periastro elongitudine del nodo ascendente mostrano il tipico andamento di precessione linea-re. Queste particolarità non si notano più nella colonna di destra dove sono plottati gli elementiorbitali del nostro asteroide creato ad hoc e posizionato nella fascia di risonanza 3:1 con Giove.Semiasse maggiore ed eccentricità variano in modo molto brusco senza avere un path prevedibile.Questo porta l’asteroide su orbite differenti creando uno svuotamento della regione di risonanza3:1.
Dato che per la simulazione non si è usato un algoritmo simplettico, non si è potuti andare oltre ai10000 anni di integrazione. Qui di seguito presenteremo però delle mappe create da A. Morbidellinel 1996 in cui ha integrato per un intervallo di tempo di circa 1010anni
19
Figura 1: Distribuzione asteroidi nella fascia principale. Si notano le risonanze con Giove e contutti gli altri pianeti.
Figura 2: Andamento dei 6 elementi orbitali in funzione del tempo di un asteroide posizionatonella risonanza 3:1
3.1 Avvicinamenti radenti-Closest Approach
Per simulare il comportamento dell’avvicinamento radente di un asteroide abbiamo consideratol’asteroide: 4581 Asclepius. Nel sito del Minor Planet Center è segnalato un suo Closest Ap-proches avvenuto il 22 Marzo 1989. L’integrazione è stata fatta a cavallo dell’istante di minimadistanza. La finestra usata copre 20 anni prima e 20 anni dopo l’evento in questione, l’integrazionetotale quindi risulta essere di un tempo totale di 40 anni. Qui di seguito sono visibili i grafici dellevariazioni degli elementi orbitali nel tempo.
20
Semiasse maggiore ed eccentricità:
Inclinazione e Argomento del Perielio:
Longitudine del nodo e Anomalia Media:
Com’ era ben prevedibile, il passaggio radente dell’asteroide 4581 ha provocato un netto cam-biamento dei suoi parametri orbitali. Per effetto fionda infatti l’asteroide è stato lanciato lontanoquindi come è ben visibile dal primo grafico, il valore del semiasse maggiore è notevolmente au-mentato. Così come è aumentata in modo netto anche la sua eccentricità. Particolarmente stranoè risutato il grafico dell’inclinazione dell’asse. Dopo un primo hump aumenta vorticosamente perpoi stabilizzarsi su volori molto alti. Argomento del Perielio e Longitudine del nodo indicanoinizialmente la tipica precessione newtoniana con una discontinuità in prossimità del passaggioradente. Ovviamente l’orbita è cambiata in modo repentino e i parametri si sono ”rimescolati”di conseguenza. In funzione al cambio repentino l’argomento del perielio ha raggiunto una suaquota. Poi dopo che i parametri si sono stabilizzati, ha continuato a precedere col suo tipico motodi precessione lagrangiana. Stessa affermazione per la longitudine del nodo in cui alla tipica pre-cessione si somma, in prossimità del passaggio radente, un cambiamento netto per poi ritornarenuovamente con un moto di precessione simile a quello iniziale.
21
Aumento dell’incertezza nei passaggi radenti
La parte di critica dell’analisi dei Closest Approaches sta nel fatto che determinare come e in chemodo avverrà il cambiamento è praticamente impossibile. Le innumerevoli perturbazioni rendenoimpossibile un calcolo preciso, in più ora vedremo che anche solo effettuando una piccolissimavariazione di un elemento orbitale, le due orbite che si vengono a creare dopo passaggio radentesono totalmente differenti. Quindi si ha molta incertezza sull’orbita finale!Piccolissime variazioni dell’orbita (o per meglio dire: le imprecisioni che noi abbiamo quandocalcoliamo un orbit) anche se molto piccole nel momento in cui l’oggetto effettua un Closest Ap-proach dopo la sua orbita può cambiare in modo notevole ! Osserviamo per esempio sempre ilcaso dell’asteroide 4581: In questo primo caso abbiamo variato di un fattore di 10−8 il valore delsemiasse maggiore ottenendo il seguente risultato:
Semiasse maggiore ed eccentricità:
Inclinazione e anomalia media:
Argomento del perielio e longitudine del nodo:
Qesti primi dati sono di facile interpretazione infatti, il primo grafico identifica l’andamento delsemiasse maggiore. C’è un brusco cambiamento e l’orbita finale va posizionarsi in modo più in-terno, l’anomalia media decresce linearmente. L’eccentricità subisce una variazione netta dopo laquale si stabilizza su di un valore. Conseguentemente cambiano anche argomento del perielio e
22
longitudine del nodo.Ora vediamo cosa succede a cambiare di un fattore di 10−8 il valore dell’eccentricità:Semiasse maggiore ed eccentricità:
Inclinazione e anomalia media:
Argomento del perielio e longitudine del nodo:
Ora vediamo cosa succede a cambiare di un fattore di 10−3 il valore dell’inclinazione:
Semiasse maggiore ed eccentricità:
Inclinazione e anomalia media:
23
Argomento del perielio e longitudine del nodo:
Ora vediamo cosa succede a cambiare di un fattore di 10−6 il valore dell’anomalia media:
Semiasse maggiore ed eccentricità:
Inclinazione e anomalia media:
Argomento del perielio e longitudine del nodo:
24
Ora vediamo cosa succede a cambiare di un fattore di 10−6 il valore dell’argomento del perie-lio:
Semiasse maggiore ed eccentricità:
Inclinazione e anomalia media:
Argomento del perielio e longitudine del nodo:
25
Ora vediamo cosa succede a cambiare di un fattore di 10−6 il valore della longitudine del nodo:
Semiasse maggiore ed eccentricità:
Inclinazione e anomalia media:
Argomento del perielio e longitudine del nodo:
26
3.1.1 99942 Aphopis
Molto interessante è anche la seguente figura in cui viene plottata l’orbita integrata su 5000 anni.In particolare notare il plot verde e blu, indicano le due orbite di Aphopis in cui è stata variata lanona cifra decimale del semiasse maggiore.
4 Utilizzo del software presente in Osservatorio
4.1 Come generare l’effemeride di un asteroidePer prima cosa generare un’effemeride collegandosi con il sito web del Minor Planet Center e sele-zionando la pagina Minor Planet & comet Ephem. Compilare il modulo di richiesta dell’effemerideinserendo l’identificazione di un asteroide (il suo nome, oppure il suo numero di catalogo se è unasteroide numerato, oppure la sua designazione provvisoria se non è ancora numerato); in questocaso chiediamo un’effemeride dell’asteroide (2) Pallas: Lasciare pure invariate tutte le opzioni alloro valore preimpostato, quindi premere il tasto Get ephemerides/HTML page. In questo modovi comparirà sullo schermo una pagina con tutte le effemeridi. Ora cerchiamo ri riproporre lestesse effemeridi ma generate dal nostro software.Cerchiamo ora di riprodurre l’effemeride. Per prima cosa procuriamoci gli elementi orbitali del-l’asteroide Pallas dal file di elementi orbitali di tutti gli asteroidi distribuito da Ted Bowell:
$ cd ∼ /bowell$ zcat astorb.dat.gz | tail -n +2 | head -1 > 2.bow$ mv 2.bow ∼/cav/work
(queste istruzioni estraggono il secondo record del file compresso astorb.dat.gz, che contienegli elementi orbitali di Pallas, e lo scrivono nel file 2.bow, che viene poi spostato nel directory∼/cav/work; naturalmente ciò presuppone che si sappia in anticipo che Pallas è il secondo aste-
27
roide numerato e che quindi i suoi elementi orbitali sono riportati nel secondo record del file; senon si ha questa informazione, occorre decomprimere il file e cercare il record che serve). Quinditorniamo al directory di lavoro e trasformiamo il file 2.bow dal formato originario in un file nelformato interno usato dal software cav, che chiameremo 2.ele:
$ cd ∼/cav/work$ trbow.x
Input file name: 2.bowOutput file name: 2.ele
Controlliamo che gli elementi orbitali contenuti nel file 2.ele siano riferiti a una data di oscu-lazione non troppo distante dalla data a cui dobbiamo generare le effemeridi (11-31 maggio 2010)
$ tr-time.xSelect one of the following input formats:1) Julian Day2) Modified Julian Day3) Calendar date (year, month, day) at 0h4) Calendar date, hour, min, sec5) Year, day within year, sec within day6) GPS week, day, second7) GPS week, second within weekSelected format = 2Modified Julian Day = 55300JD = 2455300.50000000MJD = 55300.00000000Date = 2010/04/14 (Apr 14, 2010) 0:00: 0.000Day within year = 104Second within day = 0.000GPS week and day = 1579 3Second within GPS week = 259200.000Modified Julian Day =al 14 aprile 2010, quindi può andar bene. Procuriamoci ora gli elementi orbitali della Terra allastessa epoca (anche se, siccome genereremo l’effemeride usando il programma efem1.x, che usail metodo analitico basato sulle formule del problema di Keplero, non sarebbe necessario che ledate di osculazione dell’asteroide e della Terra siano le stesse). Estraiamo un insieme di condizioniorbitali per tutti i pianeti del Sistema Solare dalle effemeridi digitali del JPL:
28
$ pian403.xTJM = 55300Output file name: pian403.55300Copiamo il file di output pian403.55300 in un nuovo file che chiamiamo t55300.ele (elementi orbi-tali della Terra a MJD=55300) e modifichiamolo con un text editor lasciando solo i primi 6 recorddi intestazione e i 5 record finali che danno gli elementi orbitali del baricentro TerraLuna:
$ cp pian403.55300 t55300.ele$ vi (o emacs, o il text editor che preferite) t55300.ele$ cat t55300.eleEQUATORIALMEANCARTESIANET2000.00055300.0000000000 1.00000000000000Terra+Luna-0.917548418762705 -0.371345775006698 -0.1609844120001596.663604901633222E-03 -1.450178023262848E-02 -6.286910308744536E-031.000 3.289005603917905E+05-99.000 -99.000
Ora possiamo procedere a generare l’effemeride:
$ efem1.xFile degli elementi orbitali della Terra: t55300.eleFile degli elementi orbitali del pianeta: 2.eleFile di uscita (effemeridi) : 2.ephFile di uscita (osservazioni simulate) : (Return)Pianeta osservatore: Terra+LunaEpoca (ET): 14 APR 2010 ore 0.000 (MJD = 55300.00000)Pianeta osservato : 2 PallasEpoca (ET): 14 APR 2010 ore 0.000 (MJD = 55300.00000)File dei tempi di osservazione (<CR> per effemeride) = (Return)Scala di tempo dell’effemeride (UT/ET) = UTTempo iniziale (Data/tempo Giuliano) = dTempo iniziale (giorno, mese, anno, ora) = 11 5 2010 0.0Numero di punti = 20Distanza tra i punti (giorni, ore) = 1 0Correzione per aberrazione da applicare(Nessuna/Planetaria/Stellare/Entrambe) = pCorrezione baricentro Terra-Luna (S/N) = sCorrezione topocentrica (S/N) = nSistema di riferimento:(eClittico/eQuatoriale) = q(Medio/true-of-Date/true-of-Epoch) = m(1950/2000/...) = 2000
A questo punto controlliamo che l’effemeride prodotta (nel file 2.eph) sia uguale o quasi a quellagenerata dal servizio web del MPC.
29
4.2 Determinazione di un’orbita con Gauss2) Determinazione preliminare di un’orbita con il metodo di Gauss Per prima cosa procuriamoci ilfile di tutte le osservazioni astrometriche di un asteroide (ad esempio,27855 Giorgilli) dal servizioMPCOBS del sito web del Minor Planet Center. Salviamo il file (a cui è assegnato automatica-mente il nome obs.txt) nel directory di lavoro ∼/cav/work. Trasformiamo il file di osservazionidal formato .txt al formato letto dai programmi e chiamiamo il file risultante 27855.oss
$ tross.xFile di ingresso: obs.txtFile di uscita : 27855.oss
Con il programma che implementa il metodo di Gauss non lavoriamo sull’insieme completo delleosservazioni (file 27855.oss) ma su un sottoinsieme selezionato in un periodo piuttosto breve (12mesi) attorno alla data di scoperta (l’osservazione a cui ufficialmente è stata assegnata la scopertadell’asteroide è identificata da un asterisco nel file obs.txt; in questo caso la data di scoperta è il4 gennaio 1995). Quindi copiamo il file 27855.oss in un nuovo file 27855s.oss, da cui cancellere-mo tutte le osservazioni che non vogliamo usare in questa fase (ad esempio, quelle precedenti aMJD=49720, 3 gennaio 1995, e quelle seguenti a MJD=49768, 20 febbraio 1995):
$ cp 27855.oss 27855s.oss$ nano 27855s.ossDopo la cancellazione, il file 27855s.oss avrà quindi questo contenuto:
Procuriamoci un insieme di elementi orbitali della Terra (in realtà del baricentro del sistema Ter-raLuna) a un’epoca vicina alla data delle osservazioni che vogliamo usare, diciamo a MJD=49750.Abbiamo già visto come si fa:
$ pian403.xTJM = 49750Output file name: pian403.49750$ cp pian403.49750 t49750.ele$ nano t49750.ele(eliminazione di tutti i pianeti tranne Terra+Luna)$ cat t49750.eleEQUATORIALMEANCARTESIANET2000.000
30
49750.0000000000 1.00000000000000Terra+Luna-0.668715021967153 0.664071515612289 0.287916427289192-1.291728012987708E-02 -1.077052147518467E-02 -4.669736607180794E-031.000 3.289005603917905E+05-99.000 -99.000Siamo ora pronti per eseguire il programma che implementa il metodo di Gauss:
$ gauss.xFile elementi orbitali pianeta osservatore : t49750.eleFile di uscita (elementi orbitali) : 27855_1.eleFile di uscita (rapporto finale) : 27855_1.logPianeta osservatore: Terra+LunaEpoca (ET): 2 FEB 1995 ore 0.000 (MJD = 49750.00000)File delle osservazioni: 27855s.ossCorrezione topocentrica (S/N) = sCorrezione baricentro Terra-Luna (S/N) = sTempi delle osservazioni disponibili (UT):1) 49720.92671 3 GEN 1995 22.2410402) 49721.94789 4 GEN 1995 22.7493603) 49721.97657 4 GEN 1995 23.4376804) 49722.93426 5 GEN 1995 22.4222405) 49722.94328 5 GEN 1995 22.6387206) 49723.88015 6 GEN 1995 21.1236007) 49723.89137 6 GEN 1995 21.3928808) 49725.83480 8 GEN 1995 20.0352009) 49725.86906 8 GEN 1995 20.85744010) 49747.92668 30 GEN 1995 22.24032011) 49747.95794 30 GEN 1995 22.990560...........20) 49768.86581 20 FEB 1995 20.77944021) 49768.90255 20 FEB 1995 21.661200Dammi il numero progressivo dei 3 punti da usare = 1 4 9Radici possibili per la distanza Sole-pianeta all’osservazione 4:r(1) = 0.778966r(2) = 0.983696r(3) = 2.363590
RADICE NO. 3Orbita preliminare: a = 2.55754; ecc = 0.07588RADICE ACCETTATA: convergenza raggiunta dopo 10 iterazioniOrbita finale: a = 2.55764; ecc = 0.07587RMS = 40.111 secondi d’arco
Vuoi riprovare scegliendo altre osservazioni? (s/n) : nMIGLIORE SOLUZIONE TROVATA:
31
Tentativo 1, radice 3; RMS = 4.011E+01 secondi d’arco(naturalmente non è detto che tutto vada bene al primo tentativo come in questo esempio).
4.3 Correzione dell’orbita con un fit ai minimi quadrati basato sullaformula analitica dei 2 corpi
Correggiamo ora la soluzione orbitale trovata con il metodo di Gauss con il metodo dei minimiquadrati, usando ancora non tutte le osservazioni disponibili, ma solo quelle contenute nel file27855s.oss:
Si noti che l’orbita di partenza riproduce bene le osservazioni solo nel periodo 49720 ≤ MJD≤ 49726 (3-9 gennaio 1995), e naturalmente passa esattamente per le tre osservazione selezionateper il metodo di Gauss, mentre l’orbita finale ha errori dell’ordine del secondo d’arco su tuttoil periodo osservativo usato. Il programma fit1.x genera anche un file fit1.trc che contiene unatraccia di come cambiano gli elementi orbitali e i residui nel corso delle varie iterazioni. Si trattadi un file accessorio che serve soprattutto per scopi di debugging (e che per questo motivo vienesovrascritto ogni volta che si usa il programma).
4.4 Correzione dell’orbita con un fit ai minimi quadrati basato sullamodello degli N-corpi
Eseguiamo ora la correzione differenziale dell’orbita usando un integratore numerico, secondo unmodello di forza a N corpi che tiene conto delle perturbazioni gravitazionali dei pianeti principalidel Sistema Solare (programma fit2.x). Usiamo le stesse condizioni iniziali dell’esempio prece-dente (file p49750.ele e 27855_2t.ele) ed eseguiamo il fit su tutte le osservazioni disponibili (file27855.oss):
$ fit2.xFile elementi orbitali pianeti perturbatori: p49750.eleFile elementi orbitali pianeta osservato : 27855_2t.eleFile di uscita (elementi orbitali corretti): 27855_3.eleFile di uscita (rapporto finale) : 27855_3.logPianeti perturbatori:1) Mercurio2) Venere3) Terra4) Luna5) Marte6) Giove7) Saturno8) Urano9) NettunoPianeta osservatore = 3
33
Pianeta osservatore: TerraEpoca (ET): 2 FEB 1995 ore 0.000 (MJD = 49750.00000)Pianeta osservato : Soluzione migliore (tent. 1; rad. 3)Epoca (ET): 2 FEB 1995 ore 0.000 (MJD = 49750.00000)File delle osservazioni: 27855.ossCorrezione topocentrica (S/N) = sCorrezione baricentro Terra-Luna (S/N) = nIterazione no. 1MJD (UT ) Dalfa() Ddelta() STAT OUT49754.851520 1.289 -1.041 58749754.856610 1.191 0.224 58749754.870100 0.874 -0.309 587.....................48293.426350 38.363 -10.882 69148293.403530 38.186 -11.062 69148293.370210 38.376 -10.973 691RMS(sel) = 8.101E+01RMS(tot) = 8.101E+01Ind Cond = 5.021E+01.....................Iterazione no. 4MJD (UT ) Dalfa() Ddelta() STAT OUT49754.851520 1.307 -0.771 58749754.856610 1.209 0.494 58749754.870100 0.892 -0.039 587.....................48293.426350 -0.028 -0.016 69148293.403530 -0.206 -0.195 69148293.370210 -0.016 -0.103 691RMS(sel) = 9.196E-01RMS(tot) = 9.196E-01Come si vede, la deviazione standard dei residui passa da 81 secondi d’arco (prima iterazione)a 0.92 secondi d’arco (ultima iterazione), valore che si può ritenere rappresentativo degli erroriosservativi.
