Scipy integrate.quad ne pas retourner les valeurs attendues

Question

je la fonction suivante que je voudrais intégrer numériquement en utilisant python,

En utilisant scipy, j'ai écrit ce code:

def voigt(a,u): 

fi = 1 
er = Cerfc(a)*np.exp(np.square(a)) 
c1 = np.exp(-np.square(u))*np.cos(2*a*u) 
c1 = c1*er #first constant term 
pis = np.sqrt(np.pi) 

c2 = 2./pis #second constant term  

integ = inter.quad(lambda x: np.exp(-(np.square(u)- 
np.square(x)))*np.sin(2*a*(u-x)), 0, u) 

print integ 
ing = c1+c2*integ[0] 

return ing 

Pour la fonction Cerfc (a), j'utilise simplement scipy.erfc pour calculer la fonction d'erreur complémentaire.

Cette fonction fonctionne donc très bien pour les valeurs faibles de u, mais les valeurs les plus grandes (au-delà de 60 ish) cassent le code et j'obtiens de très petits nombres. Par exemple, si j'entre a = 0.01 et u = 200, le résultat est 1.134335928072937e-40, où la vraie réponse est: 1.410526851411200e-007

En plus de cela, le retour d'erreur scipy pour le calcul quad est sur un ordre similaire à la réponse. Je suis vraiment perplexe ici et j'apprécierais vraiment l'aide.

Ceci est pour un devoir de devoirs mais c'est un cours de physique. Ce calcul n'est donc qu'une étape dans une question plus large en physique. Vous ne serez pas me contribuerez triche si vous me aider :)

Answer 1

Selon l'article de wikipedia Voigt profile, le Voigt functions U (x, t) et V (x, t) peut être exprimée en termes de Faddeeva function complexe w (z):

U(x,t) + i*V(x,t) = sqrt(pi/(4*t))*w(i*z) 

la fonction Voigt H (a, u) peuvent être exprimées en termes de U (x, t) comme

H(a,u) = U(u/a, 1/(4*a**2))/(a*sqrt(pi)) 

(voir également le DLMF section on Voigt functions.)

scipy a une implémentation de la fonction Faddeeva dans scipy.special.wofz. En utilisant cette , voici une implémentation des fonctions Voigt:

from __future__ import division 

import numpy as np 
from scipy.special import wofz 


_SQRTPI = np.sqrt(np.pi) 
_SQRTPI2 = _SQRTPI/2 

def voigtuv(x, t): 
    """ 
    Voigt functions U(x,t) and V(x,t). 

    The return value is U(x,t) + 1j*V(x,t). 
    """ 
    sqrtt = np.sqrt(t) 
    z = (1j + x)/(2*sqrtt)      
    w = wofz(z) * _SQRTPI2/sqrtt 
    return w 

def voigth(a, u): 
    """ 
    Voigt function H(a, u). 
    """ 
    x = u/a 
    t = 1/(4*a**2) 
    voigtU = voigtuv(x, t).real 
    h = voigtU/(a*_SQRTPI) 
    return h 

Vous avez dit que vous savez que la valeur de H (a, u) est 1.410526851411200e-007 quand a = 0,01 et u = 200. Nous pouvons vérifier:

In [109]: voigth(0.01, 200) 
Out[109]: 1.41052685142231e-07 

---

Ce qui précède ne répond pas à la question de savoir pourquoi votre code ne fonctionne pas lorsque u est grande. Pour utiliser quad avec succès, c'est toujours une bonne idée d'avoir une bonne compréhension de votre integrand. Dans votre cas, lorsque u est grand, seul un très petit intervalle proche de x = u apporte une contribution significative à l'intégrale. quad ne le détecte pas, il manque donc une grande partie de l'intégrale et renvoie une valeur trop petite. Un moyen de résoudre ce problème est d'utiliser l'argument points de quad avec un point très proche du point de fin de l'intervalle.Par exemple, j'ai changé l'appel de quad à:

integ = inter.quad(lambda x: np.exp(-(np.square(u)-np.square(x))) * np.sin(2*a*(u-x)), 
        0, u, points=[0.999*u]) 

Avec ce changement, voici ce que votre fonction retourne pour voigt(0.01, 200):

In [191]: voigt(0.01, 200) 
Out[191]: 1.4105268514252487e-07 

Je n'ai pas une justification rigoureuse de la valeur 0.999*u; C'est juste un point assez proche de la fin de l'intervalle pour donner une réponse raisonnable pour environ 200 ou plus. Une étude plus approfondie de l'intégrande pourrait vous donner un meilleur choix. (Par exemple, vous pouvez trouver une expression analytique pour l'emplacement du maximum de la integrand? Si oui, ce serait beaucoup mieux que 0.999*u.)

Vous pouvez également essayer de peaufiner les valeurs de epsabs et epsrel, mais mes quelques expériences, en ajoutant l'argument points a fait le plus grand impact.