Проблема ортогональности в полиномах легенд легенд Scipy

Question

Недавно я наткнулся на любопытный вопрос относительно scipy.special.legendre() (scipy documentation). Полиномы legendre должны быть попарно ортогональными. Однако, когда я вычисляю их в диапазоне 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 

Скалярное произведение первого многочлена с самим собой (как и следовало ожидать), равной один, и половина других результатов практически равны нулю, но есть некоторые значения в порядка 10e-3, и я понятия не имею, почему. Я также пробовал функцию scipy.special.eval_legendre(n, x) - тот же результат: - \

Является ли это ошибкой в ​​функции scipy.special.legendre()? Или я делаю что-то неправильно? Ищу конструктивные откликается :-)

приветствий, Маркус

Answer 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, что намного лучше!