Nel XIX secolo Gauss, studiando la teoria degli errori, si imbattè nella seguente funzione:

    con   0 e costanti reali arbitrarie,

oggi di fondamentale importanza in statistica e in calcolo delle probabilità, nota come funzione gaussiana o funzione normale o, per la sua particolare forma del grafico (figura 1), curva a campana e le cui salienti caratteristiche sono elencate di seguito:

• è definita in tutto l'asse reale;

• è sempre positiva;

• è crescente per x < , decrescente per x >  e assume un massimo per x = ;

• tende a 0 sia per x che tende a - che per x che tende a +;

• è dotata di derivate di ordine comunque elevato, ossia, come dicono i matematici, è "di classe C-infinito".

• il grafico della funzione è simmetrico rispetto alla retta x = ;

• possiede due flessi rispettivamente nei punti  -  e  + .

Figura 1. Grafico di N(x).

Dopo questa breve presentazione della funzione gaussiana, adesso veniamo al nocciolo della domanda. Nella teoria delle funzioni a variabile reale sussiste un importante teorema di esistenza delle funzioni primitive, secondo il quale ogni funzione f(x) continua in un intervallo I, ammette primitiva in I, ossia una funzione F(x), tale che D F(x) = f(x). Poiché N(x) è continua su tutto l'asse reale, ammette sicuramente una primitiva su tutto l'asse reale. Come moltissimi teoremi di esistenza che si trovano in matematica, anche questo ha un grave difetto, non indica in che modo la primitiva può essere calcolata, quindi ci troviamo di fronte alla seguente situazione: sappiamo che esiste una funzione G(x), primitiva di N(x), ma non siamo in grado di calcolarla. In realtà ciò non è del tutto esatto, perchè il modo di calcolare la primitiva in questo teorema di esistenza viene fornito ma la strada che indica non è molto praticabile, infatti considerata una funzione f(x) continua in [a, b], la funzione F(x) :=  è una sua primitiva in [a, b]. Questo è un integrale definito che si calcola usando il famoso teorema di Torricelli, altrimenti detto teorema fondamentale del calcolo integrale: considerata una funzione f(x) continua in [a, b], se F(x) è una sua primitiva, allora . Come è noto, rappresenta l'area sottesa dalla funzione f(x) nell'intervallo [a, b]. Per calcolare la primitiva occorre saper calcolare l'integrale definito e per calcolare quest'ultimo occorre saper calcolare la primitiva, accade un po' come il gatto che si morde la coda. In effetti però la primitiva di una funzione si calcola con altri metodi, ad esempio procedendo con l'integrazione per parti. Comunque eccettuati rari e semplici casi, non siamo mai in grado di calcolare la primitiva di una data funzione f(x), a differenza di come avviene per la derivata; data una funzione siamo quasi sempre in grado di trovare la sua derivata, usando le varie regole di derivazione.

Quadratura numerica

      Nei casi in cui non siamo in grado di calcolare analiticamente la primitiva della funzione f(x), si deve ricorrere al calcolo dell'integrale con tecniche di approssimazione numerica le quali essenzialmente consistono nell'approssimare, nell'intervallo [a, b], la funzione f(x) con un'opportuna famiglia di funzioni {g_n(x) | n  } facilmente integrabili analiticamente in [a, b] e che al crescere di n si "avvicinano" sempre di più a f(x), ossia, detto in modo un po' più formale, comunque scegliamo un numero  > 0, esiste un numero naturale n, per il quale accade che |f(x) - g_n(x)| <  per ogni x  [a, b], come illustrato nella figura 2.

Figura 2. Approssimazione uniforme.

Un esempio di questa categoria di funzioni può essere fornita da una opportuna famiglia di polinomi approssimanti f(x), come per esempio i polinomi interpolanti di Lagrange. Procedendo in questa maniera si vanno a calcolare gli integrali . A seconda delle proprietà di f(x) e di come si sceglie la successione (g_n(x))_n, la successione ()_n può convergere, oppure no, al valore e se converge lo può fare più o meno lentamente.

      Tra le innumerevoli funzioni di cui non siamo in grado di trovare la primitiva in modo analitico vi è anche la funzione N(x); per l'importanza fondamentale che riveste in statistica e nel calcolo delle probabilità, una sua primitiva, che viene indicata con erf(x), è stata calcolata numericamente per molti valori della variabile x e tabulata. Con l'avvento dei personal computer molti software matematici calcolano erf(x) con metodi di approssimazione numerica.

      Più precisamente, erf(x) è la primitiva di N(x) nel caso in cui sia  = 1 e  = 0, ossia nel caso si tratti della cosiddetta funzione normale standard. Da questa primitiva, mediante facili trasformazioni, si può risalire alle primitive di tutte le altre funzioni N(x).

