Problema su una formula di quadratura gaussiana

[.::Luisa::.]

Ciao, avrei bisogno di un aiuto.
Devo implementare un codice in matlab per l'approssimazione di un integrale mediante la formula di quadratura di Gauss-Tchebychev.
L'approssimazione che faccio è la seguente:

$int_(-1)^(1)1/sqrt(1-x^2)f(x)d(x)approxpi/n*sum_{i=1}^{n}f(cos((2*i-1)/(2*n)*pi))$

Quindi utilizza i nodi di tchebychev.
Il programma lavora in questo modo:
-crea i primi 10 nodi ed effettua la prima approssimazione;
-parte un ciclo while in cui incremento i nodi di 5 (cioè creo altri 5 nodi) e viene effettuata una nuova approssimazione.
Il ciclo termina quando viene soddisfatta la tolleranza di errore richiesta, oppure quando viene raggiunto il numero massimo di nodi che l'utente decide di utilizzare.

Il mio problema è il seguente:
in che modo devo calcolare l'errore commesso ad ogni passo?
Avevo pensato di effettuare una differenza in modulo di due approssimazioni successive, ma il risultato che ricevo è lontano dalla soluzione vera dell'integrale.
E mi rendo conto che non è corretto fare in questo modo.
Quindi come fare?
Grazie in anticipo.

[G.D.5]

Non so rispondere al quesito matematico ma posso indicarti come scrivere il circa uguale:
\$\sim\$:$\sim$
\$\approx\$:$\approx$

[.::Luisa::.]

"WiZaRd":Non so rispondere al quesito matematico ma posso indicarti come scrivere il circa uguale:
\$\sim\$:$\sim$
\$\approx\$:$\approx$

Grazie lo stesso.

[.::Luisa::.]

So che l'errore di questa formula è dato da:

$(2*pi)/(2^(2*n)*(2*n)!)*f^(2*n)(xi)$

ma non so in che modo fare la stima nel programma.

[Sk_Anonymous]

Cara studentessa
ti sarai resa conto che qui in quanto a 'chiacchera' molti se la cavano egregiamente, quanto a 'sostanza' meglio non fare commenti ...

Oggi purtroppo sono un poco 'preso' ma ti prometto che domani esaminerò con cura il tuo problema e [sempre che nel frattempo qualcuno non si sia mosso... mai dire mai! ...] vedrò di darti una risposta...

cordiali saluti

lupo grigio



... chè perder tempo a chi più sa più spiace... Dante Alighieri, Divina Commedia, Purgatorio, III, 78

[.::Luisa::.]

"lupo grigio":Cara studentessa
ti sarai resa conto che qui in quanto a 'chiacchera' molti se la cavano egregiamente, quanto a 'sostanza' meglio non fare commenti ...

Oggi purtroppo sono un poco 'preso' ma ti prometto che domani esaminerò con cura il tuo problema e [sempre che nel frattempo qualcuno non si sia mosso... mai dire mai! ...] vedrò di darti una risposta...

Ti ringrazio e attenderò.

[Sk_Anonymous]

Come ho promesso ieri di fare eccomi qui!… Per prima cosa mi scuso con ‘studentessa’ per non essere intervenuto prima ma il fatto è che da una parte sono preso col lavoro, dall’altra il ‘clima’ nei miei confronti esistente in questi giorni sul forum non è di quelli che agevolano [diciamo così…] le ‘analisi impegnative’ …

L’espressione generale di una formula di integrazione ‘alla Gauss’ è la seguente…

$int_a^b w(x)*f(x)*dx$ ~ $sum_(i=1)^n a_i*f(x_i)$ (1)

… in cui…

a) $w(x)$ è una ‘funzione peso’ che caratterizza la ‘famiglia di appartenenza’ della funzione da integrare
b) $n$ è il ‘numero dei punti’ sui quali è effettuata l’integrazione
c) le $a_i$, $i=1,2,…,n$ sono i ‘coefficienti’ della formula di integrazione
d) le $x_i$, $i=1,2,…,n$ sono i ‘punti’ in cui occorre fornire il valore di $f(x)$
e) $a$ e $b$ sono gli estremi inferiore e superiore di integrazione

