Polinomios de Legendre

Ecuación diferencial ( 1 x 2 ) d 2 y d x 2 2x dy dx +n(n+1)y=0n=1,2,3...
Fórmula de Rodríguez P n (x)= 1 2 n n! d n d x n [ ( x 2 1 ) n ]
Relación de recurrencia (n+1) P n+1 (x)=(2n+1)x P n (x)n P n1 (x)
Ortogonalidad 1 1 P m (x) P n (x)dx= 2 2n+1 δ mn

Los primeros polinomios de Legendre son.

P 0 (x)=1 P 1 (x)=x P 2 (x)= 1 2 ( 3 x 2 1 ) P 3 (x)= 1 2 ( 5 x 3 3x ) P 4 (x)= 1 8 ( 35 x 4 30 x 2 +3 ) P 5 (x)= 1 8 ( 63 x 5 70 x 3 +15x ) P 6 (x)= 1 16 ( 231 x 6 315 x 4 +105 x 2 5 ) P 7 (x)= 1 16 ( 429 x 7 693 x 5 +315 x 3 35x ) .....

Mediante del comando diff de MATLAB, obtenemos la expresión de P4(x) y P5(x) comprobando las relaciones de ortogonolidad. Comprobamos también, que la integral

1 1 P n (x)dx=0

>> syms x;
>> n=4;
>> P4=diff((x^2-1)^n,n)/(2^n*factorial(n));
>> P4=simplify(P4)
P4 =(35*x^4)/8 - (15*x^2)/4 + 3/8
>> n=5;
>> P5=diff((x^2-1)^n,n)/(2^n*factorial(n));
>> P5=simplify(P5)
P5 =(x*(63*x^4 - 70*x^2 + 15))/8

>> int(P4*P5,x,-1,1)
ans =0
>> int(P5*P5,x,-1,1)
ans =2/11

>> int(P4,x,-1,1)
ans =0
>> int(P5,x,-1,1)
ans =0

Representamos gráficamente los polinomios de P1(x) a P5(x) utilizando el comando ezplot

syms x;
hold on
for n=1:5
    y=diff((x^2-1)^n,n)/(2^n*factorial(n));
    ezplot(y,[-1,1]);
end
hold off
grid on
xlabel('x')
ylabel('P_n(x)')
title('Polinomios de Legendre')

Creamos la función legendre_p para obtener los coeficientes de los polinomios de Legendre, utilizando la relación de recurrencia

function p=legendre_p(n)
    p1=1;    
    p2=[1,0]; 
    if n==0
        p=p1; %P0
    elseif n==1
        p=p2; %P1
    else
        for i=2:n
            p=((2*(i-1)+1)*[p2,0]-(i-1)*[0,0,p1])/i;
            p1=p2;
            p2=p;
        end
    end
end

Generamos los polinomios de Legendre y comprobamos sus propiedades, etc. La función poly2sym convierte el vector de los coeficientes de un polinomio en un polinomio simbólico p(x).

>> syms x;
>> P4=poly2sym(legendre_p(4))
P4 =(35*x^4)/8 - (15*x^2)/4 + 3/8
>> P5=poly2sym(legendre_p(5))
P5 =(63*x^5)/8 - (35*x^3)/4 + (15*x)/8
>> int(P4*P5,x,-1,1)
ans =0
>> int(P5*P5,x,-1,1)
ans =2/11

Creamos la función recursiva legendre_r para obtener el valor de los polinomios cuando se proporciona el valor de x utilizando la relación de recurrencia, partiendo de P0(x) =0 y P1(x) =x

function res=legendre_r(n,x)
    if n==0
        res=ones(1,length(x));
    elseif n==1
        res=x;
    else       
        res=((2*(n-1)+1)*x.*legendre_r(n-1,x)-(n-1)*legendre_r(n-2,x))/n;
    end
end

Las llamadas a las funciones legendre_r y legendre_p, para calcular el valor del polinomio n=4 para x=0.5, P4(0.5)

>> x=0.5;
>> legendre_r(4,x)
ans =   -0.2891
>> polyval(legendre_p(4),x)
ans =   -0.2891

producen los mismos resultados.

