L’Integrazione Numerica di Equazioni Differenziali: 50 anni fa l’ uomo ha messo piede sulla Luna

Nel giorno in cui ho iniziato a scrivere questo articolo ricorreva l’anniversario della prima esplorazione da parte dei cosmonauti americani Neil Armstrong, Michael Collins, Buzz Aldrin della nostra Luna. Anche se non ho una memoria diretta di questi eventi, le missioni delle progetto Apollo evocano in tutti noi una forte emozione rappresentando un momento unico ed epico nella storia della conquista dello spazio. Nessun altro uomo ha messo di nuovo piede sulla Luna dopo l’ultima missione Apollo 17 nel 1972, per cui i recenti annunci della NASA di nuove esplorazioni umane del nostro satellite rende l’anniversario ancora piu’ eccitante.

L’ immagine di copertina di questo articolo è stata creata da mio figlio per una sua ricerca scolastica sulle missioni Apollo ed è un collage d’immagini ottenute usando i programmi Google Earth e Sketchup. Nella figura si confrontano le dimensioni del razzo Saturno V con quelle della meravigliosa cattedrale di Lincoln in Gran Bretagna per dare un’idea dell’enorme grandezza del vettore spaziale. Il pensiero che l’uomo sia riuscito nel giro pochissimi anni in uno sforzo tecnologico e scientifico immenso a costruire questa cattedrale volante della tecnologia mi ha suscitato un senso di forte curiosità e non ho saputo trattenermi nel spiluccare tra la miriade di documenti disponibili sul sito della NASA sul programma Apollo.

In un recente articolo ho descritto il metodo di Eulero per la soluzione delle equazioni differenziali ordinarie usando come motivazione la versione romanzata nel film Hidden Figures ispirato dalla storia di Katherine Goble e delle sue due colleghe e dal loro contributo come “human computers” al programma spaziale della NASA. A tale scopo ho presentato come esempio un programma per il calcolo delle orbite dei pianeti del sistema solare con il metodo di Eulero. Le posizioni e le velocità iniziali dei pianeti alla data del 1 Luglio 1969 sono state ottenute usando il servizio Horizon online Ephemeris System curato dal Solar System Dynamics Group presso il Jet Propulsion Laboratory a Pasadena negli Stati Uniti d’America. Nell’usare questo servizio, mi sono incuriosito sul modo in cui questi dati sono stati calcolati e così ho iniziato a esplorare la ricchissima documentazione disponibile negli archivi online dell’ente spaziale americano. Ecco cosa ho scoperto.

Le posizioni più accurate dei pianeti e del nostro satellite a disposizione della NASA all’inizio dell’esplorazione spaziale erano quelle ottenute dai calcoli pionieristici di Wallace John Eckert (1902-1971), D. Brouwer, and G.M. Clemence. Nella loro pubblicazione del 1951, le interazioni dei primi cinque pianeti esterni (Plutone, Urano, Nettuno, Saturno e Giove) furono calcolate rigorosamente usando un sistema di 15 equazioni differenziali del secondo ordine (i calcoli furono effettuati rispetto al un sistema di riferimento eliocentrico). Le equazioni furono risolte numericamente usando intervalli di 40 giorni dal 1653 fino al 2060. Agli albori del calcolo scientifico essa rappresentò una notevole impresa computazionale tanto che Brouwer e Clemence nel loro saggio del 1961 [1] hanno evidenziato che “la simulazione richiese il calcolo di 5 milioni di moltiplicazioni e divisioni e 7 milioni di addizioni e sottrazioni di numeri con 14 cifre decimali (traduzione in italiano)“. Occorre ricordare che W.J. Eckert fu, tra le altre cose, un pioniere nell’uso dei primi calcolatori elettronici nel calcolo scientifico. Fu anche uno dei primi ad applicare nella soluzione di problemi scientifici il lettore/scrittori di schede perforate. Questo sistema di memorizzazione è stato impiegato per decine di anni come metodo d’inserimento dati e programmi nei primi calcolatori elettronici. Un completo resoconto del contributo di Eckert all’uso del calcolo scientifico può essere letto nel sito dell’ Università della Columbia in USA dedicato alla storia del calcolo elettronico [2].

