Get Polynomial X at Y? (Python 3.10, NumPy)

Question

I'm attempting to calculate all possible real X-values at a certain Y-value from a polynomial given in descending coefficent order, in Python 3.10. I want the resulting X-values to be provided to me in a list.

I've tried using the roots() function of the numpy library, as shown in one of the answers to this post, however it does not appear to work:

import numpy as np
import matplotlib.pyplot as plt

def main():
    coeffs = np.array([1, 2, 2])
    y = 1.5

    polyDataX = np.linspace(-2, 0)
    polyDataY = np.empty(shape = len(polyDataX), dtype = float)

    for i in range(len(polyDataX)):
        polyDataY[i] = coeffs[0] * pow(polyDataX[i], 2) + coeffs[1] * polyDataX[i] + coeffs[2]

    coeffs[-1] -= y
    x = np.roots(coeffs).tolist()

    plt.axhline(y, color = "orange")
    plt.plot(polyDataX, polyDataY, color = "blue")
    plt.title("X = " + str(x))
    plt.show()
    plt.close()
    plt.clf()

if (__name__ == "__main__"):
    main()

In my example above, I have the coefficents of my polynomial stored in the local variable coeffs, in descending order. I then attempt to gather all the X-values at the Y-value of 0.5, stored within the x and y local variables respectivelly. I then display the gathered X-values as the title of the shown plot.

The script above results in the following plot:

With the X-values being shown as [-2.0, 0.0], instead of the correct:

What is the proper way to get all real X-values of a polynomial at a certain Y-value in Python?

Thanks for reading my post, any guidance is appreciated.

score 3 · Answer 1 · edited Jul 15 '23 at 10:04

Your code is correct and you are using the root function correctly. There is only a small issue with datatypes. You can see the error with a debugger or by printing coeff after coeffs[-1] -= y. At your line coeffs = np.array([1, 2, 2]), numpy will create the array with an integer type and therefore your subtraction, which should result in 0.5 actually results in 0, so your coefficients become [1, 2, 0]. For those, the correct result is calculated. So we need to simply set the dtype to float:

import numpy as np
import matplotlib.pyplot as plt

def main():
    coeffs = np.array([1, 2, 2], dtype=float)   #<-----------
    y = 1.5

    polyDataX = np.linspace(-2, 0)
    polyDataY = np.empty(shape = len(polyDataX), dtype = float)

    for i in range(len(polyDataX)):
        polyDataY[i] = coeffs[0] * pow(polyDataX[i], 2) + coeffs[1] * polyDataX[i] + coeffs[2]

    coeffs[-1] -= y
    x = np.roots(coeffs).tolist()

    plt.axhline(y, color = "orange")
    plt.plot(polyDataX, polyDataY, color = "blue")
    plt.title("X = " + str(x))
    plt.show()
    plt.close()
    plt.clf()

if __name__ == "__main__":
    main()

Thanks for sharing this solution, it worked. Just out of curiousity, what does the `coeffs[-1] -= y` line actually do? Does that shift the entire polynomial to make it so the `0` alligns with the value `y`, so we can calculate the roots? If so, why are we just changing the coefficents? — Runsva, Jul 03 '23 at 20:26
Solving the equation f(x) = y is equivalent to solving f(x) - y = 0, i.e., finding the roots. That's what's happening here when you subtract y. — meferne, Jul 03 '23 at 22:12

jared · Accepted Answer · 2023-07-03T22:53:12.427

You should be making use of the numpy.polynomial.Polynomial class that was added in numpy v1.4 (more information here). With that class, you can create a polynomial object. To find your solution, you can subtract y from the Polynomial object and then call the roots method. Another nice feature is that you can directly call the object, which you can use to compute polyDataY.

Just note that the Polynomial class expects the coefficients to be given backward from np.roots, i.e. a quadratic of x^2 + 2x + 2 should have the coefficients (2, 2, 1). To keep things consistent with what you gave, I just pass the reversed coeffs.

import numpy as np
from numpy.polynomial import Polynomial
import matplotlib.pyplot as plt

plt.close("all")

coeffs = np.array([1, 2, 2])
poly = Polynomial(coeffs[::-1])

polyDataX = np.linspace(-2, 0)
polyDataY = poly(polyDataX)

y = 1.5
x = (poly - y).roots()
plt.axhline(y, color = "orange")
plt.plot(polyDataX, polyDataY, color = "blue")
plt.title("X = " + str(x))
plt.show()

score 1 · Answer 3 · answered Jul 03 '23 at 10:50

1

The values in coeffs are int32 and the result after coeffs[-1] -= y is 0 not 0.5.

Just change coeffs = np.array([1, 2, 2]) for coeffs = np.array([1., 2., 2.])

answered Jul 03 '23 at 10:50

Andy R

144
1
6

Get Polynomial X at Y? (Python 3.10, NumPy)

3 Answers3

Linked