Integraz. numerica 01

tony19
Questo topic è una propaggine del topic originale "semplice integrale"
originato qui su "matematicamente/università" il 28/02/2005 - 15:35:31 da "bracardi" col quesito

come si fa a fare il seguente integrale? "int(senx/x*dx)"

Dopo la conclusione che l'int. da 0 all'inf. = pi/2, il topic si è orientato su modi di calcolarlo numericamente, formando un ramo che abbiamo potato e trasferito in questo topic "tutto suo".

perchè trapiantarlo? così, solo per un esperimento di pulizia.


tony

Risposte
g.schgor1
Mi aggiungo al coro (anche se sono stonato
e lo risolvo numericamente).

Il risultato e' ca. 3,55. (Mi sarei aspettato PI)

tony19
- il msg originale in "semplice integrale" era del 4/3/05, 03:07:38 -

si parla del doppio dell'integrale da o a +inf:
quote:
... Il risultato e' ca. 3,55. (Mi sarei aspettato PI) [g.schgor]

e, con le tue stesse espressioni ma con altri mezzi (un loop in basic), viene circa PI (~3,13)
come mai?

tony

g.schgor1
Gia'. Me lo sono chiesto anch'io.
E nel frattempo ho concluso che la
particolare funzione ad oscillazioni
smorzate, non gradisce l'integrazione
rettangolare.
Ho quindi riprovato con queste varianti:
- data la simmetria della funzione, calcolo
solo la parte positiva di x (raddoppiando
poi alla fine)
- aumento il numero di iterazioni (N)
- ma soprattuttio uso l'algoritmo di
integrazione trapezioidale.
Il risultato ora torna.

tony19
- il msg originale in "semplice integrale" era del 5/3/05, 01:15:58 -

dicendo che mi era venuto correttamente circa PI "con le tue stesse espressioni",
avevo tralasciato un particolare forse non irrilevante:
nella mia routine, per evitar le noie di una divisione per zero, avevo
"messo un tappo" nel buco che la funzione ha per x=0, in questo modo:

IF xn = 0 THEN yn = 1 ELSE yn = SIN(xn) / xn

invece mi sembra che nella tua prova ad integrazione rettangolare tu lasci calcolare sin(0)/0;
come si comporta il tuo software con quei valori?
forse è qui l'inghippo, la causa del tuo abnorme 3,55 invece di 3,14;
prova a modificare quella routine in questo senso.

inoltre tu dici

quote:
... E nel frattempo ho concluso che la
particolare funzione ad oscillazioni
smorzate, non gradisce l'integrazione
rettangolare. ... [g.schgor]

su questo non sono affatto d'accordo: con le mie prove i risultati dei due metodi differiscono, ma non di molto;

in più c'è da dire che nella tua prima prova, applicando il tuo metodo a rettangoli ad una funzione pari e quindi simmetrica rispetto all'asse y, gli errori in eccesso che si verificavano a destra su un tratto di curva decrescente si compensavano largamente con quelli in difetto sul simmetrico tratto crescente (e viceversa), e il risultato era "di ottima qualità" (uno schizzo ti convince immediatamente).

con la funzione asimmetrica, ovviamente, le differenze tra i due metodi aumentano, ma, ripeto, non di molto.

son curioso di sentire il risultato del "tappo" su x=0 nella tua prima integrazione, tra -30 e + 30.
(se non viene come dovrebbe, il mistero si infittisce!)
tony

