Metodi numerici di integrazione¶

Al termine del capitolo precedente siamo incappati in un integrale impossibile da risolvere con le tecniche di calcolo alla nostra portata e per questo ci siamo avvalsi del risultato fornito da un programma adatto. Possiamo evitare l’uso di un software specifico. Con un linguaggio di programmazione possiamo implementare algoritmi che simulano l’attività di questi software. Otterremo risultati approssimati, interessanti per i metodi che utilizzati, più che per i valori in sè. Come già si è visto per il calcolo differenziale, questi metodi si dicono metodi numerici.

Il metodo dei rettangoli¶

Il primo metodo che studiamo è il più semplice, ma anche il meno preciso. L’intervallo di integrazione viene suddiviso in un numero finito di parti, che quindi sono piccole ma non infinitesime.

Data la funzione continua $f:[a,b] \to \mathbf{R}$ , costruiamo una successione di $n$ punti in modo che sia $a=x_0<x_1<x_2<...<x_n=b$ . In ogni sottointervallo $[x_k,x_{k+1}]$ prendiamo un numero $\overline{x}_k$ e approssimiamo l’integrale $\int_a^bf(x)dx$ con la somma $\sum_0^{n-1}f(\overline{x}_k)\Delta x_k$ .

Si potrebbe, per semplificare, distanziare ugualmente i punti. Allora, per ogni $\Delta x_k$ abbiamo $\Delta x_k=h=\frac{b-a}{n}$ e $x_k=a+kh$ , con $k=0 ... n$ . Se prendiamo come altezza di ogni rettangolo il valore della funzione nell’estremo destro della base, $\overline{x}_k=x_{k+1},\ f(x_k)=y_k.\ hy_k$ è l’area di ogni rettangolo e la somma che approssima l’integrale diventa:

$\sum_0^{n-1}f(\overline{x}_k)\Delta x_k= \sum_0^{n-1}f(x_{k+1})h= h\sum_0^{n-1}y_{k+1}=h\sum_1^ny_k.$

Ne risulta quindi l’uguaglianza approssimata:

$\int_a^bf(x)dx\cong h\sum_1^ny_k.$

Questa è la formula alla quale facciamo riferimento per costruire l’algoritmo che svolgerà i primi calcoli.

inizio

leggi 
;
per  esegui
inizio
;
;
fine
scrivi 

fine

Due esempi¶

$\int_0^3x^2dx$ .

Calcoliamo con questo metodo un integrale facile: $S=\int_0^3x^2dx$ , il cui valore esatto è $9$ .

$&n=10,\quad S=10,395\\ &n=100,\quad S=9,13545\\ &n=1000,\quad S=9,0135045\\ &n=10000,\quad S=9,00135005$

Man mano che cresce il numero degli intervalli, e l’algoritmo prevede un maggior numero di cicli, il risultato approssima sempre meglio quello esatto.

$\int_0^x \sqrt{1+cos^2 x}\, dx$ .

L’integrale impossibile del capitolo precedente, cioè $\int_0^x \sqrt{1+cos^2 x}\, dx$ viene calcolato con la stessa precisione prevedendo $1000$ cicli.

Il metodo dei trapezi¶

Come indica il nome, questo metodo si basa sulla sostituzione dei rettangoli con i trapezi. Il grafico della funzione viene approssimato da una spezzata, che è il grafico di una funzione approssimante $\phi$ , una funzione lineare a tratti. In pratica il valore di $\int_a^bf(x)dx$ si approssima con il valore esatto di $\int_a^b\phi(x)dx$ .

$\int_a^bf(x)dx\cong\int_a^b\phi(x)dx.$

Manteniamo la stessa suddivisione in intervalli adottata con il metodo dei rettangoli. Data $y_k~=~f(x_k)$ , definiamo la funzione $\phi$ come la funzione che ha per grafico la spezzata che unisce i punti $P_k(x_k,y_k),\ k=0 ... n$ .

L’area di ogni trapezio è data da $\frac{y_{k-1}+y_k}{2}\Delta x_k$ , quindi

$\frac{y_{k-1}+y_k}{2}\Delta x_k=&\int_{x_{k-1}}^{x_k}\phi(x)dx.\\ \int_a^bf(x)dx \cong \int_a^b\phi(x)dx=&\sum_1^n\int_{x_{k-1}}^{x_k}\phi(x)dx.$

Possiamo quindi scrivere la formula approssimata

$\int_a^bf(x)dx\cong \frac{1}{2}\sum_1^n (y_{k-1}+y_k)\Delta x_k.$

Nella suddivisione in progressione aritmetica, ogni sottointervallo ha ampiezza costante: $\Delta x_k=h$ e quindi si ha:

