Qualcuno può spiegare perché scipy.integrate.quad fornisce risultati diversi per intervalli ugualmente lunghi integrando sin(X)?

StackOverflow https://stackoverflow.com/questions/581186

Domanda

Sto cercando di integrare numericamente una funzione arbitraria (nota quando codice) nel mio programma usando metodi di integrazione numerica.Sto utilizzando Python 2.5.2 insieme al pacchetto di integrazione numerica di SciPy.Per farmi un'idea, ho deciso di provare a integrare sin(x) e ho osservato questo comportamento-

>>> from math import pi
>>> from scipy.integrate import quad
>>> from math import sin
>>> def integrand(x):
...     return sin(x)
... 
>>> quad(integrand, -pi, pi)
(0.0, 4.3998892617846002e-14)
>>> quad(integrand, 0, 2*pi)
(2.2579473462709165e-16, 4.3998892617846002e-14)

Trovo strano questo comportamento perché...
1.Nell'integrazione ordinaria, l'integrazione sull'intero ciclo dà zero.
2.Nell'integrazione numerica, questo (1) non è necessariamente il caso, perché potresti essere approssimato nell'area totale sotto la curva.

In ogni caso, supponendo che 1 sia vero o presupponendo che 2 sia vero, trovo che il comportamento sia incoerente.Entrambe le integrazioni (da -pi a pi greco e da 0 a 2*pi) dovrebbero restituire 0,0 (il primo valore nella tupla è il risultato e il secondo è l'errore) oppure restituire 2,257...

Qualcuno può spiegare perché questo sta accadendo?Si tratta davvero di un'incoerenza?Qualcuno può anche dirmi se mi manca qualcosa di veramente basilare sui metodi numerici?

In ogni caso, nella mia applicazione finale, intendo utilizzare il metodo sopra descritto per trovare la lunghezza dell'arco di una funzione.Se qualcuno ha esperienza in quest'area, per favore consigliami la migliore politica per farlo in Python.

Modificare
Nota
Ho già i primi valori differenziali in tutti i punti dell'intervallo memorizzati in un array.
L'errore attuale è tollerabile.
Nota finale

Ho letto Wikipedia su questo.Come ha sottolineato Dimitry, integrerò sqrt(1+diff(f(x), x)^2) per ottenere la lunghezza dell'arco.Quello che volevo chiedere era: esiste un'approssimazione migliore/migliore pratica (?)/un modo più veloce per farlo.Se è necessario più contesto, lo pubblicherò separatamente/pubblicherò il contesto qui, come desideri.

È stato utile?

Soluzione

IL quad function è una funzione di una vecchia libreria Fortran.Funziona a giudicare dalla planarità e dalla pendenza della funzione che sta integrando come trattare la dimensione del passo che utilizza per l'integrazione numerica al fine di massimizzare l'efficienza.Ciò significa che potresti ottenere risposte leggermente diverse da una regione all'altra anche se analiticamente sono le stesse.

Senza dubbio entrambe le integrazioni dovrebbero restituire zero.Restituire qualcosa che è 1/(10 trilioni) è abbastanza vicino allo zero!Le leggere differenze sono dovute al percorso quad si sta ribaltando sin e modificandone le dimensioni del passo.Per l'attività pianificata, quad sarà tutto ciò di cui hai bisogno.

MODIFICARE:Per quello che stai facendo, credo quad è ok.È veloce e abbastanza preciso.La mia affermazione finale è di usarlo con fiducia a meno che non trovi qualcosa che è andato davvero storto.Se non restituisce una risposta priva di senso, probabilmente funziona perfettamente.Nessun problema.

Altri suggerimenti

Penso che probabilmente sia la precisione della macchina poiché entrambe le risposte sono effettivamente zero.

Se vuoi una risposta dalla bocca del cavallo, pubblicherei questa domanda su forum di discussione scipy

Direi che un numero O(10^-14) è effettivamente zero.Qual è la tua tolleranza?

Potrebbe essere che l'algoritmo alla base del quad non sia il migliore.Potresti provare un altro metodo di integrazione e vedere se questo migliora le cose.Una Runge-Kutta del 5° ordine può essere una tecnica molto bella per scopi generali.

Potrebbe essere proprio la natura dei numeri in virgola mobile: "Ciò che ogni informatico dovrebbe sapere sull'aritmetica in virgola mobile" .

Questo output mi sembra corretto poiché qui hai una stima dell'errore assoluto.Il valore integrale di sin(x) dovrebbe infatti avere valore zero per l'intero periodo (qualsiasi intervallo di lunghezza 2*pi) sia nell'integrazione ordinaria che numerica e i risultati sono vicini a quel valore.
Per valutare la lunghezza dell'arco dovresti calcolare l'integrale per la funzione sqrt(1+diff(f(x), x)^2), dove diff(f(x), x) è la derivata di f(x).Guarda anche Lunghezza dell'arco

0.0 == 2.3e-16 (absolute error tolerance 4.4e-14)

Entrambe le risposte sono le stesse e corretto cioè zero entro la tolleranza data.

La differenza deriva dal fatto che sin(x)=-sin(-x) esattamente anche in precisione finita.Mentre la precisione finita fornisce solo sin(x)~sin(x+2*pi) approssimativamente.Certo sarebbe bello se quad fosse abbastanza intelligente da capirlo, ma in realtà non ha modo di sapere a priori che l'integrale sui due intervalli forniti sia equivalente o che il primo sia un risultato migliore.

Autorizzato sotto: CC-BY-SA insieme a attribuzione
Non affiliato a StackOverflow
scroll top