Laboratorio di Calcolo. AA 2019/20.

Esercitazione n. 4

 

Richiami sulla sintassi

Il valore dell' attributo kind corrispondente ad una data precisione/estensione può variare da sistema operativo a sistema operativo e da compilatore a compilatore. Il valore numerico dell' attributo kind non è di nessuna importanza o utilità in sé. Serve solo a stabilire, attraverso un unico numero, una corrispondenza biunivoca tra richieste di precisione e/o estensione e codifica binaria dei tipi dati (quanti bit e con che significato).

Il modo corretto e consigliato di usare il meccanismo offerto dall' attributo kind è di definire preventivamente, per ogni richiesta di diverso kind, una costante simbolica intera da associare a quel parametro kind, attraverso l' uso delle due funzioni di libreria

SELECTED_INT_KIND(d) , che interroga il compilatore sul kind da associare al tipo integer, in modo da garantire almeno d cifre decimali, ovvero in grado di rappresentare esattamente tutti gli interi strettamente minori in valore assoluto di 10d.

SELECTED_REAL_KIND(d) oppure SELECTED_REAL_KIND(d,e), che interrogano il compilatore sul kind da associare al tipo real, in modo da garantire, nel primo caso, almeno d cifre decimali di precisione nella mantissa e nel secondo, oltre a d cifre decimali di precisione nella mantissa, anche un intervallo tra le potenze di 10 compreso tra -10e e 10e.

Per i real le rappresentazioni standard a 32,64,80 e 128 bit corrispondono a richieste di precisione con numero di cifre decimali fino a (incluso) 6, 15, 18 e 33.

p.es:

integer, parameter :: ik=selected_int_kind(12)
integer, parameter :: rk=selected_real_ind(5)
integer(kind=ik) :: big_int
real(kind=rk) :: x,y,z

in cui si richiede al compilatore un tipo dati per la variabile intera big_int in grado di rappresentare intero compreso tra -10-12 e 1012 (estremi esclusi). E un tipo dati per le variabili x,y, e z, in grado di rappresentare reali con almeno 5 cifre decimali di accuratezza nella mantissa

Il parametro kind utilizzabile nelle dichiarazioni di tipo dati, va anche utilizzato per le costanti:

12349871238_ik
3.1234_rk
22.234e-2_rk

