Ajuste de datos a un combinación lineal de polinomios de Chebyshev

 

http://www.sc.ehu.es/sbweb/fisica3/especial/aprox_2/aprox_2.html

Ajuste de datos a un combinación lineal de polinomios de Chebyshev

En la página titulada Ajuste de datos se aplicó el procedimiento de mínimos cuadrados para aproximar un polinomio de grado n a un conjunto de pares de datos. En esta página, vamos a aproximar los pares de datos a una combinación lineal de polinomios de Chebyshev

En la página titulada Ajuste de datos. Regresión, aproximamos un polinomio de grado n, a un conjunto de m pares de datos (x_j, y_j) de modo que n<m. Sea el polinomio

P(x)=a₁xⁿ+a₂x^n-1+...a_nx+a_n+1

Se calcula la cantidad

$S=\sum _{j=1}^m{\left(P\left(x_j\right)-y_j\right)}^2=\sum _{j=1}^m{\left(a_1{x}_j^n+a_2{x}_j^{n-1}+a_n{x}_j+a_{n+1}-y_j\right)}^2$

Para obtener los valores de los coeficientes del polinomio aproximador se tienen que determinar los valores de los coeficientes a₁, a₂, a₃, ...a_n, a_n+1 de forma que la cantidad S tome un valor mínimo.

En la página titulada Ajuste de datos (no lineal), aproximamos una combinación lineal de las funciones: lnx, cosx y e^x, a un conjunto de m pares de datos

En esta página, el polinomio P_n(x) es una combinación lineal de n funciones base φ_i(x)

$P_n\left(x\right)=a_1{\phi }_1\left(x\right)+a_2{\phi }_2\left(x\right)+\mathrm{...}+a_n{\phi }_n\left(x\right)$

Aplicando el método de los mínimos cuadrados, queremos que la suma S sea mínima

$S={\sum _{j=1}^m\left(P_n\left(x_j\right)-y_j\right)}^2=\sum _{j=1}^m{\left(a_1{\phi }_1\left(x_j\right)+a_2{\phi }_2\left(x_j\right)+\mathrm{...}{a}_n{\phi }_n\left(x_j\right)-y_j\right)}^2$

La condición de mínimo es

$\frac{\partial S}{\partial a_i}=0i=\mathrm{1,2,...}n$

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

$\begin{array}{c}\frac{1}{2}\frac{\partial S}{\partial a_1}=0\Rightarrow \sum _{j=1}^m{\phi }_1\left(x_j\right)\left(a_1{\phi }_1\left(x_j\right)+a_2{\phi }_2\left(x_j\right)+\mathrm{...}+a_n{\phi }_n\left(x_j\right)-y_j\right)=0\\ \frac{1}{2}\frac{\partial S}{\partial a_2}=0\Rightarrow \sum _{j=1}^m{\phi }_2\left(x_j\right)\left(a_1{\phi }_1\left(x_j\right)+a_2{\phi }_2\left(x_j\right)+\mathrm{...}+a_n{\phi }_n\left(x_j\right)-y_j\right)=0\\ \mathrm{...}\\ \frac{1}{2}\frac{\partial S}{\partial a_n}=0\Rightarrow \sum _{j=1}^m{\phi }_n\left(x_j\right)\left(a_1{\phi }_1\left(x_j\right)+a_2{\phi }_2\left(x_j\right)+\mathrm{...}+a_n{\phi }_n\left(x_j\right)-y_j\right)=0\end{array}$

o bien,

$\begin{array}{c}{a}_1\left(\sum _{j=1}^m{\phi }_1^2\left(x_j\right)\right)+a_2\left(\sum _{j=1}^m{\phi }_1\left(x_j\right){\phi }_2\left(x_j\right)\right)+\mathrm{...}+a_n\left(\sum _{j=1}^m{\phi }_1\left(x_j\right){\phi }_n\left(x_j\right)\right)=\left(\sum _{j=1}^m{\phi }_1\left(x_j\right)y_j\right)\\ a_1\left(\sum _{j=1}^m{\phi }_2\left(x_j\right){\phi }_1\left(x_j\right)\right)+a_2\left(\sum _{j=1}^m{\phi }_2^2\left(x_j\right)\right)+\mathrm{...}+a_n\left(\sum _{j=1}^m{\phi }_2\left(x_j\right){\phi }_n\left(x_j\right)\right)=\left(\sum _{j=1}^m{\phi }_2\left(x_j\right)y_j\right)\\ \mathrm{...}\\ a_1\left(\sum _{j=1}^m{\phi }_n\left(x_j\right){\phi }_1\left(x_j\right)\right)+a_2\left(\sum _{j=1}^m{\phi }_n\left(x_j\right){\phi }_2\left(x_j\right)\right)+\mathrm{...}+a_n\left(\sum _{j=1}^m{\phi }_n^2\left(x_j\right)\right)=\left(\sum _{j=1}^m{\phi }_n\left(x_j\right)y_j\right)\end{array}$

En forma matricial, A·a=b. Donde A es una matriz simétrica por lo que solamente tendremos que calcular los elementos pertenecientes al triángulo superior, a es el vector de las incógnitas y b el de los términos independientes.

Ajuste de los datos a una combinación lineal de polinomios de Chebyshev

Polinomios de Chebyshev en el intervalo [0,1]

Queremos ajustar estos pares de datos a un polinomio P₄(x) combinación lineal de los cuatro primeros polinomios de Chebyshev

T₀(x)=1
T₁(x)=x
T₂(x)=2x²-1
T₃(x)=4x³-3x

Si las abscisas de los datos están comprendidas en el intevalo [a,b], en vez del intervalo [-1,1] hacemos el cambio de variable

$y=\frac{x-\frac{1}{2}(b+a)}{\frac{1}{2}(b-a)}$

Cuando x varía de a a b, y varía de -1 a 1.

Las abscisas de los datos están en el intervalo [0,1]. Modificamos los polinomios sustituyendo x por 2x-1. Los primeros polinomios de Chebyshev en el intervalo [0,1] son

φ₁(x)=1
φ₂(x)=2x-1
φ₃(x)=2(2x-1)²-1=8x²-8x+1
φ₄(x)=4(2x-1)³-3(2x-1)=32x³-48x²+18x-1

>> syms x;
>> f=subs(2*x^2-1,x,2*x-1);
>> collect(f)
ans =8*x^2 - 8*x + 1
>> f=subs(4*x^3-3*x,x,2*x-1);
>> collect(f)
ans =32*x^3 - 48*x^2 + 18*x - 1

Ajuste de los datos

Supongamos que hemos obtenido de alguna experiencia los siguientes pares de datos:

x	0	0.05	0.10	0.15	0.20	0.25	0.30	0.35	0.40	0.45	0.50
y	0.486	0.866	0.944	1.144	1.103	1.202	1.166	1.191	1.124	1.095	1.122
x	0.55	0.60	0.65	0.70	0.75	0.80	0.85	0.90	0.95	1.00
y	1.102	1.099	1.017	1.111	1.117	1.152	1.265	1.380	1.575	1.875

Primer procedimiento

Calculamos los elementos del triángulo superior de la matriz simétrica A

x=0:0.05:1;
y=[0.486,0.866,0.944,1.144,1.103,1.202,1.166,1.191,1.124,1.095,1.122,
1.102,1.099,1.017,1.111,1.117,1.152,1.265,1.380,1.575,1.857];
phi={@(x)1,@(x)2*x-1,@(x)8*x.^2-8*x+1,@(x)32*x.^3-48*x.^2+18*x-1};
A=zeros(length(phi));
%calcula la parte superior de la matriz simétrica A
for i=1:length(phi)
    for j=i:length(phi)
        A(i,j)=sum(phi{i}(x).*phi{j}(x));
    end
end
A(1,1)=length(x);

En la tercera línea hemos definido el vector phi formado por las cuatro funciones base, [φ₁(x), φ₂(x), φ₃(x), φ₄(x)] primeros polinomios de Chebyshev en el intervalo [0,1]

El cálculo del primer elemento de la matriz no se efectúa del mismo modo que el resto de los elementos.

$A_{11}=\sum _{j=1}^m{\phi }_1^2\left(x_j\right)$

Como se ha explicado en la página Sentencias iterativas, si la función f es constante el resultado de pasarle un vector x no es un vector, sino el valor constante. El primer elemento del vector phi es φ₁(x) que vale 1

Una vez que hemos calculado los elementos de la diagonal y los que están por encima de éstos y los hemos guardado en la matriz A. Creamos la matriz simétrica A mediante el siguiente línea de código

A=A+triu(A,1)';

Calculamos el vector b de los términos independientes y despejamos el vector a de las incógnitas utilizando el operador división por la izquierda \

for i=1:length(phi)
    b(i)=sum(y.*phi{i}(x));
end
a=A\b;

Definimos la función f(x) combinación lineal de las cuatro funciones base, f(x)=a₁φ₁(x)+a₂φ₂(x)+a₃φ₃(x)+a₄φ₄(x)

phi={@(x)1,@(x)2*x-1,@(x)8*x.^2-8*x+1,@(x)32*x.^3-48*x.^2+18*x-1};
f=@(x) a(1)*phi{1}(x);
for i=2:length(phi)
    f=@(x)(f(x)+a(i)*phi{i}(x));
end

Representamos los datos experimentales y la función que mejor ajusta f(x) en el intervalo [0,1]

hold on
plot(x,y,'bo', 'markersize',4,'markeredgecolor','b','markerfacecolor','b');
f=@(x) a(1)*phi{1}(x)+a(2)*phi{2}(x)+a(3)*phi{3}(x)+a(4)*phi{4}(x);
fplot(f,[x(1),x(end)],'r')
hold off
grid on
xlabel('x')
ylabel('y')
title('Regresión Chebyshev')

Juntamos todas las porciones de código en el mismo script

x=0:0.05:1;
y=[0.486,0.866,0.944,1.144,1.103,1.202,1.166,1.191,1.124,1.095,
1.122,1.102,1.099,1.017,1.111,1.117,1.152,1.265,1.380,1.575,1.857];
phi={@(x)1,@(x)2*x-1,@(x)8*x.^2-8*x+1,@(x)32*x.^3-48*x.^2+18*x-1};
A=zeros(length(phi));
b=zeros(length(phi),1);
%calcula la parte superior de la matriz simétrica A
for i=1:length(phi)
    for j=i:length(phi)
        A(i,j)=sum(phi{i}(x).*phi{j}(x));
    end
end
A(1,1)=length(x);
%completa la matriz simétrica
A=A+triu(A,1)';
for i=1:length(phi)
    b(i)=sum(y.*phi{i}(x));
end
a=A\b;
f=@(x) a(1)*phi{1}(x);
for i=2:length(phi)
    f=@(x)(f(x)+a(i)*phi{i}(x));
end
hold on
plot(x,y,'bo', 'markersize',4,'markeredgecolor','b','markerfacecolor','b');
f=@(x) a(1)*phi{1}(x)+a(2)*phi{2}(x)+a(3)*phi{3}(x)+a(4)*phi{4}(x);
fplot(f,[x(1),x(end)],'r')
hold off
grid on
xlabel('x')
ylabel('y')
title('Regresión Chebyshev')

>> a
a =1.160969479033553
   0.393514467988152
   0.046849832090107
   0.239646175715970

Segundo procedimiento

Con las funciones a ajustar φ₁(x), φ₂(x), φ₃(x) y φ₄(x) y el vector x de las abscisas de los pares de datos, creamos la matriz V. Despejamos el vector columna a de las incógnitas utilizando el operador división por la izquierda \ en la relación (V^T·V)a=V^T·Y. Siendo Y el vector columna de las ordenadas

$V=\left(\begin{array}{cccc}{\phi }_1\left(x_1\right)& {\phi }_2\left(x_1\right)& {\phi }_3\left(x_1\right)& {\phi }_4\left(x_1\right)\\ {\phi }_1\left(x_2\right)& {\phi }_2\left(x_2\right)& {\phi }_3\left(x_2\right)& {\phi }_4\left(x_2\right)\\ \mathrm{...}& \mathrm{...}& \mathrm{...}& \mathrm{...}\\ {\phi }_1\left(x_m\right)& {\phi }_2\left(x_m\right)& {\phi }_3\left(x_m\right)& {\phi }_4\left(x_m\right)\end{array}\right)$

x=0:0.05:1;
y=[0.486,0.866,0.944,1.144,1.103,1.202,1.166,1.191,1.124,1.095,1.122,1.102,
1.099,1.017,1.111,1.117,1.152,1.265,1.380,1.575,1.857];
phi={@(x)1,@(x)2*x-1,@(x)8*x.^2-8*x+1,@(x)32*x.^3-48*x.^2+18*x-1};
V=zeros(length(x),length(phi));
for i=1:length(phi)
    V(:,i)=phi{i}(x);
end
a=(V'*V)\(V'*y');
f=@(x) a(1)*phi{1}(x);
for i=2:length(phi)
    f=@(x)(f(x)+a(i)*phi{i}(x));
end
hold on
plot(x,y,'bo','markersize',4,'markeredgecolor','b','markerfacecolor','b');
f=@(x) a(1)*phi{1}(x)+a(2)*phi{2}(x)+a(3)*phi{3}(x)+a(4)*phi{4}(x);
fplot(f,[x(1),x(end)],'r')
hold off
grid on
xlabel('x')
ylabel('y')
title('Regresión Chebyshev')

Obtenemos el mismo resultado ajustando los datos experimentales a un polinomio de tercer grado P(x)=a₁x³+a₂x²+a₃x+a₄.

Como se afirma en el documento que se cita en las referencias, la ventaja del primer procedimiento (combinación lineal de los polinomios de Chebyshev) es la estabilidad de las soluciones frente a una perturbación en los valores del vector b de los términos independientes

Referencias

Los datos de este ejemplo, han sido tomados del documento titulado LEAST SQUARES DATA FITTING Table 2,