Laboratorio di Calcolo. AA 2019/20.

Esercitazione n. 8


Si ricorda che tutte le procedure devono essere implementate come procedure di modulo.


  1. L' algoritmo  di Euclide  per la ricerca del massimo comun divisore  (MCD)  tra due numeri naturali può  essere descritto in questo modo:
    dati due numeri naturali, a e b non entrambi nulli, se b è zero allora  il MCD è a, altrimenti il risultato è dato dal MCD tra  b e  il resto della divisione tra a e b.
    Implementare questo algoritmo ricorsivo mediante una function ricorsiva  (da dichiarare come  RECURSIVE FUNCTION ...). Esiste una funzione intrinseca MOD(x,y)  il cui valore è il resto della divisione tra x e y. Nell' implementazione ipotizzare che la function non venga mai chiamata con tutti e due gli argomenti pari a zero.


  2. Il programma che integra le equazioni del moto del problema gravitazionale dei 2 corpi della scorsa esercitazione dovrebbe essere in grado, con poche modifiche, di risolvere anche il problema degli N-corpi in interazione gravitazionale. Si può quindi aggiungere Giove al sistema Terra-Sole per studiarne l' effetto sull'orbita terrestre. Per semplicità,   assumiamo che l' orbita di Giove sia esattamente complanare con quella terrestre e prendiamo come condizioni iniziali quelle in cui sia la Terra sia Giove sono all' afelio ed allineati dalla stessa parte: Sole nell' origine, Terra sull' asse x a x=152.1 109 m  e Giove a x=816.62 109 m con velocità nulla per il Sole, vy=29.29 103  m s-1 per la Terra e vy=12.44 103  m s-1 per Giove. Una domanda a cui si potrebbe cercare di rispondere è: come la presenza di Giove modifica la posizione del perielio della terra dopo molte orbite?  La massa di Giove è: 1898.6 1024 kg.
    Suggerimento: subito dopo la lettura delle velocità, modificare le velocità iniziali sottraendo alla velocità iniziale di ciascun corpo la velocità del centro di massa (da calcolare). In tal modo il sistema di riferimento sarà quello in cui il centro di massa è fermo. Può  anche esser utile scrivere su file la posizione della Terra relativa al Sole,  oltre a quella riferita all'    origine del sistema inerziale.
    Un modo per avere un' idea qualitativa delle modifiche dell'orbita (che non sono enormi) è  di graficare solo la prima e l'ultima delle orbite terrestri su un periodo di una trentina di anni terrestri (non molto più  di 10000 giorni). Un modo semplice per ottenere il grafico in gnuplot è  mediante il seguente comando:

    set siz ratio -1
    plot [-1.55e11:1.55e11] [-1.55e11:1.55e11] 'fort.1' u 6:7 every ::1::365 w l , 'fort.1' u 6:7 every ::(30000-365)::30000 w l


    dove si suppone di aver effettuato 30000 passi com time step di 1 giorno. Si è  assunto che le coordinate dell'orbita terrestre siano nelle colonne 6 e 7  e che l' intervallo di integrazione temporale dt sia pari a un giorno terrestre (86400 s). Eventualmente modificare questi parametri se non adeguati al run effettuato.

    Per chi volesse invece considerare Venere, al posto della Terra, la sua massa è 4.8685 10 24 kg, la distanza all' afelio 108.94 109 m, velocità all' afelio: 34.79  103 m s-1.
    Si consiglia di cancellare i file fort.... dopo aver finito.


  3. Definiamo un tipo derivato punto,  costituito da  due campi reali corrispondenti alle coordinate cartesiane di un punto nel piano, ed un tipo derivato triangolo, costituito da tre campi di tipo punto. Sulla base degli esempi, forniti nel seguente modulo, per la  somma di due punti e per l' inversione nel piano rispetto all' origine,  apportare le seguenti modifiche:
    module geom
    implicit none

    type :: punto
    real ::x,y
    end type punto

    type :: triangolo
    type(punto) :: vertice_A
    type(punto) :: vertice_B
    type(punto) :: vertice_C
    end type triangolo

    contains

    function sum_points(p,q) result(r)
    type(punto), intent(in) :: p,q
    type(punto) :: r
    r%x = p%x + q%x
    r%y = p%y + q%y
    end function sum_points

    function invert(p) result(r)
    type(punto), intent(in) :: p
    type(punto) :: r
    r%x = -p%x
    r%y = -p%y
    end function invert

    end module geom


    program piano
    use geom
    implicit none
    type(punto) :: A,B,C
    type(triangolo) :: T1, T2

    print*,'coordinate 2D primo vertice'
    read*,A
    print*,'coordinate 2D secondo vertice'
    read*,B
    print*,'coordinate 2D terzo vertice'
    read*,C
    T1 = triangolo(A,B,C) ! il nome del tipo dati "costruisce" una variabile
    ! dello stesso tipo a partire dai 3 campi
    print*,' vertice A: ',T1%vertice_A
    print*,' vertice B: ',T1%vertice_B
    print*,' vertice C: ',T1%vertice_C


    T2%vertice_A = invert(T1%vertice_A)
    T2%vertice_B = invert(T1%vertice_B)
    T2%vertice_C = invert(T1%vertice_C)

    print*,'coppie di coordinate dei vertici del triangolo T2', T2

    end program piano


    I triangoli potranno essere visualizzati mediante l' istruzione gnuplot
    p 'triangolo.dat' index 0 w l,'' index 1 w l
    Index 0 e index 1 fanno riferimento al primo ed al secondo set di dati del file. Set diversi sono separati da due linee vuote.


  4. Scrivere un programma che implementi l' algoritmo di bisezione per la determinazione della radice quadrata di un numero reale z  mediante una subroutine cui passare come parametri di input: gli estremi dell' intervallo iniziale e il parametro di tolleranza che controlla il criterio di convergenza scelto. La funzione di cui vogliamo cercare lo zero positivo come funzione di x (a zeta fissato) è  x2-z. Cercare di implementare l' algoritmo nel modo più semplice possibile senza però comprometterne la generalità (in particolare evitare implementazioni mediante funzioni ricorsive). Cercare, nell' implementazione, di seguire il più vicino possibile la descrizione dell' algoritmo. Esplorare la dipendenza del numero di iterazioni  dalla richiesta di accuratezza della soluzione,  decidendo quale criterio di stop sia più  opportuno  controllare la precisione assoluta o relativa della soluzione. Provare a calcolare (confrontando con il risultato di sqrt) la radice quadrata di:  4, 3, 1.0e6, 1.01e6 e 1e-20.

  5. Scrivere un programma che  implementi l' algoritmo di Newton per la determinazione della radice quadrata di un numero reale z (la funzione da azzerare è la stessa dell' esercizio precedente) mediante una subroutine. La derivata che appare nella formula di Newton va calcolata  analiticamente (analiticamente qui va inteso come il contrario di numericamente). Esplorare la dipendenza del numero di iterazioni  dalla richiesta di accuratezza della soluzione e  usando come criterio di stop la richiesta che i risultati (per il punto di zero)  di due iterazioni successive differiscano  meno di una soglia stabilita dall' utente (perché possiamo usare questa condizione ?). Si suggerisce di far scrivere il risultato di ogni approssimazione successiva  e confrontare i risultati con quelli ottenuti mediante l'  algoritmo di bisezione.  Provare anche a cercare uno zero della funzione  intrinseca sin(x) partendo dal punto x=3.


