2013-10-06 15 views
5

Das scheint einfach, aber ich kann es nicht ganz herausfinden. Ich habe eine Kurve berechnet aus x, y Daten. Dann habe ich eine Linie. Ich möchte die x, y-Werte für den Schnittpunkt der beiden finden.Den Schnittpunkt einer Kurve von Polyfit finden

Hier ist, was ich bisher habe. Es ist sehr verwirrend und gibt nicht das richtige Ergebnis. Ich kann den Graphen betrachten und den x-Wert der Schnittmenge finden und den korrekten y-Wert berechnen. Ich möchte diesen menschlichen Schritt entfernen.

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

Irgendwelche Ideen, um dies zum Funktionieren zu bringen?

Dank

Antwort

5

Ich glaube nicht vollständig verstehe ich, was Sie in Ihrem Code zu tun versuchen, aber was beschrieben Sie in englischer Sprache können Sie mit

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

Eine Warnung erfolgen , 'np.polynomial.polynomial.polyfit' gibt die Koeffizienten' [A, B, C] 'an' A + Bx + Cx^2 + ... 'zurück, was die umgekehrte Reihenfolge ist von dem was' np.polyfit' (was Sie hatten ursprünglich verwendet, @ user2843767) gibt zurück: '... + Ax^2 + Bx + C'. Nicht sicher, wer diese Entscheidung getroffen hat, nimm einfach nicht die erste Ausgabe und verwende sie in 'np.poly1d' oder np.polyval, es sei denn, du verwendest auch' np.polyfit'. – askewchan

+0

In der Tat eine faire Warnung. Es gibt keine Verwarnungswarnung, und es gibt vielleicht nie, aber die Dokumente [sind klar] (http://docs.scipy.org/doc/numpy/reference/routines.polynomials.html), die den Weg für neuen Code gehen ist das Polynom-Paket, nicht das ältere poly1d. – Jaime

+0

Ja, und zum Glück hat das neue (er) Paket auch die Standard-Bestellung. Danke, dass Sie auf diesen Link hingewiesen haben, ich werde jedoch nur das Polynom-Paket empfehlen. – askewchan