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)

f(x) c 0 2 + j=1 m c j T j (x)

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

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

Sea m<n. Expresamos

f( x k )= c 0 2 + j=1 m c j T j ( x k )

Empleando las relaciones de ortogonalidad, calculamos los valores de los coeficientes cj, de un modo similar a las series de Fourier.

k=1 n f( x k ) T i ( x k ) = c 0 2 k=1 n T i ( x k ) + j=1 m c j ( k=1 n T j ( x k ) T i ( x k ) ) k=1 n f( x k ) T i ( x k ) = c i n 2 c i = 2 n k=1 n f( x k ) T i ( x k )

Cuando i=0, T0(x)=1

k=1 n f( x k )= c 0 2 n c 0 = 2 n k=1 n f( x k )

Teniendo en cuenta que Tj(x)=cos() y que xk=cosθk, con θk=(2k-1)π/(2n)

T j ( x k )=cos( j(2k1)π 2n )

Los coeficientes cj con j=1,...m, se calculan mediante la fórmula

c j = 2 n k=1 n f( cos( (2k1)π 2n ) )cos( j(2k1)π 2n )

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.

y= 1 1+25 x 2

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

y = x 1 2 ( b + a ) 1 2 ( b a )

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

y= exp( x 2 /2 ) 2π

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 y= x 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

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

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

( n m )=( n nm )

>> nchoosek(7,3)
ans =    35
>> nchoosek(7,4)
ans =    35

Supongamos que n=4, un número par

( e iθ + e iθ ) 4 = e i4θ +( 4 1 ) e i3θ e iθ +( 4 2 ) e i2θ e i2θ +( 4 3 ) e iθ e i3θ + e i4θ ( e iθ + e iθ ) 4 =( e i4θ e i4θ )+( 4 1 )( e i2θ e i2θ )+( 4 2 ) 2 4 ( cosθ ) 4 =2cos( 4θ )+2( 4 1 )cos( 2θ )+( 4 2 ) ( cosθ ) 4 = 1 2 3 ( cos( 4θ )+( 4 1 )cos( 2θ )+ 1 2 ( 4 2 ) ) x 4 = 1 2 3 ( T 4 +4 T 2 +3 T 0 )

Supongamos que n=5, un número impar

( e iθ + e iθ ) 5 = e i5θ +( 5 1 ) e i4θ e iθ +( 5 2 ) e i3θ e i2θ +( 5 3 ) e i2θ e i3θ + ( 5 4 ) e iθ e i4θ + e i5θ ( e iθ + e iθ ) 5 =( e i5θ e i5θ )+( 5 1 )( e i3θ e i3θ )+( 5 2 )( e iθ e iθ ) 2 5 ( cosθ ) 5 =2cos( 5θ )+2( 5 1 )cos( 3θ )+( 5 2 )cos( θ ) ( cosθ ) 5 = 1 2 4 ( cos( 5θ )+( 5 1 )cos( 3θ )+( 5 2 )cos( θ ) ) x 5 = 1 16 ( T 5 +5 T 3 +10 T 1 )

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

f(x)=f(a)+ f ' (a)(xa)+ f '' (a) 2! (xa) 2 +...+ f (n) (a) n! (xa) n +...

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

e x 1+x+ x 2 2 + x 3 6 + x 4 24 + x 5 120 + x 6 720 +...

>> 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

e x T 0 + T 1 + 1 2 · 1 2 ( T 0 + T 2 )+ 1 6 · 1 4 ( 3 T 1 + T 3 )+ 1 24 · 1 8 ( 3 T 0 +4 T 2 + T 4 )+ 1 120 · 1 16 ( 10 T 1 +5 T 3 + T 5 )+ 1 720 · 1 32 ( 10 T 0 +15 T 2 +6 T 4 + T 6 )= c 1 T 6 + c 2 T 5 + c 3 T 4 + c 4 T 3 + c 5 T 2 + c 6 T 1 + c 7 T 0

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')