Calcolo numerico/Integrazione automatica

Descrizione del problema

Problema: approssimare $I(f)={\int }_a^{b}f(x)dx$ con una tolleranza assoluta ${\epsilon }_a$ o relativa ${\epsilon }_r$ assegnata.

Viene generata una sequenza di valori $I_k$ con $k=1,\dots ,n$, e con ${\epsilon }_k=I(f)-I_k(f)$ stima dell'errore commesso. Ci si arresta al passo $k$ tale che

$|I(f)-I_k(f)|\le max\{{\epsilon }_a,{\epsilon }_R\ast |I_k(f)|\}$

con $I_k(f)$ stima ragionevole dell'integrale.

Si possono scegliere due strategie.

1. Strategia non adattiva: nella procedura non adattiva $I_k(f)$ viene calcolato con una legge di distribuzione dei nodi fissata. Al passo $k$ consideriamo $m$ nodi equispaziati, e nel passo successivo $k+1$ ne consideriamo $2m$. Utilizziamo le formule di Newton-Coves composite.
2. Strategia adattiva: la distribuzione dei nodi non è fissata a priori, gli intervalli vengono partizionati a seconda delle esigenze e integrati separatamente.

Stima dell'errore a priori nella strategia non adattiva

Per stimare $|I(f)-I_k(f)|$, usiamo un criterio basato sull'estrapolazione di Richardson.

In base alle stime di errore a priori per le formule composite e $n$ pari, si ha

$I(f)-I_{n,m}(f)=(b-a)/(n+2)!\ast \frac{M_n}{n^{n+3}}{H}^{n+2}{f}^{(n+2)}({\xi }_n)$

Ad esempio, nel caso di Simpson $n=2$ e si ha:

$I(f)-I_{2,m}(f)=(b-a)/4!\frac{M_2}{2^5}{H}^4{f}^{(4)}({\xi }_m)={\delta }_m({\xi }_m)\ast H^4\text{formula 1}$

Riscrivo la formula dell'errore con $2m$ sottointervalli:

$I(f)-I_{2,2m}(f)=(b-a)/4!\frac{M_n}{2^5}\ast (H/2{)}^4{f}^{(4)}({\xi }_{2m})={\delta }_{2m}({\xi }_{2m})\ast \frac{H^4}{2^4}\text{formula 2}$

Sottraggo membro a membro le formule 1 e 2, e tenendo conto che

${\delta }_m({\xi }_m)\simeq {\delta }_{2m}({\xi }_{2m})=\delta$

segue che

$I_{2,m}-I_{2,2m}\simeq \delta H^4-\delta H^4/2^4=\delta H^4\ast 15/16$

cioè

$\delta H^4\simeq 16/15\ast |I_{2,2m}(f)-I_{2,m}(f)|$

e sostituendo l'espressione trovata per $\delta H^4$ nella formula 2

$I(f)-I_{2,2m}(f)\simeq 1/16\ast 16/15(I_{2,2m}(f)-I_{2,m}(f))=1/15(I_{2,m}(f)-I_{2,2m}(f))$

Possiamo stimare l'errore commesso con la formula di Simpson su $2m$ sottointervalli sulla base del calcolo di due formule di quadratura approssimata. Si stima che l'errore assoluto, quando si passa a $2m$ sottointervalli, si riduca di 15.

Stima dell'errore a posteriori nella strategia adattiva

Questa strategia è utile quando la funzione ha delle singolarità. Consideriamo un sottointervallo qualsiasi $[\alpha ,\beta ]$ dell'intervallo su cui si vuole integrare, e si pone

$I(f(\alpha ,\beta ))={\int }_{\alpha }^{\beta }f$

Poniamo $H_0=(\beta -\alpha )/2$. Per il metodo di Simpson

$s_n(f(\alpha ,\beta ))=H_0/3(f(\alpha )+4f(\alpha +H_0)+f(\beta ))$

e per l'errore si ha