$\int_a^bf(x)dx\cong &\frac{h}{2}\sum_1^n (y_{k-1}+y_k)=\\ =&\frac{h}{2}[(y_0+y_1)+(y_1+y_2)+(y_2+y_3)+ ... +(y_{n-1}+y_n)]=\\ =&\frac{h}{2}(y_0+2y_1+2y_2+2y_3+ ... +2y_{n-1}+y_n)=\\ =&\frac{h}{2}\left[y_0+2\sum_1^{n-1}y_k+y_n \right].$

L’algoritmo allora è:

inizio

leggi $(a,b,n)$ ;

$h\leftarrow\frac{b-a}{n}$ ;

per $k\leftarrow 0\ ..\ n$ esegui

$y_k\leftarrow f(a+kh)$ ;

$S \leftarrow\ 0$ ;

per $k\leftarrow 1\ ..\ n-1$ esegui

$S\leftarrow S+y_k$ ;

$S\leftarrow y_0+2S+y_n$ ;

$I \leftarrow \frac{h}{2}S$ ;

scrivi $(I)$ ;

fine

Se si testa il funzionamento dell’algoritmo sulla funzione $y=x^2$ si scopre che per $n=10$ , cioè ricoprendo la superficie con $10$ trapezi, si ottiene quasi la stessa precisione ottenuta con $1000$ rettangoli.

$&n=10,\quad S=9,045\\ &n=100,\quad S=9,00045\\ &n=1000,\quad S=9,0000045\\$

Il segno dell’area¶

Osservando il disegno della pagina precedente, si può notare come il ricoprimento a trapezi sia più efficace di quello fatto a rettangoli.

C’è un dettaglio da discutere nel caso che il tratto obliquo intersechi l’asse $x$ : il segno dell’area risulta positivo o negativo?

Supponiamo di avere $y_{k-1}$ negativo e in valore assoluto minore rispetto a $y_k$ positivo: $y_{k-1}~<~0,\ y_{k}>0,\ | y_{k-1}|<y_k$ . Al numeratore della formula dell’area del trapezio c’è la somma $y_{k-1}+y_k$ , che vale in questo caso $y_{k}- | y_{k-1}|$ (ma gli altri casi sono analoghi). Fissiamo il punto $Q(x_k,y_{k-1}+y_k$ ), quindi di ordinata corrispondente alla somma delle basi del trapezio. Uniamo $Q$ con $A(x_{k-1},0)$ . Si formano così 4 triangoli $a,c$ uguali a coppie, e un quadrilatero $b$ . Nel calcolare l’area del trapezio, i triangoli $c$ si annullano e la parte restante $a+b$ , sottesa al segmento, equivale al triangolo $a+b$ , con ipotenusa $AQ$ . Se $|y_k|<|y_{k-1}|$ lo stesso disegno si ripropone nel semipiano negativo, e anche i casi in cui $y_{k-1}>0,\ y_{k}<0$ sono analoghi. In conclusione, la formula $\frac{y_{k-1}+y_k}{2}\Delta x_k$ produce l’area effettiva del trapezio, tenendo conto correttamente di tutte le parti, negative e positive.

La formula di Simpson¶

Con il ricoprimento a trapezi la funzione è stata trasformata in una funzione a dominio discreto e per ogni coppia di suoi punti al posto della funzione originale abbiamo usato una funzione lineare. Si è trattato quindi di una interpolazione lineare.

Procediamo ora sulla stessa falsariga, isolando non coppie ma terne di punti consecutivi e considerando, al posto della funzione, l’arco di parabola per quei punti. Il grafico che passa per tre punti è infatti espressione di un polinomio al massimo di secondo grado (al massimo di primo grado se i tre punti sono allineati, di grado zero se sono allineati in orizzontale). Se dobbiamo considerare gruppi di tre punti, gli intervallini saranno necessariamente in numero pari, quindi $n=2m$ . Limitiamoci anche questa volta a intervallini di uguale ampiezza $h=\frac{b-a}{2m}$ .

Procediamo con i primi tre punti $x_0,x_1,x_2$ e poi cercheremo di generalizzare il risultato. Il fatto importante, però, non è tanto esprimere il polinomio, quanto il suo integrale in $[x_0,x_2]$ . Stiamo infatti cercando di sostituire la funzione $f$ con una funzione $\psi$ che coincide con l’unico polinomio passante per i tre punti consecutivi, che isolano una coppia di intervallini. E deve valere che l’integrale approssimato di $f$ è l’integrale esatto di $\psi$ , il quale a sua volta è la somma di una coppia di integrali:

$\int_a^b\psi(x)dx=\sum_1^m\int_{x_{2k-2}}^{x_{2k}}\psi(x)dx$

Il primo di questi integrali è

$\int_{x_0}^{x_2}\psi(x)dx.$