Metodo di Bézout

      Considerata una funzione f(x) continua nell'intervallo [a, b], scelti n + 1 punti equispaziati, a = x₀, x₁, ..., x_n = b, detti nodi e ponendo y_i = f(x_i), l'idea è quella di sostituire alla funzione f(x) all'interno di ogni intervallino [x_i-1, x_i] un polinomio di primo grado, ossia una retta o, per meglio dire, con un segmento avente per estremi i punti di coordinate (x_i-1, y_i-1) e (x_i, y_i), formando così una strana figura che prende il nome di trapezoide, la cui forma è mostrata in figura 3; l'unione di questi n segmenti rappresenta la funzione approssimante g_n(x).

Figura 3. Approssimazione alla Bézout.

Come si vede dalla figura, si vengono a formare n trapezi e l'area del generico trapezio i-esimo è (x_i - x_i-1)(y_i + y_i-1) / 2, quindi l'area del trapezoide è

tenendo presente che (x_i - x_i-1) = (b - a) / n. Indicando con E_n l'errore che si commette approssimando con S_n allora si ha che  = S_n + E_n e, se E_n tende a 0 al crescere di n, accade che S_n tende a . In particolare, se f è dotata di derivata seconda continua allora si dimostra che E_n = f''(c) (b - a) / 12n², dove c è un punto interno all'intervallo [a, b] a noi non noto (se fosse noto sapremmo il valore esatto dell'integrale) e indicando con K₁ la quantità K₁ := f''(c) (b - a) / 12, costante rispetto a n, abbiamo che E_n = K₁ / n². In questo caso quindi E_n si avvicina a 0 all'aumentare del numero del numero dei nodi, e pertanto S_n si avvicina a .

      Ovviamente noi possiamo calcolare solo un numero finito di questi S_n, bisogna trovare un criterio per fermarci a un opportuno indice. Un sistema per fare ciò è quello di fissare a priori un numero opportunamente piccolo  > 0 e calcolare le differenze |S₁ - S₀|, |S₂ - S₁|, ..., |S_k - S_k-1|, ... fino a quando non si giunge a un indice n per il quale risulti |S_n - S_n-1| < ; S_n è il valore che viene scelto per approssimare l'integrale.

Metodo di Cavalieri-Simpson

      Si considerano 2n + 1 punti equispaziati, a = x₀, x₁, ..., x_2n = b di modo che [a, b] risulti suddiviso in n intervalli del tipo [x_2k, x_2k+2] con k = 0, 1, ..., n ossia [a, b] = [x₀, x₂] U [x₂, x₄] U ... U [x_2n-2, x_2n] e in ognuno di questi si sostituisce la funzione f(x) con un polinomio di secondo grado, ossia con una parabola, che indichiamo con P_k, avente l'asse parallelo all'asse y e passante per i punti del piano di coordinate (x_2k, y_2k), (x_2k+1, y_2k+1), (x_2k+2, y_2k+2) e in questo caso la funzione approssimante g_n(x) è l'unione di tutti questi polinomi P_k. Si calcola l'area S'_n, somma di tutte le aree sottese dalle parabole P_k nell'intervallo [x_2k, x_2k+2]. Ciascuna di queste aree vale (y_2k + 4y_2k+1 + y_2k+2)(b - a) / 6 e quindi , per cui  = S'_n + E'_n dove E'_n è l'errore e se esso tende a 0 allora S'_n tende a .

      Se f è dotata di derivata seconda quarta, l'errore che si commette è E'_n = f⁽⁴⁾(c) (b - a)⁵ / 2880n⁴, dove c è un punto interno all'intervallo [a, b] a noi non noto, e indicando con K₂ la quantità K₂ := f⁽⁴⁾(c) (b - a)⁵ / 2880, costante rispetto a n, si ha allora che E'_n = K₂ / n⁴. In questo caso quindi E'_n si avvicina a 0 all'aumentare del numero del numero dei nodi, e pertanto S'_n si avvicina a .

      Abbiamo visto che, in opportune ipotesi sulla funzione f(x), i due procedimenti sono entrambi convergenti, ma confrontando i due errori:

si nota che E'_n tende molto più velocemente a zero rispetto a E_n al crescere di n; quindi se f(x) è una funzione dotata di derivata quarta continua è conveniente usare il procedimento di Cavalieri-Simpson perchè si raggiunge il risultato voluto in un numero minore di passi facendoci risparmiare tempo.

      Come nel caso delle formule di Bézout, si considera la successione S'₁, S'₂, ..., S'_n, ...: se essa converge allora si fissa un  > 0 molto piccolo, si calcolano le differenze |S'₁ - S'₀|, |S'₂ - S'₁|, ..., |S'_k - S'_k-1|, ... e si considera come risultato approssimato di la quantità S'_n per la quale risulta |S'_n - S'_n-1| < .

Finora ci siamo occupati degli integrali propri, ossia integrali i cui estremi a e b siano finiti. Che cosa succede se uno degli estremi è infinito, ossia se si vuole calcolare , dove ad esempio f(x) può essere la funzione gaussiana?

      Le formule usate fino a questo momento non vanno più bene, ma si devono considerare formule alquanto complicate che non è il caso neanche da accennare. Comunque la cosa importante da sottolineare è che in questo caso si sostituisce f(x) con una successione di funzioni più facilmente integrabili.

Impaginato da Gino Favero