$I(f(\alpha ,\beta ))-s(f(\alpha ,\beta ))=-H^5/90f^{(4)}(\xi ),\xi \in [\alpha ,\beta ],\text{formula 1}$

Consideriamo Simpson su 2 sottointervalli, e definiamo

$s_2=s(f(\alpha ,\alpha +\beta /2))+s(f(\alpha +\beta /2,\beta )),$

$I(f(\alpha ,\beta ))-s_2(f(\alpha ,\beta ))=-\frac{H^5}{90}\ast [f^{(4)}(\xi )+f^{(4)}(\eta )],\xi \in [\alpha ,\alpha +\beta /2],\eta \in [\alpha +\beta /2,\beta ]$

Se assumiamo che le due derivate quarte siano circa uguali, si può scrivere

$I(f(\alpha ,\beta ))-s_2(f(\alpha ,\beta ))\simeq -H^5/90\ast 2\ast f^{(4)}(\xi )$

e siccome $H=H_0/2$ si ha

$I(f(\alpha ,\beta ))-s_2(f(\alpha ,\beta ))\simeq -H_0^5/32\ast 1/90\ast 2\ast f^{(4)}(\xi )$

$=-1/16H_0^5/90f^{(4)}(\xi )\text{formula 2}$

e abbiamo un fattore 16 di riduzione dell'errore.

Sottraiamo membro a membro le formule di errore 1 e 2.

$s_2(f(\alpha ,\beta ))-s(f(\alpha ,\beta ))=-15/16H_0^5/90\ast f^{(4)}(\xi )$

$H_0^5/90f^{(4)}(\xi )=[s(f(\alpha ,\beta ))-s_2(f(\alpha ,\beta ))]\ast 16/15$

e sostituendo nella formula 2:

$I(f(\alpha ,\beta ))-s_2(f(\alpha ,\beta ))=-1/16\ast 16/15\ast (s(f(\alpha ,\beta ))-s_2(f(\alpha ,\beta )))$

L'errore si stima come

$\frac{-s(f(\alpha ,\beta ))-s_2(f(\alpha ,\beta ))}{15}$

Questa è una stima a posteriori dell'errore applicando Simpson quando lo applichiamo al generico intervallo $(\alpha ,\beta )$.

Test d'arresto

Nella pratica, si considera il seguente test