Nota: In particolare, se volessi una variabile con una certa precisione rk, e la inizializzassi senza mettere _rk in fondo al valore scelto, il programma leggerebbe il numero con la precisione di default prima di assegnarlo alla variabile, vanificando così un mio eventuale aumento di precisione ( dunque, se avessi dichiarato una variabile x con kind rk, per avere effettivamente la precisione richiesta all'atto di inizializzazione dovrei scrivere per esempio x=3.1_rk e non x=3.1)



  1. Modificare il programma per il calcolo dell' esponenziale della esercitazione n. 3 (es. n.2) in modo da lavorare con precisione maggiore rispetto a quella di default. Verificare quanti termini in più sono necessari per ottenere il risultato migliore.

  2. Scrivere un programma che calcoli  i valori della funzione definita dall' espressione Fortran: 
    sqrt(2/pi/x)*((3/x**2-1)*sin(x)-3/x*cos(x))
    Dove 
    pi
    rappresenta la costante pi greco ( π ), il cui valore può  esser  calcolato nel programma mediante le funzioni trigonometriche inverse (p.es.  pi=acos(-1.0) ).  Trascrivere  la formula data sopra  in notazione matematica usuale su un foglio e verificare che il limite nell' origine sia zero (ma l' x nella radice quadrata è al numeratore o al denominatore ? perché ?). Il modo più comodo (ed istruttivo)  per verificare il limite è di sostituire ai termini in seno e coseno i primi  termini delle loro espansioni in serie di Taylor.       
    Valutare la funzione in 100 punti equispaziati nell' intervallo aperto inferiormente (0, 0.1]  (escludendo quindi lo zero) e scrivere i valori di ascisse ed ordinate in un file.  Può essere utile riadattare a questo scopo il programma sviluppato per il punto 9 dell' esercitazione n.2. Utilizzare la precisione (kind) di default per i real. Discutere se, e in che parte dell' intervallo, ci sono effetti di arrotondamento (roundoff) importanti. (Può essere utile graficare la funzione mediante gnuplot). Rifare il grafico usando precisione immediatamente maggiore (verificando l' intervallo utilizzato da gnuplot per le ordinate).  Con la precisione maggiore, modificare poi l' intervallo su cui si calcolano i punti in (0, 0.001]. Ricordare che anche le costanti hanno un kind e che il risultato di ogni funzione intrinseca ha lo stesso kind dell' argomento. Un modo  per evitare del tutto le cancellazioni sottrattive a piccoli x  è  suggerito nell'ultimo degli esercizi da svolgere autonomamente alla fine della pagina.


  3. Lavorando in precisione di default, calcolare numericamente le somme  ∑i=1N 1/i,  e  ∑i=N1 1/i  (cambia solo l' ordine di somma)   in funzione di N, per diversi valori di N  da 10000  a 1000000 ,  (non tutti, ma a salti di 10000) scrivendo i risultati su file. Se i risultati differiscono quale sarà il più accurato ?

    Come possiamo verificarlo?



Ulteriori esercizi (possono essere svolti individualmente). Cercare di svolgerli ed eventualmente discuterli sul Forum (a partire da venerdì) sarebbe molto utile. Si ricorda che la partecipazione attiva al Forum dà un bonus per l' esame.

  1. Le seguenti tre formule sono algebricamente equivalenti e ci si potrebbe aspettare che numericamente convergano (per N che tende a infinito) allo stesso valore (somma della serie).  Verificare cosa succede implementandole in un programma Fortran (lavorando con la precisione di default) e  individuare la  forma che si  comporta meglio numericamente. Provare a darne una spiegazione.
    1. ∑n=12N   (-1)n  n/(1+n)
    2. ( ∑n=1N  2n/(2n+1)  )   -  (  ∑n=1N  (2n-1)/(2n))
    3. ∑n=1N    1/(2n(2n+1))



  2. Compilare il programma seguente 
    program abc
    implicit none
    integer,parameter :: rk=selected_real_kind(6)
    real(kind=rk) :: x=1.0_rk,y=1.0_rk,z
    integer ::i
    do i=1,100000000
      z=x*(y+i)
    end do
    end program abc

    ed eseguirlo col comando
    time ./a.out

    Delle tre righe di output, quella che inizia con user dà il tempo di CPU utilizzato dal processore per eseguire le istruzioni macchina corrispondenti al programma Fortran. Prendere nota del risultato e aumentare progressivamente la precisione richiesta. Cosa succede ai tempi di esecuzione ?


  3. Usare il programma che risolve l'  equazione di secondo grado per trovare le due radici dell'  equazione x2 + 80000 x + 1 = 0.  La soluzione corrispondente al segno positivo davanti alla radice quadrata del discriminante  corrisponde ad un caso di cancellazione sottrattiva. Valutare l'  errore percentuale confrontandone il risultato  con quello ottenuto usando la precisione doppia (accuratezza di 15 cifre significative)  e poi,  tornati alla precisione di default, con quello ottenuto  dalla seguente formula algebricamente equivalente a quella usuale: 
    - 2 c / ( b + sqrt(discr) )

  4. Per evitare le cancellazioni sottrattive per piccoli valori di  x nella funzione dell' esercizio N.4,  si può procedere utilizzando, per i valori di x tra 0 e una soglia  x,  i primi termini dell' espansione in serie di Taylor  della funzione che moltiplica la radice quadrata. Il punto di raccordo tra formula basata sul' espansione di Taylor e la forma completa deve essere determinato in modo da assicurare una transizione liscia tra una forma funzionale e l' altra.

  5. Costrutto di salto condizionato SELECT-CASE.
    In alternativa a una serie di IF-THEN/ELSE-IF-THEN...ELSE-ENDIF   il Fortran mette a disposizione un diverso meccanismo per realizzare salti a sezioni diverse del codice in modo diretto  senza passare per una cascata di verifiche di condizioni successive (il caso degli IF...).
    La sintassi è schematicamente la seguente:

    SELECT CASE (espressione a  valori interi o  sequenze di caratteri (stringhe))

    CASE (lista di possibili valori (costanti) interi o stringhe separati  da virgola)
            istruzioni da eseguire qualora espressione abbia valori presenti in questo CASE
    CASE (altra lista )
            ...
    CASE (ancora un'  altra...)
            ...
    CASE DEFAULT
            istruzioni da eseguire qualora espressione assuma un valore non presente nei CASE  espliciti
    END SELECT

    Nei vari CASE potranno comparire anche intervalli di valori indicati  mediante gli estremi separati da due punti ':' ). Es:  CASE(1:8) corrisponde a tutti i valori interi tra 1 e 8  inclusi;  CASE ('fabbro':)  corrisponde a tutte  le stringhe che seguono o  sono uguali, nell'  ordine alfabetico, a quella costituita dalla stringa 'fabbro'.
    Scrivere un programma che legga un mese nella forma numerica e cioè un intero da 1 (gennaio) a 12 (dicembre) e dia su schermo il numero dei giorni del mese (per anni non bisestili) utilizzando l' istruzione  SELECT CASE (espressione)... CASE (selettore)...END SELECT. In pratica, a fronte dell' input
    3
    il programma deve mostrare su schermo il numero 31 (i giorni di marzo).

  6. Scrivere un programma che determini il più piccolo numero positivo denormalizzato per ciascuno dei kind possibili per i numeri di tipo REAL.