Nel caso delle formule di Gauss-Tchebychev si ha…

$w(x)=1/sqrt(1-x^2)$ (2)

… e le $x_i$ sono gli zeri del polinomio di Tchebychev di grado $n$…

$T_n(x)= cos (n*cos^(-1) x)$ (3)

Ad onta delle apparenze si tratta di un polinomio algebrico di grado $n$ 3e i suoi zeri sono…

$x_i= cos [pi*(2i-1)/(2n)]$ (4)

Tutte le $a_i$ hanno lo stesso valore [grande vantaggio di questo tipo di formula ‘alla Gauss’…] ed è…

$a_i=pi/n$ (5)

Infine gli estremi sono rispettivamente $a=-1$ e $b=1$. In definitiva quindi la (1) diviene…

$int_(-1)^1 f(x)/sqrt(1-x^2)*dx$ ~ $pi/n* sum_(i=1)^n f(x_i)$ (6)

Si può dimostrare che il cosiddetto ‘errore di troncamento’ è…

$e=pi/(2^(2n-1)*(2n!))*f^((2n))(xi)$ (7)

… ove $-1
Per stimare l’errore che si commette ‘nel caso peggiore’ ovviamente la (7) è un preciso riferimento. Dal momento che compare la derivata della $f(x)$ di ordine $2n$ è immediato osservare [altro grande vantaggio delle formule ‘alla Gauss’…] che la (6) fornisce il risultato ‘esatto’ dell’integrale se $f(x)$ è un polinomio di grado $k<=2n-1$. In questo caso la scelta del valore di $n$ sarebbe immediata. In altri casi ovviamente non è detto che tale scelta sia così semplice, specie quando le derivate di $f(x)$ in $-1
$f(x)= ln (1-x)$ (8)

Per poter rispondere in maniera più dettagliata a studentessa mi occorrerebbe sapere alcuni dettagli sulle possibili $f(x)$ ‘candidate’ e precisamente…

a) $f(x)$ è ricavata sperimentalmente?…
b) la scelta di $f(x)$ è in qualche modo ‘vincolata’ oppure $f(x)$ può essere qualsiasi
c) la scelta delle formule di Gauss-Tchebychev è dovuta al fatto che in tutte le funzioni da integrare compare al denominatore il termine $sqrt(1-x^2)$?…

cordiali saluti

lupo grigio



… chè perder tempo a chi più sa più spiace… Dante Alighieri, Divina Commedia, Purgatorio, III, 78

[.::Luisa::.]

"lupo grigio":

Per poter rispondere in maniera più dettagliata a studentessa mi occorrerebbe sapere alcuni dettagli sulle possibili $f(x)$ ‘candidate’ e precisamente…

a) $f(x)$ è ricavata sperimentalmente?…
b) la scelta di $f(x)$ è in qualche modo ‘vincolata’ oppure $f(x)$ può essere qualsiasi
c) la scelta delle formule di Gauss-Tchebychev è dovuta al fatto che in tutte le funzioni da integrare compare al denominatore il termine $sqrt(1-x^2)$?…

Innanzitutto ti ringrazio per l'ampia spiegazione iniziale, è davvero molto precisa e interessante.
Adesso rispondo alle tue domande:

a) non è ricavata sperimentalmente, ma ci è stata assegnata.
b) $f(x)$ può essere qualsiasi
c) il progetto che devo realizzare ha lo scopo non solo di creare codici che approssimino l'integrale con le formule di Gauss-Tchebychev, ma di creare anche un codice che approssimi lo stesso integrale con la formula di quadratura di Cavalieri-Simpson. Quindi data una qualsiasi funzione e un qualsiasi intervallo di integrazione, sta a me modificare il tutto in maniera opportuna e ricondurre l'integrale alla forma (6) da te indicata.

[david_e1]

Nello caso di Gauss-Chebischev (perchè Tchebischev?!? Magari parliamo di due formule diverse... boh!), c'è una formula globale per l'errore con la norma $H^s$ della funzione:

$ | \int_{-1}^1 f(x) dx - I_N f | \leq C N^{-s} || f ||_{H_w^s(-1,1)} $

dove $H^s_w(-1,1)$ è lo spazio di Sobolev con regolarità $s$ e con peso $w$ il peso di Chebischev. Ora se $n$ è un numero piccolo e $N$ un numero grande (si matematicamente è una porcheria, ma nn mi viene in mente altro):

$ | int_{-1}^1 f(x) dx - I_n f | \leq | I_n f - I_N f | + | int_{-1}^1 f(x) dx - I_N f | \approx | I_n f - I_N f | \leq C n^{-s} || f ||_{H_w^s(-1,1)} $

ovvero potresti porre stimare $ C n^{-s} || f ||_{H_w^s(-1,1)} $ dal basso usando $ | I_n f - I_N f|$ e, a quel punto, usare questa stima per valutare, nei passi successivi $ C k^{-s} || f ||_{H_w^s(-1,1)} $ e trovare l'errore su $I_k$... L'unica fregatura è che così ti si richiede di conoscere la regolarità delle funzioni, quanto meno una stima dal basso della regolarità.

Altro proprio non mi viene in mente come stimatore, magari si può pensare di approssimare $f^{(s)}$ per via derivata spettrale alla Legendre, ma mi sembra un po' troppo complicato così. (Anche se sicuramente ha più senso del mio delirio di prima! )

[.::Luisa::.]

"david_e":

dove $H^s_w(-1,1)$ è lo spazio di Sobolev con regolarità $s$ e con peso $w$ il peso di Chebischev.

Beh, io ancora non ho studiato lo spazio di Sobolev Quindi non posso utilizzare questo tipo di stima.

[Sk_Anonymous]

Cara studentessa
innanzitutto grazie per l'apprezzamento! ...

Poi vorrei formulare meglio la domanda da me espressa nel punto c). Le formule di Gauss-Tchebychev sono una sorta di 'versione specializzata' delle formule 'alla Gauss' e il loro impiego presuppone due cose...

a) l'intervallo di integrazione deve essere $-1entrambi finiti

b) la funzione da integrare contiene al denominatore il termine $sqrt(1-x^2)$

Per capire le conseguenze pratiche del punto b) supponiamo si debba calcolare il seguente integrale...

$int_(-1)^1 cos x*dx$ (1)

Ovviamente tale integrale è del tutto elementare. Supponiamo però di volerlo risolvere con le formule di Gauss-Tchebychev. In tal caso la funzione $f(x)$ da utilizzare nelle formule sarebbe la seguente...

$f(x)= sqrt(1-x^2)*cos x$ (2)

... tutte le derivate della quale sono illimitate in $-1>x>1$. In altre parole vi è il 'rischio' di ottenere in (1) un valore completamente 'sballato'...

cordiali saluti

lupo grigio



... chè perder tempo a chi più sa più spiace... Dante Alighieri, Divina Commedia, Purgatorio, III, 78

[david_e1]

Non serve conoscere la teoria degli spazi di Sobolev. Semplicamente quel simbolo $|| f ||_{H_w^s(-1,1)}$ interpretalo così. Prima si definisce la norma $L^2(-1,1)$ con peso $w$. Nel caso di Chebischev:

$ || f ||_{L_w^2(-1,1)} = \int_{-1}^1 f^2(x) \frac{1}{(1-x^2)^{1/2}} dx $

poi definisci la norma $H^s$ (se $s$ è intero) come:

$ || f ||_{H_w^s(-1,1)} = \sum_{i=1}^s || f^{(i)} ||_{L_w^2(-1,1)} $

indicando con $f^{(i)}$ la derivata $i$-sima di $f$.

In pratica questa formula non è altro che la versione "globale" dell'errore di troncamento. Niente di speciale.

[.::Luisa::.]

"lupo grigio":
In tal caso la funzione $f(x)$ da utilizzare nelle formule sarebbe la seguente...

