scipy의 legendre 다항식의 직교성 문제

Question

나는 최근 scipy.special.legendre() (scipy documentation)과 관련된 호기심 문제를 접했습니다. legendre 다항식은 쌍으로 직교해야합니다. 그러나, 범위 x=[-1,1] 이상의 그들을 계산하고 내가 항상 제로 또는 0에 가까운 값을 얻을하지 않는 정도의 두 다항식의 스칼라 제품을 빌드 할 때. 함수 동작을 잘못 해석합니까? 의 I는 르장 드르 다항식의 특정 쌍 스칼라 제품을 생산하는 간단한 예제를 작성한 다음에 :

하지만 계산의 경우 : 하나의 다항식의 음모가 실제로 잘 보이는

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

스칼라 제품은 수동 - ... 다른 학위 elementwise 두 르장 드르 다항식을 곱 (500 정상화위한 것입니다)을 요약

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

... I 출력으로서 다음 값을 얻을 :

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, 훨씬 더이다!