Нахождение пересечения кривой из полифита

Это кажется простым, но я не могу понять это. У меня есть кривая, рассчитанная по данным 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
person user2843767    schedule 06.10.2013    source источник

Ответы (1)


arrow_upward
9
arrow_downward

Я не думаю, что полностью понимаю, что вы пытаетесь сделать в своем коде, но то, что вы описали на английском языке, можно сделать с помощью

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

введите здесь описание изображения

person Jaime    schedule 07.10.2013
comment
Справедливое предупреждение, np.polynomial.polynomial.polyfit возвращает коэффициенты от [A, B, C] до A + Bx + Cx^2 + ..., что является порядком, противоположным тому, что возвращает np.polyfit (то, что вы изначально использовали, @user2843767): ... + Ax^2 + Bx + C. Не уверен, кто принял это решение, просто не берите первый вывод и не используйте его в np.poly1d или np.polyval, если только вы не используете np.polyfit. - person askewchan; 07.10.2013
comment
Справедливое предупреждение. Предупреждения об устаревании нет и, возможно, никогда не будет, но документы ясно, что для нового кода лучше использовать пакет Polynomial, а не старый poly1d. - person Jaime; 07.10.2013
comment
Да, и, к счастью, новый пакет также имеет более стандартный порядок. Спасибо за указание на эту ссылку, я обязательно посоветую только полиномиальный пакет. - person askewchan; 07.10.2013
comment
Извините за путаницу. Также извините за беспорядочный код, я перепробовал с десяток вещей, и я думаю, что из этого была куча затяжных строк. Код, который я представил, действительно работает, мне просто нужно не использовать одни и те же значения x для обеих строк, изменение xs=np.r_[w,w] на xs=np.r_[w,a] исправляет это. Но этот метод медленный. Я попробую полиномиальный пакет. Спасибо. - person user2843767; 08.10.2013