python - Scipy integrate.quad not returning the expected values -

i have following function numerically integrate using python,

using scipy, have written 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 cerfc(a) function, use scipy.erfc calculate complimentary error function.

so function works low values of u, larges u values (beyond 60 ish) breaks code , ger very small numbers. example, if enter = 0.01 , u = 200, result 1.134335928072937e-40, true answer is: 1.410526851411200e−007

in addition this, error scipy return quad calculation on similar order answer. i'm stumped here , appreciate help.

this homework assignment it's physics course. calculation 1 step in broader question in physics. not helping me cheat if me :)

according wikipedia article voigt profile, voigt functions u(x,t) , v(x,t) may expressed in terms of 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 expressed in terms of u(x,t) as

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

(also see dlmf section on voigt functions.)

scipy has implementation of faddeeva function in scipy.special.wofz. using that, here's implementation of voigt functions:

from __future__ import division  import numpy np scipy.special import wofz   _sqrtpi = np.sqrt(np.pi) _sqrtpi2 = _sqrtpi/2  def voigtuv(x, t):     """     voigt functions u(x,t) , v(x,t).      return value 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 know value of h(a,u) 1.410526851411200e−007 when a=0.01 , u=200. can check:

in [109]: voigth(0.01, 200) out[109]: 1.41052685142231e-07 

---

the above doesn't answer question of why code doesn't work when u large. use quad successfully, idea have understanding of integrand. in case, when u large, small interval near x = u makes significant contribution integral. quad doesn't detect this, misses big part of integral , returns value small.

one way fix use points argument of quad point close end point of interval. example, changed 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 change, here's function returns voigt(0.01, 200):

in [191]: voigt(0.01, 200) out[191]: 1.4105268514252487e-07 

i don't have rigorous justification value 0.999*u; point close enough end of interval give reasonable answer u around 200 or so. further investigation of integrand give better choice. (for example, can find analytical expression location of maximum of integrand? if so, better 0.999*u.)

you try tweaking values of epsabs , epsrel, in few experiments, adding points argument made biggest impact.