La funzione $\psi(x)$ , come abbiamo detto è un polinomio al massimo di secondo grado: $P(x)~=~Ax^2+Bx+C$ , del quale sappiamo che $P(x_0)=y_0,\ P(x_1)=y_1,\ P(x_2)=y_2$ , con $A,B,C$ da ricavare.

Per facilitare la ricerca di questi tre coefficienti (e poiché, in fondo, ci interessa più l’integrale del polinomio), ricordiamoci che gli intervallini hanno uguale ampiezza $h$ : potremo allora traslare i due intervallini: $x_0=-h,\ x_1=0,\ x_2=h$ , mantenendo, ovviamente, uguali valori per gli $y$ corrispondenti. Il vantaggio di questa trasformazione è che ora dobbiamo calcolare

$\int_{x_0}^{x_2}\psi(x)dx=\int_{-h}^{h}P(x)dx=\int_{-h}^{h}(Ax^2+Bx+C)dx= \int_{-h}^{h}Ax^2dx+\int_{-h}^{h}Bxdx+\int_{-h}^{h}Cdx.$

Poiché ora abbiamo l’intervallo di integrazione simmetrico, le funzioni dispari hanno integrale nullo mentre le funzioni pari hanno integrale doppio rispetto all’integrale su metà intervallo:

$\int_{-h}^{h}(Ax^2+Bx+C)dx= 2 \int_0^{h}(Ax^2+C)dx =2\left[\frac{Ah^3}{3}+Cx\right]_0^h=\frac{2}{3}Ah^3+2Ch.$

Abbiamo così scoperto che il calcolo di $B$ non è importante. Poi, da $P(0)=C$ segue $C=y_1$ . Per il calcolo di $A$ :

$&P(-h)=Ah^2-Bh+C=y_0\\ &P(h)=Ah^2+Bh+C=y_2\\ &P(-h)+P(h)=2Ah^2+2C=y_0+y_2\\ &2Ah^2=y_0-2y_1+y_2\\ &\int_{-h}^{h}(Ax^2+Bx+C)dx=\frac{2}{3}Ah^3+2Ch=\frac{h}{3}\left[2Ah^2+6C \right]= \frac{h}{3}(y_0+4y_1+y_2).$

Abbiamo quindi trovato l’integrale dell’arco di parabola per i primi tre punti, in funzione dell’ampiezza $h$ degli intervallini e delle ordinate dei tre punti. La cosa si estende facilmente a tutte le terne di punti successive, così:

$\int_a^bf(x)dx\cong &\frac{h}{3}(y_0+4y_1+y_2)+\frac{h}{3}(y_2+4y_3+y_4)+ ...=\\ =&\frac{h}{3}(y_0+4y_1+2y_2+4y_3+2y_4+ ... +2y_{n-2}+4y_{n-1}+y_n).$

Tranne che per il primo e per l’ultimo termine, sui termini di indice pari si raddoppia la somma, mentre sui termini di indice dispari si quadruplica. Ricordando che $n=2m$ , riscriviamo e sintetizziamo:

$\int_a^bf(x)dx\cong &\frac{h}{3}(y_0+4y_1+2y_2+ ... +2y_{2m-2}+4y_{2m-1}+y_{2m})=\\ &=\frac{h}{3}\left[y_0+4\sum_1^m y_{2k-1}+2\sum_1^{m-1}y_{2k}+y_{2m} \right].$

Questa formula è detta di Simpson, o di integrazione parabolica. L’algoritmo allora è

inizio

leggi $(a,b,m)$ ;

$h\leftarrow\frac{b-a}{2m}$ ;

per $k\leftarrow 0\ ..\ 2m$ esegui

$y_k\leftarrow f(a+kh)$ ;

$S_1 \leftarrow\ 0$ ;

per $k\leftarrow 1\ ..\ m$ esegui

$S_1\leftarrow S_1+y_{2k-1}$ ;

$S_2 \leftarrow\ 0$ ;

per $k\leftarrow 1\ ..\ m-1$ esegui

$S_2\leftarrow S_2+y_{2k}$ ;

$S\leftarrow y_0+4S_1+2S_2+y_{2m}$ ;

$I \leftarrow \frac{h}{3}S$ ;

scrivi $(I)$ ;

fine

La precisione della formula di Simpson¶

Calcolare con questo algoritmo l’integrale di una funzione di secondo grado non ha molto senso, perché l’interpolazione della parabola, fatta con una parabola, non può che dare il risultato esatto.

È meno ovvio che risulti esatto anche l’integrale di una cubica. Infatti, se si impiega l’algoritmo per calcolare $\int_0^2 x^3 dx$ anche solo con due intervalli, cioè $m=1$ , si ottiene il risultato esatto, cioè $4$ .

Vediamolo in dettaglio. Con $m=1$ ci sono due intervalli e i punti sono $x_0=0, x_1=1, x_2~=~2$ . I valori corrispondenti di $y$ sono: $y_0=0, y_1=1, y_2=8$ . Poiché $h=1$ , la formula di Simpson dà:

$\frac{1}{3} (0+4+8)=4.$

Il polinomio interpolatore $P(x)=Ax^2+Bx+C$ , al massimo di secondo grado, è caratterizzato da $P(0)=0, P(1)=1, P(2)=2$ . Date queste condizioni, abbiamo

$&P(0)=0\ \to C=0,\ P(x)=Ax^2+Bx\\ &P(1)=1\ \to A+B=1,\quad P(2)=8\ \to 4A+2B=8 \to 2A+B=4\\ & A=3, B=1-A=-2\\ &P(x)=3x^2-2x.$

I grafici della funzione e del polinomio sono ovviamente diversi, ma gli integrali coincidono:

$\int_0^2(3x^2-2x)dx=[x^3-x^2]_0^2=8-4=4.$

La cosa notevole è che la formula di Simpson fornisce il risultato esatto per tutti i polinomi di terzo grado. Infatti in questi polinomi l’unica parte da discutere è il termine in $x^3$ , perchè i termini di grado inferiore coincidono con la loro interpolazione, come abbiamo già notato. Discutiamo quindi solo di

$\int_a^b x^3dx=\left[\frac{x^4}{4} \right]_a^b=\frac{b^4-a^4}{4}.$

Con la formula di Simpson, con $m=1$ , abbiamo:

$& x_0=a,\quad x_1=\frac{a+b}{2},\quad x_2=b.\\ &y_0=a^3,\quad y_1=\frac{(a+b)^3}{8},\quad y_2=b^3.\\ &h=\frac{a+b}{2}.\\ &\frac{h}{3}(y_0+4y_1+y_2)=\frac{b-a}{6}\left[ a^3+\frac{(a+b)^3}{2}+b^3\right]= \frac{b^4-a^4}{4}$

come si verifica facilmente sviluppando la parentesi quadrata. Il risultato è anche indipendente dal numero di intervalli, perché se vale per una coppia di intervalli, vale anche per ogni altro numero di coppie. La ragione di questa esattezza risiede nel fatto, dimostrabile, che l’imprecisione della formula di Simpson dipende dai valori della derivata quarta della funzione, che vale zero nel nostro polinomio.

I tre metodi a confronto¶

Vale la pena di testare i tre metodi descritti, per calcolare un integrale di cui non conosciamo la primitiva: $\int_0^2e^{-x^2}dx$ .

Ecco la tabella che registra i risultati di prove diverse.

intervalli	rettangoli	trapezi	Simpson
10	0,71460477	0,7462108	0,74682495
20	...	...	0,74682418
40	...	...	0,74682414
1000	0,74650801	0,74682407	...
10000	...	0,74682413	...
1000000	0,74682382	...	...

Il risultato che si ottiene con Derive, un software dedicato, è 0,74682413. Derive usa un numero di suddivisioni dell’ordine del miliardo e per questo posiamo prendere il suo risultato come il più preciso. Si può notare che il metodo di Simpson raggiunge l’esattezza fino alla settima cifra decimale con $40$ intervalli ( $m=20$ ), una precisione uguagliata dal metodo dei trapezi con $10000$ intervalli.

Il metodo Montecarlo¶

Il metodo Montecarlo ha un’origine e un’impostazione del tutto diverse dai precedenti. Proviene infatti da studi di statistica e il suo nome ricorda il gioco d’azzardo e i numeri che escono a caso. Applichiamo i concetti di questo metodo al calcolo degli integrali.

Abbiamo una funzione $f$ continua e positiva su $[a,b]$ e conosciamo anche un maggiorante della funzione, cioè un numero $M$ tale che $f(x)\le M$ per ogni $x\in [a,b]$ . Il trapezoide della funzione è certamente contenuto nel rettangolo di base l’intervallo $[a,b]$ e altezza $M$ .

Per capirci, possiamo immaginare che il trapezoide sia la vasca di una piscina e la parte restante del rettangolo sia il bordo piastrellato della vasca. Gettando casualmente sassi sul rettangolo, la probabilità che cadano nell’acqua è legata all’area del trapezoide. Per esempio, se $85$ % dei sassi cade in acqua, diremo che l’area della vasca è $85$ % dell’area del rettangolo. Quindi, se cade in acqua la frazione $m/n$ dei lanci, questa è la nostra stima dell’area del trapezoide:

$\int_a^bf(x)dx=\frac{m}{n}(b-a)M.$

