# Plotting Voigt Function in Python

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

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

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

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


• 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! Apr 28 '17 at 13:17

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

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.