La función MATLAB legendre dipone de varias posibilidades, véase en la ayuda, comando >> help legendre. Llamamos a esta función para calcular el valor del polinomio P4(x) para un vector x de abscisas comprendidas entre -1 y +1. En la primera fila de la matriz que devuelve la función legendre, con la opción 'sch', se encuentran los valores de dicho polinomio P4(x), como podemos comprobar en el código que sigue. Alternativamente, podemos utilizar nuestras propias funciones legendre_p y legendre_r que hemos definido anteriormente.

>> x=-1:0.25:1;
>> p=legendre(4,x,'sch');
>> p(1,:) %primera fila de la matriz que devuelve legendre 
ans=1.0000 -0.3501 -0.2891 0.1577 0.3750 0.1577 -0.2891 -0.3501 1.0000
>> polyval(legendre_p(4),x)
ans=1.0000 -0.3501 -0.2891 0.1577 0.3750 0.1577 -0.2891 -0.3501 1.0000
>> legendre_r(4,x)
ans=1.0000 -0.3501 -0.2891 0.1577 0.3750 0.1577 -0.2891 -0.3501 1.0000

Representamos gráficamente los polinomios de P1(x) a P6(x) empleando el comando fplot. Utilizamos indistintamente las funciones legrendre_p que devuelve los coeficientes del polinomio en combinación con polyval, o bien, la función recursiva legendre_r

for n=1:6
    %f=@(x) polyval(legendre_p(n),x);
    f=@(x) legendre_r(n,x);
    fplot(f,[-1,1]);
end
xlabel('x')
ylabel('P_n(x)')
title('Polinomios de Legendre')
grid on
hold off

En la gráfica vemos que

Pn(1)=1

Pn(-1)=(-1)n

Si n es impar Pn(x)=-Pn(-x), de modo que Pn(0)=0

Comprobamos numéricamente mediante la función integral que el área bajo Pn(x) es nula, para n>1

>> for n=1:6
    integral(@(x)legendre_r(n,x),-1,1)
end
ans =  -8.3267e-17
ans =   4.1633e-17
ans =     0
ans =  -6.9389e-18
ans =   3.4694e-17
ans =  -9.3675e-17

Aproximación de una función mediante los polinomios de Legendre

Aproximamos la función f(x) mediante la suma de polinomios de Legendre multiplicados por un coeficiente cn. Las propiedades de ortogonalidad de los polinomios nos permiten calcular cada uno de los coeficientes, como vamos a ver a continuación

f(x)= n=0 c n P n (x) 1 1 f(x) P m (x)dx= n=0 c n 1 1 P n (x) P m (x) dx 1 1 f(x) P m (x)dx= 2 2m+1 c m c m = 2m+1 2 1 1 f(x) P m (x)dx

Vamos a aproximar la función f(x)=sgn(x) (función signo de MATLAB) de la figura a la combinación lineal de los N+1 primeros polinomios de Legendre

Para calcular cada uno de los coeficientes cn (n=0...N) tenemos que definir la función integrando, f(xPn(x) que denominaremos f_legendre. Utilizamos la función integral para calcular numéricamente la integral definida entre los límites -1 y 1.

function res=f_legendre(n,x)
    res=sign(x).*polyval(legendre_p(n),x);
end

sign es una función de MATLAB que devuelve el signo de x.

Primero, calculamos los coeficientes cn y después la aproximación a la función f(x)

f(x) n=0 N c n P n (x)

N=10;
c=zeros(1,N+1);
for n=0:N
    c(n+1)=integral(@(x)f_legendre(n,x),-1,1)*(n+1/2);
end

x=-1:0.02:1;
y=zeros(1,length(x));
for n=0:N
    y=y+c(n+1)*polyval(legendre_p(n),x);
end
hold on
line([-1,0],[-1,-1])
line([0,1],[1,1])
plot(x,y,'r')
hold off
xlabel('x')
ylabel('f(x)')
title('Aproximación mediante polinomios de Legendre')
grid on

En color azul se dibuja la función f(x) y en color rojo la aproximación mediante la combinación lineal de los N+1=11 primeros polinomios de Legendre.