求解所有六个根的位置 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)
我用的是牛顿法,所以想求出六次多项式的所有六个根的位置,基本上就是函数为零的点。
我用下面的这段代码在我的图表上找到了粗略的值,但我想输出所有六个根的那些位置。我正在考虑使用 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)