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 多项式相乘并将它们相加（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()函数中的错误吗？还是我做错了什么？我正在寻找建设性的回应:-)

干杯，马库斯

4

1 回答 1

1

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

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

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。

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

[[  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，这要好得多！

于 2016-09-17T00:02:52.670 回答