Python integrazione con infinita

Question

Python integrazione con infinita

voti

0

A volte, ho soluzione sbagliata quando integro con infinite confini in Python. Ecco un semplice esempio per illustrare la mia confusione:

from scipy import *
import numpy as np
from scipy.integrate import quad

def integrand1(x):
    output = exp(-(x-1.0)**2.0)
    return output

def integrand2(x):
    output = exp(-(x-100.0)**2.0)
    return output

solution1 = quad(integrand1,-np.inf,np.inf)
print solution1

solution1 = quad(integrand2,-np.inf,np.inf)
print solution2

mentre l'uscita è:

(1.7724538509055159, 3.668332157626072e-11)
(0.0, 0.0)

Non capisco il motivo per cui il secondo integrale è sbagliato, mentre il primo evitare di essere sbagliata. Sarà grande per dirmi alcuni trucchi per gestire infinita in Python.

python

È pubblicato 03/05/2017 alle 09:05
dall'utente phchen

In altre lingue...

1 risposte

themiurge · Accepted Answer · 2017-05-03 12:29

Non v'è nulla di intrinsecamente sbagliato con il codice. I risultati ottenuti sono dovuti alla quadalgoritmo essendo un metodo approssimato cui accuratezza, da quello che ho raccolto facendo alcune prove, notevolmente dipende da dove si trova il punto medio dell'intervallo di integrazione, rispetto all'intervallo asse x dove l'integrando è significativamente diverso da 0.

Il punto medio dell'intervallo di integrazione in caso di un (-inf,+inf)intervallo è sempre 0 (si veda il commento al codice Fortran rilevante qui , a partire dalla riga 238), e (purtroppo) non può essere configurato. La vostra integrand2funzione è centrato su x = 100, che è troppo lontano dal quadpunto medio di un algoritmo per essere sufficientemente accurate.

Sarebbe bello per essere in grado di specificare il punto medio in caso di integrazione tra -infe +inf, ma la buona notizia è che è possibile implementare da soli con decoratori funzione. In primo luogo è necessario un wrapper per la funzione integranda, al fine di spostare arbitrariamente sull'asse x:

def shift_integrand(integrand, offset):
    def dec(x):
        return integrand(x - offset)
    return dec

Questo genera una nuova funzione in base a qualsiasi integrando ti piace, basta spostare lungo l'asse x in base al offsetparametro. Quindi, se si fa qualcosa di simile (con i tuoi integrand1e integrand2funzioni):

new_integrand1 = shift_integrand(integrand1, -1.0)
print new_integrand1(0.0)
new_integrand2 = shift_integrand(integrand2, -100.0)
print new_integrand2(0.0)

si ottiene:

1.0
1.0

Ora è necessario un altro wrapper per la quadfunzione in modo da essere in grado di passare un nuovo punto centrale:

def my_quad(func, a, b, midpoint=0.0, **kwargs):
    if midpoint != 0.0:
        func = shift_integrand(func, -midpoint)
    return quad(func, a, b, **kwargs)

Infine, sapendo dove i vostri integrandi sono centrate, si può chiamare così:

solution2 = my_quad(integrand2, -np.inf, np.inf, midpoint=100.0)
print solution2

Che produce:

(1.772453850905516, 1.4202639944499085e-08)