Computational Science Stack Exchange is a question and answer site for scientists using computers to solve scientific problems. It's 100% free, no registration required.

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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.

share|improve this question

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))
share|improve this answer
    
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:
        return integrate.quad(f, 0, x)
share|improve this answer
    
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 at 19:11

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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