我正在编写一个程序,通过 Legendre-Gauss 求积求解积分。 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 阶的方式运行正交以获得高精度,但是这个根发现会导致失败。有任何想法吗?
最佳答案
由于 np.roots 依赖于文档所述的“查找伴随矩阵的特征值”,您可能会遇到错误传播问题,导致根上的虚部非零。也许您可以使用 np.real 函数丢弃虚部。
您可以尝试使用根的泰勒近似计算根的不同方法:
https://math.stackexchange.com/questions/12160/roots-of-legendre-polynomial
关于python - 在python中找到勒让德多项式的根,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/11794937/