Orthogonalitätsproblem in Scipys legendre Polynomen

Question

Ich bin vor kurzem auf ein merkwürdiges Problem bezüglich der scipy.special.legendre() (scipy documentation) gestoßen. Die Legenden-Polynome sollten paarweise orthogonal sein. Wenn ich sie jedoch über einen Bereich x=[-1,1] berechne und das Skalarprodukt von zwei Polynomen unterschiedlichen Grades erstelle, werde ich nicht immer Null oder Werte nahe Null erhalten. Interpretiere ich das Funktionsverhalten falsch? Im Folgenden möchte ich ein kurzes Beispiel geschrieben haben, die das Skalarprodukt bestimmte Paare von Legendre-Polynome erzeugt:

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

Die Handlung der einzelnen Polynome fein tatsächlich aussieht:

Aber wenn ich berechnen manuell die Skalarprodukte - multiplizieren zwei Legendre-Polynome verschiedenen Grade element und sUM sie (die 500 für die Normalisierung ist) ...

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

... I erhalten Sie die folgenden Werte als Ausgabe:

0vs0: +1.000000e+00 
0vs1: -5.906386e-17 
0vs2: +2.004008e-03 
0vs3: -9.903189e-17 
0vs4: +2.013360e-03 
0vs5: -1.367795e-16 

Das Skalarprodukt des ersten Polynom mit sich selbst gleich ein (wie zu erwarten), und die Hälfte der anderen Ergebnisse sind fast Null, aber es gibt einige Werte in der Bestellung von 10e-3 und ich habe keine Ahnung warum. Ich versuchte auch die scipy.special.eval_legendre(n, x) Funktion - das gleiche Ergebnis: - \

Ist das ein Fehler in der scipy.special.legendre() Funktion? Oder mache ich etwas falsch? Ich bin auf der Suche nach konstruktiven antwortet :-)

prost, Markus

Answer 1

Wie andere haben gesagt, Sie gehen, um etwas Fehler zu bekommen, da Sie ein in exakte Integral Leistung erbringt.

Aber Sie können den Fehler reduzieren, indem Sie das Integral so gut Sie können. In Ihrem Fall können Sie Ihre Stichprobenpunkte noch verbessern, um das Integral genauer zu machen. Verwenden Sie bei der Probenahme den "Mittelpunkt" der Intervalle anstelle der Kanten:

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 

Dies gibt eine ganze Menge Verbesserung! Wenn ich die alte Sampling-Methode verwenden:

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) 

Das gibt:

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

Fehler Begriffe sind ~ 10^-3.

Aber mit dem "Mittelpunkt Abtastschema", erhalte ich:

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

Fehler sind jetzt ~ 10^-5 oder sogar 10^-6, was viel besser ist!