How to find the points of intersection of a line and multiple curves in Python?

Question

I have data represented in the figure.

enter image description here

The curves were extrapolated and I have a line whose equation is known. The equation of curves are unknown. Now, how do I find the points of intersection of this line with each of the curves?

The reproducible code:

import numpy as np
import matplotlib.pyplot as plt
from scipy import interpolate


x = np.array([[0.12, 0.11, 0.1, 0.09, 0.08],
              [0.13, 0.12, 0.11, 0.1, 0.09],
              [0.15, 0.14, 0.12, 0.11, 0.1],
              [0.17, 0.15, 0.14, 0.12, 0.11],
              [0.19, 0.17, 0.16, 0.14, 0.12],
              [0.22, 0.19, 0.17, 0.15, 0.13],
              [0.24, 0.22, 0.19, 0.16, 0.14],
              [0.27, 0.24, 0.21, 0.18, 0.15],
              [0.29, 0.26, 0.22, 0.19, 0.16]])

y = np.array([[71.64, 78.52, 84.91, 89.35, 97.58],
              [66.28, 73.67, 79.87, 85.36, 93.24],
              [61.48, 69.31, 75.36, 81.87, 89.35],
              [57.61, 65.75, 71.7, 79.1, 86.13],
              [55.12, 63.34, 69.32, 77.29, 83.88],
              [54.58, 62.54, 68.7, 76.72, 82.92],
              [56.58, 63.87, 70.3, 77.69, 83.53],
              [61.67, 67.79, 74.41, 80.43, 85.86],
              [70.08, 74.62, 80.93, 85.06, 89.84]])


x1 = np.linspace(0, 0.4, 100)
y1 = -100 * x1 + 100
plt.figure(figsize = (5.15,5.15))
plt.subplot(111)
for i in range(5):
    x_val = np.linspace(x[0, i] - 0.05, x[-1, i] + 0.05, 100)
    x_int = np.interp(x_val, x[:, i], y[:, i])
    poly = np.polyfit(x[:, i], y[:, i], deg=2)
    y_int  = np.polyval(poly, x_val)
    plt.plot(x[:, i], y[:, i], linestyle = '', marker = 'o')
    plt.plot(x_val, y_int, linestyle = ':', linewidth = 0.25, color =  'black')
plt.plot(x1, y1, linestyle = '-.', linewidth = 0.5, color =  'black')


plt.xlabel('X')
plt.ylabel('Y')
plt.show()

jme · Accepted Answer · 2015-04-28 17:27:20Z

17

We do know the equations of the curves. They are of the form a*x**2 + b*x + c, where a,b, and c are the elements of the vector returned by np.polyfit. Then we just need to find the roots of a quadratic equation in order to find the intersections:

def quadratic_intersections(p, q):
    """Given two quadratics p and q, determines the points of intersection"""
    x = np.roots(np.asarray(p) - np.asarray(q))
    y = np.polyval(p, x)
    return x, y

The above function isn't super robust: there doesn't need to be a real root, and it doesn't really check for that. You're free to make it better.

Anyways, we give quadratic_intersections two quadratics, and it returns the two points of intersection. Putting it into your code, we have:

x1 = np.linspace(0, 0.4, 100)
y1 = -100 * x1 + 100
plt.figure(figsize = (7,7))
plt.subplot(111)

plt.plot(x1, y1, linestyle = '-.', linewidth = 0.5, color =  'black')
for i in range(5):
    x_val = np.linspace(x[0, i] - 0.05, x[-1, i] + 0.05, 100)
    poly = np.polyfit(x[:, i], y[:, i], deg=2)
    y_int  = np.polyval(poly, x_val)
    plt.plot(x[:, i], y[:, i], linestyle = '', marker = 'o')
    plt.plot(x_val, y_int, linestyle = ':', linewidth = 0.25, color =  'black')
    ix = quadratic_intersections(poly, [0, -100, 100])
    plt.scatter(*ix, marker='x', color='black', s=40, linewidth=2)

plt.xlabel('X')
plt.ylabel('Y')
plt.xlim([0,.35])
plt.ylim([40,110])
plt.show()

This makes the following figure:

enter image description here

Now, if you don't know that the functions you are dealing with are polynomials, you can use the optimization tools in scipy.optimize to find the root. For example:

import scipy.optimize
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.style

matplotlib.style.use('fivethirtyeight')
%matplotlib inline

f = lambda x: np.cos(x) - x
g = np.sin

h = lambda x: f(x) - g(x)

x = np.linspace(0,3,100)
plt.plot(x, f(x), zorder=1)
plt.plot(x, g(x), zorder=1)

x_int = scipy.optimize.fsolve(h, 1.0)
y_int = f(x_int)

plt.scatter(x_int, y_int, marker='x', s=150, zorder=2, 
            linewidth=2, color='black')

plt.xlim([0,3])
plt.ylim([-4,2])

Which plots:

enter image description here

edited Apr 28, 2015 at 17:27

answered Apr 27, 2015 at 20:23

jme

20.8k6 gold badges44 silver badges40 bronze badges

Sign up to request clarification or add additional context in comments.

8 Comments

Tom Kurushingal Over a year ago

This is great, but in case I have an unknown function, what can be done?

user554319 Over a year ago

can you share your mpl config file? I find this plot very beautiful ^^

Gareth Rees Over a year ago

Might be better to return separate coordinate arrays from quadratic_intersections, for example x = np.roots(np.asarray(p) - q); y = np.polyval(p, x); return x, y.

jme Over a year ago

@Juanlu001 I'm just using the fivethirtyeight style file that comes with recent versions of matplotlib. That is, import matplotlib.style; matplot.style.use('fivethirtyeight').

jme Over a year ago

@nxkryptor What do you mean by "unknown function"? Which function is unknown? What is known about it?

|

Collectives™ on Stack Overflow

How to find the points of intersection of a line and multiple curves in Python?

1 Answer 1

8 Comments

Your Answer

Linked

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

8 Comments

Your Answer

Sign up or log in

Post as a guest

Linked

Related