Lezione 7 : quadratura numerica

In questa lezione implementeremo alcuni algoritmi di quadratura numerica, cioè algoritmi per il calcolo di integrali definiti di funzioni continue in un intervallo chiuso e limitato. I metodi numerici che studieremo in questa sessione si possono rendere necessari in varie occasioni: quando non sappiamo valutare analiticamente l'integrale in esame, quando la soluzione analitica è molto complicata ed il calcolo numerico è molto più semplice, oppure quando la funzione è conosciuta in un numeri finito di punti.

ESERCIZIO 7.0 - Integrazione con il metodo midpoint a numero di passi fissato:

In questo primo esercizio proviamo ad implementare un codice per il calcolo dell'integrale della funzione $x sin(x)$ tra $[0,\pi /2 ]$ con il metodo midpoint.

1. Per prima cosa costruiamo un programma di test che calcoli l'integrale utilizzando un numero di intervalli fissato e passato da riga di comando.
2. Per controllare la precisione ottenuta con un numero di passi fissato proviamo a stampare una tabella con la differenza tra il risultato numerico ed il valore esatto (ottenuto analiticamente) in funzione del numero di passi (o della lunghezza del passo h). In aggiunta alla tabella si può rappresentare l'andamento dell'errore in funzione della lunghezza del passo h con un TGraph di ROOT (una traccia su come visualizzare i dati usando un TGraph di ROOT si può trovare qui).

Il metodo midpoint

Ricordiamo che in questo metodo l'approssimazione dell'integrale è definita dalla formula

$$\int_{a}^{b} f(x) dx \approx ( f(x_{0}) + f(x_{1}) + ... + f(x_{N-1}) ) h$$

dove

$$h = \frac{b-a}{N}$$

e

$$x_{k} = a + \left( k + \frac{1}{2} \right) h \qquad k = 0, 1.... N-1$$

Questa formula fornisce un'accuratezza dell'integrale di $O(h^{2})$. Notate che questo metodo non richiede il calcolo della funzione negli estremi di integrazione.

Cenni sull'implementazione :

L'algoritmo può essere implementato con uno schema del tipo classe madre astratta e classe derivata che abbiamo già utilizzato nella lezione precedente. Si potrebbe pensare ad una organizzazione di questo tipo:

• Prepariamo una classe astratta Integral che rappresenta un generico algoritmo di integrazione. Il suo costruttore prende in ingresso gli estremi di integrazione e controlla il segno.
• La classe implementa (tra le varie cose) un metodo virtuale puro
  

  virtual double Integra(unsigned int nstep, const FunzioneBase & f ) = 0 ;
  

  che accetti in input il numero di passi (step) con cui si desidera effettuare il calcolo e una referenza (f) alla funzione da integrare.
• Il metodo di integrazione midpoint viene concretamente implementato in una classe derivata Midpoint sovrascrivendo il metodo Integra.

Qui vediamo l'esempio di un possibile header file (.h).

#ifndef __INTEGRAL_H__
#define __INTEGRAL_H__

#include "Funzioni.h"

using namespace std;                                                        

// base class

class Integral {

 public:

  Integral (double a, double b) {
    checkInterval (a,b);
    m_nstep = 0;
    m_h = 0 ;
    m_sum = 0 ;
    m_integral = 0 ;
  } ;

  virtual double Integra(unsigned int nstep, const FunzioneBase &) = 0 ;

 protected:

  void checkInterval( double a, double b ) {
    m_a = min(a,b);
    m_b = max(a,b);
    if ( a > b ) m_sign = -1;
    else m_sign = 1;
  }

  unsigned int m_nstep;
  double m_a, m_b;
  double m_sum, m_integral, m_h;
  int m_sign;

};

// derived class

class Midpoint : public Integral {

 public:

  Midpoint (double a, double b) : Integral(a,b) { ; } ;

  double Integra(unsigned int nstep, const FunzioneBase &f) override {

    if ( nstep <= 0 ) { cout << "Error, number of steps is negative" << endl; exit(-1) ; } ;

    m_nstep = nstep;
    m_h = (m_b-m_a)/m_nstep;

    m_sum = 0.;
    for (unsigned int i=0; i<m_nstep; i++ ) {
      m_sum += f.Eval( m_a  + (i+0.5)*m_h );
    }
    m_integral = m_sign*m_sum*m_h;
    return m_integral;
  };

};

#endif // __INTEGRAL_H__                                   

Per verificare il funzionamento dell'algoritmo di integrazione appena scritto possiamo utilizzare un semplice programma di questo tipo:

#include "Integral.h"
#include "Funzioni.h"

