plotting discontinuous functions

Question

I need help plotting the Heaviside function:

Real analysis often involves constructing bizarre functions which are intuitively correct, but ultimately wrong. See the great book Counterexamples in Analysis.

---

I wrote a numpy script to plot the function $\displaystyle \sum_{n=1}^\infty \frac{1}{n^2} H\left(x - \frac{1}{n}\right)$ to illustrate some concepts in real analysis. This is in a response to another question on math.stackexchange:

• Demonstrating Continuity Properties of $f(x) = \sum_{n=1}^\infty [a_n H(x - x_n)]$

Here was is example with accompanying plot:

B = [ [1,n] for n in range(1,50) ]

N = 1000
y = np.zeros(N)
count = 1
for b in B:
    x = np.zeros(N)
    x[ np.arange(b[0]*1.0*N/b[1], N).astype(int)]=count**-2
    y += x
    count += 1

plt.plot(np.arange(0,1,1.0/N),y)
plt.axis("Equal")

** How can I remove the vertical lines and really show the discontinuity? **

Answer 1

Getting rid of the vertical lines in this plot would involve finding all the discontinuities and then plotting each segment in a separate call to plt.plot.

Another, less elegant, solution is to plot points instead of lines. Since the points are very close together (high value of N), they will look like they form a line. To do this, replace plt.plot(np.arange(0,1,1.0/N), y) with plt.plot(np.arange(0,1,1.0/N), y, '.k') (the .tells plot to use dots, k is for black). This is what you'll get:

Answer 2

As long as you know the exact positions of the discontinuities, you just have to set the jump positions to nan in x, y or both. You can set this manually in the desired positions or use some criteria - for example, you can use the np.diff function to calculate the difference between contiguous positions in an array. For the simple case of the Heaviside function:

def heaviside(x):
    """See http://stackoverflow.com/a/15122658/554319"""
    y = 0.5 * (np.sign(x) + 1)
    y[np.diff(y) >= 0.5] = np.nan
    return y

x = np.linspace(-1, 1)
plt.plot(x, heaviside(x))
plt.ylim(-1, 2)

Using this in your specific example:

def U(x, n):
    return sum([heaviside(x - 1 / i) / i ** 2 for i in range(1, n)])

N = 1000
x = np.linspace(0, 1, N)
plt.plot(x, U(x, 50))

Answer 3

You can ' hold ' the plot and plot each component separately. Loop over the different functions and plot.