Ovviamente, dovremo avere lanciato in modo omogeneo, senza concentrare i lanci su un’area piuttosto che su un’altra. I lanci devono essere quindi casuali e distribuiti in modo omogeneo su tutto il rettangolo. Quindi le coordinate del punto in cui arriva ogni lancio devono essere calcolate da una funzione opportuna, che fornisce numeri uniformemente distribuiti e indipendenti, come si richiede a numeri correttamente casuali. Uniformemente distribuiti significa che in ogni porzione dell’intervallo cade più o meno la stessa quantità di numeri. Indipendenti significa che ogni numero non è legato al numero generato in precedenza. Se queste caratteristiche sono rispettate, allora il caso viene simulato abbastanza bene. Esiste un algoritmo apposito, implementato nei software, detto generatore di numeri pseudocasuali, che garantisce un buon rispetto di queste regole. I numeri di questo generatore sono compresi fra $0$ e $1$ e si ripetono all’incirca dopo milioni di volte, in modo che, generato un numero, è impossibile prevedere il successivo. La funzione dei sofware che implementa la generazione di questi numeri è detta random e noi la richiameremo con il comando rnd.

Dunque, dobbiamo generare casualmente punti nel rettangolo $[a,b]\times[0,m]$ , con il generatore di numeri pseudocasuali rnd, che fornisce numeri fra $0$ e $1$ . Avremo quindi coppie di numeri da trasformare, in modo che risultino appartenere all’intervallo voluto. La funzione che trasforma questi numeri è una semplice retta. Infatti, perché un numero casuale $\mathit{rnd} \in [0,1]$ sia trasformato in un numero dell’intervallo $x\in[a,b]$ , basta che sia

$X=(b-a)\mathit{rnd}+a.$

Analogamente, dovrà essere $y=(M-0)\mathit{rnd}\to y=M\cdot \mathit{rnd}$ .

Per costruire l’algoritmo che conta quanti punti casuali appartengono all’area richiesta abbiamo bisogno di controllare se la coordinata $Y$ dei punti è nel trapezoide, cioè se ha valore inferiore a quello corrispondente della funzione: $Y<f(x)$ .

L’algoritmo riceve in ingresso il numero $n$ di lanci, cioè di valori casuali da generare, e gli intervalli entro cui trasformarli. Poi incrementa un contatore $m$ ogni volta che il numero $y$ è inferiore a $f(x)$ e infine produce la frazione $m/n$ che, moltiplicata per l’area del rettangolo, fornisce la stima dell’integrale.

inizio

leggi ;
per  esegui
inizio
;
;
se  allora ;
fine;
;
scrivi ;

fine

Un test per il metodo Montecarlo¶

Proviamo il metodo sul calcolo di $\int_0^3 x^2dx$ e confrontiamolo con i risultati dei metodi precedenti e con il risultato esatto che conosciamo, cioè $9$ . Come maggiorante useremo $M=9$ .

lanci	Montecarlo
10	8,1
100	10,8
1000	8,937
10000	9,0315
100000	8,99532
1000000	8,967537

Osserviamo che il metodo è impreciso e non è detto che aumentare il numero dei lanci migliori la precisione. Si tratta di un metodo statistico, cioè che produce risultati in media vicini al valore esatto. Infatti lo stesso gruppo di lanci non dà, in genere, lo stesso risultato. Bisogna quindi ripetere i gruppi di lanci e fare la media. Per esempio si trova che ripetendo dieci volte $100$ lanci, il risultato medio è $9,1098$ . In conclusione, il metodo Montecarlo non è da usare per un integrale definito, cosa che d’altronde ci era nota anche in relazione agli altri metodi. In effetti il metodo Montecarlo eccelle in altre situazioni.

Possiamo provare ad applicarlo all’integrale $\int_0^1 e^{-x^2}dx$ Con $M=1$ e $100000$ lanci un risultato è $0,73678$ , che è un’ottima approssimazione del risultato a $5$ decimali ottenuto con Derive.

Un problema di stima¶

Due cilindri si intersecano perpendicolarmente. Occorre una stima del volume comune ai due cilindri.

Il volume cercato è proporzionale a $r^3$ , con $r$ raggio dei cilindri. Per semplificare poniamo $r=1$ . Immaginiamo che gli assi dei due cilindri coincidano con gli assi cartesiani $x,y$ . Lungo l’asse $x$ si sviluppa il volume interno al cilindro orizzontale, caratterizzato dalla disequazione $y^2+z^2<1$ . Per il cilindro verticale c’è una disequazione analoga: $x^2+z^2<1$ . La regione comune ai due cilindri è contenuta in un volume cubico di lato $2$ e quindi volume $8$ .

Ora fingiamo che il volume cubico sia un panettone e di volerci conficcare dei canditi a caso. Dopo avere inserito $n$ canditi in posizione casuale, rileveremo quanti di essi sono nel volume comune ai due cilindri. La frazione $\frac{m}{n}8$ sarà la stima richiesta.

