viernes, 20 de noviembre de 2020

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. Regresión, aproximamos un polinomio de grado n, a un conjunto de m pares de datos (xj, yj) de modo que n<m. Sea el polinomio

P(x)=a1xn+a2xn-1+...anx+an+1

Se calcula la cantidad

S=j=1m(P(xj)yj)2=j=1m(a1xjn+a2xjn1+anxj+an+1yj)2

Para obtener los valores de los coeficientes del polinomio aproximador se tienen que determinar los valores de los coeficientes a1, a2, a3, ...an, an+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 ex, a un conjunto de m pares de datos

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

Pn(x)=a1φ1(x)+a2φ2(x)+...+anφn(x)

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

S=j=1m(Pn(xj)yj)2=j=1m(a1φ1(xj)+a2φ2(xj)+...anφn(xj)yj)2

La condición de mínimo es

Sai=0i=1,2,...n

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

12Sa1=0j=1mφ1(xj)(a1φ1(xj)+a2φ2(xj)+...+anφn(xj)yj)=012Sa2=0j=1mφ2(xj)(a1φ1(xj)+a2φ2(xj)+...+anφn(xj)yj)=0...12San=0j=1mφn(xj)(a1φ1(xj)+a2φ2(xj)+...+anφn(xj)yj)=0

o bien,

a1(j=1mφ12(xj))+a2(j=1mφ1(xj)φ2(xj))+...+an(j=1mφ1(xj)φn(xj))=(j=1mφ1(xj)yj)a1(j=1mφ2(xj)φ1(xj))+a2(j=1mφ22(xj))+...+an(j=1mφ2(xj)φn(xj))=(j=1mφ2(xj)yj)...a1(j=1mφn(xj)φ1(xj))+a2(j=1mφn(xj)φ2(xj))+...+an(j=1mφn2(xj))=(j=1mφn(xj)yj)

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 P4(x) combinación lineal de los cuatro primeros polinomios de Chebyshev

T0(x)=1
T1(x)=x
T2(x)=2x2-1
T3(x)=4x3-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=x12(b+a)12(ba)

Cuando x varía de a a by 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

φ1(x)=1
φ2(x)=2x-1
φ3(x)=2(2x-1)2-1=8x2-8x+1
φ4(x)=4(2x-1)3-3(2x-1)=32x3-48x2+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:

x00.050.100.150.200.250.300.350.400.450.50
y0.4860.8660.9441.1441.1031.2021.1661.1911.1241.0951.122
x0.550.600.650.700.750.800.850.900.951.00
y1.1021.0991.0171.1111.1171.1521.2651.3801.5751.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, [φ1(x), φ2(x), φ3(x), φ4(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.

A11=j=1mφ12(xj)

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 φ1(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(xcombinación lineal de las cuatro funciones base, f(x)=a1φ1(x)+a2φ2(x)+a3φ3(x)+a4φ4(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 φ1(x), φ2(x), φ3(x) y φ4(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 (VT·V)a=VT·Y. Siendo Y el vector columna de las ordenadas

V=(φ1(x1)φ2(x1)φ3(x1)φ4(x1)φ1(x2)φ2(x2)φ3(x2)φ4(x2)............φ1(xm)φ2(xm)φ3(xm)φ4(xm))

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)=a1x3+a2x2+a3x+a4.

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,



No hay comentarios:

Publicar un comentario