Integrazione automatica

Descrizione del problema[modifica | modifica wikitesto]

Problema: Approssimare $I(f)=\int _{a}^{b}f(x)\,dx$ $I(f)=\int _{a}^{b}f(x)\,dx$ con una tolleranza assoluta $\varepsilon _{a}$ $\varepsilon _{a}$ o relativa $\varepsilon _{r}$ $\varepsilon _{r}$ assegnata.

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

|I(f)-I_{k}(f)|\leq \max\{\varepsilon _{a},\varepsilon _{R}*|I_{k}(f)|\}

con

I_{k}(f)

stima ragionevole dell'integrale.
Si possono scegliere due strategie:

Strategia non adattiva: Nella procedura non adattiva $I_{k}(f)$ $I_{k}(f)$ viene calcolato con una legge di distribuzione dei nodi fissata. Al passo $k$ $k$ cnsidero $m$ $m$ nodi equispaziati, e nel passo successivo $k+1$ $k+1$ ne considero $2m$ $2m$ . Utilizzo le formule di Newton-Coves composite.
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[modifica | modifica wikitesto]

Per stimare $|I(f)-I_{k}(f)|$ $|I(f)-I_k(f)|$ , uso un criterio basato sull'estrapolazione di Richardson.

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

I(f)-I_{n,m}(f) = (b-a)/(n+2)!*\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})*H^{4}\quad {\hbox{formula 1}}

Riscrivo la formula dell'errore con

2m

sottointervalli:

I(f)-I_{2,2m}(f)=(b-a)/4!{\frac {M_{n}}{2^{5}}}*(H/2)^{4}f^{(4)}(\xi _{2m})=\delta _{2m}(\xi _{2m})*{\frac {H^{4}}{2^{4}}}\quad {\hbox{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}*15/16

cioè

\delta H^{4}\simeq 16/15*|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*16/15 (I_{2,2m}(f)-I_{2,m}(f)) = 1/15 (I_{2,m}(f)-I_{2,2m}(f))

Posso 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[modifica | modifica wikitesto]

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

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

Pongo

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 ],\quad {\hbox{formula 1}}

Considero Simson su 2 sottointervalli, e definisco

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}}*[f^{(4)}(\xi )+f^{(4)}(\eta )],\quad \xi \in [\alpha ,\alpha +\beta /2],\eta \in [\alpha +\beta /2,\beta ]

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

I(f(\alpha ,\beta ))-s_{2}(f(\alpha ,\beta ))\simeq -H^{5}/90*2*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*1/90*2*f^{(4)}(\xi )

=-1/16H_{0}^{5}/90f^{(4)}(\xi )\quad {\hbox{formula 2}}

e ho un fattore

16

di riduzione dell'errore.

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

s_{2}(f(\alpha ,\beta ))-s(f(\alpha ,\beta ))=-15/16H_{0}^{5}/90*f^{(4)}(\xi )

H_{0}^{5}/90f^{(4)}(\xi )=[s(f(\alpha ,\beta ))-s_{2}(f(\alpha ,\beta ))]*16/15

e sostituendo nella formula 2:

I(f(\alpha ,\beta ))-s_{2}(f(\alpha ,\beta ))=-1/16*16/15*(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 applico al generico intervallo

(\alpha ,\beta )

.

Test d'arresto[modifica | modifica wikitesto]

Nella pratica, si considera il seguente test

|I(f(\alpha ,\beta ))-s_{2}(f(\alpha ,\beta ))|\leq {\frac {|s_{2}(f(\alpha ,\beta ))-s(f(\alpha ,\beta )|}{10}}

dove chiamo

\varepsilon (f(\alpha ,\beta ))

la quantità al numeratore.
Per ottenere la stima d'errore sull'intero intervallo

[a,b]

, ci si arresta quando

|\varepsilon (f(\alpha ,\beta ))|/10\leq \varepsilon *(\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\leq \sum _{(\alpha _{i},\beta _{i})}\varepsilon /(b-a)*|\beta _{i}-\alpha _{i}|\leq \varepsilon

perché la somma delle lunghezze degli intervallini è la lunghezza di

(a,b)

.

Algoritmo dell'integrazione automatica adattiva[modifica | modifica wikitesto]

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

|E(f(\alpha ,\beta ))/10|\leq \varepsilon (\beta -\alpha )/(b-a)

è verificato.
Chiamo

N

l'intervallo di integrazione non ancora esaminato.

Fisso i tre intervalli:
$N=\emptyset$
$N=\emptyset$
$A=[a,b]$
$A=[a,b]$
$S=\emptyset$
$S=\emptyset$
Al generico istante dell'algoritmo chiamo
$J_{s}=\int _{s}f\simeq \int _{a}^{\alpha }f$
$J_{s}=\int _{s}f\simeq \int _{a}^{\alpha }f$
$J_{A}=\int _{\alpha }^{\beta }f$
$J_{A}=\int _{\alpha }^{\beta }f$
Suppongo di aver calcolato $J(\alpha ,\beta )(f)$ $J(\alpha ,\beta )(f)$ . Se il test d'arresto è verificato, aggiungo a $J_{s}(f)$ $J_{s}(f)$ il nuovo pezzo di integrale, altrimenti divido $J_{a}(f)$ $J_{a}(f)$ in due parti.
${\rm {{if}({\frac {|E(f(\alpha ,\beta ))}{10}}|\leq \varepsilon {\frac {\beta -\alpha }{b-a}})}}$
${\rm {{if}({\frac {|E(f(\alpha ,\beta ))}{10}}|\leq \varepsilon {\frac {\beta -\alpha }{b-a}})}}$
$S=S\cup A$
$S=S\cup A$
$A=N$
$A=N$
$\alpha =\beta$
$\alpha =\beta$
$\beta =b$
$\beta =b$
$N=\emptyset$
$N=\emptyset$
${\rm {else}}$
${\rm {else}}$
$\alpha '=\alpha +\beta /2$
$\alpha '=\alpha +\beta /2$
$A=(\alpha ,\alpha ')$
$A=(\alpha ,\alpha ')$
$\alpha =\alpha$
$\alpha =\alpha$
$\beta =\alpha '$
$\beta =\alpha '$
$N=N\cup (\alpha ',\beta )$
$N=N\cup (\alpha ',\beta )$
${\rm {end}}$
${\rm {end}}$
Sul calcolatore si controlla anche che $(\alpha ,\alpha ')$ $(\alpha ,\alpha ')$ non diventi eccessivamente piccolo.

Function quad[modifica | modifica wikitesto]

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

Considero il punto medio di $[a,b]$ $[a,b]$ , pongo $h=(b-a)/C$ $h=(b-a)/C$ , e considero 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$
$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={\bar {q}}_{1}+{\bar {q}}_{2}+{\bar {q}}_{3}$
$q={\bar {q}}_{1}+{\bar {q}}_{2}+{\bar {q}}_{3}$
dove ${\bar {q}}_{1},{\bar {q}}_{2},{\bar {q}}_{3}$ ${\bar {q}}_{1},{\bar {q}}_{2},{\bar {q}}_{3}$ si ottengono applicando Simpson rispettivamente a $x_{1},x_{2},x_{3}$ $x_{1},x_{2},x_{3}$ , a $x_{3},x_{4},x_{5}$ $x_{3},x_{4},x_{5}$ e a $x_{5},x_{6},x_{7}$ $x_{5},x_{6},x_{7}$ .
Pongo ${\bar {q}}=\int _{\bar {a}}^{\bar {b}}f$ ${\bar {q}}=\int _{\bar {a}}^{\bar {b}}f$ . Divido $({\bar {a}},{\bar {b}})$ $({\bar {a}},{\bar {b}})$ esattamente a metà, considerandone il punto medio $c$ $c$ .Pongo $d={\bar {a}}+c/2$ $d={\bar {a}}+c/2$ $e=c/2+{\bar {b}}$ $e=c/2+{\bar {b}}$ , $h={\bar {b}}-{\bar {a}}$ $h={\bar {b}}-{\bar {a}}$ .
Per poter applicare il test d'arresto, calcolo
${\bar {q}}(s_{1})=h/6[f({\bar {a}})+4f(c)+f({\bar {b}})]$
${\bar {q}}(s_{1})=h/6[f({\bar {a}})+4f(c)+f({\bar {b}})]$
e raddoppiando il numero dei nodi:
${\bar {q}}(s_{2})=H/12[f({\bar {a}})+4f(d)+2f(c)+4f(e)+f({\bar {b}})]$
${\bar {q}}(s_{2})=H/12[f({\bar {a}})+4f(d)+2f(c)+4f(e)+f({\bar {b}})]$
e pongo poi ${\bar {q}}-{\bar {q}}(s_{2})=({\bar {q}}(s_{2})-{\bar {q}}(s_{1}))/15$ ${\bar {q}}-{\bar {q}}(s_{2})=({\bar {q}}(s_{2})-{\bar {q}}(s_{1}))/15$ (è la stima dell'errore tra il passo 0 e il passo $2s$ $2s$ ).
Se $|{\bar {q}}-{\bar {q}}(s_{2})|<\varepsilon$ $|{\bar {q}}-{\bar {q}}(s_{2})|<\varepsilon$ , il metodo si arresta, altrimenti $({\bar {a}},{\bar {b}})$ $({\bar {a}},{\bar {b}})$ viene spezzato in due parti.
${\bar {q}}=\int _{\bar {a}}^{c}f+\int _{c}^{\bar {b}}f$
${\bar {q}}=\int _{\bar {a}}^{c}f+\int _{c}^{\bar {b}}f$
e si applica ricorsivamente la function quad ai due sottointervalli ottenuti.