Plot integral function with scipy and matplotlib

Question

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):
    return integrate.quad(f,0,x)

plot(X,F(X))

I also tried to vectorize the function F without success.

Answer 1

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):
        y,err = integrate.quad(f,0,val)
        res[i]=y
    return res

plt.plot(X,F(X))

Answer 2

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:
        return integrate.quad(f, 0, x)

Answer 3

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

Answer 4

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)