在Python中寻找Legendre多项式的根

Question

在Python中寻找Legendre多项式的根

4

我正在编写一个通过Legendre-Gauss积分求解的程序。对于n阶积分法算法，需要在某个时候找到n阶Legendre多项式Pn（x）的根，并将它们指定到Absc数组中。Pn是一个n阶多项式，在[-1,1]区间上具有n个独立实数根。我想能够计算出这些根，而不仅仅是从某个库中导入它们。我已经可以创建一个给出多项式系数的数组，我称之为PCoeff。为了找到这些根，我尝试过

 Absc = numpy.roots(PCoeff)

这段代码可以在n小于等于40的情况下正常工作，但是如果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阶的积分以达到高精度，但这个根查找会导致失败。你有什么建议吗？

- johndmalcolm

2

你是否了解numpy.polynomial.legendre模块？它包含了legroots方法，用于查找根。你可以在这里进行检查：http://docs.scipy.org/doc/numpy/reference/routines.polynomials.legendre.html - Teudimundo

1

能够自己计算是一件值得骄傲的事情，但这个模块对我所需的功能非常好用。比我之前做的要好得多。非常感谢。 - johndmalcolm

2个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Francesco Gramano · Answer 1

您可以在SymPy文档中找到一个简单的示例实现，以解决您的问题：

>>> for n in range(5):
...     nprint(polyroots(taylor(lambda x: legendre(n, x), 0, n)[::-1]))
...
[]
[0.0]
[-0.57735, 0.57735]
[-0.774597, 0.0, 0.774597]
[-0.861136, -0.339981, 0.339981, 0.861136]

正如之前的答案所建议的那样，该示例使用了多项式的泰勒展开。

- Nicolas Barbey · Answer 2

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

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

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