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.

score 8 · Answer 1 · answered Jan 20 '16 at 19:35

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

Nice solution! You also could speed up f(x) by catching the ZeroDivisionError and avoiding the if clause. — Christoph Wehmeyer, Jan 20 '16 at 20:02
@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... — GertVdE, May 04 '17 at 14:24

Christoph Wehmeyer · Answer 2 · 2016-01-20T19:44:46.750

5

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)

edited Jan 20 '16 at 19:44

answered Jan 20 '16 at 18:53

Christoph Wehmeyer

308
1
2
6

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? – student Jan 20 '16 at 19:11

score 4 · Answer 3 · answered Jan 17 '21 at 08:40

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

score 2 · Answer 4 · answered Jan 17 '21 at 13:44

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)

Plot integral function with scipy and matplotlib

4 Answers4

Linked