Laboratorio di Calcolo. AA 2018/19.

Esercitazione n. 10

1. 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.


2. Scrivere un programma che integri la funzione e^-x2 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).



3. 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)
   

1. 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...


2. 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.



3. 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.