# plotting discontinuous functions

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:

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? **

• I edited my answer to show how to plot your specific function using np.nans. Would you mind reviewing the existing answers and selecting a good one? Thanks! – astrojuanlu May 4 '14 at 14:39

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:

• Can you also get the filled-in and empty circles respectively? See the image in the question: there needs to be a filled circle at (0, 0.5), and empty ones at (0,0) and (0,1). – ShreevatsaR Jan 17 '15 at 16:45

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))


• Can you also get the filled-in and empty circles respectively? See the image in the question: there needs to be a filled circle at (0, 0.5), and empty ones at (0,0) and (0,1). – ShreevatsaR Jan 17 '15 at 16:45

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