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
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
de modo que la integral
sea mínima. Para calcular el mínimo, se deriva respecto a a1, a2...an+1
Obtenemos el siguiente sistema de n ecuaciones con n incógnitas
o bien,
Teniendo en cuenta la ortogonalidad de las funciones φk(x). El sistema de ecuaciones se reduce considerablemente
donde
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
Estos polinomios son ortogonales respecto al peso w(x)=1
Aproximamos la función f(x) en el intervalo [-1,1] a la combinación lineal Pn(x) de polinomios de Legendre
Los coeficientes ak se calculan del siguiente modo
Ejemplo 1
Aproximamos la función f(x)=exp(x) al polinomio P3(x)
Integrando por partes, o con la ayuda de Math Symbolic, calculamos los coeficientes a0, a1, a2 y a3
El polinomio P3(x) es
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')
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=legendreP(k, x); a=int(p*exp(x),x,-1,1)*(2*k+1)/2; pol=pol+p*a; end hold on fplot(@(x) exp(x),[-1,1]); fplot(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)
Calculamos los coeficientes a0, a1, a2, a3, a4 y a5
El polinomio P5(x) es
syms x; n=5; pol=1/2; for k=1:n p=legendreP(k,x); 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]) fplot(pol,[-1,1],'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=legendreP(k,x); a=int(p*abs(x),x,-1,1)*(2*k+1)/2; pol=pol+p*a; end hold on fplot(abs(x),[-1,1]); fplot(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)
Los polinomios son ortogonales respecto al peso
Aproximamos la función f(x) en el intervalo [-1,1] a la combinación lineal Pn(x) de polinomios de Chebyshev.
Los coeficientes ak se calculan mediante la fórmula
Ejemplo 1
Aproximamos una semicircunferencia de radio r=1 a una combinación lineal P2(x) de polinomios de Chebyshev
La función que describe la semicircunferencia de radio r=1, es
Los coeficientes ak se calculan del siguiente modo
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)
Calculamos los coeficientes a0, a1, a2, a3, a4 y a5
A partir de la definición de polinomio de Chebyshev Tk(θ)=cos(kθ) y haciendo el cambio de variable x=cos(θ) los coeficientes ak se calculan del siguiente modo
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
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*chebyshevT(k,x); end hold on line([0,1],[1,1]); line([-1,0],[0,0]) fplot(pol,[-1,1],'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,
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(kθ) y haciendo el cambio de variable x=cos(θ) los coeficientes ak se escriben
Llevamos a cabo la integración numérica, mediante la función MATLAB,
>> 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+chebyshevT(k,x)*a; end hold on fplot(exp(x),[-1,1]); fplot(pol,[-1,1]); legend('exponencial','Chebyshev','location','northwest') hold off grid on xlabel('x') ylabel('y') title('Aproximación de e^x')