Corso:Calcolo Numerico I1/Formule di quadratura/Integrazione automatica

Descrizione del problema

Problema: Approssimare I(f) = ∫_a^b f(x) dx con una tolleranza assoluta εₐ
o relativa ε_r assegnata.

Viene generata una sequenza di valori Iₖ con k=1,…,n, e con εₖ ₌ I(f)-Iₖ₍f)
stima dell'errore commesso.

Ci si arresta al passo k tale che

|I(f)-I_k(f)| \le \max {ε_a, ε_R*|I_k(f)| }

con Iₖ₍f) stima ragionevole dell'integrale.

Si possono scegliere due strategie:

- Strategia non adattiva: Nella procedura non adattiva Iₖ₍f) viene
  calcolato con una legge di distribuzione dei nodi fissata. Al passo k
  cnsiderom nodi equispaziati, e nel passo successivo k+1 ne considero 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

Per stimare |I(f)-Iₖ₍f)|, uso 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)!*\frac{Mₙ}{n^{n+3}} H^{n+2} f^{(n+2)}(ξₙ₎

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

I(f)-I_{2,m}(f) = (b-a)/4! \frac{M₂}{2⁵} H^4 f^{(4)}(ξ_m) = δ_m(ξ_m)*H^4
\hbox{formula 1}

Riscrivo la formula dell'errore con 2m sottointervalli:

I(f)-I_{2,2m}(f) = (b-a)/4! \frac{Mₙ}{2⁵} *(H/2)^4 f^{(4)}(ξ_{2m}) =
δ_{2m}(ξ_{2m})*\frac{H^{4}}{2⁴} \hbox{formula 2}

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

δ_m(ξ_m) \simeq δ_{2m}(ξ_{2m}) = δ

segue che

I_{2,m}-I_{2,2m} \simeq δH^4-δH^4/2^4 = δH⁴*15/16

cioè

δH^4 \simeq 16/15*|I_{2,2m}(f)-I_{2,m}(f)|

e sostituendo l'espressione trovata per δH⁴ 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

Questa strategia è utile quando la funzione ha delle singolarità.

Considero un sottointervallo qualsiasi [α,β] dell'intervallo su cui si
vuole integrare, e si pone

I(f(α,β)) = ∫_α^βf

Pongo H_0 = (β-α)/2. Per il metodo di Simpson

s_n(f(α,β)) = H_0/3(f(α)+4 f(α+H_0)+f(β))

e per l'errore si ha

I(f(α,β))-s(f(α,β)) = -H^5/90 f^{(4)}(ξ), ξ\in [α,β], \hbox{formula 1}

Considero Simson su 2 sottointervalli, e definisco

s_2 = s(f(α,α+β/2))+s(f(α+β/2,β)),I(f(α,β))-s_2(f(α,β)) = -\frac{H⁵}{90}
*[f^{(4)}(ξ)+f^{(4)}(η)], ξ\in [α,α+β/2], η\in [α+β/2,β]

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

I(f(α,β))-s_2(f(α,β)) \simeq -H^5/90*2*f^{(4)}(ξ)

e siccome H=H₀/2 si ha

I(f(α,β))-s_2(f(α,β)) \simeq -H_0^5/32*1/90*2*f^{(4)}(ξ)= -1/16 H_0^5/90
f^{(4)}(ξ) \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(α,β))-s(f(α,β)) = -15/16 H_0^5/90*f^{(4)}(ξ)H_0^5/90 f^{(4)}(ξ) =
[s(f(α,β))-s_2(f(α,β))]*16/15

e sostituendo nella formula 2:

I(f(α,β))-s_2(f(α,β)) = -1/16*16/15*(s(f(α,β))-s_2(f(α,β)))

L'errore si stima come

\frac{-s(f(α,β))-s_2(f(α,β))}{15}

Questa è una stima a posteriori dell'errore applicando Simpson quando lo
applico al generico intervallo (α,β).

Test d'arresto

Nella pratica, si considera il seguente test

|I(f(α,β))-s_2(f(α,β))| \le \frac{|s_2(f(α,β))-s(f(α,β)|}{10}

dove chiamo ε(f(α,β)) la quantità al numeratore.

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

|ε(f(α,β))|/10 \le ε*(β-α)/(b-a)

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

|I(f(α,β))-Σ_{(α_i,β_i) t.c. ∪_i (α_i,β_i) = (a,b)} s_2(f(α_i,β_i))\simeq
Σ_{(α_i,β_i)} |E_f(α_i,β_i)|/10 \le Σ_{(α_i,β_i)} ε/(b-a)*|β_i-α_i| \le ε

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. Chiamo S l'intervallo di integrazione già esaminato, in cui il
test

|E(f(α,β))/10| \le ε(β-α)/(b-a)

è verificato.

Chiamo N l'intervallo di integrazione non ancora esaminato.

- Fisso i tre intervalli:N = \emptyset A = [a,b]S = \emptyset
- Al generico istante dell'algoritmo chiamoJ_s = ∫_s f \simeq ∫_a^αfJ_A =
  ∫_α^βfSuppongo di aver calcolato J(α,β)(f). Se il test d'arresto è
  verificato, aggiungo a Jₛ₍f) il nuovo pezzo di integrale, altrimenti
  divido Jₐ₍f) in due parti.if(\frac{|E(f(α,β))}{10} | \le ε\frac{β-α}{b-a}
  )S = S ∪AA = Nα= ββ= bN = \emptyset elseα' = α+β/2A = (α,α')α= αβ= α'N =
  N ∪(α',β)end
- Sul calcolatore si controlla anche che (α,α') non diventi eccessivamente
  piccolo.

Function quad

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

- Considero il punto medio di [a,b], pongo 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₆ ₌ b-h, x₇ ₌ bLa 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₃dove \bar q_1,\bar q_2, \bar
  q₃ si ottengono applicando Simpson rispettivamente a x₁,x₂,x₃, a x₃,x₄,x₅
  e a x₅,x₆,x₇.
- Pongo \bar q = ∫_{\bar a}^{\bar b} f. Divido (\bar a, \bar b) esattamente
  a metà, considerandone il punto medio c.Pongo d = \bar a+c/2e = c/2+\bar
  b, h = \bar b-\bar a.
- Per poter applicare il test d'arresto, calcolo\bar q(s_1) = h/6 [f(\bar
  a)+4 f(c)+f(\bar b)]e raddoppiando il numero dei nodi:\bar q(s_2) = H/12
  [f(\bar a)+4 f(d)+2 f(c)+4 f(e)+f(\bar b)]e pongo poi \bar q-\bar q(s_2)
  = (\bar q(s_2)-\bar q(s₁₎₎/15 (è la stima dell'errore tra il passo 0 e il
  passo 2s).
- Se |\bar q-\bar q(s_2)|<ε, il metodo si arresta, altrimenti (\bar a,\bar
  b) viene spezzato in due parti.\bar q = ∫_{\bar a}^c f+∫_c^{\bar b} fe si
  applica ricorsivamente la function quad ai due sottointervalli ottenuti.

Ho 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.