Introduction

When two curves are known only through sampled x and y values, the intersection is not usually present as an exact data point. High-precision intersection finding therefore becomes a numerical problem: build continuous approximations, locate sign changes, and refine each candidate root carefully.

Turn Sampled Data Into Continuous Functions

A reliable workflow has three steps:

• interpolate each curve
• subtract one interpolated function from the other
• solve for the roots of the difference

If the samples are smooth and ordered by x, cubic splines are a strong default. They give a continuous function and continuous first derivatives without forcing a single global polynomial over the whole interval.

1import numpy as np
2from scipy.interpolate import CubicSpline
3from scipy.optimize import brentq
4
5x1 = np.array([0.0, 1.0, 2.0, 3.0, 4.0])
6y1 = np.array([1.0, 2.2, 2.8, 3.1, 3.0])
7
8x2 = np.array([0.0, 1.0, 2.0, 3.0, 4.0])
9y2 = np.array([3.2, 2.9, 2.4, 2.1, 1.8])
10
11curve1 = CubicSpline(x1, y1)
12curve2 = CubicSpline(x2, y2)
13
14def difference(x):
15    return curve1(x) - curve2(x)

The intersection points are exactly the x values where difference(x) equals zero.

Use Bracketing Before Root Finding

High precision does not come from guessing a root at random. It comes from first identifying a small interval where the sign changes. A sign change means the difference moved from positive to negative or the other way around.

1search_grid = np.linspace(max(x1.min(), x2.min()), min(x1.max(), x2.max()), 400)
2values = difference(search_grid)
3
4brackets = []
5for left, right, y_left, y_right in zip(
6    search_grid[:-1], search_grid[1:], values[:-1], values[1:]
7):
8    if y_left == 0:
9        brackets.append((left, left))
10    elif y_left * y_right < 0:
11        brackets.append((left, right))

Once you have brackets, scipy.optimize.brentq is a good choice. Brent's method is fast, robust, and designed for bracketed roots.

1roots = []
2for left, right in brackets:
3    if left == right:
4        root = left
5    else:
6        root = brentq(difference, left, right, xtol=1e-12, rtol=1e-12)
7    roots.append(root)
8
9for root in roots:
10    print(f"x = {root:.12f}, y = {curve1(root):.12f}")

This combination of spline interpolation and Brent root finding is usually accurate enough for scientific and engineering workflows built around double precision.

Handling Multiple Intersections

Two curves may cross several times. The search grid matters because it determines how many sign changes you can detect. If the grid is too coarse, narrow crossings may be missed.

A practical approach is to choose a grid density based on how quickly the curves vary. If the underlying signal oscillates rapidly, increase the number of search points or segment the domain more carefully.

It also helps to deduplicate very close roots, especially when the grid catches the same crossing more than once.

1unique_roots = []
2for root in sorted(roots):
3    if not unique_roots or abs(root - unique_roots[-1]) > 1e-9:
4        unique_roots.append(root)
5
6print(unique_roots)

Choose the Right Interpolator

Spline interpolation is not always the right tool. If the data is noisy, a smoothing fit may be better than an exact interpolant. Exact interpolation through noisy measurements can create wiggles that introduce fake intersections.

For noisy data, consider these alternatives:

• use a smoothing spline
• fit a parametric model if the physics is known
• pre-filter the data before solving for intersections

The method should match the meaning of the samples. If the data is measured from instruments, “high precision” should not imply pretending the measurements are more certain than they are.

A Complete Reusable Function

Wrapping the workflow into a function makes it easier to test.

1def find_intersections(x1, y1, x2, y2, grid_size=1000):
2    curve1 = CubicSpline(x1, y1)
3    curve2 = CubicSpline(x2, y2)
4
5    def diff(x):
6        return curve1(x) - curve2(x)
7
8    start = max(np.min(x1), np.min(x2))
9    end = min(np.max(x1), np.max(x2))
10    grid = np.linspace(start, end, grid_size)
11    values = diff(grid)
12
13    roots = []
14    for left, right, y_left, y_right in zip(grid[:-1], grid[1:], values[:-1], values[1:]):
15        if y_left * y_right < 0:
16            roots.append(brentq(diff, left, right, xtol=1e-12, rtol=1e-12))
17
18    return [(root, curve1(root)) for root in roots]

Because the function accepts raw arrays, it is easy to plug into a notebook or a batch analysis pipeline.

Common Pitfalls

A common mistake is interpolating unsorted data. Most interpolation routines expect x values in increasing order, so sort the samples first if needed.

Another problem is using a very high-order polynomial fit on the full interval. That often creates oscillations near the edges and shifts the apparent intersection points.

Many failed results also come from ignoring the overlap interval. You can only search for an intersection where both curves are defined.

Finally, do not confuse numerical precision with measurement accuracy. A root reported to twelve decimal places may still be physically uncertain if the sampled data is noisy.

Summary

• Build continuous functions from sampled data before looking for intersections.
• Use sign-change detection to bracket candidate roots.
• 'CubicSpline plus brentq is a robust high-precision combination.'
• Increase grid density when curves vary rapidly or cross many times.
• Match the interpolation method to the noise level and meaning of the data.