求解所有六个根的位置 PYTHON

Solve for the positions of all six roots PYTHON

我用的是牛顿法,所以想求出六次多项式的所有六个根的位置,基本上就是函数为零的点。

我用下面的这段代码在我的图表上找到了粗略的值,但我想输出所有六个根的那些位置。我正在考虑使用 x 作为数组来输入值以找到那些位置但不确定。我现在使用 1.0 来定位粗略值。这里有什么建议吗??

def P(x):
        return 924*x**6 - 2772*x**5 + 3150*x**4 - 1680*x**3 + 420*x**2 - 42*x + 1
def dPdx(x):
        return 5544*x**5 - 13860*x**4 + 12600*x**3 - 5040*x**2 + 840*x - 42

accuracy = 1**-10
x = 1.0                                       
xlast = float("inf")

while np.abs(x - xlast) > accuracy:
        xlast = x
        x = xlast - P(xlast)/dPdx(xlast)
print(x)

p_points = []
x_points = np.linspace(0, 1, 100)
y_points = np.zeros(len(x_points))
for i in range(len(x_points)):
        y_points[i] = P(x_points[i])
p_points.append(P(x_points))


plt.plot(x_points,y_points)
plt.savefig("roots.png")
plt.show()

传统的方法是使用紧缩来分解出已经找到的根。如果您想避免对系数数组进行操作,则必须将根除掉。

找到 z[1],...,z[k] 作为根近似值,形成

g(x)=(x-z[1])*(x-z[2])*...*(x-z[k])

并将牛顿法应用于 h(x)=f(x)/g(x)h'(x)=f'/g-fg'/g^2。在牛顿迭代中,这给出了

xnext = x - f(x)/( f'(x) - f(x)*g'(x)/g(x) )

幸好商g'/g有一个简单的形式

g'(x)/g(x) = 1/(x-z[1])+1/(x-z[2])+...+1/(x-z[k])

因此,稍微修改一下牛顿步,您就可以避免再次找到相同的根。


这一切仍然保持迭代真实。要获得复数根,请使用复数开始迭代。


概念验证,将 eps=1e-8j 添加到 g'(x)/g(x) 允许迭代变得复杂而不妨碍实际值。解决等效问题 0=exp(-eps*x)*f(x)/g(x)

import numpy as np
import matplotlib.pyplot as plt

def P(x):
        return 924*x**6 - 2772*x**5 + 3150*x**4 - 1680*x**3 + 420*x**2 - 42*x + 1
def dPdx(x):
        return 5544*x**5 - 13860*x**4 + 12600*x**3 - 5040*x**2 + 840*x - 42



accuracy = 1e-10

roots = []
for k in range(6):
    x = 1.0                                       
    xlast = float("inf")

    x_points = np.linspace(0.0, 1.0, 200)
    y_points = P(x_points)
    for rt in roots:
        y_points /= (x_points - rt)
    y_points = np.array([ max(-1.0,min(1.0,np.real(y))) for y in y_points ])

    plt.plot(x_points,y_points,x_points,0*y_points)
    plt.show()
    while np.abs(x - xlast) > accuracy:
        xlast = x
        corr = 1e-8j
        for rt in roots:
            corr += 1/(xlast-rt)
        Px = P(xlast)
        dPx = dPdx(xlast)
        x = xlast - Px/(dPx - Px*corr)
    print(x)
    roots.append(x)