inizio

leggi ;
per  esegui
inizio
;
;
;
se  AND  allora ;
fine;
;
scrivi ;

fine

Il controllo sulla posizione di ogni candito avviene attraverso la verifica delle due disequazioni. L’intervallo entro cui far cadere i numeri pseudocasuali è $[-1,1]$ e quindi ogni numero si ottiene da $(b-a)rnd-1=2rnd-1$ .

L’algoritmo genera un risultato che si assesta attorno a $V=5,33 r^3$ , con un milione di canditi. Il metodo Montecarlo ci consente di stimare rapidamente un integrale non facile, senza impegnarci in calcoli: è questo il suo vantaggio principale.

La soluzione esatta¶

Il problema precedente ha una soluzione esatta, che proviamo a calcolare. Visto da $z$ , l’incrocio dei due cilindri avviene lungo segmenti inclinati di $45^\circ$ rispetto agli assi dei cilindri (in realtà, dall’esterno non si vedono segmenti, ma curve).

La parte comune ai due cilindri può essere sezionata con piani perpendicolari all’asse $z$ , formando scaglie quadrate. A distanza $z$ dall’origine, una scaglia ha il lato che misura $2\sqrt{1-z^2}$ e ha quindi l’area di $4(1-z^2)$ . Ogni scaglia ha spessore infinitesimo $dz$ ed è indistinguibile da un prisma a base quadrata, per cui l’integrale esatto si calcola così:

$4\int_{-1}^1(1-z^2)dz=8\int_0^1(1-z^2)dz=8\left[z-\frac{z^3}{3} \right]_0^1 =\frac{16}{3}.$

Il risultato corrisponde pienamente a quanto calcolato con il metodo Montecarlo e ne conferma la potenza.

La densità dei numeri casuali¶

Negli esempi precedenti ci siamo avvalsi del generatore di numeri pseudocasuali rnd, implementato nei linguaggi di programmazione, e abbiamo visto che è sufficientemente affidabile per i nostri scopi. Abbiamo visto anche che i numeri rnd, generati nell’intervallo $[0,1]$ , possono essere trasformati in altri numeri $y$ nell’intervallo desiderato $[a,b]$ , grazie alla funzione di trasferimento $y=(b-a)rnd+a$ , il cui grafico è una retta. L’intervallo $[a,b]$ viene allora riempito di numeri casuali in modo omogeneo, cioè senza che si addensino preferibilmente in una zona dell’intervallo.

Tuttavia proprio questa condizione di uniformità, in genere desiderata, potrebbe in qualche caso non essere preferibile. Si vorrebbe quindi poter guidare la distribuzione di questi numeri nell’intervallo $[a,b]$ , in modo che alcune zone dell’intervallo siano preferite rispetto ad altre e quindi vi sia una distribuzione disomogenea, più concentrata in alcune zone rispetto ad altre.

Per esempio, con una funzione come quella in grafico nel disegno, i numeri rnd vengono generati uniformemente nell’intervallo $[0,1]$ , ma si concentrano agli estremi dell’intervallo di arrivo. La loro concentrazione è guidata dall’andamento della curva, e quindi dalla derivata della funzione in grafico.

In generale, ogni punto $x$ dell’intervallo di arrivo $[a,b]$ è legato alla probabilità che siano generati numeri $rnd\le x$ . Per esempio la probabilità legata a $b$ è la certezza, cioè $1$ , mentre la probabilità legata a $a$ è così piccola da poter essere considerata nulla. Chiamiamo $F$ la funzione che rappresenta questa probabilità. Si tratta di una funzione non decrescente. $F$ ha caratteristiche analoghe alla massa di una barra, nulla all’inizio e pari a quella di tutta la barra, all’estremo opposto. Lavoreremo su funzioni di probabilità crescenti e continue, per comodità di ragionamento.

Nel caso della barra avevamo definita una densità (di massa). Densità media, che è costante per la barra omogenea, e densità puntuale, quando la concentrazione della massa varia da punto a punto. Lo stesso possiamo dire della densità di probabilità. La densità di probabilità media in un tratto è il rapporto fra la probabilità che in quel tratto cadano i valori (generati casualmente e trasformati) e la lunghezza del tratto. Per un tratto $[x,x+\Delta x]$ , sarà $F(x+\Delta x)-F(x)~=~\Delta F(x)$ e la densità media sarà: $\frac{\Delta F(x)}{\Delta x}$ .

La densità di probabilità puntuale sarà quindi

$f(x)=F'(x)$

ed ha la caratteristica $f(x)>0$ . Per il fatto di determinare, grazie al suo andamento, la distribuzione dei valori casuali nell’intervallo $[a,b]$ , $F$ è chiamata funzione di ripartizione.