#include <cmath>
#include <iostream>
#include <iomanip>

using namespace std;

int main (int argc, char** argv ) {

  if ( argc!=2 ) {
    cout << "Usage: " << argv[0] << " <nstep>" << endl;
    return -1;
  }

  unsigned int nstep = atoi(argv[1]);

  xsinx f;

  Midpoint myInt(0,M_PI/2.) ;

  double I = myInt.Integra(nstep, f); 

  cout << " Passi= " << setw(20) << nstep << " I= " << setw(20)<< I << endl;

  return 0;
}

Si può poi aggiungere in coda o in un altro codice lo studio dell'andamento dell'errore in funzione del passo $h$ o del numero di step sfruttando un frammento di codice di questo tipo :

  // ............

  TGraph g_errore;

  nstep = 4;
  double Iv = 1. ;
  double I = 0, err = 0 , h = 0 ;

  for ( unsigned int k = 0 ; k < 10; k++ ) {
    h = (b-a)/nstep; 
    I = myInt.Integra(nstep, f); 
    err = fabs(I-Iv);
    cout << " Passi= " << setw(20) << nstep << " Errore= " << setw(20)<< err << endl;
    g_errore.SetPoint(k, h, err );
    nstep *= 2;
  }

  // ............

ESERCIZIO 7.1 - Integrazione con metodo di Simpson a numero di passi fissato (da consegnare):

Implementare l'integrazione con il metodo di Simpson con un numero di passi fissato. Si può utilizzare lo stesso schema dell'esercizio precedente costruendo una classe derivata Simpson . Come nell'esercizio 7.0 fornire il valore dell'integrale calcolato con numero di passi inserito da linea di comando e stampare una tabella (o costruire un TGraph di ROOT) con la precisione raggiunta in funzione del numero di passi.

Il metodo Simpson

