Aproximación de una función a polinomios ortogonales

Polinomios ortogonales

Un conjunto de funciones φ0(x), φ1(x),....φn(x) continuas en [a,b] se denominan ortogonales, si

a b w(x)· φ j (x) · φ i (x)dx={ 0siij C i sii=j

donde Ci es un número positivo, i,j=0,1,2...n, y w(x) se denomina función peso

Vamos a ajustar una función f(x) a una combinación lineal de n polinomios ortogonales φk(x) en el intervalo [a,b].

Dada una función f(x) continua en [a,b], encontraremos el polinomio Pn(x) calculando los coeficientes a0, a1...an

P n (x)= a 0 φ 0 (x)+ a 1 φ 1 (x)+...+ a n φ n (x)

de modo que la integral

E= a b w(x) ( P n (x)f(x) ) 2 dx

sea mínima. Para calcular el mínimo, se deriva respecto a a1, a2...an+1

E a i =0i=1,2,...n+1

Obtenemos el siguiente sistema de n ecuaciones con n incógnitas

1 2 E a 0 =0 a b w(x)· φ 0 (x)( a 0 φ 0 (x)+ a 1 φ 1 (x)+...+ a n φ n (x)f(x) ) dx=0 1 2 E a 1 =0 a b w(x)· φ 1 (x)( a 0 φ 0 (x)+ a 1 φ 1 (x)+...+ a n φ n (x)f(x) ) dx=0 ... 1 2 E a n =0 a b w(x)· φ n (x)( a 0 φ 0 (x)+ a 1 φ 1 (x)+...+ a n φ n (x)f(x) ) dx=0

o bien,

a 0 a b w(x)· φ 0 (x)· φ 0 (x)dx + a 1 a b w(x)· φ 0 (x) φ 1 (x)dx +...+ a n a b w(x)· φ 0 (x) φ n (x)dx = a b w(x)· φ 0 (x)f(x)dx a 0 a b w(x)· φ 1 (x) φ 0 (x)dx + a 1 a b w(x)· φ 1 (x) φ 1 (x)dx +...+ a n a b w(x)· φ 1 (x) φ n (x)dx = a b w(x)· φ 1 (x)f(x)dx ... a 0 a b w(x)· φ n (x) φ 0 (x)dx + a 1 a b w(x)· φ n (x) φ 1 (x)dx +...+ a n a b w(x)· φ n (x) φ n (x)dx = a b w(x)· φ n (x)f(x)dx

Teniendo en cuenta la ortogonalidad de las funciones φk(x). El sistema de ecuaciones se reduce considerablemente

a 0 C 0 = a b w(x)· φ 0 (x)f(x)dx a 1 C 1 = a b w(x)· φ 1 (x)f(x)dx ... a n C n = a b w(x)· φ n (x)f(x)dx

donde

C k = a b w(x)· φ k 2 (x)dx k=0,1,2...n

Consideraremos dos conjuntos de polinomios ortogonales: los polinomios de Legendre y los polinomios de Chebyshev

Aproximación a polinomios de Legendre

Los primeros polinomios de Legendre son

φ 0 (x)=1 φ 1 (x)=x φ 2 (x)= 1 2 ( 3 x 2 1 ) φ 3 (x)= 1 2 ( 5 x 3 3x ) φ 4 (x)= 1 8 ( 35 x 4 30 x 2 +3 ) φ 5 (x)= 1 8 ( 63 x 5 70 x 3 +15x )

Estos polinomios son ortogonales respecto al peso w(x)=1

1 1 φ i (x) φ j (x)dx ={ 0ij C j = 2 2j+1 i=j

Aproximamos la función f(x) en el intervalo [-1,1] a la combinación lineal Pn(x) de polinomios de Legendre

P n (x)= a 0 φ 0 (x)+ a 1 φ 1 (x)+...+ a n φ n (x)

Los coeficientes ak se calculan del siguiente modo

a k = 1 C k 1 1 f(x) φ k (x)dx C k = 2 2k+1 k=0,1,2...n

Ejemplo 1

Aproximamos la función f(x)=exp(x) al polinomio P3(x)

P 3 (x)= a 0 φ 0 (x)+ a 1 φ 1 (x)+ a 2 φ 2 (x)+ a 3 φ 3 (x)

Integrando por partes, o con la ayuda de Math Symbolic, calculamos los coeficientes a0, a1, a2 y a3

a 0 = 1 2 1 1 e x dx= 1 2 ( e 1 e ) a 1 = 3 2 1 1 x e x dx= 3 e a 2 = 5 2 1 1 1 2 ( 3 x 2 1 ) e x dx= 5 2 ( e 7 e ) a 3 = 7 2 1 1 1 2 ( 5 x 3 3x ) e x dx= 1 2 ( 259 e 35e )

El polinomio P3(x) es

P 3 (x)= 1 2 ( e 1 e ) φ 0 (x)+ 3 e φ 1 (x)+ 5 2 ( e 7 e ) φ 2 (x)+ 1 2 ( 259 e 35e ) φ 3 (x)

e=exp(1);
p=@(x) (e-1/e)/2+3*x/e+(5*(e-7/e)/2)*(3*x^2-1)/2+(259/e-35*e)*(5*x^3-3*x)/4;
hold on
fplot(@exp,[-1,1]);
fplot(p,[-1,1])
hold off
legend('exponencial','Legendre','location','northwest')
grid on
xlabel('x')
ylabel('y')
title('Aproximando exp(x) a un polinomio')

En la página Polinomios de Legendre, elaboramos una función denominada legendre_p, que nos devuelve los coeficientes de un polinomio Pn(x) de Legendre. Mediante la función MATLAB poly2sym convertimos los coeficientes en dicho plonomio en un polinomio simbólico de grado n

Elaboramos un script que nos permite aproximar la función ex en el intervalo [-1,1] a una combinación Pn(x) de polinomios de Legendre. Sea n=3

syms x;
n=3;
pol=int('exp(x)',x,-1,1)/2;
for k=1:n
    p=poly2sym(legendre_p(k));
    a=int(p*exp(x),x,-1,1)*(2*k+1)/2;
    pol=pol+p*a;
end
hold on
ezplot(exp(x),[-1,1]);
ezplot(pol,[-1,1]);
legend('exponencial','Legendre','location','northwest')
hold off
grid on
xlabel('x')
ylabel('y')
title('Aproximación de e^x')

El polinomio P3(x) ordenado en potencias de x es

>> collect(pol)
ans =((1295*exp(-1))/4 - (175*exp(1))/4)*x^3 + ((15*exp(1))/4 -
 (105*exp(-1))/4)*x^2 + ((105*exp(1))/4 - (765*exp(-1))/4)*x
 + (33*exp(-1))/4 - (3*exp(1))/4

Ejemplo 2

Aproximamos la función escalón f(x) al polinomio P5(x)

f(x)={ 01<x<0 10<x<1

Calculamos los coeficientes a0, a1, a2, a3, a4 y a5

a 0 = 1 2 0 1 1·dx = 1 2 a 1 = 3 2 0 1 x·dx = 3 4 a 2 = 5 2 0 1 1 2 ( 3 x 2 1 )·dx =0 a 3 = 7 2 0 1 1 2 ( 5 x 3 3x )·dx = 7 16 a 4 = 9 2 0 1 1 8 ( 35 x 4 30 x 2 +3 )·dx =0 a 5 = 11 2 0 1 1 8 ( 63 x 5 70 x 3 +15x )·dx = 11 32

El polinomio P5(x) es

P 5 (x)= 1 2 φ 0 (x)+ 3 4 φ 1 (x) 7 16 φ 3 (x)+ 11 32 φ 5 (x)

syms x;
n=5;
pol=1/2;
for k=1:n
    p=poly2sym(legendre_p(k));
    a=int(p,x,0,1)*(2*k+1)/2;
    pol=pol+p*a;
end
hold on
line([0,1],[1,1]);
line([-1,0],[0,0])
fg=ezplot(pol,[-1,1]);
set(fg,'color','r')
hold off
grid on
xlabel('x')
ylabel('y')
title('Aproximación de f(x)')

>> collect(pol)
ans =(693*x^5)/256 - (525*x^3)/128 + (525*x)/256 + 1/2

Ejemplo 3

Aproximamos la función f(x)=|x| al polinomio P4(x) en el intervalo [-1,1]

syms x;
n=4;
pol=int('abs(x)',x,-1,1)/2;
for k=1:n
    p=poly2sym(legendre_p(k));
    a=int(p*abs(x),x,-1,1)*(2*k+1)/2;
    pol=pol+p*a;
end
hold on
ezplot(abs(x),[-1,1]);
ezplot(pol,[-1,1]);
legend('|x|','Legendre','location','north')
hold off
grid on
xlabel('x')
ylabel('y')
title('Aproximación de |x|')

El polinomio P4(x) ordenado en potencias de x es

>> collect(pol)
ans =- (105*x^4)/128 + (105*x^2)/64 + 15/128

Aproximación a polinomios de Chebyshev

Las funciones base φk(x), son los polinomios de Chebyshev Tk(x)

T 0 (x)=1 T 1 (x)=x T 2 (x)=2 x 2 1 T 3 (x)=4 x 3 3x T 4 (x)=8 x 4 8 x 2 +1 T 5 (x)=16 x 5 20 x 3 +5x

Los polinomios son ortogonales respecto al peso

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

Aproximamos la función f(x) en el intervalo [-1,1] a la combinación lineal Pn(x) de polinomios de Chebyshev.

P n (x)= a 0 T 0 (x)+ a 1 T 1 (x)+...+ a n T n (x)

Los coeficientes ak se calculan mediante la fórmula

a k = 1 C k 1 1 f(x) 1 x 2 T k (x)dx C 0 =π C k = π 2 k=1,2...n

Ejemplo 1

Aproximamos una semicircunferencia de radio r=1 a una combinación lineal P2(x) de polinomios de Chebyshev

P 2 (x)= a 0 T 0 (x)+ a 1 T 1 (x)+ a 2 T 2 (x)

La función que describe la semicircunferencia de radio r=1, es

f(x)= 1 x 2

Los coeficientes ak se calculan del siguiente modo

a 0 = 1 π 1 1 dx= 2 π a 1 = 2 π 1 1 xdx=0 a 2 = 2 π 1 1 ( 2 x 2 1 )dx= 4 3π P 2 (x)= 2 π 4 3π ( 2 x 2 1 )

p=@(x) 2/pi-4*(2*x^2-1)/(3*pi);
f=@(x) sqrt(1-x^2);
hold on
fplot(f,[-1,1]);
fplot(p,[-1,1])
hold off
legend('semicircunferencia','Chebyshev','location','northwest')
grid on
axis equal
xlabel('x')
ylabel('y')
title('Aproximando f(x) a un polinomio')

El lector interesado puede calcular a mano o con ayuda de Math Symbolic los coeficientes a3=0, a4..., para obtener mejores aproximaciones

Ejemplo 2

Aproximamos la función escalón f(x) al polinomio P5(x)

f(x)={ 01<x<0 10<x<1

Calculamos los coeficientes a0, a1, a2, a3, a4 y a5

a k = 1 C k 0 1 T k (x) 1 x 2 dx

A partir de la definición de polinomio de Chebyshev Tk(θ)=cos() y haciendo el cambio de variable x=cos(θ) los coeficientes ak se calculan del siguiente modo

a 0 = 1 π 0 π/2 dθ = 1 2 a k = 2 π 0 π/2 cos( kθ )dθ

syms x;
for k=1:5
    f=cos(k*x);
    int(f,x,0,pi/2)*2/pi  
end
ans =2/pi
ans =0
ans =-2/(3*pi)
ans =0
ans =2/(5*pi)

El polinomio P5(x) es

P 5 (x)= 1 2 φ 0 (x)+ 2 π φ 1 (x) 2 3π φ 3 (x)+ 2 5π φ 5 (x)

syms x;
n=5;
pol=1/2;
for k=1:n
    f=cos(k*x);
    a=int(f,x,0,pi/2)*2/pi;
    pol=pol+a*poly2sym(chebypoly(k));   
end
hold on
line([0,1],[1,1]);
line([-1,0],[0,0])
fg=ezplot(pol,[-1,1]);
set(fg,'color','r')
hold off
grid on
xlabel('x')
ylabel('y')
title('Aproximación de f(x)')

>> collect(pol)
ans =(32*x^5)/(5*pi) - (32*x^3)/(3*pi) + (6*x)/pi + 1/2

Ejemplo 3

Aproximamos la función f(x)=exp(x) en el intervalo [-1,1], mediante la combinación lineal P3(x) de polinomios de Chebyshev,

P 3 (x)= a 0 T 0 (x)+ a 1 T 1 (x)+ a 2 T 2 (x)+ a 3 T 3 (x)

Para calcular los coeficientes ak, observamos que el integrando diverge en los extremos del intervalo [-1,1]. La integración numérica en términos de la variable x no parece la más apropiada. A partir de la definición de polinomio de Chebyshev Tk(θ)=cos() y haciendo el cambio de variable x=cos(θ) los coeficientes ak se escriben

a 0 = 1 π 0 π exp( cosθ )·dθ a k = 2 π 0 π exp( cosθ )·cos( kθ )·dθ k=1,2...n

Llevamos a cabo la integración numérica, mediante la función MATLAB, integral

>> f=@(x) exp(cos(x));
>> integral(f,0,pi)/pi
ans =    1.2661
>> f=@(x) exp(cos(x)).*cos(x);
>> 2*integral(f,0,pi)/pi
ans =    1.1303
>> f=@(x) exp(cos(x)).*cos(2*x);
>> 2*integral(f,0,pi)/pi
ans =    0.2715
>> f=@(x) exp(cos(x)).*cos(3*x);
>> 2*integral(f,0,pi)/pi
ans =    0.0443

El polinomio aproximador es
P3(x)=1.2661·T0(x)+1.1303·T1(x)+0.2715·T2(x)+0.0443·T3(x)

P3(x)=1.2661+1.1303·x+0.2715·(2x2-1)+0.0443·(4x3-3x)

p=@(x) 1.2661+1.1303*x+0.2715*(2*x^2-1)+0.0443*(4*x^3-3*x);
hold on
fplot(@exp,[-1,1]);
fplot(p,[-1,1])
hold off
legend('exponencial','Chebyshev','location','northwest')
grid on
xlabel('x')
ylabel('y')
title('Aproximando e^x a un polinomio')

Elaboramos un script que nos permita aproximar la función ex a Pn(x) en el intervalo [-1,1]. Sea por ejemplo, n=3

n=3;
f=@(x) exp(cos(x));
pol=integral(f,0,pi)/pi; %a_0
for k=1:n
    f=@(x) exp(cos(x)).*cos(k*x);
    a=integral(f,0,pi)*2/pi;   %a_k
    pol=pol+poly2sym(chebypoly(k))*a;
end
hold on
ezplot(exp(x),[-1,1]);
ezplot(pol,[-1,1]);
legend('exponencial','Chebyshev','location','northwest')
hold off
grid on
xlabel('x')
ylabel('y')
title('Aproximación de e^x')

Hemos definido la función chebypoly en la página Polinomios de Chebyshev, calcula y devuelve los coeficientes del polinomio Tk(x). La función MATLAB poly2sym convierte los coeficientes en un polinomio simbólico