Intersection points of curves
We often need to compute the points where two lines intersect or where a function intersects the x-axis (zero crossing). Here, we provide a Python script that does precisely that. Computes the points where two curves (lines, time series, etc.) intersect.
Hence, we provide a simple Python code that can compute the points $(x, y)$ where two curves, $ f(t) $ and g(t), intersect. The first step, and the easy one, is to take two points, one from each curve, $ (t_1, Y_f) $, and $ (t_1, Y_g) $, compute their distance, and check if it is smaller than an $ \epsilon $, where $ \epsilon $ is a tiny number close to zero. If that’s the case, then we are done and have found the intersection point being $ (t_1, Y_f) $. However, if the distance is larger than $ \epsilon $, we have to check for any crossing point. We first consider a new function $ h(t) = f(t) - g(t $, and we evaluate that function on the points $ (t_0, Y_f^0) $, $ (t_0, Y_g^0) $, $ (t_1, Y_f^1) $, and $ (t_1, Y_g^1) $. Thus, we obtain $ \delta_0 = Y_f^0 - Y_g^0 $ and $ \delta_1 = Y_f^1 - Y_g^1 $. If the product $ \delta_0 \delta_1 < 0 $, then from Bolzano’s theorem, we know that the two lines must intersect, and the intersection point $ (t^{\ast}, y^{\ast}) $ is:
$$ (t^{\ast}, y^{\ast}) = \Big( \frac{-t_0 \delta_1 + t_1 \delta_0}{\Delta}, \frac{Y_f^0 Y_g^1 - Y_f^1 Y_g^0}{\Delta} \Big), $$ where $ \Delta = Y_f^0 - Y_g^0 - Y_f^1 + Y_g^1 $. The coordinates of the point $ (t^{\ast}, y^{\ast}) $ are derived by computing the following determinant(s) (which gives us the intersection point):
When one runs the Python code provided in the above snippet, one obtains the results shown in the Figure below.
import numpy as np
import matplotlib.pylab as plt
import matplotlib.style as style
style.use('tableau-colorblind10')
params = {'font.size': 14,
}
plt.rcParams.update(params)
eps = 1e-12
def findCrossings(X, Y, t=None):
if len(X) != len(Y):
raise ValueError
if t is None:
t = np.arange(len(X))
crossings = []
deltaOld = 0
tPrev, xOld, yOld = None, None, None
for tCurr, xNew, yNew in zip(t, X, Y):
deltaNew = yNew - xNew
# When the difference of two points is virtually zero then accept that
# as a crossing point (x, y)
if np.abs(deltaNew) < eps:
crossings.append((tCurr, xNew))
# Otherwise check if the product of the differences is negative and then
# from Bolzano's theorem we know that there should be a crossing. To find
# the intersection we compute the determinant as it is defined in the text
# above and compute the crossing point (x, y)
elif deltaNew * deltaOld < 0:
denom = deltaOld - deltaNew
crossings.append(((-deltaNew * tPrev + deltaOld * tCurr) / denom,
(xNew * yOld - yNew * xOld) / denom))
tPrev, xOld, yOld, deltaOld = tCurr, xNew, yNew, deltaNew
return crossings
def nicePlot(t, x, y, crossing_pts, ax):
""" A simple plotting function
"""
ax.plot(t, x, 'b-', alpha=0.7, lw=2)
ax.plot(t, y, 'k-', alpha=0.7, lw=2)
ax.plot(*zip(*crossing_pts),
ls="",
color='orange',
marker='o',
alpha=1.,
ms=10)
ax.set_xticks([])
ax.set_yticks([])
if __name__ == "__main__":
N = 20
fig = plt.figure(figsize=(19, 5))
ax = fig.add_subplot(131)
# Two curves that do not intesect thus no crossing points
t = np.linspace(-np.pi, np.pi, N)
x = np.sin(t * .08+1.4)*np.random.uniform(0.5, 0.9) + 1
y = -np.cos(t * .07+.1)*np.random.uniform(0.7, 1.0) + 1
crossing_pts = findCrossings(x, y, t=t)
nicePlot(t, x, y, crossing_pts, ax)
# Two time series intesecting at multiple points
ax = fig.add_subplot(132)
t = np.arange(N)
x = np.random.normal(0, 1, (N,))
y = np.random.normal(0, 1, (N,))
crossing_pts = findCrossings(x, y, t=t)
nicePlot(t, x, y, crossing_pts, ax)
# x^2 intersecting at two points with the line x
ax = fig.add_subplot(133)
t = np.linspace(-1, 1, N)
x = t**2
y = t
crossing_pts = findCrossings(x, y, t=t)
nicePlot(t, x, y, crossing_pts, ax)
plt.show()
When one runs the Python code provided in the above snippet then they obtain the results shown in the Figure below.
Figure 1. Three exampels of identifying crossing points between two curves. The orange discs indicate the crossing points.