Nel metodo di integrazione alla Simpson, la funzione integranda è approssimata, negli intervalli $[x_{k},x_{k+2}]$ (dove $k$ è un intero pari e $x_{k}=a+kh$) con un polinomio di secondo grado i cui nodi sono nei punti $(x_{k},f(x_{k}))$, $(x_{k+1},f(x_{k+1})$ e $(x_{k+2},f(x_{k+2}))$. Esso fornisce una valutazione dell'integrale con una precisione pari a $h^{4}$. La sua applicazione richiede che il numero di passi sia pari e la formula che approssima l'integrale è :

$$\int_{a}^{b} f(x) dx \approx \left( \frac{1}{3}f(x_{0}) + \frac{4}{3} f(x_{1}) + \frac{2}{3} f(x_{2}) + \frac{4}{3} f(x_{3}) + ... + \frac{2}{3} f(x_{N-2}) + \frac{4}{3} f(x_{N-1}) + \frac{1}{3} f(x_{N}) \right) h$$

dove

$$h = \frac{b-a}{N}$$

e

$$x_{k} = a + k h \qquad k = 0, 1.... N$$

ESERCIZIO 7.2 - Integrazione con metodo dei trapezi a precisione fissata (da consegnare):

Concludiamo l'esercitazione implementando l'integrazione della funzione $x sin(x)$ tra $[0,\pi /2]$ con il metodo dei trapezi. In quest'ultimo esercizio proviamo a riflettere sull'uso di un algoritmo di integrazione numerica a precisione fissata invece che a numero di passi fissato. Negli esercizi precedenti il calcolo a numero di passi fissato non ci da alcuna indicazione sulla qualità del risultato : integrare a 10 step è sufficiente ? L'idea di base che vogliamo sviluppare è che l'utente fornisca una precisione desiderata e l'algoritmo debba essere in grado di aumentare automaticamente il numero di passi fino a raggiungere la precisione richiesta sul valore dell'integrale. L'algoritmo dovrà accettare in input il valore della precisione e raddoppiare il numero di passi finchè l'errore (stimato runtime, si veda sotto) non diventa inferiore alla precisione impostata.

1. Come nei casi precedenti si può costruire una classe dedicata per l'implementazione del metodo dei trapezi
2. Possiamo implementare come nei casi precedenti un metodo

   double Integra(unsigned int nstep, const FunzioneBase & f) override {...}
   

   che contenga il calcolo con il metodo dei trapezi a numero di passi fissato.
3. Il punto che ci interessa maggiormente è provare a realizzare un metodo dei trapezi che lavori a precisione fissata, quindi qualcosa del tipo

   double Integra(double prec, const FunzioneBase &f ) override {...}
   

   L'algoritmo dovrà raddoppiare il numero di passi finché l'errore non diventa minore della richiesta dell'utente (prec). La struttura di base sarè quindi un ciclo while nel quale la condizione di uscita è rappresentata dal controllo sull'errore runtime.

Un algoritmo a precisione fissata si può implementare anche per il metodo midpoint e per il metodo di Simpson ma nel caso dei trapezoidi si presta ad una implementazione particolarmente efficiente (si veda la nota sotto).

Il metodo dei trapezi

Il metodo dei trapezi rappresenta la prima formula chiusa di Newton-Cotes e si basa su una interpolazione polinomiale di ordine 1 sui due punti che delimitano l'intervallo. La formula che approssima l'integrale è :

$$\int_{a}^{b} f(x) dx \approx \left( \frac{1}{2}f(x_{0}) + f(x_{1}) + ... + f(x_{N-1}) + \frac{1}{2} f(x_{N}) \right) h$$

dove

$$h = \frac{b-a}{N}$$

e

$$x_{k} = a + k h \qquad k = 0, 1.... N$$

Stima runtime dell'errore

Come possiamo fare a stimare l'errore che stiamo commettendo nel calcolo di un integrale se non conosciamo il valore vero $I$ dell'integrale ? Partendo dalla conoscenza dell'andamento teorico dell'errore in funzione del passo $h$ possiamo trovare un modo semplice per la stima dell'errore che stiamo commettendo. Nel caso della regola dei trapezi l'adamento dell'errore è $err= k h^{2}$. Calcolando l'integrale $I_{N}$ utilizzando un passo $h$ e successivamente l'intgrale $I_{2N}$ con un passo $h/2$ possiamo scrivere il seguente sistema di equazioni: $$I - I_{N} = k h^{2}$$ $$I - I_{2N} = k (h/2)^{2}$$ Sottraendo per esempio la prima equazione alla seconda, non è difficile ricavare che una stima dell'errore pari a $$err(N) = \frac{4}{3} fabs(I_{2N} - I_{N})$$

Cenni sull'implementazione

La condizione di uscita dell'algoritmo è basata sul confronto del risultato di due passaggi successivi. Se la differenza della stima dell'integrale tra due passaggi consecutivi è più piccola della precisione richiesta allora l'algoritmo si ferma. Questo poichè l'integrale calcolato con una partizione più fine è sempre una stima migliore dell'approssimazione con la partizione originaria. L'algoritmo dei trapezoidi si presta ad una implementazione ottimizzata descritta qui di seguito. Nel costruttore calcoliamo la prima approssimazione:

$$\Sigma_{0} = ( f(a) + f(b) ) / 2$$ $$I_{0} = \Sigma_{0} * (b-a)$$

Al primo passaggio dell'algoritmo suddividiamo l'intervallo in due:

$$\Sigma_{1} = \Sigma_{0} + f[x_{11}]$$ $$I_{1} = \Sigma_{1} * (b-a)/2$$

Al secondo passaggio otteniamo :

$$\Sigma_{2} = \Sigma_{1} + f[x_{21}] + f[x_{22}]$$ $$I_{2} = \Sigma_{2} * (b-a)/4$$

Al terzo pasaggio otteniamo :

$$\Sigma_{3} = \Sigma_{2} + f[x_{31}] + f[x_{32}] + f[x_{33}] + f[x_{34}]$$ $$I_{3} = \Sigma_{3} * (b-a)/8$$

e così secondo lo schema in figura ( $L=(b-a)$ ) :

I valori dell'ultima approssimazione dell'integrale e dell'ultima somma calcolata sono memorizzati all'interno dell'oggetto. In questo modo, se viene richiesto di ricalcolare l'integrale, non è necessario ricominciare da capo.

ESERCIZIO 7.3 - Integrazione di una funzione Gaussiana ( facoltativo ):

Come esercizio facoltativo vi proponiamo il calcolo di un integrale non risolubile analiticamente. Assumendo che la distribuzione limite per la misura di una grandezza $x$ sia rappresentata da una gaussiana con valore medio $\mu$ e deviazione standard $\sigma$, proviamo a calcolare la probabilità che una misura di $x$ cada entro $\pm t \sigma$ dal valore medio per $t$ che va da 0 a 5.

1. Provate a costruire una funzione Gaussiana generica. I parametri $\mu$ e $\sigma$ della funzione possono essere passati nel costruttore.
2. Integrate la funzione Gaussiana con il metodo dei Trapezi a precisione fissata in un intervallo $[ \mu - t \sigma, \mu + t \sigma]$ con $t$ variabile da 0 a 5
3. Produrre un grafico dell'integrale ottenuto in funzione del numero di $\sigma$. Vi torna ?

Risultato atteso

Il risultato non dovrebbe sorprendere :