g.schgor1
A seguito delle osservazioni di tony, ho riesaminato
i programmi di integrazione numerica della fuzione in
questione. In effetti tony ha ragione.
Riprogrammando le 2 procedure, non ho trovato differenze
nei risultati (=3.136 in entrambi i casi, con un errore quindi
di circa il 2 per mille).
Non so proprio ricostruire da dove sia saltato fuori quel risultato
fasullo di 3.55 da me indicato nel primo caso (e che mi ha
indotto a provare con il secondo).
Me ne dispiace e mi scuso (pur se e’ stata un’occasione di
citare il metodo d’integrazione trapezia)
Anche la seconda osservazione e’ giusta, ma non e’ la causa
dell’errore precedente: per x=0, Mathcad calcola y=0 anziche’
=1, quindi va corretto assegnando y(3000)=1, ma questo
comporterebbe un errore dell’1/100 in meno (=3.126 , non 3.55).
Rimane il fatto che integrando da zero la sola parte positiva e poi
raddoppiando, si migliora la precisione.

Ringraziando quindi tony dell’attenta verifica, approfitto
dell’occasione per chiedergli cosa ne pensa lui dell’argomento
citato in “Matematica e Calcolatore”.
La mia intenzione era di fare un’indagine sull’impatto del
calcolatore nello studio della matematica (stato di diffusione
delle conoscenze sulle possibilita’ d’impiego e suo effettivo uso),
ma finora non ho raccolto granche’ e non capisco se sia utile
continuare.


G.Schgör

tony19
- il msg originale in "semplice integrale" era del 9/3/05, 00:00:30 -
- qui è stato largamente abbreviato -


grazie, g. schgoer, di aver rifatto la prova!
il mistero sul 3.55 si è infittito e, a questo punto, rimarrà tale.
peccato ! il mettere a fuoco un errore di programma (mio o di qualcun altro) è stata per me una passione (e un lavoro) di una vita.