Possiamo anche ragionare al contrario, cioè dire che la funzione di ripartizione relativa al punto $x\in [a,b]$ si ottiene da

$F(x)=\int_a^xf(t)dt$

con $F(a)=0, F(b)=1$ , e quindi

$\int_a^bf(x)dx=1.$

Una funzione $f$ definita su un intervallo $[a,b]$ continua, positiva e con l’integrale sull’intervallo che vale $1$ , è una densità di probabilità. La probabilità che la variabile casuale assuma valori compresi nell’intervallo $[x_1,x_2]$ può essere data in due modi equivalenti: o integrando la densità di probabilità oppure calcolando la differenza fra i valori della funzione di ripartizione. Vale, ancora una volta, il teorema fondamentale del calcolo integrale:

$\int_{x_1}^{x_2}f(x)dx=F(x_2)-F(x_1).$

La funzione di ripartizione $F$ ha per dominio l’intervallo $[a,b]$ e le sue differenze $\Delta F(x)=F(x+\Delta x)-F(x)$ sono proporzionali al numero di valori che cadono nell’intervallo $[x, x+\Delta x]$ , cioè sono proporzionali a $\int_x^{x+\Delta x}f(t)dt$ , la probabilità che la variabile casuale assuma valori in quell’intervallo.

Il procedimento che determina quale sia la funzione di trasformazione $F$ necessaria a determinare una certa distribuzione di valori in $[a,b]$ è il seguente:

Si parte da una funzione f che esprime la densità. Essa deve essere: $f:[a,b]\to \mathbf{R}$ , continua, positiva e con $\int_a^bf(x)dx=1$ . Se manca quest’ultima condizione, per esempio in una diversa funzione $\phi(x)$ , si può operare in modo da ottenere $f$ in questo modo: $f(x)=\frac{\phi(x)}{\int_a^b\phi(x)dx}$ .
Si considera la funzione di ripartizione $F(x)=\int_a^xf(t)dt$ che risulta automaticamente continua, invertibile e con $F(a)=0, F(b)=1$ .
Si ricava la sua inversa $F^{-1}:[0,1]\to [a,b]$ e la si applica ai valori pseudocasuali $F^{-1}(rnd)$ .

Esempi¶

Una funzione di prova¶

La funzione $\phi(x)=\frac{1}{1+x^2}$ nell’intervallo $[-1,1]$ è pari e ha un massimo assoluto per $x=0$ , è continua e positiva ma non è una densità di probabilità, perché

$\int_{-1}^1\frac{dx}{1+x^2}=2\int_0^1\frac{dx}{1+x^2}=2[\arctan x]_0^1=\frac{\pi}{2}.$

Allora trasformiamo $\phi(x)$ in una funzione di densità: $f(x)=\frac{2}{\pi(1+x^2)}$ . La funzione di ripartizione è allora:

$F(x)=\int_{-1}^x\frac{2}{\pi(1+t^2)}dt=\frac{2}{\pi}[\arctan]_{-1}^x= \frac{2}{\pi}\left(\arctan x +\frac{\pi}{4} \right)= \frac{2}{\pi}\arctan x+\frac{1}{2}.$

Per avere $F^{-1}$ , dobbiamo invertire la funzione $F$

$&y=\frac{2}{\pi}\arctan x+\frac{1}{2}\\ &\frac{2}{\pi}\arctan x=y-\frac{1}{2}\\ &\arctan x=\frac{\pi}{2}\left(y- \frac{1}{2}\right)\\ &x=F^{-1}(y)=\tan\left[\frac{\pi}{2}\left(y-\frac{1}{2} \right) \right].$

Poiché la densità di probabilità non ha un picco molto accentuato, la funzione $F$ è abbastanza simile a una retta. Cerchiamo allora di accentuare le differenze e controlliamo i risultati, che otterremo al calcolatore con Derive.

Variando opportunamente...¶

Questa volta $\phi(x)=\frac{1}{1+25x^2}$ . L’integrale, calcolato con Derive, risulta $\frac{2}{5}\arctan 5$ e quindi usiamo, come funzione densità, $f(x)=\frac{5}{2\arctan 5(1+25 x^2)}$ .

Da qui, sempre con l’aiuto di Derive, abbiamo: $F(x)=\frac{\arctan 5x}{2\arctan 5}+\frac{1}{2}$ . Seguono i grafici:

Questa volta il picco di $f$ è accentuato, quindi la funzione $F'$ di distribuzione è sufficientemente ondulata. Abbiamo perciò una funzione di trasferimento $F^{-1}$ adatta a concentrare i valori (pseudocasuali trasformati) in una zona ristretta di $[a,b]$ , in questo caso la zona centrale.