$f(x)= sqrt(1-x^2)*cos x$ (2)

... tutte le derivate della quale sono illimitate in $-1>x>1$. In altre parole vi è il 'rischio' di ottenere in (1) un valore completamente 'sballato'...

Si, è una cosa che ho notato anche io, ma non saprei in che modo risolvere il problema.

[.::Luisa::.]

"david_e":
$ || f ||_{H_w^s(-1,1)} = \sum_{i=1}^s || f^{(i)} ||_{L_w^2(-1,1)} $

indicando con $f^{(i)}$ la derivata $i$-sima di $f$.

In pratica questa formula non è altro che la versione "globale" dell'errore di troncamento. Niente di speciale.

Ma dato che devo scrivere un programma, non posso inserire il calcolo delle derivate.
Mi serve qualcosa di numerico.
Ad esempio, utilizzando la formula di Cavalieri-Simpson so che la stima dell'errore è data dalla differenza in modulo di due approssimazioni successive diviso 15. Mi servirebbe qualcosa del genere anche per questo problema.

[david_e1]

No no no. Gauss-Chebishev non è un metodo di data-fitting, lo scopo non è approssimare una funzione e le sue derivate puntualmente, le formule di stima locale dell'errore non vanno bene (derivata illimitata etc...), bisogna usare delle stime globali.

Se noti il mio metodo non presuppone il calcolo delle derivate di $f$ etc..., ma approssimo la grandezza $C n^{-s} || f ||_{H_w^s}$ e la uso per stimare l'errore. Quello che serve è una stima di regolarità di $f$, ad esempio se $f$ è continua puoi prendere $s=1$ nei conti. Se $f$ è più regolare puoi andare su $s$ maggiori.

Altrimenti puoi usare metodi di derivazione numerica per stimare $f^{(i)}$, basta che non usi la derivata pseudo-spettrale di Gauss-Chebishev, perchè altrimenti è il serpente che si morde la coda...

PS: Se vuoi qualche cosa di più fine ti conviene cercare in letteratura: esistono libri interi solo sulle formule di Guass stracolmi di ricette più o meno euristiche per stimare l'errore.

*** EDIT ***
O era il gatto che si morde la coda?

[.::Luisa::.]

"david_e":Quello che serve è una stima di regolarità di $f$, ad esempio se $f$ è continua puoi prendere $s=1$ nei conti. Se $f$ è più regolare puoi andare su $s$ maggiori.

Il mio deve essere un codice che valga per una qualsiasi funzione scelta dall'utente, quindi non posso sapere a priori se tale funzione sia regolare.

[david_e1]

Un minimo di regolarità la devi avere comunque... altrimenti il metodo non converge nemmeno (ad esempio non puoi mica integrare numericamente la funzione di Dirichlet, quella che vale 1 sui razionali e zero altrove.).

Io ho sempre visto usare questa formula in contesti nei quali si conosce la regolarità (Gauss-Chebishev si usava di brutto nei codici di calcolo spettrali per problemi fluido-dinamici), per cui la formula che conosco io sfrutta questa conoscenza... altre formule non ne conosco (ma ne esistono di certo)...

[.::Luisa::.]

"david_e":Un minimo di regolarità la devi avere comunque... altrimenti il metodo non converge nemmeno (ad esempio non puoi mica integrare numericamente la funzione di Dirichlet, quella che vale 1 sui razionali e zero altrove.).

Beh, sicuramente non posso integrare funzioni impossibili.

[Sk_Anonymous]

Cara studentessa
se ho ben capito quello che devi realizzare è un programma general purpose che, dati una certa $f(x)$ e due numeri reali $a$ e $b$, calcoli…

$int_a^b f(x)*dx$ (1)

Volendo utilizzare delle formule ‘alla Gauss’ e limitandosi alle ipotesi seguenti…

a) $a$ e $b$ sono entrambi finiti
b) la funzione $f(x)$ è limitata, ossia esiste un numero $F$ tale che per ogni $a
… la scelta migliore l’impiego delle formule di Gauss-Legendre. In tal caso è $w(x)=1$ e i limiti di integrazione valgono $a=-1$ e $b=1$ , col che la formula di integrazione approssimata si riduce a…