comunque rimane qualcosa da sgranocchiare sull'argomento e la aggiungerei qui, ma ho l'impressione che tutto questo ramo "numerico" della discussione esuli dal tema principale del thread.
quindi proporrei, se sei d'accordo, di spostare in blocco tutti i nostri msg sull'integrazione numerica, togliendoli di qui (dove lascerei un semplice rimando) e mettendoli in un nuovo thread che intitolerei ad es. "integraz. numerica 01" (e che comincerebbe con un brevissimo msg introduttivo per spiegare da dove proviene e di che si tratta.
- - - - -
lo spostamento è stato fatto.
sgranocchierò al più presto

tony

g.schgor1
Bene, tony, adesso lo spostamento e' stato fatto.
Il mio tentativo di richiamare l'attenzione sulle possibilita'
d'impiego del calcolatore nella matematica (di cui
l'integrazione numerica e' piccola parte) rimane.
Non mi hai ancora risposto cosa ne pensi tu dell'argomento
(se vale o no la pena d'insistere).
Segnalo, con l'occasione, il nuovo topic "Massa e molla", in cui
tento di illustrare una simulazione di comportamento dinamico.

G.Schgör

tony19
caro g.schgoer, cerco di portare avanti l'esercitazione che avevamo impostata e riprendo i miei appunti del mese scorso;
comincio dall'aspetto della precisione del calcolo, che mi sta assai a cuore, come si vede da altri miei post sul forum.
e, forse, possiamo interessare qualche lettore (anche se rimane silenzioso): ho l'impressione che tra i liceali e i da_poco_ex_liceali ci sia una certa passione per il calcolo fatto a macchina, non ancora attutita dal fascino per la "teoria pura" che tanto avvince molti dei più grandicelli.

sono incuriosito dalla differenza tra il tuo risultato (=3.136) e quello che ottengo io; mi piacerebbe arrivare a risultati comuni, o per lo meno a giustificarne le differenze.
definisco le condizioni della prova:

la funzione: y=sin(x)/x per x#0; y=1 per x=0
l'intervallo d'integrazione:
x0=-30; xN=+30; dx=0,01; Nintervalli=6000

calcolo dell'integrale da x0 ad xN col metodo dei rettangoli
( gli N rettangoli sono: y(x0)*dx, ..., y(xN-1)*dx ):

ambiente software: variabili floating point "single precision" (4 byte)

il risult. "teorico" calcolato con derive a 24 cifre di precisione è
3.13351308006070222196744 (ci possiamo fidare?)

sottolineo, per il lettore distratto, che NON stiamo cercando di calcolare l'integrale da -inf a +inf;
quindi NON stiamo cercando di ottenere Pi;
stiamo calcolando da -30 a +30 (e, en passant, si può vedere (diagrammando) che tronchiamo la funz. circa a metà di semionde negative)

mio risultato con una procedurina in basic:
I = 3.133513
tuo risultato con mathcad: I = 3.136
possibile che tra il mio Basic e il tuo Mathcad ci sia una diff. di quasi 3 su 3000, l'1 per mille? troppo grossa!

la differenza è dovuta anche al fatto che tu riporti (ovviamente arrotondando) solo 3 decimali?
se sì, sciorina tutte le tue cifre; risolviamo questo aspetto, andiamo poi (se il tuo software lo consente) a "double precision"
ed eventualmente a "long double" (o "extended") fino ad esser soddisfatti della precisione dei nostri calcoli.

poi ho in mente altri aspetti da approfondire:
1 - il calcolo + raddoppio della sola metà destra (stesso dx oppure stesso N. di intervalli?)
2 - i calcoli rifatti col metodo dei trapezi.
3 - idem con Simpson
4 - ci sarebbe anche un certo Weddle ...
ma tutti questi trattiamoli dopo; facciamo una cosa per volta.

grazie a chiunque voglia contribuire

tony

g.schgor1
Caro tony, ti trasmetto il calcolo
numerico dell'integrale svolto in Mathcad,
con le cifre decimali da te richieste.

Faccio pero' notare che la precisione dipende
dall'entita' dell'ultimo termine, quindi le
cifre dopo la terza decimale non hanno alcun
significato.(Ho usato N=8000 termini)
Sara' la mia mentalita' troppo da ingegnere,
ma la precisione oltre il millesimo non mi
appassiona, anche se capisco l'importanza di
controllare l'efficienza dei vari metodi di
calcolo.
Sono cmq a disposizione per eventuali ulteriori
precisazioni

G.Schgör

tony19
grazie, g.sgoer, della pronta risposta al mio pigro intervento
però con la tua prova da -40 a +40 e 8000 intervalli, hai usato condizioni diverse da quelle che avevo chiesto, e che riporto qui
quote:

definisco le condizioni della prova:

la funzione: y=sin(x)/x per x#0; y=1 per x=0
l'intervallo d'integrazione:
x0=-30; xN=+30; dx=0,01; Nintervalli=6000

calcolo dell'integrale da x0 ad xN col metodo dei rettangoli
( gli N rettangoli sono: y(x0)*dx, ..., y(xN-1)*dx )

il risult. "teorico" calcolato con derive a 24 cifre di precisione è
3.13351308006070222196744 (ci possiamo fidare almeno delle prime 18?)

mio risultato con una procedurina in basic, in "single precision":
I = 3.133513
tuo risultato con mathcad: I = 3.136


come faccio a comparare il tuo risult. da 40+40 col mio da 30+30?
sii cortese, prova con i miei dati.

inoltre, scusa, non capisco il tuo
quote:
Faccio pero' notare che la precisione dipende
dall'entita' dell'ultimo termine, quindi le
cifre dopo la terza decimale non hanno alcun
significato.

questo varrebbe se cercassimo di calcolare l'integrale da -inf a +inf e troncassimo a +/- 30
ma, come dicevo,
"... NON stiamo cercando di calcolare l'integrale da -inf a +inf;
quindi NON stiamo cercando di ottenere Pi

mi aspetterei che due strumenti basati sullo stesso hardware dessero risultati moolto più simili

su quanto segue, poi,
quote:

Sara' la mia mentalita' troppo da ingegnere,
ma la precisione oltre il millesimo non mi
appassiona, anche se capisco l'importanza di
controllare l'efficienza dei vari metodi di
calcolo.

aggiungo che non solo trovo importante "tarare" uno strumento e definirne la "classe di precisione"; così, anche se il risultato finale di un calcolo mi serve con tre sole cifre significative, posso valutare se non sia affetto da errori che han pesantemente inquinato i risultati intermedi.

Inoltre, non tutti gli ingegneri vanno con la stessa precisione.
un teodolite "serio" arrivava ai miei tempi ad una approssimaz. di +/- 0,2 secondi d'arco, che su 90 gradi corrispondono a 0,2/3600/90 ~= 6*10^-7 !
Il manualetto (tavole) dei logaritmi e delle funz. trigonometriche a 5 decimali usato al liceo fu considerato insufficiente alle esercitazioni di topografia, e ci fecero usare il Bruhns a 7 decimali (*) (anche se è meglio non si sappia di come arronzavo io le misure con la stadia in quei piacevoli pomeriggi di primavera ! [:D])

grazie, alla prossima.
tony

(*)tutta carta da macero, oggi

g.schgor1
Caro tony, anch'io non ha nostalgia di quei
"pomeriggi di primavera" ed ho usato il calcolatore
fin dai suoi inizi.
E oggi mi batto perche' i giovani lo utilizzino di piu'
(ma non solo per "chattare").

Ecco qui il calcolo in Mathcad come mi hai richiesto:

N=6000 n=1...N dx=0.01 x(n)=n*dx-30
y(0)=0 y(n)=sin(x(n))/x(n) y(3000)=1
iy(0)=0 iy(n)=iy(n-1)+y(n)*dx
iy(N)=3.133513184053061

(ovviamente iy(N) e' il risultato finale dell'integrazione
e mi sembra corrispondere al tuo).

Credo quindi che possiamo fidarci dei calcolatori.
(e' una battuta, naturalmente!). Ciao.

G.Schgör

rocco.g1
ciao,
non ho letto i vari post con le ipotesi per la soluzione, tuttavia il mio prof. di analisi I mi disse che un tale integrale non è calcolabile...

cioè mi posso sbagliare... ma credo sia proprio quello...

GIOVANNI IL CHIMICO
Forse che non è calcolabile analiticamente...

tony19
x rocco.g (e, di rimbalzo, il chimico):
l'aspetto dell'integrabilità di questa funzione è stato esaurientemente trattato nel topic originale, di cui questo, come dicevo nel prologo, è una propaggine.

era:
"semplice integrale" originato qui su "matematicamente/università" il 28/02/2005 - 15:35:31 da "bracardi" col quesito

'come si fa a fare il seguente integrale? "int(senx/x*dx)" '

contiene, tra l'altro, una elegante soluzione di goblyn, basata sulla trasformata di Fourier;
appare in lista con l'ultima data: 10/03/2005.

consiglio di vederlo, ed eventualmente di prosegure là con questioni teoriche.

qui stiamo semplicemente facendo (come dice il titolo) una esercitazione di calcolo numerico, e la funzione su cui lavoriamo è solo per caso il sin(x)/x

x g.schgoer:
- ora saresti d'accordo di rimuovere da quel topic la coda di messaggi che abbiamo traslocato qui? a un tuo cenno provvedo a togliere i miei, in blocco.
- grazie del calcolo; ho ancora qualche domanda; proseguo domani

tony

tony19
Grazie, g.schgoer
con la tua prova hai mostrato che il tuo strumento ha una precisione almeno uguale a quella del mio basic in precisione "double".
Il mio risultato infatti, con quella precisione, è 3.13351 31840 5307 (differisce dal tuo per un 10^-15).
non ottengo apprezzabile vantaggio (in questo problema) dalla precisione "extended" (quella con variabili da 10 byte, 80 bit, coerente col formato hardware del co-processore numerico x87).
Butto quindi alle ortiche la precisione "single" in tutti gli esperimenti che seguiranno.

Ricordo ai passanti che la differenza (1*10^-7) tra questo nostro risultato e quello "campione" 3.13351 30800 60702 22196 744 calcolato con "derive" a 24 cifre
è dovuta all'imprecisione intrinseca nel metodo (rettangoli) e nel passo scelto (0,01), e non ad errori di macchina.

Ora, dopo aver verificato la taratura dei nostri strumenti, passo alla fase successiva di questa esercitazione;
vorrei approfondire una cosa di cui non sono convinto:
dicevi che "integrando da zero la sola parte positiva e poi raddoppiando, si migliora la precisione."
Non mi convince, anzi direi il contrario, per le considerazioni di simmetria fatte tempo fa.
E propongo la seguente prova:

stessa funzione: y=sin(x)/x per x#0; y=1 per x=0
intervallo d'integrazione dimezzato: x0=0; xN=+30, data la simmetria della funzione
stesso passo: dx=0.01
Nintervalli: dimezzato, =3000

calcolo del DOPPIO dell 'integrale da x0 a xN col metodo dei rettangoli
(gli N rettangoli sono: y(x0)*dx, ... , y(xN-1)*dx )

ambiente software: variabili floating point "double precision" (8 byte)
stesso risultato "campione", calcolato con derive a 24 cifre di precisione:
3.13351 30800 60702 22196 744

il mio risultato:
3.14384 25279 2776
diff. rspetto al campione: 1*10^-2
decisamente assai peggiore dell'integrale tra -30 e + 30 (che dava diff=1*10^-7) !!

con ciò ho portato acqua al mio mulino;
a te, o a qualche altro lettore, risulta qualcosa di diverso?

tony

p.s. come già ho accennato, esaurito questo capitolo passerò all'integrazione a trapezi.

g.schgor1
Scusa, con l'integrazione numerica se raddoppi il
risultato devi escludere il valore dello zero e
poi aggiungerlo una sola volta. Hai fatto cosi'?

tony19
no di sicuro!
coerentemente con la definizione che avevo data

"calcolo del DOPPIO dell 'integrale da x0 a xN col metodo dei rettangoli
(gli N rettangoli sono: y(x0)*dx, ... , y(xN-1)*dx )"


il primo intervallo è tutta la striscia tra x0 e x1, l'ultimo quella tra xN-1 e xN.

(come avevo fatto anche per l'esercizio precedente)

tony

g.schgor1
Gia'. Ma cosi' sommi 2 volte x(0)*dx!
(Anch'io l'ho fatto, ma io non cerco un'estrema precisione)

G.Schgör

tony19
quote:
Gia'. Ma cosi' sommi 2 volte x(0)*dx!
(Anch'io l'ho fatto, ma io non cerco un'estrema precisione) [g.schgor]

ottimo!,
ma allora cambiamo le carte in tavola all'integrale a rettangoli come lo definivo sopra:
I=y(x0)*dx+y(x1)*dx+...+y(xN-1)*dx
lì il termine y(xN) non viene affatto considerato;
se uso solo metà del termine in x0, mi tocca aggiungere metà del termine in xN (e in questo modo i rettangoli son quelli "a cavallo" dei punti scelti)

E' un modo lecito (e più preciso) di specificare l'integrale, basta mettersi d'accordo: e, se lo si sceglie, va applicato anche all'esercizio precedente, togliendo metà del rettangolo in x=-30 e aggiungendo metà di quello in x=+30 (che, vedi caso coincidono)

basta essed'accordo sul metodo; io mi sono attenuto a quello che avevo specificato e che credo sia quello generalmente adottato:
i segmenti sono x0_x1, x1_x2, xN-1_xN.

butto lì una bomba, che tenevo in serbo per il post immediatamente successivo:
il metodo a rettangoli "a cavallo" equivale al metodo a trapezi!

tony

g.schgor1
"il metodo a rettangoli "a cavallo" equivale al metodo a trapezi!"

Si (ma solo se consideri una funzione lineare nell'intervallo dx)!

PS Non so editare il "quote". Come si fa?

G.Schgör

Rispondi
Per rispondere a questa discussione devi prima effettuare il login.