Newsgroup: it.scienza.fisica
Titolo: Errore sull'integrazione numerica

Chiedo a voi in quanto ho già ricevuto risposte alquanto discordanti...
Mi trovo a dover calcolare l'integrale di una serie di dati sperimentali,
ovvero l'area sottesa dalla spezzata che congiunge questi punti (tipicamente
la forma e' gaussianoide ma non posso fare troppo affidamento su questa
affermazione). I punti hanno un certo errore noto sull'ordinata e un errore
sostanzialmente trascurabile sull'ascissa.
Fin qui nessun problema, trapezi o Simpson e si fa...
Ma l'errore sul valore di integrale ottenuto?!? Pensavo di trovare l'errore
come differenza tra area minima e massima (ovvero sottraendo o aggiungendo a
tutti i punti l'errore sulla y). Ma mi e' stato giustamente obiettato che
compierei una sovrastima notevole dell'errore (un po' come quando si prova a
fare un fit di retta con la massima e minima pendenza). Allora mi e' stato
proposto di fare un fit di una polinomiale, ma a questo punto dovrei
decidere a priori di che grado.
Non c'e' una terza via?
Grazie
Lorenzo

 decidere a priori di che grado.
Hai provato semplicemente ad applicare la formuletta della propagazione 
dell'errore?  Se non hai altre informazioni a-priori, in genere, e' 
quella la formula piu' appropriata anche se sovrastima l'errore.
Non credo che hai modo di giustificare un ottimismo maggiore!
Se utilizzi uno schema semplice di integrazione, sospetto che la 
propagazione degli errori coincida con la differenza che citi.
Daniele Fua'
Uni. Milano-Bicocca

 Ma l'errore sul valore di integrale ottenuto?!?
Ni sembra semplice...
Non entro sull'adeguatezza o meno del metodo d'integrazione, ma
qualunque sia il metodo che userai, sara' sempre del tipo
I = \sum a_k y_k
dove gli a_k sono coeffcienti caratteristici del metodo, e y_k le
ordinate (dati sperimentali).
Se puoi assumere gli errori sulle y_k come indipendenti, avrai, dalla
propagazione degli errori:
Var(I) = \sum a_k^2 Var(y_k).
------------------------------
Elio Fabri
Dip. di Fisica - Univ. di Pisa
------------------------------

 Chiedo a voi in quanto ho gią ricevuto risposte alquanto discordanti...

Io userei le spline cubiche per approssimare l'integrale.
Francesco Daddi
--------------------------------
Inviato via 
http://arianna.libero.it/usenet/

Non c'e' una terza via?
Ci sarà di sicuro. : )
saluti,
Peltio
[1] Ma povero Cavalieri: lui la formula di quadratura l'ha scoperta almeno
un secolo prima di Simpson!

 decidere a priori di che grado.
Come al solito ognuno dice la sua ma sono d'accordo con una citazione di 
Shakespeare che dice qualcosa come "tanto casino per nulla".
Ripeto e Fabri (a parte un quadrato che non mi "quadra") mi sembra 
d'accordo con me: se non sai (o vuoi ignorare) altre informazioni sul 
legame tra il risultato finale e le misure, la  propagazione degli 
errori e' il metodo piu' corretto per stimare l'errore.  Sovrastima 
proprio perche' ignori altre informazioni.
Nel tuo caso, se usi un metodo lineare di integrazione (cosa molto 
ragionevole), la propagazione dell'errore da' un errore che e' la META' 
(chiedo venia per quanto scritto prima) di quella differenza che avevi 
gia' ipotizzato fin dall'inizio.
Se utilizzi dei fits, vuol dire che hai altre informazioni e allora e' 
tutta un'altra storia...
Daniele Fua'
Uni. Milano-Bicocca

 tutta un'altra storia...
D'accordo. In quel caso nel fit ci saranno dei parametri, e dal
procedimento di fit si ricavera' la matrice di covarianza dei
parametri.
Poi si puo' calcolare l'integrale della funzione, che contiene ancora
i parametri, e si puo' propagare l'errore sull'integrale (tenendo conto
che i parametri *non sono* indipendenti: la matrice di cov. non sara'
diagonale, di regola).
------------------------------
Elio Fabri
Dip. di Fisica - Univ. di Pisa
------------------------------
X-Moz

        I[f] - Tn[f] = (b-a)/12 h^2 f''[chi]
