Scipy integrate.quad not returning the expected va

I have the following function which I would like to numerically integrate using python,

Using scipy, I have written this 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

For the Cerfc(a) function, I just use scipy.erfc to calculate the complimentary error function.

So this function works really well for low values of u, however larges u values (beyond 60 ish) breaks the code and I ger very very small numbers. For example, if I enter a = 0.01 and u = 200, the result is 1.134335928072937e-40, where the true answer is: 1.410526851411200e−007

In addition to this, the error scipy return for the quad calculation is on a similar order to the answer. I'm really stumped here and would really appreciate help.

This is for a homework assignment but it's a physics course. So this calculation is just one step in a broader question in physics. You will not be helping me cheat if you help me :)

According to the wikipedia article Voigt profile, the Voigt functions U(x,t) and V(x,t) may be expressed in terms of the complex Faddeeva function w(z):

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

The Voigt function H(a,u) can be expressed in terms of U(x,t) as

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

(Also see the DLMF section on Voigt functions.)

scipy has an implementation of the Faddeeva function in scipy.special.wofz. Using that, here's an implementation of the Voigt functions:

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

You said that you know that value of H(a,u) is 1.410526851411200e−007 when a=0.01 and u=200. We can check:

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

The above doesn't answer the question of why your code doesn't work when u is large. To use quad successfully, it is always a good idea to have a good understanding of your integrand. In your case, when u is large, only a very small interval near x = u makes a significant contribution to the integral. quad doesn't detect this, so it misses a big part of the integral and returns a value that is too small.

One way to fix this is to use the points argument of quad with a point that is very close to the end point of the interval. For example, I changed the call of quad to:

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

With that change, here's what your function returns for voigt(0.01, 200):

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

I don't have a rigorous justification for the value 0.999*u; that is just a point close enough to the end of the interval to give a reasonable answer for u around 200 or so. Further investigation of the integrand could give you a better choice. (For example, can you find an analytical expression for the location of the maximum of the integrand? If so, that would be much better than 0.999*u.)

You could also try tweaking the values of epsabs and epsrel, but in my few experiments, adding the points argument made the biggest impact.