从polyfit找到曲线的交点

6

这似乎很简单,但我还是无法完全理解。我有一个由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))

有什么办法可以让这个工作起来吗?
谢谢。
1个回答

9

我不确定我完全理解你在代码中尝试做什么,但是根据你用英语描述的内容可以使用以下方式实现:

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')

enter image description here


1
提醒一下,np.polynomial.polynomial.polyfit 返回系数 [A, B, C],对应的多项式为 A + Bx + Cx^2 + ...,与 np.polyfit(你最初使用的函数,@user2843767)返回的顺序相反:... + Ax^2 + Bx + C。不确定是谁做出了这个决定,只是不要直接使用第一个输出并在 np.poly1dnp.polyval 中使用,除非你也使用了 np.polyfit - askewchan
确实需要提醒。虽然没有弃用警告,也可能永远不会有,但文档清楚表明,对于新代码而言,使用Polynomial包而不是旧的poly1d是正确的选择。请参考文档:http://docs.scipy.org/doc/numpy/reference/routines.polynomials.html - Jaime
是的,而且幸运的是新的包也有更标准的排序。谢谢你指出那个链接,我会确保只建议使用多项式包。 - askewchan
抱歉造成困惑。同时也很抱歉代码看起来有些混乱,我尝试了很多次,可能还有一些残留的代码行。实际上,我提交的代码是可以工作的,我只需要不要在两条线中使用相同的x值,将xs=np.r_[w,w]更改为xs=np.r_[w,a]即可解决问题。但这种方法速度较慢。我将尝试使用Polynomial包。谢谢。 - user2843767

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接