2013-10-06 14 views
5

Wydaje się to proste, ale nie mogę tego dokładnie rozgryźć. Mam krzywą obliczoną na podstawie danych x, y. Potem mam linię. Chcę znaleźć wartości x, y dla miejsca, w którym te dwie przecinają się.Znaleźć przecięcie krzywej z polyfit

Oto, co mam do tej pory. Jest to bardzo mylące i nie daje prawidłowego wyniku. Mogę spojrzeć na wykres i znaleźć wartość przecięcia x oraz obliczyć prawidłową wartość y. Chciałbym usunąć ten ludzki krok.

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

Jakieś pomysły, aby to zadziałało?

Dzięki

Odpowiedz

5

Nie sądzę w pełni zrozumieć, co chce zrobić w kodzie, ale to, co opisane w języku angielskim można zrobić z

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

Fair ostrzegawczego , 'np.polynomial.polynomial.polyfit' zwraca współczynniki' [A, B, C] 'na' A + Bx + Cx^2 + ... ', które jest odwrotną kolejnością od tego, co' np.polyfit' (co pierwotnie użyłeś, @ user2843767) zwraca: '... + Ax^2 + Bx + C'. Nie jestem pewien, kto podjął taką decyzję, po prostu nie przyjmuj pierwszego wyniku i używaj go w 'np.poly1d' lub np.polyval, chyba że użyjesz' np.polyfit'. – askewchan

+0

Rzeczywiście, uczciwe ostrzeżenie. Nie ma ostrzeżenia o wycofaniu i może nigdy nie być, ale dokumenty [są jasne] (http://docs.scipy.org/doc/numpy/reference/routines.polynomials.html), że sposób, aby przejść do nowego kodu jest pakietem wielomianowym, a nie starszym poly1d. – Jaime

+0

Tak, i na szczęście nowy (er) pakiet ma również bardziej standardowe zamówienie. Dzięki za wskazanie tego linku, na pewno doradzę tylko pakiet wielomianowy. – askewchan