$F^{-1}(y)=\frac{1}{5\tan\left[(2y-1)\arctan\frac{1}{5}-\pi y \right]}.$

è la funzione che raccoglie i valori rnd in $[0,1]$ e li trasferisce in $[-1,1]$ , più densi intorno allo zero e più radi agli estremi.

Generaiamo mille valori $y$ con il solito algoritmo che contiene $F^{-1}(rnd)$ e mettiamo in grafico i punti $(x,y)$ corrispondenti. Comprimendo il più possibile l’asse orizzontale, si ottiene come risultato della distribuzione un’immagine piuttosto evidente.

Riassunto¶

Ci sono metodi approssimati per calcolare un integrale, utili quando non si è in grado di esprimere la primitiva. Si tratta di metodi numerici, cioè di algoritmi realizzati attraverso un linguaggio di programmazione. Sono efficaci, ma la loro efficienza varia a seconda del metodo.
Il metodo dei rettangoli approssima l’integrale realizzando un ricoprimento fatto a rettangoli dell’area da calcolare, detta trapezoide. Se i rettangoli hanno tutti uguale base $h$ , la formula di riferimento è:

$\int_a^bf(x)dx\cong h\sum_1^ny_k.$

Il metodo dei trapezi realizza un’approssimazione migliore. I rettangoli sono sostituiti da trapezi, quindi il grafico di $f$ è approssimato da una spezzata, cioè dal grafico di una funzione a dominio discreto $\phi$ . La formula di riferimento per realizzare l’algoritmo, nel caso di $n$ trapezi di altezza $h$ costante, è:

$\int_a^bf(x)dx\cong \frac{h}{2}\left[y_0+2\sum_1^{n-1}y_k+y_n \right].$
L’approssimazione polinomiale è il metodo che permette di sostituire tratti del grafico di $f$ con archi del grafico di polinomi. Nel caso di archi di parabola, il metodo si dice di integrazione parabolica e la formula di riferimento si chiama formula di Simpson:

$\int_a^bf(x)dx\cong \frac{h}{3}\left[y_0+4\sum_1^m y_{2k-1}+2\sum_1^{m-1}y_{2k}+y_{2m} \right].$

L’algoritmo realtivo fornisce l’integrale esatto dei polinomi fino al terzo grado e la sua efficienza è notevolmente superiore agli altri metodi descritti.
Il metodo Montecarlo è un metodo statistico. La superficie di cui si cerca l’area viene pensata come un bersaglio da colpire con lanci casuali. Ogni lancio individua un punto del piano entro la superficie o fuori di essa. Alla termine del gruppo di lanci si fa il conteggio di quanti sono andati a bersaglio, sul totale. Il risultato è proporzionale al valore dell’integrale. L’essenza del metodo è la generazione di numeri casuali, che avviene attraverso il software: la funzione che li esprimne è random (rnd). La formula di riferimento è:

$\int_a^bf(x)dx=\frac{m}{n}(b-a)M.$

Il metodo Montecarlo non è di grande precisione, ma è semplice ed efficiente. Viene usato per avere una rapida stima di integrali difficili.
Esistono software specializzati per la matematica che risolvono con grande efficienza anche integrali definiti e indefiniti. La loro efficienza è di molti ordini di grandezza superiore a quelli realizzzati dai nostri algoritmi. Uno di questi è Derive.
La generazione di numeri casuali viene simulata dai software: per questo si chiamano più correttamente pseudocasuali. Sono numeri nell’intervallo $[0,1]$ , le cui sequenze iniziano a ripetersi dopo milioni di numeri. La simulazione della casualità è affidabile, perchè ogni numero è indipendente dai precedenti e le sequenze generate sono omogenee, cioè ricoprono l’intervallo con uguale densità.
L’intervallo di integrazione $[a,b]$ è in genere diverso dall’intervallo $[0,1]$ dei numeri rnd. La funzione che trasforma un rnd in modo che appartenga ad $[a,b]$ è $y=(b-a)rnd+a$ .
Si può modificare questa funzione, quando si vuole una particolare densità in $[a,b]$ , in modo che sia disomogeneo secondo una certa legge. Il procedimento è spiegato in dettaglio nel Par.10.6.

Metodi numerici di integrazione¶

Il metodo dei rettangoli¶

Due esempi¶

Il metodo dei trapezi¶

Il segno dell’area¶

La formula di Simpson¶

La precisione della formula di Simpson¶

I tre metodi a confronto¶

Il metodo Montecarlo¶

Un test per il metodo Montecarlo¶

Un problema di stima¶

La soluzione esatta¶

La densità dei numeri casuali¶

Esempi¶

Una funzione di prova¶

Variando opportunamente...¶

Riassunto¶

Esercizi¶