在python中找到勒让德多项式的根 [英] Finding roots of Legendre polynomial in python

问题描述

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

I'm writing a program that solves an integral by Legendre-Gauss quadrature. The algorithm for nth-order quadrature requires, at one point, finding the roots of the nth-order Legendre polynomial, Pn(x), assigning them to the array Absc (for 'abscissa'). Pn is an nth order polynomial with n independent real roots on the interval [-1,1]. I'd like to be able to compute the roots, instead of just importing them from some library. I am able to create an array that gives the polynomial coefficients, which I call PCoeff. To find the roots I have tried

 Absc = numpy.roots(PCoeff)

这最多可以工作到大约 n = 40，但除此之外它开始失败，当它真的不应该时给出复杂的根.我也试过使用

This works up to about n = 40, but beyond that it starts to fail, giving complex roots when it really shouldn't be. I've also tried defining the polynomial using

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

但这会产生相同的问题，大概是因为它使用相同的 numpy 求根算法.

but this gives the same problems, presumably because it uses the same numpy root-finding algorithm.

另一种看起来很有前景的方法是使用 scipy.optimize.fsolve(Pn, x0)，其中 x0 是我猜测的根处的 n 元素数组.问题在于，根据我的 x0 选择，此方法可能会多次给出一个特定的根来代替其他根.我试过将 x0 填充为 [-1,1]

Another method, which seemed promising was using scipy.optimize.fsolve(Pn, x0), where x0 is an n-element array of my guess at the roots. The problem with this is that, depending on my x0 choices, this method may give one particular root more than once in the place of other roots. I've tried populating x0 as equidistant points on [-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 出现错误

but once I get to n = 5, fsolve gives some roots repeated and neglects others. I have also tried using the results of numpy.roots as x0. However, in the problem area where np.roots gives complex values, these cause an error in fsolve

TypeError: array cannot be safely cast to required type

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

I saw online that there is a scipy.optimize.roots() routine that could work, but it is not in the scipy library on my computer. Updating is a hassle since I don't have permission to download things on this computer.

为了高精度，我希望能够以 64 阶运行正交，但此根本发现会导致失败.有什么想法吗?

I'd like to be able to run the quadrature at 64th order for high accuracy, but this root finding causes failure. Any ideas?

推荐答案

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

As np.roots rely on "finding eigenvalues of the Companion matrix" as stated by the documentation you may be hit by an error propagation issue resulting in non-zero imaginary part on the roots. Maybe you could just discard the imaginary part using the np.real function.

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

You could try a different way to compute the roots using taylor approximation of the roots :

https://math.stackexchange.com/questions/12160/roots-of-legendre-多项式

这篇关于在python中找到勒让德多项式的根的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！