Esercizi supplementari.

  1. Una function può  anche ritornare un risultato di tipo  array. In tal caso,  il risultato dovrà  avere l' attributo dimension. Se l' array ritornato dalla function deve avere le stesse dimensioni di un argomento, queste possono essere determinate facendo ricorso alla function intrinseca  size

    Data la function
    function raddoppia(x) result(dop)
    real,dimension(:),intent(in) ::x
    real,dimension(size(x,1)) ::dop
    dop= 2*x
    end function raddoppia

    Utilizzarla all' interno di un programma che legga l' array in input da tastiera e scriva l' array risultato dell' elaborazione della function su schermo.
    Verificare cosa succede se si chiama   la function precedente con un argomento che non sia un array (Es. raddoppia(3.0) ).
    Per poter scrivere fuction che si comportino come le function intrinseche del Fortran (che funzionano indifferentemente con input di scalari, array di rango 1,2,3 etc) si deve ricorrere alle cosiddette funzioni elementali.

    Trasformare raddoppia in elemental function e utilizzarla in un program in presenza di argomento di tipo scalare e uno di tipo array.


  2. Si può tentare la soluzione numerica dell' equazione di terzo grado x3 + 4 x -2 =0  con diversi metodi iterativi, ottenendo in vario modo dall'  equazione una successione del tipo  xn+1 = g(xn).  Manipolando algebricamente i termini dell' equazione di terzo grado  si possono facilmente ottenere le seguenti  iterazioni:
    xn+1=(2-xn3)/4    (ottenuta isolando il termine 4x e dividendo per 4)
    xn+1=(2-4xn)/xn2 (ottenuta isolando il termine x3 e dividendo per x2 )

    xn+1=2 - 3 xn - xn3 (ottenuta scrivendo 4x come x+3x e isolando il termine in x )
    Inoltre, dal metodo di Newton se ne può ottenere una quarta:
    xn+1=(2 xn3 + 2)/(3 xn2 + 4)
    Scrivere un programma che implementi queste formule e scriva i risultati per ciascuna formula su un file diverso. Verificare (approssimativamente, per tentativi) quale ha il massimo bacino di attrazione (il bacino di attrazione di un punto fisso è costituito dall' insieme dei punti iniziali che danno luogo ad una successione convergente a quel punto fisso). Verificare che le instabilità sono legate alla presenza di  regioni in cui |g'(x)|>1,  dove g(x)  è  il membro a destra di ciascuna formula.Si suggerisce di iterare ogni formula per non più di 100 passi, facendo scrivere ad ogni passo il valore della variabile x e della derivata prima della funzione (g'). Non è  richiesto di implementare criteri di stop ma solo di scrivere i valori ad ogni passo.
  3. Può essere utile utilizzare la seguente funzione che verifica se il valore di una variabile real, di kind=rk, definito nel modulo prec, è "normale" ovvero diverso da Inf  oppure NaN.
    module prec 
    implicit none
    integer,parameter :: rk=selected_real_kind(12)
    end module prec

    module ieee_test
    use prec
    implicit none

    contains

    function is_normal(x) result(ans)
    real(kind=rk),intent(in) :: x
    logical :: ans
    ans=.true.
    if(abs(x)>huge(x)) ans=.false.
    if(x/=x) ans=.false.
    end function is_normal

    end module ieee_test



  4. Dati inziali accurati per i corpi del sistema solare sono ottenibili dal sito NASA http://ssd.jpl.nasa.gov/horizons.cgi. Sulla Web interface si consiglia di selezionare Ephemeris Type: VECTOR TABLE, Coordinate Origin "@0" (centro di massa del sistema solare),output units (nelle table settings) KM-S (le lunghezze sono in km e velocità in km/s).
    Provare a vedere che forma ha l' orbita della luna nel sistema di riferimento centrato sul Sole.
    Le coordinate iniziali (data 10/5/2010) ricavate dal sito NASA sono (copiarle con tutte le cifre, senza tentare arrotondamenti):

          x               y             z
    -8.5914189E+07 -6.9058925E+08 -1.101059E+07    sole
    -1.0086797E+11 -1.1314793E+11 -8.425213E+06    terra
    -1.0047512E+11 -1.1322453E+11  2.409495E+06    luna

    le velocità iniziali (in m/s) sono

          vx              vy            vz
    1.96391722E+03 -2.40343250E+02 -2.497495E+02   sole
    2.36617978E+04 -2.02444140E+04 -2.497575E+02   terra
    2.38176782E+04 -1.92814345E+04 -2.161803E+02   luna

    La massa della luna è 0.073 1024 kg. 

  5. Si può verificare che la traiettoria della Terra rispetto al Sole ottenuta integrando le equazioni del moto nel problema dei due corpi è  un' ellisse.
    Per dimostrarlo numericamente  si può procedere in modi diversi. Ad esempio, verificando  che la somma delle distanze della Terra da due punti, di cui uno è  il centro del Sole, è costante. Conviene lavorare nel sistema di coordinate centrato sul Sole,  usando la differenza tra il vettore posizione della Terra e  quello del Sole. Uno dei due punti cercati  (fuochi) coincide quindi con l' origine. Il secondo, facendo partire la Terra dall' afelio sull' asse x positivo, deve stare sull' asse x ad un valore positivo pari alla differenza tra distanza dell' afelio  e quella del perielio.  Usando l'  afelio come punto iniziale,   la  distanza del perielio sarà determinabile dalla coordinata dell' altro punto in cui la traiettoria interseca l' asse x (dopo metà  anno). Mostrando la distanza da uno dei fuochi e la somma delle distanze sullo stesso grafico si dovrebbe poter verificare che la somma è  costante sulla scala delle variazioni della prima. Accertarsi però di lavorare nel sistea di riferimento in cui il centro di massa è  in quiete.