Problema su una formula di quadratura gaussiana
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.
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.
Risposte
Per avere un’idea della ‘efficacia’ dell’impiego delle formule di Gauss-Tchebishev vale la pena di provare il seguente test. Si vuol calcolare l’integrale…
$int_(-1)^(1) dx$ (1)
… utilizzando le formule di Gauss-Tchebychev e variando il numero di punti $n$. Dire che l’integrale (1) è ‘elementare’ è chiaramente una presa in giro e quello che si vuole verificare è la precisione con la quale si ottiene, su un computer che lavora a 32 bit, il risultato ‘noto’. Volendo applicare le formule di Gauss-Tchebychev occorre naturalmente moltiplicare la funzione da integrare per il fattore peso…
$w(x)=sqrt(1-x^2)$ (2)
… funzione che ha il non trascurabile difetto di avere le tutte le derivate illimitate in $-1
$n=2$, $int = 2.214414146908$
$n=4$, $int = 2.05234430595$
$n=8$, $int = 2.0129090856$
$n=16$, $int = 2.00321637817$
$n=32$, $int = 2.00080341631$
$n=64$, $int = 2.00020081173$
$n=128$, $int = 2.00005020029$
$n=256$, $int = 2.00001254991$
$n=512$, $int = 2.00000313747$
$n=1024$, $int = 2.00000078437$
$n=2048$, $int = 2.00000019609$
$n=4096$, $int = 2.00000004902$
$n=8192$, $int = 2.00000001226$
Diciamo che l’errore scende all’incirca come $1/(n^2)$ per cui se si vuole una precisione di $12$ cifre $n$ dovrebbe essere preso dell’ordine di $130.000$
… Anche se il computer riuscisse a fare un lavoro del genere, sarebbe in ogni caso difficile stimare l’errore di arrotondamento accumulato in un numero così elevato di moltiplicazioni…
Domanda a studentessa: sicura che la scelta delle formule di Gauss-Tchebychev sia ‘azzeccata'?
...
cordiali saluti
lupo grigio
... chè perder tempo a chi più sa più spiace... Dante Alighieri, Divina Commedia, Purgatorio, III, 78
$int_(-1)^(1) dx$ (1)
… utilizzando le formule di Gauss-Tchebychev e variando il numero di punti $n$. Dire che l’integrale (1) è ‘elementare’ è chiaramente una presa in giro e quello che si vuole verificare è la precisione con la quale si ottiene, su un computer che lavora a 32 bit, il risultato ‘noto’. Volendo applicare le formule di Gauss-Tchebychev occorre naturalmente moltiplicare la funzione da integrare per il fattore peso…
$w(x)=sqrt(1-x^2)$ (2)
… funzione che ha il non trascurabile difetto di avere le tutte le derivate illimitate in $-1
$n=2$, $int = 2.214414146908$
$n=4$, $int = 2.05234430595$
$n=8$, $int = 2.0129090856$
$n=16$, $int = 2.00321637817$
$n=32$, $int = 2.00080341631$
$n=64$, $int = 2.00020081173$
$n=128$, $int = 2.00005020029$
$n=256$, $int = 2.00001254991$
$n=512$, $int = 2.00000313747$
$n=1024$, $int = 2.00000078437$
$n=2048$, $int = 2.00000019609$
$n=4096$, $int = 2.00000004902$
$n=8192$, $int = 2.00000001226$
Diciamo che l’errore scende all’incirca come $1/(n^2)$ per cui se si vuole una precisione di $12$ cifre $n$ dovrebbe essere preso dell’ordine di $130.000$
… Anche se il computer riuscisse a fare un lavoro del genere, sarebbe in ogni caso difficile stimare l’errore di arrotondamento accumulato in un numero così elevato di moltiplicazioni… Domanda a studentessa: sicura che la scelta delle formule di Gauss-Tchebychev sia ‘azzeccata'?
... cordiali saluti
lupo grigio
... chè perder tempo a chi più sa più spiace... Dante Alighieri, Divina Commedia, Purgatorio, III, 78
Eseguiamo ora un altro test, scegliendo un integrale leggermente ‘meno elementare’ del precedente, vale a dire…
$int_(-1)^1 ln (1+x)*dx= 2*(ln 2-1)= -.61370563888…$
Ecco alcuni risultati ottenuti con vari valori di $n$…
$n=2$, $int = -.769892945536$
$n=4$, $int = -.692300760348$
$n=8$, $int = -.642347443302$
$n=16$, $int = -.623090388535$
$n=32$, $int = -.616607917421$
$n=64$, $int = -.614570321544$
$n=128$, $int = -.613956598957$
$n=256$, $int = -.613777077293$
$n=512$, $int = -.613725673168$
$n=1024$, $int = -.61371119113$
$n=2048$, $int = -613707162863$
$n=4096$, $int = -.613706053856$
$n=8192$, $int = -.613705751119$
E’ visibile l’effetto della singolarità della funzione logaritmo per $x=-1$. Mentre in precedenza con $8192$ punti si ottenevano otto cifre decimali esatte qui dobbiamo accontentarci di sei. L’errore in questo caso dipende esclusivamente dalle ‘limitazioni intrinseche’ dell’algoritmo utilizzato, dal momento che l’errore di arrotondamento deve essere all’incirca lo stesso nei due casi. Si potrebbero certamente testare altri tipi di funzioni per dare un parere definitivo, anche se a me pare che i risultati ottenuti sinora siano abbastanza inequivocabili…
cordiali saluti
lupo grigio
... chè perder tempo a chi più sa più spiace... Dante Alighieri, Divina Commedia, Purgatorio, III, 78
$int_(-1)^1 ln (1+x)*dx= 2*(ln 2-1)= -.61370563888…$
Ecco alcuni risultati ottenuti con vari valori di $n$…
$n=2$, $int = -.769892945536$
$n=4$, $int = -.692300760348$
$n=8$, $int = -.642347443302$
$n=16$, $int = -.623090388535$
$n=32$, $int = -.616607917421$
$n=64$, $int = -.614570321544$
$n=128$, $int = -.613956598957$
$n=256$, $int = -.613777077293$
$n=512$, $int = -.613725673168$
$n=1024$, $int = -.61371119113$
$n=2048$, $int = -613707162863$
$n=4096$, $int = -.613706053856$
$n=8192$, $int = -.613705751119$
E’ visibile l’effetto della singolarità della funzione logaritmo per $x=-1$. Mentre in precedenza con $8192$ punti si ottenevano otto cifre decimali esatte qui dobbiamo accontentarci di sei. L’errore in questo caso dipende esclusivamente dalle ‘limitazioni intrinseche’ dell’algoritmo utilizzato, dal momento che l’errore di arrotondamento deve essere all’incirca lo stesso nei due casi. Si potrebbero certamente testare altri tipi di funzioni per dare un parere definitivo, anche se a me pare che i risultati ottenuti sinora siano abbastanza inequivocabili…
cordiali saluti
lupo grigio
... chè perder tempo a chi più sa più spiace... Dante Alighieri, Divina Commedia, Purgatorio, III, 78
"lupo grigio":
Domanda a studentessa: sicura che la scelta delle formule di Gauss-Tchebychev sia ‘azzeccata'?
Beh, la scelta non l'ho fatta io, ma il mio prof.
Infatti la traccia dice:
Sviluppo dei codici matlab per il calcolo di un integrale definito mediante la formula di Gauss- Tchebychev.
Questa "dimostrazione" sull'efficacia delle formule di Gauss-Chebishev (a me risulta senza la "T" davanti) mi sembra molto "unfair".
L'utilità nelle formule di integrazione Gaussiana non è la precisione assoluta numerica con ogni tipo di funzioni, per quello un Cavalieri-Simpson composito è l'ideale, come compromesso fra onere di calcolo e precisione numerica. La bellezza delle formule di integrazione Gaussiana è che la precisione è legata alla norma dell'energia delle funzioni che si integrano, questo conferisce a questo metodo delle particolari proprietà di convergenza. Infatti si lavora su spazi in cui le EDP trovano la loro formulazione più naturale e l'errore numerico è dello stesso tipo di quello indotto dall'approssimazione dei termini noti e dello spazio funzionale, tanto è vero che nell'applicazione ai metodi spettrali l'uso di formule di quadratura gaussiana, invece dell'integrazione esatta, non modifica l'ordine di convergenza teorico del metodo.
Tutto questo per sottolineare l'importanza di questa famiglia di metodi, che funziona molto bene se si lavora con funzioni regolari, nel senso degli spazi di Sobolev (con derivate limitate è il massimo, ma basta in effetti che la derivata distribuzionale sia quadrato-integrabile e la convergenza è ottima).
Infine io credo che un programma C++ scritto con criterio (non certo un codice Matlab con un ciclo for dentro!
) possa fare integrazione di Gauss-Chebishev tranquillamente su qualche milione di nodi in pochi secondi. Basta stare molto attenti a che la funzione integranda sia compilata inline all'interno dell'integratore dal compilatore... magie della meta-programmazione template!
L'utilità nelle formule di integrazione Gaussiana non è la precisione assoluta numerica con ogni tipo di funzioni, per quello un Cavalieri-Simpson composito è l'ideale, come compromesso fra onere di calcolo e precisione numerica. La bellezza delle formule di integrazione Gaussiana è che la precisione è legata alla norma dell'energia delle funzioni che si integrano, questo conferisce a questo metodo delle particolari proprietà di convergenza. Infatti si lavora su spazi in cui le EDP trovano la loro formulazione più naturale e l'errore numerico è dello stesso tipo di quello indotto dall'approssimazione dei termini noti e dello spazio funzionale, tanto è vero che nell'applicazione ai metodi spettrali l'uso di formule di quadratura gaussiana, invece dell'integrazione esatta, non modifica l'ordine di convergenza teorico del metodo.
Tutto questo per sottolineare l'importanza di questa famiglia di metodi, che funziona molto bene se si lavora con funzioni regolari, nel senso degli spazi di Sobolev (con derivate limitate è il massimo, ma basta in effetti che la derivata distribuzionale sia quadrato-integrabile e la convergenza è ottima).
Infine io credo che un programma C++ scritto con criterio (non certo un codice Matlab con un ciclo for dentro!
) possa fare integrazione di Gauss-Chebishev tranquillamente su qualche milione di nodi in pochi secondi. Basta stare molto attenti a che la funzione integranda sia compilata inline all'interno dell'integratore dal compilatore... magie della meta-programmazione template!
"david_e":
Infine io credo che un programma C++ scritto con criterio (non certo un codice Matlab con un ciclo for dentro!
Io credo che per essermi stato assegnato questo tipo di problema, vuol dire che è possibile creare anche codici matlab per la risoluzione.
Certo che è possibile farlo in Matlab, ma senza farlo con un ciclo for, ma sfruttando le operazioni vettoriali del Matlab. Fare cicli in un linguaggio interpretato è molto molto molto inefficiente. Bisogna vettorializzare il tutto. Ad esempio memorizzando tutti i nodi in un vettore e valutando la funzione sull'intero vettore, in questo modo si usano le routine intrinseche di Maltab che sono implementate in C e FORTRAN.
*** EDIT ***
Certamente il codice Matlab risulterà comunque molto più lento di un codice C++ scritto bene.
*** EDIT2 ***
Non è perversione usare i vettori per fare ogni operazione in Matlab: il programma risulta anche 100 volte più veloce alla fine!
*** EDIT ***
Certamente il codice Matlab risulterà comunque molto più lento di un codice C++ scritto bene.
*** EDIT2 ***
Non è perversione usare i vettori per fare ogni operazione in Matlab: il programma risulta anche 100 volte più veloce alla fine!
"david_e":
Certo che è possibile farlo in Matlab, ma senza farlo con un ciclo for, ma sfruttando le operazioni vettoriali del Matlab. Fare cicli in un linguaggio interpretato è molto molto molto inefficiente. Bisogna vettorializzare il tutto. Ad esempio memorizzando tutti i nodi in un vettore e valutando la funzione sull'intero vettore, in questo modo si usano le routine intrinseche di Maltab che sono implementate in C e FORTRAN.
Io avevo pensato di utilizzare un ciclo while che si fermasse solo nel caso in cui:
-l'errore sia minore della tolleranza richiesta dall'utente;
-il numero di nodi utilizzato fino a quel punto supera il numero di nodi che l'utente ha deciso di utilizzare.
Quindi qualcosa del tipo
While (abs(err>TOLL)) & (n
dove
- err = errore sull'approssimazione
- TOLL = tolleranza di errore
- n = numero di nodi utilizzati
- max = massimo numero di nodi da utilizzare
No il ciclo while non lo vedo problematico, perchè i cicli da evitare sono quelli del tipo:
in questo caso il Maltab perderà un sacco di tempo a interpretare il codice for e pochissimo tempo a fare la somma: in pratica il 99.999% del tempo lo passa sull'istruzione for ... end, non sulla parte di somma. Invece un codice del tipo:
non spreca tantissime risorse, perchè il vero onere è dentro nell'inverti matrice. L'unica cosa da evitare è mettere un ciclo for che somma le valutazioni della funzione nei nodi, ma bisognarebbe fare una cosa di questo tipo:
il fatto di usare la funzione intrinseca "sum" (così mi pare si chiami), accelera tantissimo il codice. Anche il fatto di chiamare la funzione una volta sola invece che n è un grande vantaggio: meno operazioni si fanno in esplicito e più ci si rifà a funzioni implementate a basso livello meglio è. Chiamare sin 1000 volte è molto ma molto ma molto più lento che chiamare sin una volta sula su un vettore di 1000 elementi e salvare 1000 risultati.
Tutto questo è contrario a ogni logica del buon programmatore C, perchè si spreca tantissima memoria (è inutile memorizzare tutti i risultati...), ma è l'unico modo per spremere ogni MHz dalla tua macchina usando Matlab.
sum=0;
for i=1:1e6
sum=sum+1/i;
end
in questo caso il Maltab perderà un sacco di tempo a interpretare il codice for e pochissimo tempo a fare la somma: in pratica il 99.999% del tempo lo passa sull'istruzione for ... end, non sulla parte di somma. Invece un codice del tipo:
for i=1:1e6
inverti una matrice 10'000x10'000
end
non spreca tantissime risorse, perchè il vero onere è dentro nell'inverti matrice. L'unica cosa da evitare è mettere un ciclo for che somma le valutazioni della funzione nei nodi, ma bisognarebbe fare una cosa di questo tipo:
while (abs(err>TOLL)) & (n<max)
nodi=qualcosa per generare un vettore di n nodi;
risultati=f(nodi); % qui dentro ci saranno n valutazioni.
output=sum(risultati);
end
il fatto di usare la funzione intrinseca "sum" (così mi pare si chiami), accelera tantissimo il codice. Anche il fatto di chiamare la funzione una volta sola invece che n è un grande vantaggio: meno operazioni si fanno in esplicito e più ci si rifà a funzioni implementate a basso livello meglio è. Chiamare sin 1000 volte è molto ma molto ma molto più lento che chiamare sin una volta sula su un vettore di 1000 elementi e salvare 1000 risultati.
Tutto questo è contrario a ogni logica del buon programmatore C, perchè si spreca tantissima memoria (è inutile memorizzare tutti i risultati...), ma è l'unico modo per spremere ogni MHz dalla tua macchina usando Matlab.
"david_e":
bisognarebbe fare una cosa di questo tipo:
while (abs(err>TOLL)) & (n<max) nodi=qualcosa per generare un vettore di n nodi; risultati=f(nodi); % qui dentro ci saranno n valutazioni. output=sum(risultati); end
Si, io infatti ho scritto un codice di questo tipo:
%inizializzazioni varie di variabili
[x]=funzione che crea n nodi
[s]=funzione che esegue le somme (quelle della formula)
appr=(($pi/n$)*somma);
n=n + incremento il numero di nodi
while ((abs(errore)>TOLL) & (n<max))
[x]=creo un certo numero di nodi pari all'incremento precedente
%esegue le operazioni precedenti fino ad appr che questa volta chiamo appr1
%a questo punto dovrei calcolare l'errore commesso sull'approssimazione e ho pensato di fare qualcosa del tipo
err=abs(appr1-appr);
(ma credo che sia sbagliato)
end
%verifica causa di uscita dal while
Basta che le somme e le valutazioni le hai fatte vettorialmente nelle funzioni e non come cicli... anche la generazione dei nodi credo si possa fare senza ricorrere a dei cicli, ma usando le funzioni di creazione automatica dei vettori. Altrimenti il codice potrebbe avere qualche difficoltà con n grandi... credo che con i for più di 10'000 non si riesca in tempi accettabili.
Ad esempio recentemente ho dovuto fare un integrale in 8 dimensioni usando Monte Carlo. Con il codice "lineare" stile C con i vari cicli for si arrivava a fare circa 1000 campionamenti e impiegava diversi minuti a dare il risultato. Con il codice vettoriale faceva 5milioni di valutazioni in pochi istanti. L'unico problema è che mi riempiva tutta la memoria RAM, quindi oltre un certo tot il metodo vettoriale torna a diventare lentissimo, dato che il computer va in swap.
Ad esempio recentemente ho dovuto fare un integrale in 8 dimensioni usando Monte Carlo. Con il codice "lineare" stile C con i vari cicli for si arrivava a fare circa 1000 campionamenti e impiegava diversi minuti a dare il risultato. Con il codice vettoriale faceva 5milioni di valutazioni in pochi istanti. L'unico problema è che mi riempiva tutta la memoria RAM, quindi oltre un certo tot il metodo vettoriale torna a diventare lentissimo, dato che il computer va in swap.
Però il mio problema ancora esiste, cioè la stima dell'errore.
Facendo ogni volta la differenza di due approssimazioni successive, il programma si ferma quando queste sono vicine fra loro, ma non è detto che siano vicine alla soluzione vera dell'integrale.
Facendo ogni volta la differenza di due approssimazioni successive, il programma si ferma quando queste sono vicine fra loro, ma non è detto che siano vicine alla soluzione vera dell'integrale.
Io purtroppo non posso aiutarti su questo punto... tutto ciò che mi veniva in mente l'ho scritto, guardarei in letteratura se fossi in te.
"studentessa":
Però il mio problema ancora esiste, cioè la stima dell'errore.
Facendo ogni volta la differenza di due approssimazioni successive, il programma si ferma quando queste sono vicine fra loro, ma non è detto che siano vicine alla soluzione vera dell'integrale.
Ma su questo non ci puoi fare niente!
E' come pretendere che sia per forza convergente una serie (a termini non negativi) di cui sai solo che il termine generale converge a zero.
Esempio banale (numerico): se la soluzione approssimata è $\log n$, questa non converge da nessuna parte. Ma $\log (n+1) - \log n$ va a zero e quindi il tuo criterio d'arresto, per quanto piccolo lo metti, sarà sempre soddisfatto prima o poi.
Se non hai altre informazioni che ti permettano di avere una stima a priori dell'errore o che ti garantiscano comuqnue qualche informazione addizionale, questa tua "preoccupazione" è ineliminabile.
Esatto: infatti le stime a priori e a posteriori sono sempre costruite postulando certe proprietà dell'oggetto che si studia. In questo caso della funzione integranda. Senza ipotesi su questa non sò cosa ci sia come stima... più che altro ci saranno molti stimatori euristici, che funzionano in molti casi, ma sui quali non avrai alcuna garanzia.
L'avevo promesso... ecco qui un'implementazione C++ fatta usando un po' di metaprogrammazione:
E questo e' l'output sul mio MacBook:
come vedete ho un tempo di esecuzione "user" di meno di un secondo, su 5milioni di nodi. Ho compilato il codice usando il gcc 4.0.1, con ottimizzazione "-march=prescott" (che non e' il top sul macbook ma e' gia' qualche cosa), "-O3" e poi ho strippato il binario usando il comando UNIX "strip".
Dato che ho usato le estensioni SGI del C++ non credo proprio che questo codice compili su un compilatore che non sia il gcc. Sotto windows forse il bloodshed C++ (che credo usi il gcc) lo compila. Altrimenti usato un sistema operativo serio: Linux o BSD o OS-X o Solaris... Basta che sia UNIX!
*** EDIT ***
Bisognerebbe disabilitare il MathML in modalita' "code"...
#include <iostream>
#include <cmath>
#include <functional>
#include <ext/functional>
#define PI 3.14159265358979
using namespace std;
using namespace __gnu_cxx;
template <class T, class F> double integrate(T const & weight, F const & f, int N) {
double xi;
double ris=0;
for(int i=1;i<N+1;++i) {
xi=double(2*i-1)/double(2*N)*PI;
xi=cos(xi);
ris+=f(xi)/weight(xi);
}
return PI/N*ris;
}
struct Weight : public unary_function<double,double> {
inline double operator()(double x) const {
return 1.0/sqrt(1-x*x);
}
};
struct Integranda : public unary_function<double,double> {
inline double operator()(double x) const {
return log(1+x);
}
};
int main(){
Weight w;
Integranda f;
cout.precision(30);
double risultato=integrate(w,f,5000000);
double errore=abs(2*(log(2)-1)-risultato);
cout << "Risultato: " << risultato << endl;
cout << "Errore: " << errore << endl;
}
E questo e' l'output sul mio MacBook:
[davide@macbook temp]\$ time ./main Risultato: -0.613705638880640669263755171414 Errore: 5.31241717283137404592707753181e-13 real 0m0.874s user 0m0.863s sys 0m0.006s
come vedete ho un tempo di esecuzione "user" di meno di un secondo, su 5milioni di nodi. Ho compilato il codice usando il gcc 4.0.1, con ottimizzazione "-march=prescott" (che non e' il top sul macbook ma e' gia' qualche cosa), "-O3" e poi ho strippato il binario usando il comando UNIX "strip".
Dato che ho usato le estensioni SGI del C++ non credo proprio che questo codice compili su un compilatore che non sia il gcc. Sotto windows forse il bloodshed C++ (che credo usi il gcc) lo compila. Altrimenti usato un sistema operativo serio: Linux o BSD o OS-X o Solaris... Basta che sia UNIX!
*** EDIT ***
Bisognerebbe disabilitare il MathML in modalita' "code"...
Accedi a tutti gli appunti
Tutor AI: studia meglio e in meno tempo