La NASA e l’ente spaziale sovietico iniziarono a raffinare le coordinate dei pianeti e, in particolare, quelle del sistema Terra-Luna fin dagli inizi degli anni ’60 in concomitanza con le prime esplorazione spaziali del nostro satellite e dei pianeti limitrofi (nel sito web curato da S. Gerovitch, si può trovare un interessante e dettagliato resoconto storico dello sviluppo del calcolo elettronica nell’unione sovietica [3]). La disponibilità di dati telemetrici e lo sviluppo dei calcolatori elettronici ha permesso di ottenere traiettorie più’ accurate necessarie per compiere le prime imprese spaziali.

In particolare, con l’acquisizione del mainframe IBM 7094, il primo super computer commerciale fornito di compilatori per poter essere programmato in FORTRAN (versione II), il primo linguaggio ad alto livello, rese possibile agli ingegneri della NASA di effettuare i calcoli necessari per le missioni Mercury, Gemini e le prime Apollo. La prima versione (la versione 3) delle Development Ephemeris (DE3) fu calcolata usando l’IBM7090/94 e distribuita in una serie di nastri magnetici contenenti le posizioni dei pianeti e del baricentro Luna-Terra calcolate dal 1950 fino all’anno 2000 [4]. In varie relazioni tecniche disponibili online ho trovato delle informazioni dettagliate sui metodi usati per effettuare questi calcoli. L’approccio numerico era basato sull’uso di metodi perturbativi per correggere le orbite kepleriane dei pianeti. La perturbazione è stata calcolata con vari metodi numerici ottimizzati per usare la meglio la limitata memoria del IBM 7094. In particolare, i calcoli usati per le effemeridi DE3 sono basati sull’integrazione delle equazioni usando lo schema, del predictor-corrector dell’ottavo ordine di Gauss-Jackson [ref], originariamente usato da Cowell and Crommelin per lo studio della cometa di Halley ([5], vedi oltre). Questo è diventato uno dei metodi più usati in problemi d’integrazione numerica per studi e applicazioni pratiche di astrodinamica e astronomia dinamica.

La serie DE19 rilasciata nel Novembre del 1967 è stata la seconda distribuzione delle DE [2]. Essa costituisce un miglioramento della serie DE3. Nel 1969, il JPL inizia a distribuire la terza serie denominata DE69. Quest’ultima è una serie speciale limitata al periodo dal 28 Ottobre 1961 fino al 23 Gennaio del 1976 che apporta una correzione sulle effemeridi della luna. La serie fu calcolata in coincidenza con le missioni Apollo, per migliorare le posizioni della Luna tenendo conto delle osservazioni radar. Queste prime serie furono calcolate usando metodi perturbativi. Questi metodi usando le orbite di N-1 corpi del sistema planetario per integrare numericamente l’orbita del rimanente corpo e raffinarne le effemeridi.

Il 15 Novembre 1967, il JPL rilascia anche la serie DE28 che comprende il periodo 1800-2000 e si basa sulla integrazione diretta delle equazioni del moto. In particolare, queste effemeridi sono state calcolate usando il programma sviluppato da J. Schubart and P. Stumpff [ref]. Nel articolo che descrive il programma, gli autori enumerano le molteplici difficoltà nei metodi perturbativi, fino ad allora comunemente usati, rispetto alla integrazione diretta delle equazioni di Newton. La lettura di queste argomentazioni ci fa capire come potesse essere difficile effettuare simulazioni numeriche con i calcolatori elettronici di quell’epoca (Siemens 2002 e IBM 7090). Per esempio, uno dei problemi più difficoltosi era la gestione di dati in entrata e uscita, la lunghezza e complessità di programmi e la disponibilità di memoria temporanea. La versione usata dal JPL (the Goddard version) è documentata in un rapporto tecnico completo con il listato completo del programma in Fortran.