$|I(f(\alpha ,\beta ))-s_2(f(\alpha ,\beta ))|\le \frac{|s_2(f(\alpha ,\beta ))-s(f(\alpha ,\beta )|}{10}$

dove chiamiamo $\epsilon (f(\alpha ,\beta ))$ la quantità al numeratore.

Per ottenere la stima d'errore sull'intero intervallo $[a,b]$, ci si arresta quando

$|\epsilon (f(\alpha ,\beta ))|/10\le \epsilon \ast (\beta -\alpha )/(b-a)$

Infatti, se questo avviene, per l'errore totale si ha:

$|I(f(\alpha ,\beta ))-\sum _{({\alpha }_i,{\beta }_i)t.c.\bigcup _i({\alpha }_i,{\beta }_i)=(a,b)}{s}_2(f({\alpha }_i,{\beta }_i))$

$\simeq \sum _{({\alpha }_i,{\beta }_i)}|E_f({\alpha }_i,{\beta }_i)|/10\le \sum _{({\alpha }_i,{\beta }_i)}\epsilon /(b-a)\ast |{\beta }_i-{\alpha }_i|\le \epsilon$

perché la somma delle lunghezze degli intervallini è la lunghezza di $(a,b)$.

Algoritmo dell'integrazione automatica adattiva

Chiamiamo $A$ l'intervallo di integrazione attivo in cui si deve calcolare l'integrale. Chiamiamo $S$ l'intervallo di integrazione già esaminato, in cui il test

$|E(f(\alpha ,\beta ))/10|\le \epsilon (\beta -\alpha )/(b-a)$

è verificato.

Chiamiamo $N$ l'intervallo di integrazione non ancora esaminato.

Fissiamo i tre intervalli:

$N=\mathrm{\varnothing }$

$A=[a,b]$

$S=\mathrm{\varnothing }$

Al generico istante dell'algoritmo chiamiamo

$J_s={\int }_{s}f\simeq {\int }_a^{\alpha }f$

$J_A={\int }_{\alpha }^{\beta }f$

Supponiamo di aver calcolato $J(\alpha ,\beta )(f)$. Se il test d'arresto è verificato, aggiungiamo a $J_s(f)$ il nuovo pezzo di integrale, altrimenti dividiamo $J_a(f)$ in due parti.

$\mathrm{i}\mathrm{f}\mathrm{(}\frac{|E(f(\alpha ,\beta ))}{10}\mathrm{|}\mathrm{\le }\mathrm{\epsilon }\frac{\beta -\alpha }{b-a}\mathrm{)}$

$S=S\cup A$

$A=N$

$\alpha =\beta$

$\beta =b$

$N=\mathrm{\varnothing }$

$\mathrm{e}\mathrm{l}\mathrm{s}\mathrm{e}$

${\alpha }^{\prime }=\alpha +\beta /2$

$A=(\alpha ,{\alpha }^{\prime })$

$\alpha =\alpha$

$\beta ={\alpha }^{\prime }$

$N=N\cup ({\alpha }^{\prime },\beta )$

$\mathrm{e}\mathrm{n}\mathrm{d}$

Sul calcolatore si controlla anche che $(\alpha ,{\alpha }^{\prime })$ non diventi eccessivamente piccolo.

Function quad

La function quad è il metodo utilizzato da Matlab per calcolare gli integrali ed è ricorsiva.

Consideriamo il punto medio di $[a,b]$, poniamo $h=(b-a)/C$, e consideriamo i sei sottointervalli ottenuti ponendo:

$x_1=a,x_2=a+h,x_3=a+2h,x_4=(a+b)/2,x_5=b-2h,x_6=b-h,x_7=b$

La distribuzione non è esattamente uniforme.

Si applica Simpson sui punti a tre a tre, e si ottiene l'approssimazione dell'integrale:

$q={\overline{q}}_1+{\overline{q}}_2+{\overline{q}}_3$

dove ${\overline{q}}_1,{\overline{q}}_2,{\overline{q}}_3$ si ottengono applicando Simpson rispettivamente a $x_1,x_2,x_3$, a $x_3,x_4,x_5$ e a $x_5,x_6,x_7$.

Poniamo $\overline{q}={\int }_{\overline{a}}^{\overline{b}}f$. Dividiamo $(\overline{a},\overline{b})$ esattamente a metà, considerandone il punto medio $c$. Poniamo $d=\overline{a}+c/2$ $e=c/2+\overline{b}$, $h=\overline{b}-\overline{a}$.

Per poter applicare il test d'arresto, calcoliamo

$\overline{q}(s_1)=h/6[f(\overline{a})+4f(c)+f(\overline{b})]$

e raddoppiando il numero dei nodi:

$\overline{q}(s_2)=H/12[f(\overline{a})+4f(d)+2f(c)+4f(e)+f(\overline{b})]$

e poniamo poi $\overline{q}-\overline{q}(s_2)=(\overline{q}(s_2)-\overline{q}(s_1))/15$ (è la stima dell'errore tra il passo 0 e il passo $2s$).

Se $|\overline{q}-\overline{q}(s_2)|<\epsilon$, il metodo si arresta, altrimenti $(\overline{a},\overline{b})$ viene spezzato in due parti.

$\overline{q}={\int }_{\overline{a}}^{c}f+{\int }_c^{\overline{b}}f$

e si applica ricorsivamente la function quad ai due sottointervalli ottenuti.

Abbiamo una chiamata ricorsiva alla funzione, l'integrale complessivo viene calcolato come somma di due integrali.

Data una funzione con una singolarità, si spezza l'integrale in quel punto e si fa la somma degli integrali separati.