# Plot integral function with scipy and matplotlib

I want to plot a numerical integral function of some function $f$ using scipy and matplotlib. How can I do this?

I tried the following but it didn't work (run with ipython %pylab):

import numpy as np
from scipy import integrate

def f(x):
return x*np.sin(1/x)

X = np.arange(-0.5,0.5,0.001)

#plot(X,f(X))

def F(x):

plot(X,F(X))


I also tried to vectorize the function F without success.

First of all, your function $x\sin(\frac{1}{x})$ is singular in $x=0$. You might want to add an if clause like this:

def f(x):
if abs(x) < 1e-10:
res = x
else:
res  = x*sin(1/x)


but this does hurt speed (masked arrays would be better).

The reason why your code doesn't work is because

• scipy.quad only accepts a single value as left and right boundary

• scipy.quad returns a tuple (val,err) containing the value of the integral and an estimate for the error, so you need to unpack it.

This code works:

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

def f(x):
if (np.abs(x)<1e-10):
res = x
else:
res = x*np.sin(1.0/x)
return res

X = np.arange(-0.5,0.5,0.001)

#plot(X,f(X))

def F(x):
res = np.zeros_like(x)
for i,val in enumerate(x):
res[i]=y
return res

plt.plot(X,F(X))

• Nice solution! You also could speed up f(x) by catching the ZeroDivisionError and avoiding the if clause. Jan 20 '16 at 20:02
• "%matplotlib inline" invalid syntax. May 4 '17 at 10:39
• @Alan Turing since the OP indicated that he used ipython, I used the %matplotlib magic to have the plot in-line. If you use vanilla python, you have to remove this line... May 4 '17 at 14:24

Replace the last line by

plot(X, [F(x)[0] for x in X])


That should do it.

Edit: you can define your function F as

def F(x):
try:
return [integrate.quad(f, 0, y) for y in x]
except TypeError:

• Thanks. That's a good workaround. Is there also an easy way to redefine the function such that the plot command works as in my post? Jan 20 '16 at 19:11
from scipy.integrate import cumtrapz
import numpy as np
import matplotlib.pyplot as plt

def f(x):
return x*np.sin(1/x) if abs(x) > 1e-10 else x

f = np.vectorize(f)

X = np.arange(-0.5, .5, 0.001)

fv = f(X)
plt.plot(fv)

F = cumtrapz(fv, x=X, initial=0)
plt.plot(F);


If you insist on using quad, a more efficient implementation would calculate the integrals over the segments of the subdivision with quad for best accuracy and then a cumulative sum for the anti-derivative value at each sample point.

def f(x):
return x*np.sin(1/x)

X = np.arange(-0.5,0.5,0.001)

DF = [ integrate.quad(f,a,b)[0] for a,b in zip(X[:-1],X[1:]) ]

F = np.cumsum(DF)

plot(X[1:],F)