Python中两个图的交集，求x值
Intersection of two graphs in Python, find the x value
设 0 <= x <= 1。我有两列 f 和 g，长度分别为 5000。现在我绘制：
plt.plot(x, f, '-')
plt.plot(x, g, '*')
我想找到曲线相交的点'x'。我不想找到 f 和 g 的交集。
我可以简单地做到这一点：
set(f) & set(g)
您可以将 np.sign 与 np.diff 和 np.argwhere 结合使用以获得线交叉点的索引（在这种情况下，点是 [ 0, 149, 331, 448, 664, 743] ):
import numpy as np
import matplotlib.pyplot as plt
x = np.arange(0, 1000)
f = np.arange(0, 1000)
g = np.sin(np.arange(0, 10, 0.01) * 2) * 1000
plt.plot(x, f, '-')
plt.plot(x, g, '-')
idx = np.argwhere(np.diff(np.sign(f - g))).flatten()
plt.plot(x[idx], f[idx], 'ro')
plt.show()
首先计算f - g和使用np.sign的相应符号。应用 np.diff 显示符号变化的所有位置（例如，线交叉）。使用 np.argwhere 为我们提供了准确的索引。
可能有多个路口，你可以通过下面的列表理解找到每个路口的(x,y)点
intersections = [(x[i], f[i]) for i,_ in enumerate(zip(f,g)) if f[i] == g[i]]
举个简单的例子
>>> x = [1,2,3,4,5]
>>> f = [2,4,6,8,10]
>>> g = [10,8,6,4,2]
>>> [(x[i], f[i]) for i,_ in enumerate(zip(f,g)) if f[i] == g[i]]
[(3, 6)]
所以这在x = 3, y = 6找到了一个交点。请注意，如果您使用 float，则两个值可能不完全相等，因此您可以使用一些公差而不是 ==。
即使f和g相交，你也不能确定f[i]== g[i]对于整数i（交点可能发生在点之间）。
你应该像
那样测试
# detect intersection by change in sign of difference
d = f - g
for i in range(len(d) - 1):
if d[i] == 0. or d[i] * d[i + 1] < 0.:
# crossover at i
x_ = x[i]
对于数组 f 和 g，我们可以简单地执行以下操作：
np.pad(np.diff(np.array(f > g).astype(int)), (1,0), 'constant', constant_values = (0,))
这将给出所有交叉点的数组。每个1是从下到上的交叉，每个-1是从上到下的交叉。
好吧，我正在为两条大小不同且 x 值不同的曲线寻找一个 matplotlib。这是我想出的：
import numpy as np
import matplotlib.pyplot as plt
import sys
fig = plt.figure()
ax = fig.add_subplot(111)
# x1 = [1,2,3,4,5,6,7,8]
# y1 = [20,100,50,120,55,240,50,25]
# x2 = [3,4,5,6,7,8,9]
# y2 = [25,200,14,67,88,44,120]
x1=[1.4,2.1,3,5.9,8,9,12,15]
y1=[2.3,3.1,1,3.9,8,9,11,9]
x2=[1,2,3,4,6,8,9,12,14]
y2=[4,12,7,1,6.3,7,5,6,11]
ax.plot(x1, y1, color='lightblue',linewidth=3, marker='s')
ax.plot(x2, y2, color='darkgreen', marker='^')
y_lists = y1[:]
y_lists.extend(y2)
y_dist = max(y_lists)/200.0
x_lists = x1[:]
x_lists.extend(x2)
x_dist = max(x_lists)/900.0
division = 1000
x_begin = min(x1[0], x2[0]) # 3
x_end = max(x1[-1], x2[-1]) # 8
points1 = [t for t in zip(x1, y1) if x_begin<=t[0]<=x_end] # [(3, 50), (4, 120), (5, 55), (6, 240), (7, 50), (8, 25)]
points2 = [t for t in zip(x2, y2) if x_begin<=t[0]<=x_end] # [(3, 25), (4, 35), (5, 14), (6, 67), (7, 88), (8, 44)]
# print points1
# print points2
x_axis = np.linspace(x_begin, x_end, division)
idx = 0
id_px1 = 0
id_px2 = 0
x1_line = []
y1_line = []
x2_line = []
y2_line = []
xpoints = len(x_axis)
intersection = []
while idx < xpoints:
# Iterate over two line segments
x = x_axis[idx]
if id_px1>-1:
if x >= points1[id_px1][0] and id_px1<len(points1)-1:
y1_line = np.linspace(points1[id_px1][1], points1[id_px1+1][1], 1000) # 1.4 1.401 1.402 etc. bis 2.1
x1_line = np.linspace(points1[id_px1][0], points1[id_px1+1][0], 1000)
id_px1 = id_px1 + 1
if id_px1 == len(points1):
x1_line = []
y1_line = []
id_px1 = -1
if id_px2>-1:
if x >= points2[id_px2][0] and id_px2<len(points2)-1:
y2_line = np.linspace(points2[id_px2][1], points2[id_px2+1][1], 1000)
x2_line = np.linspace(points2[id_px2][0], points2[id_px2+1][0], 1000)
id_px2 = id_px2 + 1
if id_px2 == len(points2):
x2_line = []
y2_line = []
id_px2 = -1
if x1_line!=[] and y1_line!=[] and x2_line!=[] and y2_line!=[]:
i = 0
while abs(x-x1_line[i])>x_dist and i < len(x1_line)-1:
i = i + 1
y1_current = y1_line[i]
j = 0
while abs(x-x2_line[j])>x_dist and j < len(x2_line)-1:
j = j + 1
y2_current = y2_line[j]
if abs(y2_current-y1_current)<y_dist and i != len(x1_line) and j != len(x2_line):
ymax = max(y1_current, y2_current)
ymin = min(y1_current, y2_current)
xmax = max(x1_line[i], x2_line[j])
xmin = min(x1_line[i], x2_line[j])
intersection.append((x, ymin+(ymax-ymin)/2))
ax.plot(x, y1_current, 'ro') # Plot the cross point
idx += 1
print "intersection points", intersection
plt.show()
交叉点可能发生在点之间。让我们探索下面的示例。
import numpy as np
import matplotlib.pyplot as plt
xs=np.arange(0, 20)
y1=np.arange(0, 20)*2
y2=np.array([1, 1.5, 3, 8, 9, 20, 23, 21, 13, 23, 18, 20, 23, 24, 31, 28, 30, 33, 37, 36])
绘制上面的 2 条曲线及其交点，使用 idx 交点之前和之后的平均坐标作为交点，所有点都更接近第一条曲线.
idx=np.argwhere(np.diff(np.sign(y1 - y2 )) != 0).reshape(-1) + 0
plt.plot(xs, y1)
plt.plot(xs, y2)
for i in range(len(idx)):
plt.plot((xs[idx[i]]+xs[idx[i]+1])/2.,(y1[idx[i]]+y1[idx[i]+1])/2., 'ro')
plt.legend(['Y1', 'Y2'])
plt.show()
使用前后平均坐标作为交点，但对于 y1 和 y2 曲线通常 更接近真实交点
plt.plot(xs, y1)
plt.plot(xs, y2)
for i in range(len(idx)):
plt.plot((xs[idx[i]]+xs[idx[i]+1])/2.,(y1[idx[i]]+y1[idx[i]+1]+y2[idx[i]]+y2[idx[i]+1])/4., 'ro')
plt.legend(['Y1', 'Y2'])
plt.show()
为了更准确地估计交叉点，我们可以使用插值法。
这是一个解决方案：
- 适用于 N 维数据
- 使用欧几里德距离而不是仅仅在 y 轴上寻找交叉点
- 对于大量数据更有效（它查询 KD-tree，应该以对数时间而不是线性时间查询）。
- 您可以更改 KD 树查询中的
distance_upper_bound 以定义足够接近的程度。
- 如果需要，您可以同时查询具有多个点的KD树。注意：如果您需要一次查询数千个点，您可以通过查询 KD-tree with another KD-tree.
来获得显着的性能提升
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy.spatial import cKDTree
from scipy import interpolate
fig = plt.figure()
ax = fig.add_axes([0, 0, 1, 1], projection='3d')
ax.axis('off')
def upsample_coords(coord_list):
# s is smoothness, set to zero
# k is degree of the spline. setting to 1 for linear spline
tck, u = interpolate.splprep(coord_list, k=1, s=0.0)
upsampled_coords = interpolate.splev(np.linspace(0, 1, 100), tck)
return upsampled_coords
# target line
x_targ = [1, 2, 3, 4, 5, 6, 7, 8]
y_targ = [20, 100, 50, 120, 55, 240, 50, 25]
z_targ = [20, 100, 50, 120, 55, 240, 50, 25]
targ_upsampled = upsample_coords([x_targ, y_targ, z_targ])
targ_coords = np.column_stack(targ_upsampled)
# KD-tree for nearest neighbor search
targ_kdtree = cKDTree(targ_coords)
# line two
x2 = [3,4,5,6,7,8,9]
y2 = [25,35,14,67,88,44,120]
z2 = [25,35,14,67,88,44,120]
l2_upsampled = upsample_coords([x2, y2, z2])
l2_coords = np.column_stack(l2_upsampled)
# plot both lines
ax.plot(x_targ, y_targ, z_targ, color='black', linewidth=0.5)
ax.plot(x2, y2, z2, color='darkgreen', linewidth=0.5)
# find intersections
for i in range(len(l2_coords)):
if i == 0: # skip first, there is no previous point
continue
distance, close_index = targ_kdtree.query(l2_coords[i], distance_upper_bound=.5)
# strangely, points infinitely far away are somehow within the upper bound
if np.isinf(distance):
continue
# plot ground truth that was activated
_x, _y, _z = targ_kdtree.data[close_index]
ax.scatter(_x, _y, _z, 'gx')
_x2, _y2, _z2 = l2_coords[i]
ax.scatter(_x2, _y2, _z2, 'rx') # Plot the cross point
plt.show()
对于正在使用或即将使用的Shapely library for geometry-related computations, getting the intersection will be much easier. You just have to construct LineString from each line and get their intersection如下：
import numpy as np
import matplotlib.pyplot as plt
from shapely.geometry import LineString
x = np.arange(0, 1000)
f = np.arange(0, 1000)
g = np.sin(np.arange(0, 10, 0.01) * 2) * 1000
plt.plot(x, f)
plt.plot(x, g)
first_line = LineString(np.column_stack((x, f)))
second_line = LineString(np.column_stack((x, g)))
intersection = first_line.intersection(second_line)
if intersection.geom_type == 'MultiPoint':
plt.plot(*LineString(intersection).xy, 'o')
elif intersection.geom_type == 'Point':
plt.plot(*intersection.xy, 'o')
要获取 x 和 y 值作为 NumPy 数组，您只需编写：
x, y = LineString(intersection).xy
# x: array('d', [0.0, 149.5724669847373, 331.02906176584617, 448.01182730277833, 664.6733061190541, 743.4822641140581])
# y: array('d', [0.0, 149.5724669847373, 331.02906176584617, 448.01182730277833, 664.6733061190541, 743.4822641140581])
或者如果一个交点只有一个点：
x, y = intersection.xy
我有一个类似的问题，但是有一个中断函数，比如正切函数。为了避免在不连续点上得到分数， 我不想考虑交叉点，我在以前使用 np.diff 和 np.sign 的解决方案上添加了一个公差参数。我将公差参数设置为两个数据点之间差异的平均值， 就足够了。
import numpy as np
import matplotlib.pyplot as plt
fig,ax = plt.subplots(nrows = 1,ncols = 2)
x = np.arange(0, 1000)
f = 2*np.arange(0, 1000)
g = np.tan(np.arange(0, 10, 0.01) * 2) * 1000
#here we set a threshold to decide if we will consider that point as a intersection
tolerance = np.abs(np.diff(f-g)).mean()
idx = np.argwhere((np.diff(np.sign(f - g)) != 0) & (np.abs(np.diff(f-g)) <= tolerance)).flatten()
#general case (tolerance = infinity)
tolerance = np.inf
idx2 = np.argwhere((np.diff(np.sign(f - g)) != 0) & (np.abs(np.diff(f-g)) <= tolerance)).flatten()
ax1,ax2 = ax
ax1.plot(x,f); ax1.plot(x,g)
ax2.plot(x,f); ax2.plot(x,g)
ax1.plot(x[idx], f[idx], 'o'); ax1.set_ylim(-3000,3000)
ax2.plot(x[idx2],f[idx2], 'o'); ax2.set_ylim(-3000,3000)
plt.show()
作为一个记录和测试的函数（该算法归功于@Matt，我只是将示例更改为更简单的内容并使用 linspace 而不是 arange 来处理 non-integers更好）：
from typing import Iterable, Tuple
import numpy as np
import doctest
def intersect(x: np.array, f: np.array, g: np.array) -> Iterable[Tuple[(int, int)]]:
"""
Finds the intersection points between `f` and `g` on the domain `x`.
Given:
- `x`: The discretized domain.
- `f`: The discretized values of the first function calculated on the
discretized domain.
- `g`: The discretized values of the second function calculated on the
discretized domain.
Returns:
An iterable containing the (x,y) points of intersection.
Test case, line-parabola intersection:
>>> x = np.linspace(0, 10, num=10000)
>>> f = 3 * x
>>> g = np.square(x)
>>> list(intersect(x, f, g))
[(0.0, 0.0), (2.999299929992999, 8.997899789978998)]
"""
idx = np.argwhere(np.diff(np.sign(f - g))).flatten()
return zip(x[idx], f[idx])
if __name__ == "__main__":
doctest.testmod()
在Python2中，只需删除类型提示。
设 0 <= x <= 1。我有两列 f 和 g，长度分别为 5000。现在我绘制：
plt.plot(x, f, '-')
plt.plot(x, g, '*')
我想找到曲线相交的点'x'。我不想找到 f 和 g 的交集。 我可以简单地做到这一点：
set(f) & set(g)
您可以将 np.sign 与 np.diff 和 np.argwhere 结合使用以获得线交叉点的索引（在这种情况下，点是 [ 0, 149, 331, 448, 664, 743] ):
import numpy as np
import matplotlib.pyplot as plt
x = np.arange(0, 1000)
f = np.arange(0, 1000)
g = np.sin(np.arange(0, 10, 0.01) * 2) * 1000
plt.plot(x, f, '-')
plt.plot(x, g, '-')
idx = np.argwhere(np.diff(np.sign(f - g))).flatten()
plt.plot(x[idx], f[idx], 'ro')
plt.show()
首先计算f - g和使用np.sign的相应符号。应用 np.diff 显示符号变化的所有位置（例如，线交叉）。使用 np.argwhere 为我们提供了准确的索引。
可能有多个路口，你可以通过下面的列表理解找到每个路口的(x,y)点
intersections = [(x[i], f[i]) for i,_ in enumerate(zip(f,g)) if f[i] == g[i]]
举个简单的例子
>>> x = [1,2,3,4,5]
>>> f = [2,4,6,8,10]
>>> g = [10,8,6,4,2]
>>> [(x[i], f[i]) for i,_ in enumerate(zip(f,g)) if f[i] == g[i]]
[(3, 6)]
所以这在x = 3, y = 6找到了一个交点。请注意，如果您使用 float，则两个值可能不完全相等，因此您可以使用一些公差而不是 ==。
即使f和g相交，你也不能确定f[i]== g[i]对于整数i（交点可能发生在点之间）。
你应该像
那样测试# detect intersection by change in sign of difference
d = f - g
for i in range(len(d) - 1):
if d[i] == 0. or d[i] * d[i + 1] < 0.:
# crossover at i
x_ = x[i]
对于数组 f 和 g，我们可以简单地执行以下操作：
np.pad(np.diff(np.array(f > g).astype(int)), (1,0), 'constant', constant_values = (0,))
这将给出所有交叉点的数组。每个1是从下到上的交叉，每个-1是从上到下的交叉。
好吧，我正在为两条大小不同且 x 值不同的曲线寻找一个 matplotlib。这是我想出的：
import numpy as np
import matplotlib.pyplot as plt
import sys
fig = plt.figure()
ax = fig.add_subplot(111)
# x1 = [1,2,3,4,5,6,7,8]
# y1 = [20,100,50,120,55,240,50,25]
# x2 = [3,4,5,6,7,8,9]
# y2 = [25,200,14,67,88,44,120]
x1=[1.4,2.1,3,5.9,8,9,12,15]
y1=[2.3,3.1,1,3.9,8,9,11,9]
x2=[1,2,3,4,6,8,9,12,14]
y2=[4,12,7,1,6.3,7,5,6,11]
ax.plot(x1, y1, color='lightblue',linewidth=3, marker='s')
ax.plot(x2, y2, color='darkgreen', marker='^')
y_lists = y1[:]
y_lists.extend(y2)
y_dist = max(y_lists)/200.0
x_lists = x1[:]
x_lists.extend(x2)
x_dist = max(x_lists)/900.0
division = 1000
x_begin = min(x1[0], x2[0]) # 3
x_end = max(x1[-1], x2[-1]) # 8
points1 = [t for t in zip(x1, y1) if x_begin<=t[0]<=x_end] # [(3, 50), (4, 120), (5, 55), (6, 240), (7, 50), (8, 25)]
points2 = [t for t in zip(x2, y2) if x_begin<=t[0]<=x_end] # [(3, 25), (4, 35), (5, 14), (6, 67), (7, 88), (8, 44)]
# print points1
# print points2
x_axis = np.linspace(x_begin, x_end, division)
idx = 0
id_px1 = 0
id_px2 = 0
x1_line = []
y1_line = []
x2_line = []
y2_line = []
xpoints = len(x_axis)
intersection = []
while idx < xpoints:
# Iterate over two line segments
x = x_axis[idx]
if id_px1>-1:
if x >= points1[id_px1][0] and id_px1<len(points1)-1:
y1_line = np.linspace(points1[id_px1][1], points1[id_px1+1][1], 1000) # 1.4 1.401 1.402 etc. bis 2.1
x1_line = np.linspace(points1[id_px1][0], points1[id_px1+1][0], 1000)
id_px1 = id_px1 + 1
if id_px1 == len(points1):
x1_line = []
y1_line = []
id_px1 = -1
if id_px2>-1:
if x >= points2[id_px2][0] and id_px2<len(points2)-1:
y2_line = np.linspace(points2[id_px2][1], points2[id_px2+1][1], 1000)
x2_line = np.linspace(points2[id_px2][0], points2[id_px2+1][0], 1000)
id_px2 = id_px2 + 1
if id_px2 == len(points2):
x2_line = []
y2_line = []
id_px2 = -1
if x1_line!=[] and y1_line!=[] and x2_line!=[] and y2_line!=[]:
i = 0
while abs(x-x1_line[i])>x_dist and i < len(x1_line)-1:
i = i + 1
y1_current = y1_line[i]
j = 0
while abs(x-x2_line[j])>x_dist and j < len(x2_line)-1:
j = j + 1
y2_current = y2_line[j]
if abs(y2_current-y1_current)<y_dist and i != len(x1_line) and j != len(x2_line):
ymax = max(y1_current, y2_current)
ymin = min(y1_current, y2_current)
xmax = max(x1_line[i], x2_line[j])
xmin = min(x1_line[i], x2_line[j])
intersection.append((x, ymin+(ymax-ymin)/2))
ax.plot(x, y1_current, 'ro') # Plot the cross point
idx += 1
print "intersection points", intersection
plt.show()
交叉点可能发生在点之间。让我们探索下面的示例。
import numpy as np
import matplotlib.pyplot as plt
xs=np.arange(0, 20)
y1=np.arange(0, 20)*2
y2=np.array([1, 1.5, 3, 8, 9, 20, 23, 21, 13, 23, 18, 20, 23, 24, 31, 28, 30, 33, 37, 36])
绘制上面的 2 条曲线及其交点，使用 idx 交点之前和之后的平均坐标作为交点，所有点都更接近第一条曲线.
idx=np.argwhere(np.diff(np.sign(y1 - y2 )) != 0).reshape(-1) + 0
plt.plot(xs, y1)
plt.plot(xs, y2)
for i in range(len(idx)):
plt.plot((xs[idx[i]]+xs[idx[i]+1])/2.,(y1[idx[i]]+y1[idx[i]+1])/2., 'ro')
plt.legend(['Y1', 'Y2'])
plt.show()
使用前后平均坐标作为交点，但对于 y1 和 y2 曲线通常 更接近真实交点
plt.plot(xs, y1)
plt.plot(xs, y2)
for i in range(len(idx)):
plt.plot((xs[idx[i]]+xs[idx[i]+1])/2.,(y1[idx[i]]+y1[idx[i]+1]+y2[idx[i]]+y2[idx[i]+1])/4., 'ro')
plt.legend(['Y1', 'Y2'])
plt.show()
为了更准确地估计交叉点，我们可以使用插值法。
这是一个解决方案：
- 适用于 N 维数据
- 使用欧几里德距离而不是仅仅在 y 轴上寻找交叉点
- 对于大量数据更有效（它查询 KD-tree，应该以对数时间而不是线性时间查询）。
- 您可以更改 KD 树查询中的
distance_upper_bound以定义足够接近的程度。 - 如果需要，您可以同时查询具有多个点的KD树。注意：如果您需要一次查询数千个点，您可以通过查询 KD-tree with another KD-tree. 来获得显着的性能提升
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy.spatial import cKDTree
from scipy import interpolate
fig = plt.figure()
ax = fig.add_axes([0, 0, 1, 1], projection='3d')
ax.axis('off')
def upsample_coords(coord_list):
# s is smoothness, set to zero
# k is degree of the spline. setting to 1 for linear spline
tck, u = interpolate.splprep(coord_list, k=1, s=0.0)
upsampled_coords = interpolate.splev(np.linspace(0, 1, 100), tck)
return upsampled_coords
# target line
x_targ = [1, 2, 3, 4, 5, 6, 7, 8]
y_targ = [20, 100, 50, 120, 55, 240, 50, 25]
z_targ = [20, 100, 50, 120, 55, 240, 50, 25]
targ_upsampled = upsample_coords([x_targ, y_targ, z_targ])
targ_coords = np.column_stack(targ_upsampled)
# KD-tree for nearest neighbor search
targ_kdtree = cKDTree(targ_coords)
# line two
x2 = [3,4,5,6,7,8,9]
y2 = [25,35,14,67,88,44,120]
z2 = [25,35,14,67,88,44,120]
l2_upsampled = upsample_coords([x2, y2, z2])
l2_coords = np.column_stack(l2_upsampled)
# plot both lines
ax.plot(x_targ, y_targ, z_targ, color='black', linewidth=0.5)
ax.plot(x2, y2, z2, color='darkgreen', linewidth=0.5)
# find intersections
for i in range(len(l2_coords)):
if i == 0: # skip first, there is no previous point
continue
distance, close_index = targ_kdtree.query(l2_coords[i], distance_upper_bound=.5)
# strangely, points infinitely far away are somehow within the upper bound
if np.isinf(distance):
continue
# plot ground truth that was activated
_x, _y, _z = targ_kdtree.data[close_index]
ax.scatter(_x, _y, _z, 'gx')
_x2, _y2, _z2 = l2_coords[i]
ax.scatter(_x2, _y2, _z2, 'rx') # Plot the cross point
plt.show()
对于正在使用或即将使用的Shapely library for geometry-related computations, getting the intersection will be much easier. You just have to construct LineString from each line and get their intersection如下：
import numpy as np
import matplotlib.pyplot as plt
from shapely.geometry import LineString
x = np.arange(0, 1000)
f = np.arange(0, 1000)
g = np.sin(np.arange(0, 10, 0.01) * 2) * 1000
plt.plot(x, f)
plt.plot(x, g)
first_line = LineString(np.column_stack((x, f)))
second_line = LineString(np.column_stack((x, g)))
intersection = first_line.intersection(second_line)
if intersection.geom_type == 'MultiPoint':
plt.plot(*LineString(intersection).xy, 'o')
elif intersection.geom_type == 'Point':
plt.plot(*intersection.xy, 'o')
要获取 x 和 y 值作为 NumPy 数组，您只需编写：
x, y = LineString(intersection).xy
# x: array('d', [0.0, 149.5724669847373, 331.02906176584617, 448.01182730277833, 664.6733061190541, 743.4822641140581])
# y: array('d', [0.0, 149.5724669847373, 331.02906176584617, 448.01182730277833, 664.6733061190541, 743.4822641140581])
或者如果一个交点只有一个点：
x, y = intersection.xy
我有一个类似的问题，但是有一个中断函数，比如正切函数。为了避免在不连续点上得到分数， 我不想考虑交叉点，我在以前使用 np.diff 和 np.sign 的解决方案上添加了一个公差参数。我将公差参数设置为两个数据点之间差异的平均值， 就足够了。
import numpy as np
import matplotlib.pyplot as plt
fig,ax = plt.subplots(nrows = 1,ncols = 2)
x = np.arange(0, 1000)
f = 2*np.arange(0, 1000)
g = np.tan(np.arange(0, 10, 0.01) * 2) * 1000
#here we set a threshold to decide if we will consider that point as a intersection
tolerance = np.abs(np.diff(f-g)).mean()
idx = np.argwhere((np.diff(np.sign(f - g)) != 0) & (np.abs(np.diff(f-g)) <= tolerance)).flatten()
#general case (tolerance = infinity)
tolerance = np.inf
idx2 = np.argwhere((np.diff(np.sign(f - g)) != 0) & (np.abs(np.diff(f-g)) <= tolerance)).flatten()
ax1,ax2 = ax
ax1.plot(x,f); ax1.plot(x,g)
ax2.plot(x,f); ax2.plot(x,g)
ax1.plot(x[idx], f[idx], 'o'); ax1.set_ylim(-3000,3000)
ax2.plot(x[idx2],f[idx2], 'o'); ax2.set_ylim(-3000,3000)
plt.show()
作为一个记录和测试的函数（该算法归功于@Matt，我只是将示例更改为更简单的内容并使用 linspace 而不是 arange 来处理 non-integers更好）：
from typing import Iterable, Tuple
import numpy as np
import doctest
def intersect(x: np.array, f: np.array, g: np.array) -> Iterable[Tuple[(int, int)]]:
"""
Finds the intersection points between `f` and `g` on the domain `x`.
Given:
- `x`: The discretized domain.
- `f`: The discretized values of the first function calculated on the
discretized domain.
- `g`: The discretized values of the second function calculated on the
discretized domain.
Returns:
An iterable containing the (x,y) points of intersection.
Test case, line-parabola intersection:
>>> x = np.linspace(0, 10, num=10000)
>>> f = 3 * x
>>> g = np.square(x)
>>> list(intersect(x, f, g))
[(0.0, 0.0), (2.999299929992999, 8.997899789978998)]
"""
idx = np.argwhere(np.diff(np.sign(f - g))).flatten()
return zip(x[idx], f[idx])
if __name__ == "__main__":
doctest.testmod()
在Python2中，只需删除类型提示。