Schubart e Stumpff usarono per l’integrazione numerica delle equazioni del moto il metodo di Adams-Störmer all’ottavo ordine [7]. Questo metodo fu inizialmente ideato da J.C. Adams e F. Bashforth, two fisici inglesi. Adams si è occupato principalmente di problemi astronomici (uno dei risultati che lo ha reso famoso è stato quello del calcolo della posizione del pianeta Nettuno prima della sua scoperta, una storia molto interessante per l’intricata controversia che ne è derivata [8]). Nel 1883, Adams e Bashforth hanno pubblicato uno studio sulla forma delle gocce pendenti in cui viene derivata sulla base di dati sperimentali una equazione differenziale che porta il loro nome [9]. Per trovare la soluzione numerica di questa equazione, i due fisici derivarono uno schema d’integrazione numerico successivamente usato anche da altri studiosi. Per questo motivo, lo schema è anche conosciuto come metodo Adams-Bashforth. Questo metodo fa parte della famiglia dei metodi lineari a passo multiplo. I metodi di soluzione numerica di equazioni differenziali si basano sul calcolo di punti successivi delle soluzioni usando valori calcolati precedentemente. Nei metodi a passo singolo, come il metodo di Eulero, considerano solo uno o pochi punti precedenti per calcolare il successivo. Altri metodi usano punti intermedi per aumentare la risoluzione. Per potere iniziare l’integrazione, i metodi multi passo di tipo lineare vengono usati un numero di posizioni precedenti a quella iniziale e della sua derivata pari all’ordine del metodo meno uno.

Il metodo di Adams-Bashforth è stato adattato e modificato da altri scienziati tra cui C. Störmer che lo ha usato nei suoi pionieristici studi sull’interazione di particelle cariche con il campo magnetico terrestre. In questa forma modificata, esso viene anche chiamato metodo di Adams-Störmer. Il metodo di Adams è stato anche usato dai due astronomi P.H. Cowell e A.C. de la Cherois Crommelin che lo hanno usato per calcolare l’orbita della cometa di Halley. Infine, nel 1960, il metodo è stato applicato per la prima volta dal fisico francese L. Verlet nelle prime simulazioni di dinamica molecolare e da allora è conosciuto nell’ambito della chimica computazionale come metodo di Verlet.

La formula del metodo Adams-Störmer per calcolare i valori delle posizioni delle coordinate a tempi successive è data dalla relazione

\vec{r}_{n+1}-2\vec{r}_{n}+\vec{r}_{n-1}=h^2 \sum_{j=0}^{k-1}\sigma_j \nabla^j \vec{a}_j

con i valori delle differenze \nabla^n a_j, calcolate usando le relazioni ricursive:

\nabla^0 \vec{a}_j =\vec{a}_j

\nabla \vec{a}_j= \vec{a}_j - \vec{a}_{j-1}

\nabla^j \vec{a}_j= \nabla^{j-1} \vec{a}_j - \nabla^{j-1} \vec{a}_{j-1}

con h il valore del passo d’integrazione e i valori di \sigma_j per j=0\dots8 sono riportati nella seguente tabella

j012345678
\sigma10\frac{1}{12}\frac{1}{12}\frac{19}{240}\frac{3}{40}\frac{863}{12096}\frac{275}{4032}\frac{33953}{518400}

Cowell e Crommelin hanno usato la formula per j=6, mentre Schubart e Stumpff nel loro programma l’espansione è estesa fino all’ottavo ordine (j=8). Il metodo produce traiettorie molto stabili anche con ampi passi temporali. Come tutti i metodi multi passo, il metodo per essere applicato richiede la conoscenza di un numero di condizioni iniziali pari al suo ordine meno uno.

