5
$\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.

$\endgroup$
6
$\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))
$\endgroup$
  • $\begingroup$ Nice solution! You also could speed up f(x) by catching the ZeroDivisionError and avoiding the if clause. $\endgroup$ – Christoph Wehmeyer Jan 20 '16 at 20:02
  • $\begingroup$ "%matplotlib inline" invalid syntax. $\endgroup$ – Henry May 4 '17 at 10:39
  • $\begingroup$ @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... $\endgroup$ – GertVdE May 4 '17 at 14:24
4
$\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)
$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – student Jan 20 '16 at 19:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.