Python 多项式根不准确
Python polynomial roots are inaccurate
我尝试实现的算法需要找到 10 次多项式的根,这是我用 sympy 创建的,它看起来像这样:
import sympy
import numpy as np
det = sympy.Poly(1.3339507303385e-16*z**10 + 6.75390067076469e-14*z**9 + 7.18791227134351e-12*z**8 + 2.27504286959352e-10*z**7 + 2.37058998324426e-8*z**6 + 1.63916629439745e-6*z**5 + 3.0608041245671e-5*z**4 + 4.83564348562906e-8*z**3 + 2.0248073519853e-5*z**2 - 4.73126177166988e-7*z + 1.1495883206077e-6)
为了求多项式的根,我使用以下代码:
coefflist = det.coeffs()
solutions = np.roots(coefflist)
print(coefflist)
[1.33395073033850e-16, 6.75390067076469e-14, 7.18791227134351e-12, 2.27504286959352e-10, 2.37058998324426e-8, 1.63916629439745e-6, 3.06080412456710e-5, 4.83564348562906e-8, 2.02480735198530e-5, -4.73126177166988e-7, 1.14958832060770e-6]
print(solutions)
[-3.70378229e+02+0.00000000e+00j -1.18366138e+02+0.00000000e+00j
2.71097137e+01+5.77011644e+01j 2.71097137e+01-5.77011644e+01j
-3.59084863e+01+1.44819591e-02j -3.59084863e+01-1.44819591e-02j
2.60969082e-03+7.73805425e-01j 2.60969082e-03-7.73805425e-01j
1.42936329e-02+2.49877948e-01j 1.42936329e-02-2.49877948e-01j]
但是,当我用根替换 z
时,假设第一个,结果不是零,而是一些数字:
print(det.subs(z,solutions[0]))
-1.80384169514123e-6
我本以为结果可能不是整数 0
,但 1e-6
非常糟糕(应该为零,对吧?)。我的代码有错误吗?这种不准确是正常的吗?任何 thoughts/suggestions 都会有所帮助。是否有更准确的方法来计算 10 次多项式的根?
不需要sympy,numpy中的方法完全够用。通过系数列表定义多项式并计算根
p=[1.33395073033850e-16, 6.75390067076469e-14, 7.18791227134351e-12, 2.27504286959352e-10, 2.37058998324426e-8, 1.63916629439745e-6, 3.06080412456710e-5, 4.83564348562906e-8, 2.02480735198530e-5, -4.73126177166988e-7, 1.14958832060770e-6]
sol= np.roots(p); sol
给出结果
array([ -3.70378229e+02 +0.00000000e+00j, -1.18366138e+02 +0.00000000e+00j,
2.71097137e+01 +5.77011644e+01j, 2.71097137e+01 -5.77011644e+01j,
-3.59084863e+01 +1.44819592e-02j, -3.59084863e+01 -1.44819592e-02j,
2.60969082e-03 +7.73805425e-01j, 2.60969082e-03 -7.73805425e-01j,
1.42936329e-02 +2.49877948e-01j, 1.42936329e-02 -2.49877948e-01j])
并计算这些近似根处的多项式
np.polyval(p,sol)
给出数组
array([ 2.28604877e-06 +0.00000000e+00j, 1.30435230e-10 +0.00000000e+00j,
1.05461854e-11 -7.56043461e-12j, 1.05461854e-11 +7.56043461e-12j,
-3.98439686e-14 +6.84489332e-17j, -3.98439686e-14 -6.84489332e-17j,
1.18584613e-20 +1.59976730e-21j, 1.18584613e-20 -1.59976730e-21j,
6.35274710e-22 +1.74700545e-21j, 6.35274710e-22 -1.74700545e-21j])
显然,求一个接近根的多项式会涉及很多灾难性的抵消,也就是说,中间项的符号很大并且抵消了,但它们的误差与它们的原始大小成正比。要估计组合误差大小,请将多项式系数替换为它们的绝对值以及评估点。
np.polyval(np.abs(p),np.abs(sol))
导致
array([ 1.81750423e+10, 8.40363409e+05,
8.08166359e+03, 8.08166359e+03,
2.44160616e+02, 2.44160616e+02,
2.50963696e-05, 2.50963696e-05,
2.65889696e-06, 2.65889696e-06])
在第一个根的情况下,尺度乘以机器常数给出了1e+10*1e-16=1e-6
的误差尺度,这意味着根处的值在double-的框架内与零一样好精度浮点。
我尝试实现的算法需要找到 10 次多项式的根,这是我用 sympy 创建的,它看起来像这样:
import sympy
import numpy as np
det = sympy.Poly(1.3339507303385e-16*z**10 + 6.75390067076469e-14*z**9 + 7.18791227134351e-12*z**8 + 2.27504286959352e-10*z**7 + 2.37058998324426e-8*z**6 + 1.63916629439745e-6*z**5 + 3.0608041245671e-5*z**4 + 4.83564348562906e-8*z**3 + 2.0248073519853e-5*z**2 - 4.73126177166988e-7*z + 1.1495883206077e-6)
为了求多项式的根,我使用以下代码:
coefflist = det.coeffs()
solutions = np.roots(coefflist)
print(coefflist)
[1.33395073033850e-16, 6.75390067076469e-14, 7.18791227134351e-12, 2.27504286959352e-10, 2.37058998324426e-8, 1.63916629439745e-6, 3.06080412456710e-5, 4.83564348562906e-8, 2.02480735198530e-5, -4.73126177166988e-7, 1.14958832060770e-6]
print(solutions)
[-3.70378229e+02+0.00000000e+00j -1.18366138e+02+0.00000000e+00j
2.71097137e+01+5.77011644e+01j 2.71097137e+01-5.77011644e+01j
-3.59084863e+01+1.44819591e-02j -3.59084863e+01-1.44819591e-02j
2.60969082e-03+7.73805425e-01j 2.60969082e-03-7.73805425e-01j
1.42936329e-02+2.49877948e-01j 1.42936329e-02-2.49877948e-01j]
但是,当我用根替换 z
时,假设第一个,结果不是零,而是一些数字:
print(det.subs(z,solutions[0]))
-1.80384169514123e-6
我本以为结果可能不是整数 0
,但 1e-6
非常糟糕(应该为零,对吧?)。我的代码有错误吗?这种不准确是正常的吗?任何 thoughts/suggestions 都会有所帮助。是否有更准确的方法来计算 10 次多项式的根?
不需要sympy,numpy中的方法完全够用。通过系数列表定义多项式并计算根
p=[1.33395073033850e-16, 6.75390067076469e-14, 7.18791227134351e-12, 2.27504286959352e-10, 2.37058998324426e-8, 1.63916629439745e-6, 3.06080412456710e-5, 4.83564348562906e-8, 2.02480735198530e-5, -4.73126177166988e-7, 1.14958832060770e-6]
sol= np.roots(p); sol
给出结果
array([ -3.70378229e+02 +0.00000000e+00j, -1.18366138e+02 +0.00000000e+00j,
2.71097137e+01 +5.77011644e+01j, 2.71097137e+01 -5.77011644e+01j,
-3.59084863e+01 +1.44819592e-02j, -3.59084863e+01 -1.44819592e-02j,
2.60969082e-03 +7.73805425e-01j, 2.60969082e-03 -7.73805425e-01j,
1.42936329e-02 +2.49877948e-01j, 1.42936329e-02 -2.49877948e-01j])
并计算这些近似根处的多项式
np.polyval(p,sol)
给出数组
array([ 2.28604877e-06 +0.00000000e+00j, 1.30435230e-10 +0.00000000e+00j,
1.05461854e-11 -7.56043461e-12j, 1.05461854e-11 +7.56043461e-12j,
-3.98439686e-14 +6.84489332e-17j, -3.98439686e-14 -6.84489332e-17j,
1.18584613e-20 +1.59976730e-21j, 1.18584613e-20 -1.59976730e-21j,
6.35274710e-22 +1.74700545e-21j, 6.35274710e-22 -1.74700545e-21j])
显然,求一个接近根的多项式会涉及很多灾难性的抵消,也就是说,中间项的符号很大并且抵消了,但它们的误差与它们的原始大小成正比。要估计组合误差大小,请将多项式系数替换为它们的绝对值以及评估点。
np.polyval(np.abs(p),np.abs(sol))
导致
array([ 1.81750423e+10, 8.40363409e+05,
8.08166359e+03, 8.08166359e+03,
2.44160616e+02, 2.44160616e+02,
2.50963696e-05, 2.50963696e-05,
2.65889696e-06, 2.65889696e-06])
在第一个根的情况下,尺度乘以机器常数给出了1e+10*1e-16=1e-6
的误差尺度,这意味着根处的值在double-的框架内与零一样好精度浮点。