Scipy中勒让德多项式的正交性问题

Question

Scipy中勒让德多项式的正交性问题

3

最近我在scipy.special.legendre() (scipy文档) 中遇到了一个奇怪的问题。Legendre多项式应该是成对正交的。但是，当我计算它们在x=[-1,1]范围内的值，并且计算不同次数的两个多项式的标量积时，我并不总是得到零或接近零的值。我是否误解了这个函数的行为？下面是一个简短的示例，可以生成某些Legendre多项式的标量积：

from __future__ import print_function, division
import numpy as np 
from scipy import special
import matplotlib.pyplot as plt

# create range for evaluation
x = np.linspace(-1,1, 500)

degrees = 6
lp_array = np.empty((degrees, len(x)))

for n in np.arange(degrees):
    LP = special.legendre(n)(x)
    # alternatively:
    # LP = special.eval_legendre(n, x)
    lp_array[n, ] = LP
    plt.plot(x, LP, label=r"$P_{}(x)$".format(n))

plt.grid()
plt.gca().set_ylim([-1.1, 1.1])
plt.legend(fontsize=9, loc="lower right")
plt.show()

单项式的图像实际上看起来不错：从0到5次的Legendre多项式

但如果我手动计算标量积——逐个元素相乘两个不同次数的Legendre多项式并将它们加起来（500用于归一化）……

for i in range(degrees):
    print("0vs{}: {:+.6e}".format(i, sum(lp_array[0]*lp_array[i])/500))

我得到以下值作为输出：

0vs0: +1.000000e+00
0vs1: -5.906386e-17
0vs2: +2.004008e-03
0vs3: -9.903189e-17
0vs4: +2.013360e-03
0vs5: -1.367795e-16

第一个多项式与自身的数量积（如预期）等于1，而另一半结果几乎为零，但是有一些值在10e-3的数量级，我不知道为什么。我还尝试了scipy.special.eval_legendre(n, x)函数，结果相同:-\ 是否这是scipy.special.legendre()函数中的一个错误？或者我做错了什么？我正在寻找建设性的回应:-) 问候，马库斯

- Markus

感谢您将图像集成进来 :-) - Markus

嗯，好的。我不太明白为什么我的方法不能正确地工作，但是把多项式乘积整合起来肯定是有道理的。感谢你提出这个观点。 - Markus

2

你的集成方法不准确。10^{-3} ~ 1/500是在使用500个点时你所期望的误差数量级。 - pv.

关于我的实际问题(在统计学上)：我有一个126维度的问题，这迫使我的勒让德多项式长度为126。因此，对于scipy.special.legendre来说，这显然太短了，无法得到适当的正交多项式。然而，这会导致在计算协方差矩阵时出现不同的方差值，具体取决于我考虑多少个勒让德多项式。有没有办法解决这个问题？ - Markus

为什么Legendre多项式只需要在126个点进行采样？样本数量与问题的维度有什么关系？（如果通过分析方法来解决问题，多项式不是从[-1,1]（无限多个点）到实数的函数吗？） - unutbu

显示剩余3条评论

1个回答

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

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

- Praveen · Accepted Answer

正如其他人所评论的那样，由于您进行了不精确的积分，因此您会得到一些错误。

但是您可以通过尽可能地进行积分来减少误差。在您的情况下，您仍然可以改进采样点以使积分更加准确。在采样时，请使用间隔的“中点”而不是边缘：

x = np.linspace(-1, 1, nx, endpoint=False)
x += 1 / nx   # I'm adding half a sampling interval
              # Equivalent to x += (x[1] - x[0]) / 2

这将会有很大的改进！如果我使用旧的采样方法：

nx = 500
x = np.linspace(-1, 1, nx)

degrees = 7
lp_array = np.empty((degrees, len(x)))

for n in np.arange(degrees):
    LP = special.eval_legendre(n, x)
    lp_array[n, :] = LP

np.set_printoptions(linewidth=120, precision=1)
prod = np.dot(lp_array, lp_array.T) / x.size
print(prod)

This gives:

[[  1.0e+00  -5.7e-17   2.0e-03  -8.5e-17   2.0e-03  -1.5e-16   2.0e-03]
 [ -5.7e-17   3.3e-01  -4.3e-17   2.0e-03  -1.0e-16   2.0e-03  -1.1e-16]
 [  2.0e-03  -4.3e-17   2.0e-01  -1.3e-16   2.0e-03  -1.0e-16   2.0e-03]
 [ -8.5e-17   2.0e-03  -1.3e-16   1.4e-01  -1.2e-16   2.0e-03  -1.0e-16]
 [  2.0e-03  -1.0e-16   2.0e-03  -1.2e-16   1.1e-01  -9.6e-17   2.0e-03]
 [ -1.5e-16   2.0e-03  -1.0e-16   2.0e-03  -9.6e-17   9.3e-02  -1.1e-16]
 [  2.0e-03  -1.1e-16   2.0e-03  -1.0e-16   2.0e-03  -1.1e-16   7.9e-02]]

错误项约为10的负3次方。

但是使用“中点抽样方案”，我得到：

[[  1.0e+00  -2.8e-17  -2.0e-06  -3.6e-18  -6.7e-06  -8.2e-17  -1.4e-05]
 [ -2.8e-17   3.3e-01  -2.8e-17  -4.7e-06  -2.7e-17  -1.1e-05  -4.1e-17]
 [ -2.0e-06  -2.8e-17   2.0e-01  -5.7e-17  -8.7e-06  -2.3e-17  -1.6e-05]
 [ -3.6e-18  -4.7e-06  -5.7e-17   1.4e-01  -2.1e-17  -1.4e-05  -5.3e-18]
 [ -6.7e-06  -2.7e-17  -8.7e-06  -2.1e-17   1.1e-01   1.1e-17  -2.1e-05]
 [ -8.2e-17  -1.1e-05  -2.3e-17  -1.4e-05   1.1e-17   9.1e-02   7.1e-18]
 [ -1.4e-05  -4.1e-17  -1.6e-05  -5.3e-18  -2.1e-05   7.1e-18   7.7e-02]]

现在的错误率已经达到了10的负5次方甚至是10的负6次方，这比之前好多了！