Polinomios de Chebyshev

Definimos un conjunto de polinomios Tn(x)=cos() de grado n, donde cosθ=x en el intervalo -1≤x≤1.

Para determinar la forma de estos polinomios recordaremos la fórmula cos(a+b)=cos(a)cos(b)-sin(a)sin(b)

{ cos( (n+1)θ )=cos( nθ )cosθsin( nθ )sinθ cos( (n1)θ )=cos( nθ )cosθ+sin( nθ )sinθ cos( (n+1)θ )+cos( (n1)θ )=2cos( nθ )cosθ T n+1 (x)+ T n1 (x)=2x T n (x)

Los primeros polinomios son:

T0(x)=cos(0)=1
T1(x)=cosθ=x;

El resto de polinomios se obtiene empleando la relación de recurrencia:

Tn+1(x)=2xTn(x)-Tn-1(x)

T2(x)=cos(2θ)=2x2-1
T3(x)=cos(3θ)=2x(2x2-1)-x=4x3-3x
T4(x)= cos(4θ)=2x(4x3-3x)-( 2x2-1)=8x4-8x2+1
T5(x)= cos(5θ)=2x(8x4-8x2+1)-(4x3-3x)=16x5-20x3+5x
T6(x)= .....

Creamos una función recursiva que nos devuelve el valor del polinomio Tn(x) de grado n para valores dados de x.

function y=chebyshev(n,x) %n=1,2,3...
    if n==0
        y=ones(1,length(x));
    elseif n==1
        y=x;
    else
        y=2*x.*chebyshev(n-1,x)-chebyshev(n-2,x);
    end
end

Otra forma equivalente de crear una función que nos devuela el valor del polinomio Tn(x) de grado n para valores dados de x es la siguiente:

function T=chebyshev(n,x) %n=0,1,2,3,...
    m=length(x);
    T=zeros(1,m);
    T1=ones(1,m);
    T2=x;
    if n==1
        T=T2;
    else
        for j=2:n
            T=2*x.*T2-T1;
            T1=T2;
            T2=T;
        end
    end
end

Mediante el siguiente script, representamos gráficamente los polinomios T1(x), T2(x), T3(x) y T4(x)

hold on
for n=1:4
    fplot(@(x) chebyshev(n,x),[-1,1], 'displayName',num2str(n))
end
title('Polinomios de Chebyshev')
legend('-DynamicLegend','location','southeast')
xlabel('x')
ylabel('y')
grid on
hold off 

Se puede obtener los polinomios de Chebyshev T0(x), T1(x), T2(x), T3(x), T4(x) para valores de x mediante el siguiente código.

n=5;    
x=-1:0.01:1; %vector de valores de x
m=length(x);
T=zeros(m,n);
T(:,1)=ones(1,m);  %T0(x)
T(:,2)=x;          %T1(x)
hold on
for j=3:n           %de T2(x) a T4(x)
    T(:,j)=2*x'.*T(:,j-1)-T(:,j-2);
end

Creamos una función que nos devuelve el vector de los coeficientes, [a1,a2,...an,an+1], (de mayor a menor grado) del polinomio Tn(x).

function p=chebypoly(n) %n es grado 0,1,2...
    p=zeros(n+1,1);
    p1=[1];
    p2=[1 0];
    if n==0
        p=p1;
     elseif n==1
       p=p2;
    else
       for i=2:n
            p=2*[p2 0]-[0 0 p1];
            p1=p2;
            p2=p;
        end   
    end       
end

La función poly2sym de Math Symbolic convierte los coeficientes del polinomio en la expresión algebraica a1xn+a2xn-1+...anx+an+1. La función ezplot representa gráficamente los polinomios

syms x;
hold on
for n=1:4
     Tk=poly2sym(chebypoly(n));
     ezplot(Tk,[-1,1]);
end
hold off
title('Polinomios de Chebyshev')
xlabel('x')
ylabel('y')
grid on

Esta misma reepresentación se puede obtener utilizando la función chebyshevT de MATLAB

hold on
fplot(chebyshevT(1:4,x),[-1,1]);
hold off
title('Polinomios de Chebyshev')
xlabel('x')
ylabel('y')
grid on

Propiedades de los polinomios de Chebyshev

Los polinomios de Chebyshev están definidos en el intervalo [-1,1].

Simetría

Los polinomios de índice par son funciones pares f(x)=f(-x) y los polinomios de índice impar son funciones impares f(x)=-f(-x)

Raíces

Obtenemos las raíces del polinomio Tn(x) de tres formas distintas:

Representamos las raíces del polinomio de grado 5, T5(x), mediante puntos de color rojo a lo largo del eje X

n=5;
k=1:n;
ang=(2*k-1)*pi/(2*n);
xk=cos(ang);
yk=sin(ang);
hold on
plot(cosd(1:180),sind(1:180),'b')
for i=1:length(xk)
    line([0,xk(i)],[0,yk(i)])
    hg=line([xk(i),xk(i)],[0,yk(i)]);
    set(hg,'linestyle','--');
end
line([-1,1],[0,0])
plot(xk,0,'ro','markersize',4,'markerfacecolor','r')
hold off
axis equal
title('Raíces')

Máximos y mínimos

Los extremos (máximo o mínimo) del polinomio Tn(x)=cos() donde cosθ=x son 1 y -1, respectivamente.

d T n (x) dx =nsin( nθ ) dθ dx = nsin( nθ ) sinθ

se producen cuando sin()=0, es decir nθ=kπ, k=1,2…n-1. En las posiciones.

x k ' =cos( kπ n )k=1,2...n1

Se excluye k=0 y k=n por que

lim θ0,π sin(nθ) sinθ =n

>> k=1:4;
>> cos(k*pi/5)
ans = 0.8090    0.3090   -0.3090   -0.8090

>> T5=poly2sym(chebypoly(5)) 
T5 =16*x^5 - 20*x^3 + 5*x
>> dT5=diff(T5) %calcula la derivada
dT5 =80*x^4 - 60*x^2 + 5
>> r=solve(dT5)
r =
   5^(1/2)/4 + 1/4
   5^(1/2)/4 - 1/4
   1/4 - 5^(1/2)/4
 - 5^(1/2)/4 - 1/4
>> y=subs(T5,x,r); %valor del polinomio T5 para cada una de las raíces r.
>> simplify(y) 
ans =
 -1
  1
 -1
  1
>> double(r)
ans =
    0.8090
    0.3090
   -0.3090
   -0.8090

Relaciones entre dos polinomiosTi(x) y Tj(x),

1 1 T i (x) T j (x) 1 x 2 dx={ 0ij πi=j=0 π 2 i=j0

>> syms x;
>> T5=poly2sym(chebypoly(5))
T5 =16*x^5 - 20*x^3 + 5*x
>> T4=poly2sym(chebypoly(4))
T4 =8*x^4 - 8*x^2 + 1
>> f=T4*T5/sqrt(1-x^2);
>> int(f,x,-1,1)
ans =0
 
>> f=T5*T5/sqrt(1-x^2);
>> int(f,x,-1,1)
ans =pi/2
 
>> f=1/sqrt(1-x^2);
>> int(f,x,-1,1)
ans =pi

Relación de ortogonalidad

Como hemos visto, las raíces del polinomio de grado n, Tn(x) son

x k =cos( (2k1)π 2n )k=1,2....n

Sean dos polinomios, Ti(x) y Tj(x), i,j<n vamos a comprobar que se cumple la siguiente relación que es importante para la aproximación de funciones f(x)

k=1 n T i ( x k ) T j ( x k )={ 0ij ni=j=0 n 2 i=j0

Esta es una relación análoga a la existente en las series de Fourier

n=7;
k=1:n;
xk=cos((2*k-1)*pi/(2*n));
y5=chebyshev(5,xk);
y4=chebyshev(4,xk);
suma=sum(y5.*y4)
suma=sum(y5.*y5)
suma=sum(y4.*y4)

En la ventana de comandos obtenemos el siguiente resultado

suma = -2.4980e-015
suma =    3.5000
suma =    3.5000

De otra forma, mediante la función chebypoly

>> n=7;
>> k=1:n;
>> xk=cos((2*k-1)*pi/(2*n));
>> T5=chebypoly(5);
>> T4=chebypoly(4);
>> suma=sum(polyval(T5,xk).*polyval(T4,xk))
suma = -1.2768e-015
>> suma=sum(polyval(T5,xk).*polyval(T5,xk))
suma =    3.5000
>> suma=sum(polyval(T4,xk).*polyval(T4,xk))
suma =    3.5000

Utilizando Math Symbolic obtenemos valores exactos

>> syms x;
>> T5=poly2sym(chebypoly(5));
>> T4=poly2sym(chebypoly(4));
>> n=sym('7');  %raices del polinomio T7
>> k=1:n;
>> xk=cos((2*k-1)*pi/(2*n));
>> suma=sum(subs(T5,x,xk).*subs(T4,x,xk))
suma =0 
>> suma=sum(subs(T5,x,xk).*subs(T5,x,xk));
>> simplify(suma)
ans =cos(pi/7) - cos((2*pi)/7) + cos((3*pi)/7) + 3
>> double(ans)
ans =    3.5000  %que es n/2
>> suma=sum(subs(T4,x,xk).*subs(T4,x,xk));
>> simplify(suma)
ans =cos((2*pi)/7) - cos(pi/7) - cos((3*pi)/7) + 4
>> double(ans)
ans =    3.5000 %que es n/2

Polinomios de Chebyshev de segunda especie

Definimos un conjunto de polinomios

U n (x)= sin( (n+1)θ ) sinθ

grado n, donde cosθ=x en el intervalo -1≤x≤1.

Para determinar la forma de estos polinomios recordaremos la fórmula sin(a+b)=sin(a)cos(b)+sin(b)cos(a)

sin( (n+2)θ )=sin( (n+1)θ )cosθ+cos( (n+1)θ )sinθ= sin( (n+1)θ )cosθ+{ cos( nθ )cosθsin( nθ )sinθ }sinθ= sin( (n+1)θ )cosθ+cos( nθ )cosθsinθsin( nθ ) sin 2 θ= sin( (n+1)θ )cosθ+cos( nθ )cosθsinθ+sin( nθ ) cos 2 θsin( nθ )

Agrupando términos

sin( (n+2)θ )+sin( nθ )=sin( (n+1)θ )cosθ+{ sin( nθ )cosθ+cos( nθ )sinθ }cosθ= =2cosθ·sin( (n+1)θ )

Dividiendo entre sinθ

sin( (n+2)θ ) sinθ + sin( nθ ) sinθ =2cosθ sin( (n+1)θ ) sinθ U n+1 (x)+ U n1 (x)=2x U n (x)

Esta es la relación de recurrencia, para utilizarla hemos de determinar lo dos primeros polinomios, U0(x) y U1(x)

Los primeros polinomios son:

U 0 (x)= sin( θ ) sinθ =1 U 1 (x)= sin( 2θ ) sinθ =2cosθ=2x

U2(x)=2x·U1(x)-U0(x)=4x2-1

Alternativamente, podríamos calcular U2(x) a partir de su definición

U 2 (x)= sin( 3θ ) sinθ = sin( 2θ )cosθ+cos( 2θ )sinθ sinθ = 2sinθ cos 2 θ+ cos 2 θsinθ sin 3 θ sinθ = 2sinθ cos 2 θ+ cos 2 θsinθ(1 cos 2 θ)sinθ sinθ =4 x 2 1

Creamos una función recursiva que nos devuelve el valor del polinomio Un(x) de grado n para valores dados de x.

function y=chebyshev_2(n,x) %n=0,1,2,3...
    if n==0
        y=ones(1,length(x));
    elseif n==1
        y=2*x;
    else
        y=2*x.*chebyshev_2(n-1,x)-chebyshev_2(n-2,x);
    end
end

Mediante el siguiente script, representamos gráficamente los polinomios U0(x), U1(x), U2(x), U3(x) y U4(x)

hold on
for n=0:4
    fplot(@(x) chebyshev_2(n,x),[-1,1], 'displayName',num2str(n))
end
title('Polinomios de Chebyshev')
legend('-DynamicLegend','location','southeast')
xlabel('x')
ylabel('y')
grid on
hold off 

Una representación gráfica similar se obtiene con la función chebyshevU de MATLAB

syms x;
hold on
fplot(chebyshevU(0:4,x),[-1,1]);
hold off
title('Polinomios de Chebyshev')
xlabel('x')
ylabel('y')
grid on

En la figura, vemos la simetría

U n (x)= ( 1 ) n U n (x)

Cuando n es un número par, 0,2,4... la función es simétrica respecto del eje Y. Cuando n es un número impar, 1,3,5...la función es antisimétrica

Cuando x=cosθ=1, θ=0, 2π,...,

U n (1)= lim θ0 sin( (n+1)θ ) sinθ =n+1

>> syms n x;
>> limit(sin((n+1)*x)/sin(x),x,0)
ans =n + 1

Teniendo en cuenta la simetría de las funciones Un(x)

U n (1)= ( 1 ) n ( n+1 )

Creamos una función que nos devuelve el vector de los coeficientes, [a1,a2,...an,an+1], (de mayor a menor grado) del polinomio Un(x).

function p=chebypoly_2(n) %n es grado 1,2...
    p=zeros(n,1);
    p0=[0];
    p1=[1];
    for i=2:n+1
        p=2*[p1,0]-[0,p0];
        p0=[0,p1];
        p1=p;
    end
end

La función poly2sym de Math Symbolic convierte los coeficientes del polinomio en la expresión algebraica a1xn+a2xn-1+...anx+an+1. La función ezplot representa gráficamente los polinomios

hold on
for n=1:4
     Un=poly2sym(chebypoly_2(n));
     ezplot(Un,[-1,1]);
end
hold off
ylim([-4,5])
title('Polinomios de Chebyshev')
xlabel('x')
ylabel('y')
grid on

Raíces

Obtenemos las raíces del polinomio Un(x) de tres formas distintas:

Relaciones entre dos polinomios Ui(x) y Uj(x),

1 1 U i (x) · U j (x) 1 x 2 dx={ 0ij π 2 i=j

>> syms x;
>>  U5=poly2sym(chebypoly_2(5))
U5 =32*x^5 - 32*x^3 + 6*x
>>  U4=poly2sym(chebypoly_2(4))
U4 =16*x^4 - 12*x^2 + 1

>>  f=U4*U5*sqrt(1-x^2);
>>  int(f,x,-1,1)
ans =0
 
>>   f=U5*U5*sqrt(1-x^2);
>>  int(f,x,-1,1)
ans =pi/2 

Como hemos visto, las raíces del polinomio de grado n, Un(x) son

x m =cos( mπ n+1 )m=1,2...n

Sean dos polinomios, Ui(x) y Uj(x), i,j<n vamos a comprobar que se cumple la siguiente relación

m=1 n U i ( x m ) U j ( x m )( 1 x m 2 )={ 0ij n+1 2 i=j

n=7;
m=1:n;
xm=cos(m*pi/(n+1)); %raíces
y5=chebyshev_2(5,xm);
y4=chebyshev_2(4,xm);
suma=sum((y5.*y4).*(1-xm.^2))
suma=sum((y5.*y5).*(1-xm.^2))
suma=sum((y4.*y4).*(1-xm.^2))
suma =  -8.8818e-16
suma =    4.0000
suma =    4.0000