esercitazione n. 9

Laboratorio di Calcolo. AA 2019/20.

Esercitazione n. 9

Scrivere un programma che integri la funzione x⁴ tra 0 e 1 mediante due sottoprogrammi (uno di tipo function e uno subroutine) che implementino rispettivamente la formula trapezoidale e quella di Simpson . La funzione integranda deve essere calcolata in una function di un argomento real il cui nome sarà passato come argomento alla procedura che valuta l' integrale. Una procedura Fortran (subroutine o function) può avere tra i suoi argomenti il nome di un'altra procedura. Questo argomento va dichiarato nel sottoprogramma mediante un' interfaccia esplicita^(*).

Verificare l' accuratezza del risultato in funzione del numero di punti, e della precisione dei real, partendo da una decina di punti fino a qualche milione. Si ricorda che la formula trapezoidale (estesa) corrisponde a sommare i valori della funzione su punti equispaziati, a distanza h, attribuendo a ciascuno un peso h, tranne gli estremi che intervengono con peso h/2. La formula di Simpson (estesa) fa invece intervenire pesi h/3,4h/3,2h/3,4h/3,2h/3,...,4h/3,h/3 e quindi presuppone un numero dispari di punti (ovvero un numero pari di intervallini). In particolare andrà verificata la dipendenza dell' accuratezza dell' integrale dal numero di punti. Si ricorda che l' errore algoritmico del metodo trapezoidale va a zero con 1/N², mentre nel caso Cavalieri-Simpson la dipendenza è 1/N⁴, per funzioni con derivata seconda e quarta limitata, rispettivamente.

(*) Esempio:
```
module prec
   integer, parameter :: rk= selected_real_kind(6)
end module prec
module funzi
   use prec
   implicit none
 contains
   function f(fun,a,b) result (r)
      real(kind=rk),intent(in) :: a,b
      real(kind=rk) :: r
      interface
         function fun(x) result(res)
           use prec
           real(kind=rk), intent(in) :: x
           real(kind=rk) :: res
         end function fun
        end interface 
     .....
     end function f

     function fun1(x) result(r1)
     ...
     end function fun1
     function fun2(x) result(r1)
     ...
     end function fun2 
end module funzi 
program pippo
use funzi
implicit none
real(kind=rk) :: q,s,a,b
...
 q=f(fun1,a,b)
 s=f(fun2,a,b)
...
```
dove fun1 e fun2 sono funzioni con la stessa interfaccia di fun.

Scrivere un programma che integri la funzione e^-x² tra -∞ e +∞ mediante la formula trapezoidale e verificare l' accuratezza del risultato in funzione del numero di punti e degli estremi di integrazione. Evidentemente non è possibile integrare la funzione su un intervallo infinito ma è possibile approssimare l' integrale sull' intervallo infinito mediante integrali su intervalli finiti di ampiezza tale da rendere arbitrariamente piccolo l' errore dovuto al troncamento. Per verifica, è utile tener presente che l' integrale in questione vale radice quadrata di pi greco. Per la determinazione ragionevole degli estremi si suggerisce di osservare preliminarmente il grafico della funzione integranda (p.es. tramite gnuplot).

