在Y点求极大多项式的X?(Python 3.10,NumPy)

kyvafyod  于 2023-08-05  发布在  Python
关注(0)|答案(1)|浏览(123)

我试图计算所有的真实的X值,在一定的Y值,一个20次多项式提供了一个系数列表的升序权力。我使用Python 3.10以及NumPy库的root函数。
(我之前发布了一个simillar question here on StackOverflow,试图通过用一个小得多的多项式生成一个更简单的例子来解决这个问题。然而,虽然该帖子的解决方案确实解决了小多项式的问题,但它并没有解决我使用大多项式时的问题。
下面是我的示例代码:

import numpy as np
import matplotlib.pyplot as plt

def main():
    # Declare 20-degree polynomial
    coeffs = np.array([0.028582380269090664, -11.07209382330515, 81.98886137780704, -231.55726577581814, 344.8963540923082, -312.0774305574022, 185.36300988588312, -75.58767327080196, 21.660357329823697, -4.380346453954649, 0.6125779086029393, -0.05524029639869697, 0.002489646479713822, 4.159253423505381e-05, -9.927184104751197e-06, 8.634292922740336e-08, 3.278890552085848e-08, -3.350925351280405e-10, -1.1510786631391544e-10, -2.3053456641329463e-13, 3.9927010006271344e-13, 8.499990550259627e-15, -1.2514130885401332e-15, -5.653183054629818e-17, 3.5580325101956145e-18, 2.634844330077943e-19, -1.2086461102288157e-20, -9.772155900613053e-22, 8.336597931160041e-23, -2.286862805703104e-24, 2.2638043801338818e-26], dtype = float)
    y = -5

    # Create polynomial data on interval (0, 17)
    polyDataX = np.linspace(0, 17)
    polyDataY = np.empty(shape = len(polyDataX), dtype = float)

    for i in range(len(polyDataX)):
        polyDataY[i] = 0

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

    # Solve for all real solutions
    coeffs[-1] -= y
    sols = np.roots(coeffs).tolist()
    i = 0

    while (i < len(sols)):
        if (sols[i].imag != 0) or (sols[i].real < 0) or (sols[i].real > 17):
            del sols[i]
            i -= 1
        else:
            sols[i] = round(sols[i].real, 2)

        i += 1

    # Plot polynomial & solutions
    plt.xlim([0, 17])
    plt.ylim([-30, 30])
    plt.plot(polyDataX, polyDataY, color = "gray")
    plt.axhline(y, color = "blue")
    
    for sol in sols:
        plt.axvline(sol, color = "red")

    plt.title("Sols: " + str(sols))
    plt.show()
    plt.close()
    plt.clf()

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

字符串
首先,我生成一个20次多项式(存储在coeffs中),并使用roots函数以获取某个Y值(存储在y中)处的所有X值(存储在sols中)。我随后过滤掉所有的虚数解,只留下在所需间隔[0, 17]处的真实的解,并绘制多项式,以及所需的Y值和找到的实数X值。我在gray中对多项式着色,在blue中对所需的Y值着色,在red中对所有过滤后的X值着色。
在上面的示例中,我尝试在-5的Y值处获取所有实X值。X值被列为[3.92, 1.34, 1.13],根据生成的图,这显然是不正确的:x1c 0d1x的数据
根据该图,对于-5的Y值,实际的X值应该是~ 11.2和~ 16
我做错什么了?为什么这里列出的roots函数的解决方案不正确?
感谢您阅读我的文章,任何指导都很感激。

yjghlzjz

yjghlzjz1#

删除大部分代码并使用内置的Polynomial支持:

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

def get_roots(poly: Polynomial, y: float = 0) -> np.ndarray:
    roots = (poly - y).roots()
    real_roots = roots.real[np.isreal(roots)]
    return real_roots[
        (poly.domain[0] <= real_roots) &
        (poly.domain[1] >= real_roots)
    ]

def plot(poly: Polynomial, real_roots: np.ndarray, y: float) -> plt.Figure:
    x = np.linspace(6, poly.domain[1], num=1_000)

    fig, ax = plt.subplots()
    ax.set_title(f"Roots of degree-{poly.degree()} polynomial")
    ax.set_ylim(-1, 1)
    ax.axhline(y, c='black', linewidth=0.5)
    ax.plot(x, poly(x))
    for root in real_roots:
        ax.annotate(
            text=f'{root:.3f}',
            xy=(root, y),
            xytext=(root - 0.5, -0.6),
            arrowprops={'arrowstyle': '-'},
        )

    return fig

def main() -> None:
    coeffs = np.array([
        2.85823803e-02, -1.10720938e+01,  8.19888614e+01, -2.31557266e+02,
        3.44896354e+02, -3.12077431e+02,  1.85363010e+02, -7.55876733e+01,
        2.16603573e+01, -4.38034645e+00,  6.12577909e-01, -5.52402964e-02,
        2.48964648e-03,  4.15925342e-05, -9.92718410e-06,  8.63429292e-08,
        3.27889055e-08, -3.35092535e-10, -1.15107866e-10, -2.30534566e-13,
        3.99270100e-13,  8.49999055e-15, -1.25141309e-15, -5.65318305e-17,
        3.55803251e-18,  2.63484433e-19, -1.20864611e-20, -9.77215590e-22,
        8.33659793e-23, -2.28686281e-24,  2.26380438e-26])
    y = -0.1
    poly = Polynomial(coeffs, domain=(0, 17))
    real_roots = get_roots(poly, y)
    plot(poly, real_roots, y)
    plt.show()

if __name__ == "__main__":
    main()

字符串


的数据

相关问题