How to deal with functions that approach infinity in NumPy?

The period of the sine function is rather higher than what is plotted, so what you’re seeing is aliasing from the difference in the sampling frequency and some multiple of the true frequency. Since one of the roots is at 0, the smallest discrepancy that happens to exist between the first few samples and a multiple of π itself scales linearly away from 0, producing a 1/x envelope.

In this example, input[5] is 5(1000/(200-1))=8π−0.007113, so the function is about −141 there, as shown. input[10] is of course 16π−0.014226, so that the function is about −70, and so on as long as the discrepancy is much smaller than π.

It’s possible for some one of the quasi-periodic sample sequences to eventually land even closer to nπ, producing a more complicated pattern like that in the second plot.

Why does one see this decreasing tendency even though the function is purely periodic?

Keep in mind that actually at every multiple of pi the function goes to infinity. And the size of the jump displayed actually only reflects the biggest value of the sampled values where the function still made sense. Therefore you get a big jump if you happen to sample a value were the function is big but not too big to be a float.

To be able to plot anything matplotlib throws away values that do not make sense. Like the np.nan you get at multiples of pi and ±np.infs you get for values very close to that. I believe what happens is that one step size away from zero you happen to get a value small enough not to be thrown away but still very large. While when you get to pi and multiples of it the largest value gets thrown away.

How can I make a plot like this consistent i.e. independent of step-size/step-amount?

You get strange behaviour around the values where your function becomes unreasonable large. Just pick a ylimit to avoid plotting those crazy large values.

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.transforms import Bbox

x = np.linspace(10**-10,50, 10**4)
plt.plot(x,1/np.sin(x))
plt.ylim((-15,15))