Problems when plotting integrated singular functions in matplotlib

Question

I want to plot the integral of an integral of (singular) function in matplotlib, but my code doesn't work. Mathematically I want this:

Code:

import numpy as np
from scipy import integrate
import matplotlib.pyplot as plt


def g(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,g(X)) ## Doesn't work

def f(x):
    res = np.zeros_like(x)
    for i,val in enumerate(x):
        y,err = integrate.quad(g,0,val)
        res[i]=y
    return res

plot(X,f(X)) ## Works

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)) ## Doesn't work

(Code is an adapted version of https://scicomp.stackexchange.com/a/21871/9417)

So I can't plot the original function g, because it says:

      5 
      6 def g(x):
----> 7     if (np.abs(x)<1e-10):
      8         res = x
      9     else:

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

And also I cannot plot the integrated integral-function as it says:

     17 def f(x):
     18     res = np.zeros_like(x)
---> 19     for i,val in enumerate(x):
     20         y,err = integrate.quad(g,0,val)
     21         res[i]=y

TypeError: 'float' object is not iterable

However plotting the first integral f works. How can I fix this? Is there a better method of plotting this with python without running into such problems?

Answer 1

There are a few errors in your code:

• When you run g(X), the parameter is an array, while in the inside of the function you treat X as a value. g(x) is not applied to every element of X, it is applied to the whole array X. This explains why The truth value of an array... error, because x is in fact the whole array X. You can solve it in 2 different ways:

  1. Mapping the function g to every value of X

     >>> y = map(X, g)  # y = np.asarray(map(X, g)) if you want a numpy array back
     >>> plot(X, y)
     

  y = map(X, g) can also be written using list comprehension as y = [g(x) for x in X].

  2. Vectorize the function g(x) to be applied to the whole array

     def g(X):
         r = X.copy() # create a copy of x
         mask = np.abs(r) < 1e-10   # create a mask with invalid values
         r[mask] *= np.sin(1.0 / r[mask])  # replace those values
         return r
     

  However, if you choose option 2, the function f(x) will fail, as it is calling g(x) with a single value, once for every value of X.

• Your function f(x) is made to work with arrays instead, as you loop over every element of X and apply g(x) to each of them. If you still want to apply g(x) to every element of X use option 1.

• Your function F(x) also works with whole arrays, as you loop through its elements. However, for every element you call f(x), which only admits lists/arrays as input (and you are giving a single number). I think F(x) is redundant as it does exactly the same as f(x).

With your definition, you can rewrite equations as follows:

def g(x):
    return x if np.abs(x) < 1e-10 else x * np.sin(1.0/x)

def f(x):
    return integrate.quad(g, 0, x)[0]

def F(x):
    return integrate.quad(f, 0, x)[0]

And then obtain the result by mapping the function for every x.

X = np.arange(0,0.5,0.01)
yg = map(g, X)
yf = map(f, X)
yF = map(F, X)

Plot results:

import matplotlib.pyplot as plt
X = np.arange(-0.5, 0.5, 0.01)
yf = map(f, X)
plt.plot(X, yf)
plt.show()

---

A Quicker version that allows limiting the number of points within the integral:

def g(x):
    return x if np.abs(x) < 1e-10 else x * np.sin(1.0/x)

def f(x, limit=50):
    return integrate.quad(g, 0, x, limit=limit)[0]

def F(x, limit=50):
    return integrate.quad(f, 0, x, args=limit, limit=limit)[0]

Then run:

X = np.arange(-0.5, 0.5, 0.01)
FY = [F(x, limit=10) for x in X]
plt.plot(X, FY)

Larger values of limit will result in better representation of the image, at the expenses of much larger running times.