2013-10-06 25 views
5

これは単純なようですが、それはわかりません。私は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)) 

これを行うには、どのような考えですか?

おかげ

答えて

5

私は、私は完全にあなたのコードでやろうとしているのか理解とは思わないが、あなたは英語で説明し、何が

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 + ... 'に返します。もともと使用していた@ user2843767)は、 '... + Ax^2 + Bx + C'を返します。誰がその決定を下したのかはっきりしないが、 'np.polyfit'も使わない限り、最初の出力を' np.poly1d'やnp.polyvalで使用しないでください。 – askewchan

+0

確かに公正な警告。廃止予定の警告はなく、決して存在しないかもしれませんが、新しいコードを作成する方法は明確です(http://docs.scipy.org/doc/numpy/reference/routines.polynomials.html)。古いpoly1dではなく、多項式パッケージです。 – Jaime

+0

はい、そして幸運なことに、新しいパッケージにはより標準的な順序付けがあります。そのリンクを指摘してくれてありがとう、私は多項式パッケージだけを助言することを確認します。 – askewchan

関連する問題