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

Question

我最近偶然发现了一个关于 scipy.special.legendre() (scipy documentation) 的奇怪问题。勒让德多项式应该成对正交。但是，当我在 x=[-1,1] 范围内计算它们并构建两个不同次数的多项式的标量积时，我并不总是得到零或接近零的值。我是否误解了函数的行为？在下面我写了一个简短的例子，它产生了某些勒让德多项式对的标量积：

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()

单个多项式的绘图实际上看起来不错：

但是如果我手动计算标量积——将两个不同次数的勒让德多项式逐元素相乘并将它们相加（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() 函数中的错误吗？或者我做错了什么？我正在寻找建设性的回应:-)

干杯，马库斯

Answer 1

正如其他人评论的那样，您会遇到一些错误，因为您正在执行不精确的积分。

但是您可以通过尽可能地进行积分来减少误差。在您的情况下，您仍然可以改进采样点以使积分更精确。采样时，使用区间的"midpoint"代替边缘：

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)

这给出：

[[  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。

但是使用 "midpoint sampling scheme"，我得到：

[[  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，好多了！

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

Orthogonality issue in scipy's legendre polynomials

python

orthogonal

scipy

polynomials