python - 在python中寻找勒让德多项式的根

Question

我正在编写一个程序，通过勒让德-高斯积分求解积分。 n 阶正交算法需要在某一时刻找到 n 阶勒让德多项式 Pn(x) 的根，并将它们分配给数组 Absc(表示“横坐标”)。 Pn 是区间 [-1,1] 上有 n 个独立实根的 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 阶运行正交以获得高精度，但这个根发现导致失败。有什么想法吗？

Answer 1

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

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

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