【问题标题】:Finding the intersection of a curve from polyfit从 polyfit 中找到曲线的交点
【发布时间】:2013-10-13 11:10:16
【问题描述】:

这看起来很简单,但我想不通。我有一条从 x,y 数据计算的曲线。然后我有一条线。我想找到两者相交处的 x、y 值。

这是我到目前为止所得到的。这非常令人困惑,并且没有给出正确的结果。我可以查看图表并找到交点 x 值并计算正确的 y 值。我想删除这个人为步骤。

import numpy as np
import matplotlib.pyplot as plt
from pylab import * 
from scipy import linalg
import sys
import scipy.interpolate as interpolate
import scipy.optimize as optimize

w = np.array([0.0, 11.11111111111111, 22.22222222222222, 33.333333333333336, 44.44444444444444, 55.55555555555556, 66.66666666666667, 77.77777777777777, 88.88888888888889, 100.0])
v = np.array([0.0, 8.333333333333332, 16.666666666666664, 25.0, 36.11111111111111, 47.22222222222222, 58.333333333333336, 72.22222222222221, 86.11111111111111, 100.0])

z = np.polyfit(w, v, 2)
print (z)
p=np.poly1d(z)
g = np.polyval(z,w)
print (g)
N=100
a=arange(N)
b=(w,v)
b=np.array(b)
c=(w,g)
c=np.array(c)
print(c)
d=-a+99
e=(a,d)
print (e)
p1=interpolate.PiecewisePolynomial(w,v[:,np.newaxis])
p2=interpolate.PiecewisePolynomial(w,d[:,np.newaxis])

def pdiff(x):
    return p1(x)-p2(x)

xs=np.r_[w,w]
xs.sort()
x_min=xs.min()
x_max=xs.max()
x_mid=xs[:-1]+np.diff(xs)/2
roots=set()
for val in x_mid:
    root,infodict,ier,mesg = optimize.fsolve(pdiff,val,full_output=True)
    # ier==1 indicates a root has been found
    if ier==1 and x_min<root<x_max:
        roots.add(root[0])
roots=list(roots)        
print(np.column_stack((roots,p1(roots),p2(roots))))

plt.plot(w,v, 'r', a, -a+99, 'b-')
plt.show()
q=input("what is the intersection value? ")
print (p(q))

有什么想法可以让它发挥作用吗?

谢谢

【问题讨论】:

    标签: python numpy matplotlib intersection curve-fitting


    【解决方案1】:

    我不认为我完全理解你在代码中试图做什么,但是你用英文描述的可以用

    from __future__ import division
    import numpy as np
    import matplotlib.pyplot as plt
    
    w = np.array([0.0, 11.11111111111111, 22.22222222222222, 33.333333333333336,
                  44.44444444444444, 55.55555555555556, 66.66666666666667,
                  77.77777777777777, 88.88888888888889, 100.0])
    v = np.array([0.0, 8.333333333333332, 16.666666666666664, 25.0,
                  36.11111111111111, 47.22222222222222, 58.333333333333336,
                  72.22222222222221, 86.11111111111111, 100.0])
    
    poly_coeff = np.polynomial.polynomial.polyfit(w, v, 2)
    poly = np.polynomial.polynomial.Polynomial(poly_coeff)
    roots = np.polynomial.polynomial.polyroots(poly_coeff - [99, -1, 0])
    
    x = np.linspace(np.min(roots) - 50, np.max(roots) + 50, num=1000)
    plt.plot(x, poly(x), 'r-')
    plt.plot(x, 99 - x, 'b-')
    for root in roots:
        plt.plot(root, 99 - root, 'ro')
    

    【讨论】:

    • 公平警告,np.polynomial.polynomial.polyfit 将系数 [A, B, C] 返回到 A + Bx + Cx^2 + ...,这与 np.polyfit(您最初使用的 @user2843767)返回的顺序相反:... + Ax^2 + Bx + C。不确定是谁做出的决定,只是不要在np.poly1d 或np.polyval 中使用第一个输出,除非你也使用np.polyfit
    • 确实是公平的警告。没有弃用警告,而且可能永远不会有,但是文档are clear 指出新代码的方法是 Polynomial 包,而不是旧的 poly1d。
    • 是的,幸运的是新(er)包也有更标准的排序。感谢您指出该链接,但我一定会只建议多项式包。
    • 很抱歉给您带来了困惑。也为混乱的代码感到抱歉,我尝试了十几件事,我认为其中有一堆挥之不去的行。我提交的代码实际上有效,我只需要不为两行使用相同的 x 值,将 xs=np.r_[w,w] 更改为 xs=np.r_[w,a] 即可修复它。但是这种方法很慢。我会尝试多项式包。谢谢。
    猜你喜欢
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    相关资源
    最近更新 更多