在python中求Legendre多项式的根

我正在编写一个程序，通过勒让德-高斯求积求解积分。 n 阶求积算法需要在某一时刻找到 n 阶勒让德多项式 Pn(x) 的根，并将它们分配给数组 Absc（表示“横坐标”）。 Pn 是一个 n 阶多项式，在区间 [-1,1] 上有 n 个独立的实根。我希望能够计算根，而不是仅仅从某个库导入它们。我能够创建一个给出多项式系数的数组，我将其称为 PCoeff。为了寻找根源我曾尝试过

 Absc = numpy.roots(PCoeff)

这在 n = 40 左右有效，但超过这个值就开始失败，在不应该的情况下给出复杂的根。我还尝试使用定义多项式

P = numpy.poly1d(PCoeff)
Absc = P.r

但这会带来同样的问题，大概是因为它使用相同的 numpy 求根算法。

另一种似乎很有前途的方法是使用 scipy.optimize.fsolve(Pn, x0)，其中 x0 是我对根的猜测的 n 元素数组。这样做的问题是，根据我的 x0 选择，此方法可能会多次给出一个特定的根来代替其他根。我尝试将 x0 填充为 [-1,1] 上的等距点

x0 = numpy.zeros(n)
step = 2./(n-1)
for i in xrange(n):
  x0[i] = -1. + i*step

但是一旦我达到 n = 5，fsolve 就会重复给出一些根并忽略其他根。我还尝试使用 numpy.roots 的结果作为 x0。然而，在 np.roots 给出复数值的问题区域中，这些会导致 fsolve 出现错误

TypeError: array cannot be safely cast to required type

我在网上看到有一个 scipy.optimize.roots() 例程可以工作，但它不在我的计算机上的 scipy 库中。更新很麻烦，因为我没有权限在这台计算机上下载东西。

我希望能够以 64 阶运行求积以获得高精度，但这种求根会导致失败。有任何想法吗？

由于 np.roots 依赖于文档中所述的“查找伴随矩阵的特征值”，因此您可能会遇到错误传播问题，导致根上的虚部不为零。也许您可以使用 np.real 函数丢弃虚部。

您可以尝试使用根的泰勒近似来尝试不同的方法来计算根：

https://math.stackexchange.com/questions/12160/roots-of-legendre-polynomial https://math.stackexchange.com/questions/12160/roots-of-legendre-polynomial