$int_(-1)^1 *f(x)*dx$ ~ $sum_(i=1)^n a_i*f(x_i)$ (1)

Dando per scontato il modo per riportare l’integrale, con opportuno cambio della variabile $x$, all’intervallo $(-1,1)$ vediamo come si ricavano le $a-i$ e le $x_i$ in questo caso. Si definisce il polinomio di Legendre di ordine $n$ come…

$P_n(x)=1/(2^n*n!) d^n/(dx^n) (x^2-1)^n$ (2)

Le $x_i$, $i=1,2,…,n$ sono gli zeri del polinomio di Legendre di ordine $n$ e le $a_i$, $i=1,2,…,n$ sono date dalla formula…

$a_i=2/(n^2)* (1-x_i^2)/(|P_(n-1)(x_i)|^2)$ (3)

Sia le $x_i$ sia le $a_i$ si possono calcolare dato $n$ oppure sono disponibili in tavole non difficili da trovare in rete. L’errore di troncamento è dato da…

$e=(f^((2n)) (xi))/(2n!) int_(-1)^1 [pi(x)]^2*dx$ (4)

... ove è...

$pi(x)= prod_(i=1)^n (x-x_i)$ (5)

Anche in questo caso dunque l’errore dipende dal valore assunto dalla derivata di ordine $2n$ di $f(x)$ in un punto dell’intervallo $(-1,1)$. Pertanto le formule di Gauss-Legendre forniscono il valore ‘esatto’ dell’integrale quando $f(x)$ è un polinomio di grado $k<=2n-1$. C’è rilevare che rispetto alle formule di Gauss- Tchebychev le $x_i$, le $a_i$ e l’errore di troncamento sono assai più difficili da calcolare. In compenso però non si corrono ‘rischi’ dovuti alla presenza dello ‘scomodo’ fattore peso $1/sqrt(1-x^2)$…

cordiali saluti

lupo grigio



... chè perder tempo a chi più sa più spiace... Dante Alighieri, Divina Commedia, Purgatorio, III, 78

[.::Luisa::.]

"lupo grigio":Cara studentessa
se ho ben capito quello che devi realizzare è un programma general purpose che, dati una certa $f(x)$ e due numeri reali $a$ e $b$, calcoli…

$int_a^b f(x)*dx$ (1)

Volendo utilizzare delle formule ‘alla Gauss’ e limitandosi alle ipotesi seguenti…

a) $a$ e $b$ sono entrambi finiti
b) la funzione $f(x)$ è limitata, ossia esiste un numero $F$ tale che per ogni $a
… la scelta migliore l’impiego delle formule di Gauss-Legendre.

Non posso scegliere a piacere quale formula usare. Devo obbligatoriamente utilizzare la formula di Gauss-Tchebychev.

[david_e1]

L'integrazione secondo Gauss - Legendre è molto più complicata da implementare, a meno di non ricorrere ad algoritmi di complessità esagerata e, comunque, anche con un'implementazione "furba", il metodo di Legendre risulta molto più oneroso computazionalmente. Infatti bisogna generare i polinomi in maniera efficiente (e non certo per via ricorsiva usando derivate o produttorie!) e bisogna calcolare in maniera numerica, ma sufficientemente accurata gli zeri di tali polinomi, usando Newton-Horner.

Il motivo per cui poi, alla fine, si usa Legendre, è proprio per il peso unitario. Nella risoluzione delle EDP il peso di Chebishev è un bel problema quando si tratta di dimostrare la buona posizione del problema variazionale. Inoltre la potenza dei computer odierni è molto superiore a quella dei calcolatori degli anni '80 (gli anni d'oro di Chebishev con i metodi spettrali) e le moderne tecniche di decomposizione del dominio o degli elementi spettrali hanno reso superfluo il ricorrere a polinomi di grado gigantesco (1024 o più), rendendo superfluo l'uso della trasformata rapida (altro cavallo di battaglia dei Chebishev).