Uff...
Mi sono perso il segno meno per strada, qui. E' chiaramente
        I[f] - Tn[f] = - (b-a)/12 h^2 f''[chi]
visto che l'area sottesa da una funzione concava viene sovrastimata e,
viceversa quella di una funzione convessa viene sottostimata.
Ho fatto due conti con un integrale banale: quello del seno tra 0 ed 1
Il risultato esatto č
        I = 0.459697694131860...
Il risultato della formula di quadratura con h=0.01 e i valori della
funzione seno esatti nei limiti della precisione macchina č
        Tn = 0.4596938633113578...
L'errore del metodo (che comprende anche l'errore di troncamento a 16 cifre,
a dire il vero) č quindi
        |I-Tn| = 3.8 x 10^-6.
La stima considerando un valore massimo per la derivata seconda pari a 1 (si
potrebbe fare di meglio) č di poco piů del doppio:
        8.3 x 10^-6
Se poi integriamo la funzione con degli errori casuali di modulo massimo
pari a 0.05 (err=0.05(rnd - 1/2)) (con rnd numero casuale compreso tra 0 e
1), uno dei risultati che si puň ottenere č
        Tnappr = 0.4431661573039723...
che differisce dalla quadratura numerica sui valori 'esatti' di f per
        |Tn-Tnappr| = 0.016528...
contro gli
        0.05.
(poco piů del triplo) previsti dalla relazione |Tn-Tnappr|<=(b-a)Max[err].
L'errore rispetto al valore esatto dell'integrale č chiaramente dominato
dall'effetto degli errori sui valori di f (per forza, li ho scelti enormi
: ) ):
        |I-Tnappr| = 0.016531...
saluti,
Peltio

 varianza!
Non mi quadrava perche' sono rimbambito... capita dopo una certa eta'.
Daniele Fua'
Uni. Milano-Bicocca

 |Tn-Tnappr|<=(b-a)Max[err].
Tanto per completezza, avresti potuto calcolare anche la varianza che
avevo dato io.


E comunque ripeto che a mio parere la domanda non era quella.
Io vedo i tre calcoli come tre diverse stime della stessa quantità, ossia
l'incertezza nel risultato dovuto alla propagazione dell'errore nei
dati (yk) nella sommatoria che rappresenta l'integrale numerico.
    Nel primo caso considero il caso peggiore possibile: il valore yapprk
dista dal valore esatto la bellezza di metà del massimo range riscontrato
in tutte le cinquine. Risolvo tutto con una sola moltiplicazione, ma ottengo
una stima molto rozza e decisamente in eccesso.
    Nel secondo caso considero un po' meno di sfortuna, ossia che il valore
yapprk disti dal valore esatto yk metà del range associato alla cinquina
k-esima. Posso avere errori minori dell'errore massimo, ma ancora li
considero che agiscono nello stesso senso.
    Nel caso della propagazione della varianza considero un errore probabile
che tiene conto delle possibili compensazioni per via della distribuzione
un po' per eccesso e un po' per difetto.
Sull'altra questione faccio un post a parte.
Mi è venuto un dubbio filosofico.
saluti,
Peltio

nessun modo di stimare f".
Io quella formula l'avevo data solo per distinguere l'errore dovuto al
metodo numerico (ossia al fatto che sto integrando un polinomio interpolante
lineare a tratti invece della funzione - qualunque essa sia) da quello
dovuto agli errori sui dati.
Ma questa tua osservazione mi fa sorgere un dubbio di carattere filosofico.
Non ho più il post dell'OP e non ricordo a che tipo di grandezza fisica (se
continua, discreta, un conteggio, o che altro) si riferisse la misura per
cui perdonami se la tua osservazione si riferiva ad altro.
Il dubbio che ho ora è questo:
    quando una serie di misure si riferisce a una grandezza fisica
    continua si può assumere che 'alle spalle' delle misure ci sia
    una funzione continua che descrive il 'valore esatto' della
    grandezza, cui poi vengono sommati gli errori di vario genere?
