在python中找到几条曲线的交点

Question

我有 4 组数据，尝试在这些数据上拟合曲线。对于每组数据，我使用如下代码来拟合曲线。

center_540=np.sort(center_540,axis=0)
y=center_540[itr,:]
x=[0.2 ,0.4, 0.6, 0.8, 1, 2, 3, 4, 5,6]
popt, pcov,info,msg, ier= curve_fit(objective, x, y,full_output=True)
# summarize the parameter values
a, b, c ,d,e,f= popt
# plot input vs output
pyplot.scatter(x, y)
# define a sequence of inputs between the smallest and largest known inputs
x_line = np.linspace(min(x), max(x), 100)
# calculate the output for the range
y_line = objective(x_line, a, b, c,d,e,f)
# create a line plot for the mapping function
pyplot.plot(x_line, y_line, '--', color='red')
pyplot.xlabel('Rate')
pyplot.ylabel('PSNR')
err = np.dot(info['fvec'], info['fvec'])
plt.title('interpolation error resolution 540p: err=%f' %err)
pyplot.show()

之后，我试图找到这些曲线的交点。为此，我使用了以下代码来显示两条曲线的交点：

find_root(0.0063 , -0.1260 , 1.0146,-4.2147,10.4717,14.9473,0.0179 , -0.3148 , 2.1038,-6.7798,11.6871,15.5829)#1080 vs 720
find_root(0.0063 , -0.1260 , 1.0146,-4.2147,10.4717,14.9473,0.0189 , -0.3332 , 2.2272,-7.1543,12.0520,15.3168)# 1080 vs 540
find_root(0.0063 , -0.1260 , 1.0146,-4.2147,10.4717,14.9473,0.0178 , -0.3181 , 2.1672,-7.0761,11.6511,15.8349)#1080 vs 360
find_root(0.0179 , -0.3148 , 2.1038,-6.7798,11.6871,15.5829,0.0189 , -0.3332 , 2.2272,-7.1543,12.0520,15.3168) # 720 vs 540
find_root( 0.0179 , -0.3148 , 2.1038,-6.7798,11.6871,15.5829,0.0178 , -0.3181 , 2.1672,-7.0761,11.6511,15.8349)#720 vs 360
find_root(0.0189 , -0.3332 , 2.2272,-7.1543,12.0520,15.3168,0.0178 , -0.3181 , 2.1672,-7.0761,11.6511,15.8349)# 540 vs 360

find_root函数中的数字是曲线拟合中多项式函数的值。 这是

find_root

函数定义：

def find_root(a1,b1,c1,d1,e1,f1,a2,b2,c2,d2,e2,f2):
    coeff=[a1-a2,b1-b2,c1-c2,d1-d2,e1-e2,f1-f2]
    print(np.roots(coeff))

但是当我使用

find_root

函数和第一个代码中的图表检查根时，我看到一些图表中不存在的根。问题是什么？

例如，在图中我们可以看到交集小于 1，但是当我使用 find_root 函数时，我得到了这些根，但在图中看不到一些真正的根

[ 7.19666881+0.j          4.11941032+2.23825508j  4.11941032-2.23825508j
  1.14334481+0.j         -0.30297219+0.j        ]
[ 7.32487047+0.j          4.13676725+2.34581369j  4.13676725-2.34581369j
  1.01965241+0.j         -0.17361292+0.j        ]
[ 7.48885498+0.j          4.24400787+2.72602922j  4.24400787-2.72602922j
  1.09680306+0.j         -0.36932595+0.j        ]
[8.56875045+0.j         4.40845674+3.11694138j 4.40845674-3.11694138j
 0.50716803+0.89895985j 0.50716803-0.89895985j]
[-47.61740723+0.j           7.2784646 +3.00665406j
   7.2784646 -3.00665406j   0.95450723+0.j
  -0.89402922+0.j        ]
[ 6.88889702+1.76503189j  6.88889702-1.76503189j -0.72213743+2.48208939j
 -0.72213743-2.48208939j  1.39375355+0.j        ]

Answer 1

我们有以下曲线：

import itertools
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize

c1080 = np.array([0.0063 , -0.1260 , 1.0146, -4.2147, 10.4717, 14.9473])
c720 = np.array([0.0179 , -0.3148 , 2.1038, -6.7798, 11.6871, 15.5829])
c540 = np.array([0.0189 , -0.3332 , 2.2272, -7.1543, 12.0520, 15.3168])
c360 = np.array([0.0178 , -0.3181 , 2.1672, -7.0761, 11.6511, 15.8349])

curves = [c1080, c720, c540, c360]
labels = ["1080p", "720p", "540p", "360p"]

我们可以使用

np.roots

（返回所有根，包括复数）来估计曲线交点：

np.roots(c1080 - c720)
# array([ 7.19666881+0.j        ,  4.11941032+2.23825508j,
#         4.11941032-2.23825508j,  1.14334481+0.j        ,  # <---- Selected root
#        -0.30297219+0.j        ])

或者

scipy.optimize.root

接近统一（只有真正的根靠近您正在寻找的交点。

def factory(c1, c2):
    def wrapped(x):
        return np.poly1d(c1 -c2)(x)
    return wrapped

optimize.root(factory(c1080, c720), x0=1.)
#    fjac: array([[-1.]])
#      fun: array([0.])
#  message: 'The solution converged.'
#     nfev: 7
#      qtf: array([-2.97761815e-13])
#        r: array([-1.40828365])
#   status: 1
#  success: True
#        x: array([1.14334481])  # <----- Selected root

现在我们可以评估所有组合：

fig, axe = plt.subplots()
for (c1, c2), (l1, l2) in zip(itertools.combinations(curves, 2), itertools.combinations(labels, 2)):
    objective = factory(c1, c2)
    roots = np.roots(c1 - c2)
    solution = optimize.root(objective, x0=1.)
    axe.plot(xlin, objective(xlin), label="%s vs %s" % (l1, l2))
    axe.scatter(solution.x, objective(solution.x))
axe.legend()
axe.grid()

我们看到 720p 与 540p 的组合不相交：

# [8.56875045+0.j         4.40845674+3.11694138j 4.40845674-3.11694138j
#  0.50716803+0.89895985j 0.50716803-0.89895985j]  # <----- No real root close to unity
#     fjac: array([[-1.]])
#      fun: array([0.15609229])
#  message: 'The iteration is not making good progress, as measured by the \n  improvement from the last ten iterations.'
#     nfev: 28
#      qtf: array([-0.15609231])
#        r: array([-0.00035679])
#   status: 5
#  success: False
#        x: array([0.69365609])   # <---- Not a root