Solving for x values of polynomial with known y in scipy / numpy

Question

I am trying to solve for the x values with a known y. I was able to get the polynomial to fit my data, and now I want to know the x value that a chosen y would land on the curve.

import numpy as np

x = [50, 25, 12.5, 6.25, 0.625, 0.0625, 0.01]
y = [0.00, 0.50, 0.68, 0.77, 0.79, 0.90, 1.00]

poly_coeffs = np.polyfit(x, y, 3)

f = np.poly1d(poly_coeffs)

I want to do 0.5 = f and solve for the x values.

I can solve this in WolframAlpha by typing:

0.5 = -9.1e-6*x^3 + 5.9e-4*x^2 - 2.5e-2*x + 9.05e-1

The real x value is ~26

Solution 1

In [1]: from numpy.polynomial import Polynomial as P

In [2]: x = [50, 25, 12.5, 6.25, 0.625, 0.0625, 0.01]

In [3]: y = [0.00, 0.50, 0.68, 0.77, 0.79, 0.90, 1.00]

In [4]: p = P.fit(x, y, 3)

In [5]: (p - .5).roots()
Out[5]: 
array([ 19.99806935-37.92449551j,  19.99806935+37.92449551j,
        25.36882693 +0.j        ])

Looks like the root you want is 25.36882693.

Solution 2

You can solve the equation f(x) - y = 0 using np.roots. Consider the function:

def solve_for_y(poly_coeffs, y):
    pc = poly_coeffs.copy()
    pc[-1] -= y
    return np.roots(pc)

Then you can use it to solve your polynomial for any y you want:

>>> print solve_for_y(poly_coeffs, 0.5)
[ 19.99806935+37.92449551j  19.99806935-37.92449551j  25.36882693 +0.j        ]
>>> print solve_for_y(poly_coeffs, 1.)
[ 40.85615395+50.1936152j  40.85615395-50.1936152j -16.34734226 +0.j       ]

Comments

[Amit Kumar Gupta]

Your method appears to modify poly_coeffs, so the second call to the method will actually give the wrong answer, i.e. it's solving for y = 1 for x in a different polynomial than the original one.

[Andrey Sobolev]

@Amit: Yes, my bad. I thought that function parameters are passed by value, and it seems to be true only for immutable types. Since np.array is mutable, it is passed by reference. I'll edit the post, thank you.