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. – Christoph Wehmeyer Jan 20 at 20:02

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: