Aproximación de funciones (I)
Aproximación de funciones
De modo análogo a las series de Fourier, expresemos la función f(x) definida en el intervalo [-1,1] como combinación lineal de m polinomios Tj(x)
Como hemos visto, las raíces del polinomio de grado n, Tn(x) son
Sea m<n. Expresamos
Empleando las relaciones de ortogonalidad, calculamos los valores de los coeficientes cj, de un modo similar a las series de Fourier.
Cuando i=0, T0(x)=1
Teniendo en cuenta que Tj(x)=cos(jθ) y que xk=cosθk, con θk=(2k-1)π/(2n)
Los coeficientes cj con j=1,...m, se calculan mediante la fórmula
El código de la función chebfit que calcula los coeficientes cj con j=0,...n-1, es el siguiente
function [c] = chebfit(n,f) y=cos((2*(1:n)-1)*pi/(2*n)); fk=f(y); c0=sum(fk)*2/n; ck=zeros(1,n-1); for j=1:n-1 coseno=cos(j*(2*(1:n)-1)*pi/(2*n)); ck(j)=sum(fk.*coseno)*2/n; end c=[c0 ck]; end
Como ejemplo, dibujamos la función f(x)=exp(x) en el intervalo de [-1,1] en color azul y aproximamos dicha función mediante polinomios de Chebyshev tomando n=3, en color rojo. Señalamos los puntos (nodos) las raíces del polinomio Tn(x).
f=@(x) exp(x); %función n=3; c=chebfit(n,f); hold on x=linspace(-1,1,100); yy=zeros(length(x),1); plot(x,f(x),'b') %dibuja la función en el intervalo [-1,1] %aproximación a la función mediante polinomios de Chebyshev yy=c(1)/2; for j=2:n yy=yy+c(j)*polyval(chebypoly(j-1),x); end plot(x,yy,'r') %puntos nodos x=cos((2*(1:n)-1)*pi/(2*n)); plot(x,f(x),'ro','markersize',4,'markerfacecolor','r') hold off grid on xlabel('x') ylabel('y') title('Aproximación Chebyshev') %utiliza Math Symbolic cheb=c(1)/2; for j=2:n cheb=cheb+c(j)*poly2sym(chebypoly(j-1)); end vpa(cheb,6)
En la ventana de comandos observamos el polinomio aproximador, cuyos coeficientes se proporcionan con 6 cifras decimales
ans =0.532042*x^2 + 1.12977*x + 1.0
Cambiamos en la primera línea del script la definición de la función anónima f para estudiar el siguiente ejemplo.
El resultado con n=7 es el siguiente
f=@(x) 1./(1+25*x.^2); %función n=7; %... el resto igual
El polinomio aproximador aparece en la ventana de comandos. Se han suprimido los coeficientes próximos a cero
ans =- 6.79168*x^6 + 12.1571*x^4 - 6.429*x^2 + 1.0
Dibujamos este polinomio con la función ezplot en el intervalo [-1,1] y obtendremos la curva en color rojo.
Aproximación de funciones definidas en el intervalo [a,b]
Si la función f(x) está definida en el intevalo [a,b], en vez del intervalo [-1,1] hacemos el cambio de variable
Cuando x varía de a a b, y varía de -1 a 1.
Modificamos la función chebfit y el script par incorporar esta posibilidad.
function [c] = chebfit_ab(a,b,n,f) y=cos((2*(1:n)-1)*pi/(2*n)); fk=f(y*(b-a)/2+(b+a)/2); c0=sum(fk)*2/n; ck=zeros(1,n-1); for j=1:n-1 coseno=cos(j*(2*(1:n)-1)*pi/(2*n)); ck(j)=sum(fk.*coseno)*2/n; end c=[c0 ck]; end
Dibujamos la función
en el intervalo de [-3,3] en color azul y aproximamos dicha función mediante polinomios de Chebyshev tomando n=5, en color rojo.
f=@(x) exp(-x.^2/2)/sqrt(2*pi); n=5; a=-3; b=3; c=chebfit_ab(a,b,n,f); hold on x=linspace(-1,1,100); xx=(a+b)/2+(b-a)*x/2; yy=zeros(length(x),1); plot(xx,f(xx),'b') yy=c(1)/2; for j=2:n yy=yy+c(j)*polyval(chebypoly(j-1),x); end plot(xx,yy,'r') %puntos nodos x=cos((2*(1:n)-1)*pi/(2*n)); xx=(a+b)/2+(b-a)*x/2; plot(xx,f(xx),'ro','markersize',4,'markerfacecolor','r') hold off grid on xlabel('x') ylabel('y') title('Aproximación Chebyshev') syms x; func=c(1)/2; for j=2:n func=func+c(j)*poly2sym(chebypoly(j-1)); end func=subs(func,x,(2*x-(a+b))/(b-a)); vpa(func,6)
El polinomio aproximador aparece en la ventana de comandos. Se han suprimido los coeficientes próximos a cero
ans =0.0105398*x^4 - 0.13397*x^2 + 0.398942
Dibujamos la función en el intervalo de [0,10] en color azul y aproximamos dicha función mediante polinomios de Chebyshev tomando n=6, en color rojo.
f=@(x) sqrt(x); n=6; a=0; b=10; c=chebfit_ab(a,b,n,f); ......
ans =0.225068*x - 0.136044*(0.2*x - 1.0)^2 + 0.0141154*(0.2*x - 1.0)^3 - 0.375426*(0.2*x - 1.0)^4 + 0.311013*(0.2*x - 1.0)^5 + 1.1024 >> vpa(expand(ans),6) ans =0.0000995242*x^5 - 0.00308879*x^4 + 0.0370076*x^3 - 0.221643*x^2 + 0.899309*x + 0.265797
La línea final es el polinomio aproximador.
Series de polinomios de Chebyshev
Los primeros polinomios de Chebyshev son
T0(x)=1
T1(x)=x
T2(x)=2x2-1
T3(x)=4x3-3x
T4(x)=8x4-8x2+1
T5(x)=16x5-20x3+5x
T6(x)= 32x6-48x4+18x2-1
....
Potencias de x
Escribimos las potencias de x, xn en términos de los polinomios de Chebyshev
Expresamos las potencias xn en términos de los polinomios de Chebyshev Tn(x),Tn-2(x), Tn-4(x)..., teniendo en cuenta la relación existente entre cosnθ y cos(n-2)θ, cos(n-4)θ... y la simetría de los números combinatorios
>> nchoosek(7,3) ans = 35 >> nchoosek(7,4) ans = 35
Supongamos que n=4, un número par
Supongamos que n=5, un número impar
Elaboramos una función, para obtener las potencias de x en términos de los polinomios de Chebyshev
function [T] = potencias_cheby(n) if rem(n,2)==0 T=zeros(1,n); for k=0:n/2-1 T(2*k+1)=nchoosek(n,k)/2^(n-1); end T0=nchoosek(n,n/2)/2^n; T=[T,T0]; else T=zeros(1,n+1); for k=0:floor(n/2) T(2*k+1)=nchoosek(n,k)/2^(n-1); end end end
Calculamos las potencias xn=c1Tn+c2Tn-1+...+cnT1+cn+1T0
for k=0:6 disp(['x^',num2str(k),'= ',num2str(potencias_cheby(k))]) end
x^0= 1 x^1= 1 0 x^2= 0.5 0 0.5 x^3= 0.25 0 0.75 0 x^4= 0.125 0 0.5 0 0.375 x^5= 0.0625 0 0.3125 0 0.625 0 x^6= 0.03125 0 0.1875 0 0.46875 0 0.3125
Desarrollo en serie en términos de los polinomios de Chebyshev
La fórmula del desarrollo en serie de Taylor es
El desarrollo en serie de log(x) cuando a=1 y tomando n=5 términos
>> syms x; >> taylor(log(x),x,'ExpansionPoint',1,'Order',5) ans =x - (x - 1)^2/2 + (x - 1)^3/3 - (x - 1)^4/4 - 1
En el caso particular de que a=0, se conoce como serie de Maclaurin. El desarrollo en serie de la función ex, tomando n=7 términos, es
>> syms x; >> taylor(exp(x),'order',7) ans =x^6/720 + x^5/120 + x^4/24 + x^3/6 + x^2/2 + x + 1
Sustituimos las potencias de x por combinaciones lineales de polinomios de Chebyshev
Elaboramos un script para aproximar la función f(x)=ex a un desarrollo en serie en términos de los polinomios de Chebyshev, calculando los coeficientes c1, c2...c6 y c7 que multiplican a los polinomios de Chebyshev, T6, T5...T1 y T0, respectivamente
syms x; n=7; %número de términos c_x=sym2poly(taylor(exp(x),'order',n)); cT=zeros(1,n); for k=1:n a=[zeros(1,k-1),potencias_cheby(n-k)]; cT=cT+c_x(k)*a; end
Partimos de la suma c1T6+c2T5+c3T4+c4T3+c5T2+c6T1+c7T0. Escribimos cada polinomio Ti(x) en función de la variable x y comprobamos que obtenemos el desarrollo en serie de Taylor original de n=7 términos
syms x; n=7; %número de términos c_x=sym2poly(taylor(exp(x),'order',n)); cT=zeros(1,n); for k=1:n a=[zeros(1,k-1),potencias_cheby(n-k)]; cT=cT+c_x(k)*a; end %comprobación Tk=0; for k=1:n Tk=Tk+cT(k)*poly2sym(chebypoly(n-k)); end
>> Tk Tk =x^6/720 + x^5/120 + x^4/24 + x^3/6 + x^2/2 + x + 1
Una vez que hemos comprobado que el script funciona correctamente con n=7 términos. El desarrollo en serie incluye infinitos términos, n=∞.
>> format long >> cT cT =0.000043402777778 0.000520833333333 0.005468750000000 0.044270833333333 0.271484375000000 1.130208333333333 1.266059027777778
Observamos que los primeros coeficientes c1=0.000043402777778 de T6, c2=0.000520833333333 de T5, c3=0.005468750000000 de T4 son próximos a cero, por tanto con n=10 términos podría ser suficiente.
syms x; n=10; %número de términos c_x=sym2poly(taylor(exp(x),'order',n)); cT=zeros(1,n); for k=1:n a=[zeros(1,k-1),potencias_cheby(n-k)]; cT=cT+c_x(k)*a; end
cT =0.000000010764578 0.000000193762401 0.000003197079613 0.000044952876984 0.000542922247024 0.005474175347222 0.044336841724537 0.271495225694444 1.130318196614583 1.266065809461806
El desarrollo en serie de la función ex en términos de los polinomios de Chebyshev es
ex∼1.266066·T0+1.130318·T1+0.271495·T2+0.044337·T3+0.005474·T4+0.000543·T5+...
Aproximación de funciones
Vamos a aproximar la función ex con la suma de los tres primeros términos del desarrollo en serie. 1.266066·T0+1.130318·T1+0.271495·T2. Sustituímos T0 por 1, T1 por x y T2 por 2x2-1. Comparamos esta aproximación con el desarrollo en serie de Taylor de tres términos. Elaboramos un script que nos permita aproximar la función ex a un número m de términos del desarrollo en serie cnT0+cn-1T1+...cn-m+1Tm-1 y la vamos a comparar con el mismo número de términos del desarrollo en serie de Taylor. Todas las representaciones gráficas se realizan en el intervalo [-1,1].
syms x; f=exp(x); n=10; %número de términos c_x=sym2poly(taylor(f,'order',n)); cT=zeros(1,n); for k=1:n a=[zeros(1,k-1),potencias_cheby(n-k)]; cT=cT+c_x(k)*a; end Tk=0; m=3; %tres términos c_8*T_2+c_9*T_1+c_10*T_0 for k=n-m+1:n Tk=Tk+cT(k)*poly2sym(chebypoly(n-k)); end hold on ezplot(f,[-1,1]) ezplot(Tk,[-1,1]) ezplot(taylor(f,'order',m),[-1,1]) hold off grid on legend('función','serie Chebyshev','Taylor','location','northwest') xlabel('x') ylabel('y') title('Desarrollo en serie')
>> vpa(Tk,6) ans =0.54299*x^2 + 1.13032*x + 0.994571 >> taylor(exp(x),'order',m) ans =x^2/2 + x + 1
Apreciamos que, la aproximación en términos de los polinomios de Chebyshev se ajusta mucho más a la función ex. Tomando m=4, no se aprecia apenas diferencia entre la función y dicho desarrollo en serie tal como se ve en la figura más abajo
>> vpa(Tk,6) ans =0.177347*x^3 + 0.54299*x^2 + 0.997308*x + 0.994571 >> taylor(exp(x),'order',m) ans =x^3/6 + x^2/2 + x + 1
Probamos ahora, la función f(x)=(1+x)-1, tomando m=6 términos del desarrollo en serie
syms x; f=(1+x)^-1; n=10; %número de términos c_x=sym2poly(taylor(f,'order',n)); cT=zeros(1,n); for k=1:n a=[zeros(1,k-1),potencias_cheby(n-k)]; cT=cT+c_x(k)*a; end Tk=0; m=6; for k=n-m+1:n Tk=Tk+cT(k)*poly2sym(chebypoly(n-k)); end hold on ezplot(f,[-1,1]) ezplot(Tk,[-1,1]) ezplot(taylor(f,'order',m),[-1,1]) hold off grid on legend('función','serie Chebyshev','Taylor') xlabel('x') ylabel('y') title('Desarrollo en serie')