scipyのlegendre多項式の直交性の問題

Question

私は最近、scipy.special.legendre()（scipy documentation）に関する興味深い問題に遭遇しました。 legendre多項式はpairwise直交する必要があります。しかし、私がそれらを範囲x=[-1,1]で計算し、異なる次数の2つの多項式のスカラー積を構築すると、常にゼロに近い値やゼロに近い値が得られるわけではありません。関数の振る舞いを誤解していますか？単一の多項式のプロットが実際に正常に見える

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

：

をしかし、私は計算した場合、次の私はルジャンドル多項式の特定のペアの内積を生産する簡単な例を、書かれているで スカラー積手動 - 異なる程度の要素ごとの2つのルジャンドルの多項式を乗算し、（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 

自体と第1の多項式のスカラー積である（予想される）、及び半他の結果の同じ一つはほぼゼロであるが、いくつかの値は、です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、はるかに優れているです！