Find intersections numerically -- Using NumPy

Intersects of 2 functions

Assuming function $f$ and $g$, scipy.optimize.fsolve would give an elegant solution:

from scipy.optimize import fsolve
def f(xy):
    x, y = xy
    return (y - x ** 2, y - np.sin(x))
fsolve(f, [1.0, 1.0])

gives array([0.87672622, 0.76864886]) is the intersect of $y=x^2$ and $y=\sin(x)$.

Intersects of 2 arrays

There is a function called numpy.intersect1d that returns the sorted, uniq values in both arrays, and reduce(np.intersect1d, list_of_arrays) gives intersect of all arrays in the list. Index can also be returned by return_indices=True option.

If the 2 arrays are data values, using np.flatnonzero(np.abs(arr1-arr2))<tol) seems to be a good solution; this method returns the indices of values in arr1 and arr2 which are closed by tol.

Another method is based on idea that there is a change in sign of $\lim_{x\to x_0^+}f(x_0)-g(x_0)$ and $\lim_{x\to x_0^-}f(x_0)-g(x_0)$, therefore, if $\mathrm{d}x$ is small enough, this method would give a very good estimation of the intersects.

def intersects(x, arr1, arr2):
    return x[np.flatnonzero(np.diff(np.sign(arr1 - arr2))]