Calcolare, usando sia il metodo dei trapezi, sia il metodo di Cavalieri-Simpson, l' integrale tra 0 e 1 della funzione x^3/2 per valori crescenti del numero di intervallini (N) fino a qualche decina di migliaia, ricavare da un grafico log-log la dipendenza da N dell' errore (si consiglia di scrivere su file log(N) e log( | intg-0.4 | ), dove intg è il valore numerico dell' integrale e 0.4 il valore analitico) e cercare di darne una spiegazione.

Non è richiesto (per svolgere l' esercizio basta una stima visuale), ma volendo andare più sul quantitativo, si potrebbe far ricavare a gnuplot la legge di potenza per l' errore, chiedendo il miglior fit dei dati mediante una retta. Le istruzioni sono:
```
b=1
a=-1
f(x)=a*x+b
fit [0:4] f(x) 'fort.1' via a,b
```
supponendo che i dati log-log siano in fort.1 (colonne 1 e 2), che il regime lineare corrisponda alle ascisse tra 0 e 4 e di non essere andati oltre qualche decina di migliaia con N.
La verifica della bontà del fit si può avere mediante:
```
p 'fort.1' ,f(x)
```

Istogrammi.
Scrivere un programma che generi N variabili casuali, ognuna delle quali sia la media aritmetica di 100 numeri casuali distribuiti uniformemente tra 0 e 1. Si ricorda che la subroutine intrinseca random_number(x), chiamata con argomento scalare o array, ritorna uno o più valori real distribuiti uniformemente tra 0 e 1.Ciascun numero così generato dovrà essere scritto in un file di nome DATI. Il file deve quindi essere chiuso dal programma e riaperto in sola lettura. Il programma deve quindi leggere i dati e costruire un istogramma delle frequenze dei dati su 100 intervalli di ampiezza 1 tra 0 e 1.
Il seguente frammento di codice gestisce la creazione dell' istogramma. Può essee utilizzato per costruire il programma completo.
```
real :: dx, xmin, xmax
real :: dato
integer,dimension(:),allocatable   :: histo
integer :: i,ind,fuori,nbin,Ndati  ! Ndati rappresenta il numero dei dati nel file 
                                   ! nbin il numero dei canali (intervallini) dell' istogramma.


! generazione di Ndati numeri casuali ottenuti sommando 100 n. casuali 
! distribuiti uniformemente tra 0 e 1 e scrittura deglistessi su file di nome DATI
! chiusura del file e riapertura in lettura

print*,' xmin, xmax, nbin:? '
read*, xmin, xmax, nbin

allocate(histo(nbin))
histo = 0
fuori = 0

dx = (xmax-xmin)/nbin

do i = 1,Ndati
  read(1,*)dato
  ind = int((dato-xmin)/dx) + 1
  if ( ind>nbin .or. ind <1)then
    fuori = fuori + 1
  else
    histo(ind) = histo(ind) + 1
  end if
end do

do i=1,nbin
  write(2,*) i, xmin+(i-1)*dx+dx/2,real(histo(i))/Ndati*nbin
end do

print*,
print*,"# fuori: ", fuori
```
L'istogramma, i cui dati saranno su file fort.2, può essee graficato con gnuplot col comando
plot [0:1] 'fort.2' u 2:3 w boxes fs solid .5
Verificare l'effetto di variare il parametro .5 tra 0 e 1. E anche l'effetto di dare il comando
set boxw 0.6 relative
prima del comando plot.

Esercizi supplementari.

Calcolare numericamente, col metodo trapezoidale, l'integrale tra 0 e x della funzione sin(10 t) e confrontarlo con il valore esatto ( 1-cos(10 x) )/10 per 101 valori di x uniformemente spaziati nell' intervallo [0,1] (inclusi gli estremi).Suggerimento: piuttosto che ricalcolare per ogni x l' integrale da 0 a x, sarebbe più efficiente aggiungere all' integrale precedentemente calcolato solo il contributo dell' ultimo intervallino...

Il programma lagrange.f90 implementa il calcolo delle interpolanti di Lagrange di vari ordini, con nodi equispaziati, di una funzione definita nella function f . Controllare come il programma implementa il calcolo delle interpolanti e verificare visivamente, graficando i risultati con gnuplot, come cresce , per gradi di interpolazione elevati, il massimo scarto tra interpolante e funzione nell' intervallo [-1,1].

Modificare la definizione dei nodi , nella function Poly, passando dai nodi equispaziati ai cosiddetti nodi di Chebyshev:
nodo=[ ( cos( real((2*i-1),k)/( 2*(g+1))*acos(-1.0_k) ),i=g+1,1,-1 ) ]
e verificare il comportamento delle interpolanti al crescere del grado, confrontandolo con il caso a nodi equispaziati.

In questo esercizio utilizzare solo variabili dei tipi dati primitivi. Un metodo semplice (anche se non molto efficiente) per minimizzare funzioni di più variabili, derivabili almeno una volta, è il cosiddetto metodo dello steepest descent (massima pendenza) a passo fisso. Nella implementazione più semplice, si consideri ad esempio la funzione da minimizzare U(x,y,z). L' algoritmo può essere descritto come segue:

     1. si parte da un punto di coordinate x_n, y_n,z_n (all' inizio, n=0) e si fissa il parametro α (vedi dopo)
     2. si calcola il gradiente della funzione nel punto   x_n, y_n,z_n, usando le formule analitiche per le tre derivate parziali
        (il gradiente è il vettore che ha per componenti le derivate parziali rispetto alle variabili indipendenti della funzione in questione).
     3. si trova un nuovo punto con la seguente formula:
        x_n+1 = x_n - α ∂U/∂x(calcolata in x_n, y_n,z_n)
        y_n+1 = y_n - α ∂U/∂y(calcolata in x_n, y_n,z_n)
        z_n+1 = z_n - α ∂U/∂z(calcolata in x_n, y_n,z_n)
     4. si assegna a x_n il valore x_n+1 , e così via per y e z, e si torna al punto 2 fino a che non si raggiunge una precisione prefissata per il minimo (vedi dopo).
Il parametro α va scelto (per tentativi) in modo che sia abbastanza grande da garantire progresso nella ricerca del minimo ma non tale da portare il sistema lontano dal minimo.

Un modo per giudicare se si è arrivati al punto fisso dell' iterazione 1-4 è di calcolare una "misura" del progresso che si ottiene tra la n e la n+1 esima iterazione. Al punto fisso della trasformazione evidentemente non ci sarà più nessuna variazione. P.es.:
      s = (x_n+1 - x_n)² + (y_n+1 - y_n)² + (z_n+1 - z_n)².
Se s decresce ad ogni passo si è scelto un α ragionevole. Se si vedono oscillazioni o divergenze di s occorre ridurre il valore di α. Occorrerà determinare la soglia per s al sotto della quale considerare l' algoritmo a convergenza.

Implementare ed applicare questo algoritmo al caso della funzione    A(x - 1)² + B(y - 1)² + C(z - 1)²
che ha minimo pari a zero nel punto x=1, y=1, z=1, nei casi   A=1, B=2, C=3    e A=1, B=20, C=300, partendo, in tutti e due i casi, dal punto   di coordinate x=3, y=4, z=3.
Verificare quale sia il massimo valore di α utilizzabile nei due casi facendo scrivere su file ad ogni iterazione: i valori delle tre variabili indipendenti, di s e della funzione da minimizzare.