Tredicesima lezione

Dopo la discussione fatta sulle equazioni differenziali, potrebbe sembrare una perdita di tempo discutere il calcolo degli integrali. Infatti sappiamo che il valore dell'integrale

è uguale al valore y(b) della soluzione di questa equazione differenziale:

per la quale sia posta la condizione iniziale y(a) = 0.

Di conseguenza, è perfettamente lecito calcolare il valore di un integrale utilizzando i metodi di risoluzione delle equazioni differenziali (il metodo di Runge-Kutta con controllo dinamico dell'incremento, che abbiamo imparato, oppure qualcuno dei metodi più sofisticati che si trovano sui libri o nelle librerie matematiche).

Questo caso relativamente più semplice presenta tuttavia (come sempre) qualche metodo di risoluzione "specializzato". Soprattutto il metodo di integrazione "Montecarlo", che vedremo, illustra una possibile applicazione di una delle tecniche più semplici e potenti dell'elaborazione numerica, sulla quale vale la pena di spendere qualche parola in più.
Ma iniziamo con ordine.

Il primo metodo che può venire in mente, per calcolare l'integrale

è sperare che la funzione non vari troppo nell'intervallo [a,b], e considerarla costante, con un valore pari a quello nel punto medio. In questo caso la prima approssimazione per il valore dell'integrale è:

L'errore che si commette possiamo valutarlo considerando lo sviluppo in serie della funzione attorno al punto medio (usiamo per maggior chiarezza il resto di Lagrange):

Si può fare di meglio, dividendo l'intervallo in N sottointervalli, ed utilizzando i rispettivi punti medi. In questo modo il valore dell'integrale si approssima con l'area del trapezoide risultante:

Ecco il valore dell'integrale calcolato con questa procedura:

Volendo raffinare ulteriormente le cose, possiamo calcolare il valore dell'integrale ad ordini di accuratezza maggiori; potremmo ancora (ma non lo facciamo) combinare i risultati con opportuni coefficienti che annullino esattamente i corrispondenti termini dello sviluppo in serie. E' possibile anche estrapolare il risultato per un numero di elementi del trapezoide tendente all'infinito, utilizzando (come avevamo mostrato nel caso delle derivate) un metodo di estrapolazione polinomiale. Per l'algoritmo risultante (integrazione di Romberg) ecco di nuovo una implementazione dell'algoritmo da Numerical Recipes in C/C++: romberg.cpp.

Come abbiamo visto, c'è un vasto ventaglio di possibilità, ed una lunghissima storia, riguardo all'integrazione delle funzioni di una variabile. Di per sé non avremmo dunque bisogno di sapere altro, se non fosse per il fatto che il problema può diventare molto più complicato, per esempio quando le dimensioni aumentano. Uno dei possibili metodi (probabilmente il più semplice) per calcolare numericamente integrali in più dimensioni, è il metodo Montecarlo. Vediamo brevemente di cosa si tratta nel caso ad una dimensione:

Supponiamo di individuare un rettangolo che contiene la funzione da integrare nell'intervallo scelto (attenzione a non "tagliare" la funzione: occorre saper valutare bene massimi e minimi). A questo punto scegliamo un certo numero di punti a caso all'interno del rettangolo, e contiamo quanti si trovano sopra, e quanti sotto la funzione. L'integrale cercato è dato dall'area del rettangolo (che conosciamo) moltiplicata per la frazione di punti contati sotto la funzione. Se i numeri casuali che utilizziamo sono casuali per davvero, possiamo andare avanti ad estrarre numeri fino a misurare l'integrale alla precisione della macchina.

Vediamo in questo come calcolare il volume di una calotta sferica con il metodo Montecarlo. L'indeterminazione che grava sulla misura del volume è regolata dalla distribuzione binomiale. La distribuzione binomiale infatti descrive la probabilità di avere n "successi" in N "prove" quando si sollecita un processo caratterizzato dalla probabilità p per un singolo successo:

Non lo dimostriamo, ma la migliore stima del valore di p ottenuta, come nel nostro caso, dal conteggio del numero di successi e di prove totali (la misura del volume è infatti p moltiplicata per il volume totale) è:

e la corrispondente indeterminazione è:

Dal programma di esempio risulta abbastanza immediato come la semplicità di impostazione del metodo si paghi con un costo computazionale piuttosto rilevante. Ovviamente, infatti, sappiamo come calcolare meglio il volume di un solido di rotazione come la calotta sferica.

Esistono tuttavia casi per i quali la complessità del fenomeno o della struttura in esame lascia come unica via d'uscita l'utilizzo di questo sistema che permette di ottenere la precisione desiderata, a patto di possedere una potenza di calcolo (o una quantità di pazienza) sufficiente.

La strategia di eseguire calcoli utilizzando valori casuali si applica ovviamente bene alla simulazione di qualunque fenomeno che possa essere predetto solo in maniera statistica. Conoscendo la distribuzione di probabilità di un certo fenomeno, è possibile "pescare" numeri casuali che seguano la distribuzione stessa, usando un metodo molto semplice e simile a quello utilizzato per l'integrazione.
Questo mostra come si possono pescare numeri che seguono una distribuzione gaussiana centrata sullo zero:

Questo metodo di generare numeri casuali che seguono una distribuzione si chiama metodo di reiezione, e si dimostra tanto più inefficiente quanto più è grande l'area utilizzata per racchiudere la distribuzione di probabilità che ci interessa (le scelte di numeri "inutili" avverranno con una probabilità proporzionale alla differenza fra tale area e quella sottesa dalla distribuzione).

Vale la pena di notare, per completezza, che esiste un altro metodo (analitico) di generare numeri casuali che seguono una distribuzione, a partire da una distribuzione di numeri casuali "piatta". Tale metodo richiede però di sapere come trovare una primitiva della distribuzione desiderata e come invertirla. Con questo metodo (detto metodo di trasformazione) si possono generare distribuzioni esponenziali, ed anche gaussiane (utilizzando uno stratagemma in due dimensioni). E' quindi quest'ultimo metodo che viene implementato nelle librerie matematiche.

Il metodo Montecarlo si può applicare anche quando non si conosce la distribuzione di probabilità del fenomeno che si desidera simulare. E' possibile infatti operare una scelta casuale fra un campione ben costruito di precedenti incarnazioni del fenomeno (per esempio provenienti da misure sperimentali).

Ma perché porre tanta importanza a metodi di simulazione su computer? La fisica non si basa su misure sperimentali?

Vero, ma nel caso di apparati sperimentali particolarmente complessi, non è possibile "eliminare" nell'apparato le fonti di errore sistematico, ed occorre in qualche modo fattorizzare la risposta dell'apparato e valutare la dipendenza della misura dalle scelte fatte durante la misura stessa. Questa è la codiddetta "misura dell'errore sistematico", e la simulazione al calcolatore dell'apparato sperimentale, fatta con metodi Montecarlo, riveste in questo caso importanza fondamentale.

Esercizio: provare a convertire in un metodo di uso generale la procedura che estrae n numeri casuali "secondo una distribuzione". La definizione della classe corrispondente potrebbe essere:

#include <functional>

template <typename DT>
class DistributedNumbers
{
  public:
    typedef DT value_type;
    typedef DT (*function_ptr_type)(DT); // Normal - nonmember function pointer
    typedef std::function<DT (DT arg)> functor_type;

    // Creator
    DistributedNumbers(functor_type &f);
    DistributedNumbers(function_ptr_type f);

    // Destructor
    ~DistributedNumbers();

    DT get_next_number();
}