Scipy integrate.quad nicht die erwarteten Werte Rückkehr

Question

Ich habe die folgende Funktion, die ich numerisch möchte integrieren Python verwenden,

scipy verwenden, ich diesen Code geschrieben haben:

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 

Für die Cerfc (a) -Funktion, verwende ich nur scipy.erfc, um die komplementäre Fehlerfunktion zu berechnen.

Also funktioniert diese Funktion wirklich gut für niedrige Werte von u, aber große Werte (über 60 ish) bricht den Code und ich ger sehr sehr kleine Zahlen. Zum Beispiel, wenn ich a = 0,01 und u = 200, ist das Ergebnis 1.134335928072937e-40, wo die wahre Antwort geben: 1.410526851411200e-007

Zusätzlich dazu, SciPy kehrt der Fehler für die Quad-Berechnung in einer ähnlichen Reihenfolge wie die Antwort. Ich bin wirklich hier ratlos und würde wirklich Hilfe zu schätzen wissen.

Dies ist für eine Hausaufgabe, aber es ist ein Physikkurs. Diese Berechnung ist also nur ein Schritt in einer größeren Frage in der Physik. Sie werden helfen, mich nicht zu betrügen, wenn Sie mir helfen :)

Answer 1

Nach dem Wikipedia-Artikel Voigt profile, die Voigt functions U (x, t) und V (x, t) kann im Hinblick auf die komplexen Faddeeva function ausgedrückt werden w (z):

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

die Voigt-Funktion H (a, u) als

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

hinsichtlich U (x, t) ausgedrückt werden (siehe auch die DLMF section on Voigt functions.)

scipy hat eine Implementierung der Faddeeva-Funktion in scipy.special.wofz. dass Verwendung, hier ist eine Implementierung der Voigt-Funktionen:

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 

Sie sagten, dass Sie diesen Wert von H kennen (a, u) 1.410526851411200e-007, wenn a = 0,01 und u = 200. Wir können überprüfen:

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

---

Die oben keine Antwort auf die Frage, warum Ihr Code funktioniert nicht, wenn u groß ist. Um quad erfolgreich zu verwenden, ist es immer eine gute Idee, ein gutes Verständnis Ihres Integranden zu haben. In Ihrem Fall, wenn u groß ist, leistet nur ein sehr kleines Intervall in der Nähe von x = u einen wesentlichen Beitrag zum Integral. quad erkennt dies nicht, daher fehlt ein großer Teil des Integrals und es wird ein zu kleiner Wert zurückgegeben.

Eine Möglichkeit, dies zu beheben, ist das points Argument von quad mit einem Punkt zu verwenden, die bis zum Endpunkt des Intervalls sehr nahe ist.Zum Beispiel habe ich den Anruf von quad geändert:

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]) 

Mit dieser Änderung ist hier, was Ihre Funktion gibt für voigt(0.01, 200):

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

Ich habe keine strenge Rechtfertigung für den Wert 0.999*u haben; das ist nur ein Punkt nahe genug am Ende des Intervalls, um eine vernünftige Antwort für u um 200 oder so zu geben. Eine weitere Untersuchung des Integranden könnte Ihnen eine bessere Wahl geben. (Zum Beispiel können Sie einen analytischen Ausdruck für die Lage des Maximums des Integra finden? Wenn ja, so viel besser sein würde als 0.999*u.)

Sie könnten auch versuchen, die Werte von epsabs und epsrel zwicken, aber in Meine wenigen Versuche, die points Argument hinzugefügt, machte den größten Einfluss.