$\begingroup$
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.
nicoguaro♦
8,1295 gold badges21 silver badges47 bronze badges
asked Jan 20, 2016 at 15:06
$\endgroup$
$\begingroup$
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))answered Jan 20, 2016 at 19:35
GertVdEGertVdE
6,0911 gold badge19 silver badges36 bronze badges
$\endgroup$
3
$\begingroup$
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)answered Jan 20, 2016 at 18:53
$\endgroup$
1
$\begingroup$
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);
answered Jan 17, 2021 at 8:40
$\endgroup$
$\begingroup$
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)answered Jan 17, 2021 at 13:44
Lutz LehmannLutz Lehmann
3,7111 gold badge12 silver badges17 bronze badges
$\endgroup$