2
$\begingroup$

I've been trying to plot the following function in Python:

$H(a,u) = \frac{a}{\pi} \int_{-\infty}^{\infty}\frac{exp(-y^2)}{a^2 + (u - y)^2}dy $

But I keep receiving the following error:

File "//anaconda/lib/python3.5/site-packages/scipy/integrate/quadpack.py", line 383, in _quad
return _quadpack._qagie(func,bound,infbounds,args,full_output,epsabs,epsrel,limit)

error: Supplied function does not return a valid float.

Here is the code I am using:

import numpy as np
from scipy import integrate
import matplotlib.pyplot as plt

a = 0.1
pi = np.pi
u = np.linspace(-100,100,200)


def integrand(a,u,y):
    top = np.exp(-y**2)
    bottom = a**2 + (u - y)**2
    return top/bottom

def func(u):
    res = np.zeros_like(u)
    for i,val in enumerate(u):
        y,err = integrate.quad(integrand, -np.inf, np.inf, args=(a,u))
        res[i] = y
    return (a/pi)*res


plt.plot(u,func(u))

u is a series of values which represent a range of frequencies. I've always struggled to get integration to work in python. Is there a good direction anyone could recommend me to go in to get it to work?

Thanks

$\endgroup$
3
$\begingroup$

Have you seen the documentation of scipy.integrate.quad? I see at least two issues with your code:

  • the first argument in the definition of integrand should be the integration variable y.
  • the argument u in integrand is a scalar, yet you are providing a vector in the args option of integrate.quad.

If all you want to do is plot your function, the following snippet is compact and clean

import numpy as np
from scipy import integrate
import matplotlib.pyplot as plt

def Voigt(a, u):
    I = integrate.quad(lambda y: np.exp(-y**2)/(a**2 + (u - y)**2),-np.inf, np.inf)[0]

    return (a/np.pi)*I

a = 0.1
u_range = np.linspace(-100,100,200)

plt.plot(u_range, [Voigt(a, u) for u in u_range])

enter image description here

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ I hadn't actually. I had been on the documentation for scipy.integrate but never thought of going to the quad part. Foolish on my behalf. Thanks for this. I never knew that you could put a for loop inside of a plot function! $\endgroup$ – Daniel Apr 28 '17 at 13:17
1
$\begingroup$

Your code is mostly correct. Just change u to val when you call integrate.quad:

import numpy as np
from scipy import integrate
import matplotlib.pyplot as plt

a = 0.1
pi = np.pi
u = np.linspace(-100,100,200)


def integrand(a,u,y):
    top = np.exp(-y**2)
    bottom = a**2 + (u - y)**2
    return top/bottom

def func(u):
    res = np.zeros_like(u)
    for i,val in enumerate(u):
        y,err = integrate.quad(integrand, -np.inf, np.inf, args=(a,val))
        res[i] = y
    return (a/pi)*res

plt.plot(u,func(u))

As mentioned by @Stelios, passing u (an array) to integrand means it will return an array, while integrate.quad expects a scalar. Here, we just use the fact that you already run over every val of the array with enumerate.

| cite | improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.