Perchè se la grandezza fisica si comporta con una certa regolarità (non ha
discontinuità, cuspidi eccetera), ossia si può ipotizzare una sua
appartenenza a Cn per un certo n sull'intervallo di integrazione, posso
utilizzare la formula
        I[f] - Tn[f] = - (b-a)/12 h^2 f''[chi]
e anche andarmi a cercare una stima asintotica in base a
        |I[f] - Tn[f]| ~ h/12  |f'[b]-f'[a]|
Ora il problema (grosso, è vero) sta nello stimare la derivata agli estremi.
Anzi no. Il problema sta nello stimare l'errore della derivata numerica che
si sa richiede a sua volta informazioni sulla derivata successiva, è molto
sensibile all'incertezza sui dati e al decrescere del passo h sotto il
valore ottimale. Per avere un'idea di massima si potrebbe vedere come cambia
il valore computato quando considero i dati agli estremi dei rispettivi
range.
ESEMPIO (lunghetto)
==================
Nel caso della funzione seno dell'esempio di integrazione si ottiene una
stima rozzissima purtroppo, in quanto l'errore in corrispondenza dello 0 si
fa parecchio sentire. Ma ho anche generato dati con una incertezza fissa e
tutt'altro che trascurabile. I valori di yapprk, dk e sk per i primi e
ultimi quattro punti sono
    0.000768, 0.012103, 0.009572
    0.018631. 0.013390, 0.013291
    0.01635, 0.013735, 0.010971
    0.033303, 0.013334, 0.010425
    ...
    0.821439, 0.011254,  0.008542
    0.830764, 0.013689, 0.010583
    0.834993, 0.008729, 0.006793
    0.839114, 0.007852, 0.005627
Uso la formula più semplice a due punti per cercare una stima della derivata
prima. Anche prendendo un passo h = 0.03 (il primo e il quarto punto), la
derivata prima all'estremo a è quasi ignota. Può spaziare tra
((0.033303 - 0.013334) - (0.000768 + 0.012103)) / 0.03= 0.24
((0.033303+0.013334) - (0.000768-0.012103)) / 0.03= 1.93
la media è 0.845 (il valore esatto è 1)
La derivata (all'indietro) all'estremo b varia tra i seguenti valori
(paura):
((0.839114 - 0.007852) - (0.821439 + 0.011254)) / 0.03 = -0.05
((0.839114 + 0.007852) - (0.821439 - 0.011254)) / 0.03 = 1.23
la media è 0.64 (il valore esatto è cos[1]= 0.54)
Del resto con h= 0.03 e d= 0.01 c'era da attendersi un errore di derivazione
_minimo_ (dovuto all'incertezza sui dati) di 0.01/0.03 pari a d/h= 1/3.
Con questi dati, per quanto poco precisi (da vergognarcisi), possiamo dare
una stima dell'errore di integrazione
     h/12  |f'[b]-f'[a]|=0.01/12 (1.23-0.24) = 0.000825 = 0.825 10^-3
contro il valore asintotico teorico di
    0.01/12 |0.54-1| = 0.000383 = 0.383 10^-3
Entrambi questi valori sono comunque eccessivamente pessimistici, visto che
la stima basata sulla massima derivata seconda (1,per eccesso) fornisce
    (b-a)/12 h^2 |f''[chi]|<=0.000008 = 0.008 10^-3
e l'errore effettivamente calcolato rispetto all'integrale con i valori
esatti yk vale:
    |I[f]-Tn[f]|=0.000004 = 0.004 10^-3
(ed è sotto tutte le stime fatte).
Questo esempio in particolare è forse un caso in cui l'errore nei dati era
un po' eccessivo (circa il 10% mi pare all'estremo b ma esagerato per valori
bassi di x). Ma in generale su funzioni più serie con errori dell'ordine
dell'1 o 2%, così come è possibile calcolare un integrale di dati empirici
dovrebbe essere possibile calcolare anche una derivata, o no?
A meno di non avere un h troppo piccolo e un errore sui dati troppo grande
una stima a quattro punti della derivata prima potrebbe fornire una stima
adatta ad essere maggiorata per determinare l'errore di integrazione.
O sto fantasticando troppo?
saluti,
Peltio