La disponibilità di calcolatori sempre più veloci e di accurati dati telemetrici ottenuti da misure ottiche, radar e da sonde spaziali, ha reso possibile l’aumento della precisione nella determinazione delle posizioni di pianeti, satelliti e altri corpi orbitanti del sistema solare. Pertanto sono state calcolate altre serie di effemeridi con accuratezza spaziale e temporale sempre più elevata [10]. La più recente serie, la DE438 è stata creata nel 2018 per migliorare la posizione di Mercurio e Giove in supporto alle nuove missioni NASA verso questi pianeti. In questi ultimi sviluppi vengono tenuto in considerazione non solo i principali pianeti con i satelliti ma anche altri oggetti come asteroidi. Vengono tenuti in conto effetti relativistici e nel caso dell’orbita della luna le accelerazioni dovute alla forma della Terra. Queste nuove serie hanno accuratezze che vanno da pochi metri (per la posizione della Luna), a poche centinaia or addirittura decine di metri per i pianeti interni del sistema solare. Basti pensare al successo della missione New Horizon che nel 2015 ha raggiunto e fotografato passando a una distanza di soli 12500 km dalla superfice del pianeta [11,12].

BIBLIOGRAFIA

  1. D. Brouwer and G.M Clemence. 3. Orbits and Masses Planets and Satellites, in Gerard P. Kuiper & Barbara M. Middlehurst, Planets and Satellites ( The Solar System Volume III ) ,University of Chicago Press 1961. 
  2. http://www.columbia.edu/cu/computinghistory/eckert.html.
  3. Computing in the Soviet Space Program:An Introduction. http://web.mit.edu/slava/space/introduction.htm
  4. Devine, C.J. JPL Development Ephemeris Number 19, technical Report 32-1181, JPL, Pasadena, Calif., Nov.15, 1967.
  5. Cowell, P.H. and De la Cherois Crommelin, A.C., 1910. Investigation of the Motion of Halley’s Comet from 1759 to 1910. Neill & Company, limited.
  6. O’Handley, D.A;, Holdridge, D.B; Melbourne, W.G.; Mulholland, D.J. JPL Development Ephemeris Number 69, Technical Report 32-1465, JPL, Pasadena, Calif. Dec. 15, 1969.
  7. J. Schubart and P. Stumpff, On an N-Body Program of High Accuracy for the Computation of Ephemerides of Minor Planets and Comets, Veröff Astron, Rechen-Inst., Heidelberg No. 18, Verlag. G. Braun, Karlsnihe, 1966.
  8. Robert W. Smith. La scoperta di Nettuno e il problema del perielio di Mercurio. http://www.treccani.it/enciclopedia/l-ottocento-astronomia-la-scoperta-di-nettuno-e-il-problema-del-perielio-di-mercurio_%28Storia-della-Scienza%29/
  9. Bashforth, F. and Adams, J. C. Theories of Capillary Action. London: Cambridge University Press, 1883.
  10. JPL Planetary and Lunar Ephemerides. https://ssd.jpl.nasa.gov/?planet_eph_export
  11. http://pluto.jhuapl.edu
  12. I segreti di Plutone e delle sue lune svelati da New Horizons. Le Scienze. http://www.lescienze.it/news/2016/03/18/news/plutone_caronte_new_horizons-3020365/

About Danilo Roccatano

I have a Doctorate in chemistry at the University of Roma “La Sapienza”. I led educational and research activities at different universities in Italy, The Netherlands, Germany and now in the UK. I am fascinated by the study of nature with theoretical models and computational. For years, my scientific research is focused on the study of molecular systems of biological interest using the technique of Molecular Dynamics simulation. I have developed a server (the link is in one of my post) for statistical analysis at the amino acid level of the effect of random mutations induced by random mutagenesis methods. I am also very active in the didactic activity in physical chemistry, computational chemistry, and molecular modeling. I have several other interests and hobbies as video/photography, robotics, computer vision, electronics, programming, microscopy, entomology, recreational mathematics and computational linguistics.
This entry was posted in Uncategorized. Bookmark the permalink.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

This site uses Akismet to reduce spam. Learn how your comment data is processed.