Ajuste de datos (no lineal)

En la página anterior, se han ajustado datos a funciones no lineales (función potencial y función exponencial) pero que se pueden transformar en funciones lineales. En esta página, vamos a estudiar algunos ejemplos de ajuste a funciones no lineales

Función y=a/x+bx

La función

y= a x +bx

aparece en el estudio del péndulo compuesto, con y=P2/(4·π2) y x es la distancia entre el centro de masas c.m. y el centro de oscilación O. P es el periodo de oscilación del péndulo compuesto cuando su amplitud es pequeña.

Dada una tabla de valores xj e yj se trata de calcular los valores de los coeficientes a y de b que mejor ajustan a los datos experimentales. El procedimiento aplicado es similar a la regresión lineal

Si (xj, yj) son las coordenadas de un dato experimental, a la abscisa xj le correspondería la ordenada y=a/xj+bxj. La diferencia es

d j = y j a x j b x j

Calcularemos los valores de los parámetros a y b que hacen que la suma

S(a,b)= 1 m d j 2 = 1 m ( y j a x j b x j ) 2

sea mínima.

S a =0 1 m 2( y j a x j b x j ) ( 1 x j )=0 ( 1 m 1 x j 2 )a+mb= 1 m y j x j S b =0 1 m 2( y j a x j b x j ) ( x j )=0 ma+( 1 m x j 2 )b= 1 m x j y j

Resolvemos el sistema de dos ecuaciones con dos incógnitas, para determinar los coeficientes a y b.

b= m( 1 m y j x j )( 1 m 1 x j 2 )( 1 m x j y j ) m 2 ( 1 m x j 2 )( 0 N1 1 x j 2 ) a= ( 1 m x j y j )( 1 m x j 2 )b m

Medimos el periodo Pi de péndulo para cada posición xi, completando una tabla con N pares de datos

x (cm) P (s)
5 2.620
10 1.936
15 1.668
20 1.568
25 1.520
30 1.512
35 1.536
40 1.576
45 1.600

Creamos un script con MATLAB para calcular los valores de a y b. Representamos los datos experimentales y la función que mejor ajusta

x=(5:5:45)/100;  
y=[2.620, 1.936, 1.668, 1.568, 1.520, 1.512, 1.536, 1.576, 1.600].^2/(4*pi^2);
m=length(x);
b=(m*sum(y./x)-sum(1./x.^2)*sum(x.*y))/(m^2-sum(x.^2)*sum(1./x.^2));
a=(sum(x.*y)-b*sum(x.^2))/m;
hold on
%representa los datos experimentales
plot(x,y,'bo','markersize',4,'markerfacecolor','b')
% función
f_ajuste=@(x) a./x+b*x;          
fplot(f_ajuste,[0.04,0.45],'r')
title('función, y=a/x+bx')
xlabel('x')
ylabel('y')
grid on
hold off 

Función, y=a+b/x2

En una experiencia de laboratorio, colocamos un sensor, a una distancia x de la fuente de luz, que vamos cambiando. Obtenemos la siguiente tabla:

Distancia x a la fuente de luz (cm) Intensidad I(lux)
0.20 171
0.25 106
0.30 73
0.35 52
0.40 39
0.45 30.5
0.50 24.5

Otros ejemplos de este ajuste de datos es la dependencia del índice de refracción de un vidrio con la longitud de onda, o del agua de una gota.

Queremos ajustar m pares de datos (xi,yi) a la función

y=a+ b x 2

de modo que la suma S sea mínima

S= j=1 m ( a+ b x j 2 y j ) 2 1 2 S a =am+b j=1 m 1 x j 2 j=1 m y j =0 1 2 S b =a j=1 m 1 x j 2 +b j=1 m 1 x j 4 j=1 m y j x j 2 =0 ( m j=1 m 1 x j 2 j=1 m 1 x j 2 j=1 m 1 x j 4 )( a b )=( j=1 m y j j=1 m y j x j 2 )

Resolvemos el sistema de dos ecuaciones con dos incógnitas. Creamos la matriz A de los coeficientes y el vector B de los términos independientes, despejamos el vector X de las incógnitas utilizando el operador división por la izquierda \. El elemento X(1) es a y el elemento X(2) del vector X es b.

x=0.2:0.05:0.5;
y=[171,106,73,52,39,30.5,24.5];
A=[length(x),sum(1./x.^2); sum(1./x.^2),sum(1./x.^4)];
B=[sum(y);sum(y./x.^2)];
X=A\B;  % X(1) es a y X(2) es b
hold on
plot(x,y,'bo','markersize',4,'markerfacecolor','b')
f_ajuste=@(x) X(1)+X(2)./x.^2;
fplot(f_ajuste,[x(1),x(end)],'r')
axis([0.19,0.51, 0,180])
grid on
xlabel('x')
ylabel('y')
title('función, y=a+b/x^2')

Función, y=a+b/x+c/x2

Queremos ajustar los siguientes pares de datos

x 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
y 6.4 15.2 15.6 15.3 15.0 14.7 14.5 14.3 14.2 14.0

a la función

y=a+ b x + c x 2

de modo que la suma S sea mínima

S= j=1 m ( a+ b x j + c x j 2 y j ) 2 1 2 S a =am+b j=1 m 1 x j +c j=1 m 1 x j 2 j=1 m y j =0 1 2 S b =a j=1 m 1 x j +b j=1 m 1 x j 2 +c j=1 m 1 x j 3 j=1 m y j x j =0 1 2 S c =a j=1 m 1 x j 2 +b j=1 m 1 x j 3 +c j=1 m 1 x j 4 j=1 m y j x j 2 =0 ( m j=1 m 1 x j j=1 m 1 x j 2 j=1 m 1 x j j=1 m 1 x j 2 j=1 m 1 x j 3 j=1 m 1 x j 2 j=1 m 1 x j 3 j=1 m 1 x j 4 )( a b c )=( j=1 m y j j=1 m y j x j j=1 m y j x j 2 )

x=0.5:0.5:5;
y=[6.4,15.2,15.6,15.3,15.0,14.7, 14.5,14.3, 14.2,14.0];
A=[length(x),sum(1./x), sum(1./x.^2); sum(1./x),sum(1./x.^2), sum(1./x.^3);
sum(1./x.^2),sum(1./x.^3), sum(1./x.^4)];
B=[sum(y);sum(y./x);sum(y./x.^2)];
X=A\B;  % X(1) es a, X(2) es b y X(3) es c
hold on
plot(x,y,'ro','markersize',4,'markerfacecolor','r')
f_ajuste=@(x) X(1)+X(2)./x+X(3)./x.^2;
fplot(f_ajuste,[x(1),x(end)],'b')
axis([0,5, 6,16])

grid on
xlabel('x')
ylabel('y')
title('función, y=a+b/x+c/x^2')

Función y=a·lnx+b·cosx+c·ex

Supongamos que queremos ajustar los datos (xj,yj), j=1,2,...m

x 1 1.5 2 2.5 3 3.5 4
y -7.3 -13.9 -24.2 -42.1 -64.7 -103.5 -180.3

a la función y=a·lnx+b·cosx+c·ex de modo que la suma S sea mínima

S= j=1 m ( aln x j +bcos x j +c e x j y j ) 2 1 2 S a =0a j=1 m ( ln x j ) 2 +b j=1 m ln x j cos x j +c j=1 m ln x j e x j = j=1 m y j ln x j 1 2 S b =0a j=1 m cos x j ln x j +b j=1 m ( cos x j ) 2 +c j=1 m cos x j e x j = j=1 m y j cos x j 1 2 S c =0a j=1 m e x j ln x j +b j=1 m e x j cos x j +c j=1 m ( e x j ) 2 = j=1 m y j e x j

Escribimos el sistema de tres ecuaciones con tres incógnitas de forma matricial

( ln x 1 ln x 2 ... ln x m cos x 1 cos x 2 ... cos x m e x 1 e x 2 ... e x m )( ln x 1 cos x 1 e x 1 ln x 2 cos x 2 e x 2 ... ... ... ln x m cos x m e x m )( a b c )= ( ln x 1 ln x 2 ... ln x m cos x 1 cos x 2 ... cos x m e x 1 e x 2 ... e x m )( y 1 y 2 ... y m )

Con las funciones a ajustar ln(x), cos(x) y exp(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 (a,b,c) 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

x=1:0.5:4;
y=[-7.3,-13.9,-24.2,-42.1,-64.7,-103.5,-180.3];
V=zeros(length(x),3);
V(:,1)=log(x);
V(:,2)=cos(x);
V(:,3)=exp(x);
a=(V'*V)\(V'*y');
hold on
plot(x,y,'ro','markersize',4,'markerfacecolor','r')
z=@(x) a(1)*log(x)+a(2)*cos(x)+a(3)*exp(x);
fplot(z,[x(1),x(end)])
xlabel('x')
ylabel('y')
grid on
title('Ajuste no lineal')
hold off

Ilustramos este procedimiento con un ejemplo más. Ajustamos los datos (xj, yj) j=1,2...m de la tabla a la función y=a+b·cos(2πt)+c·sin(2πt)+d·cos(4πt)

x 0 0.125 0.25 0.375 0.5 0.625 0.75 0.875
y 2.2 2.8 6.1 3.9 0.0 -1.1 0.6 1.1

Creamos la matriz V

V=( 1 cos( 2π x 1 ) sin( 2π x 1 ) cos( 4π x 1 ) 1 cos( 2π x 2 ) sin( 2π x 2 ) cos( 4π x 2 ) ... ... ... ... 1 cos( 2π x m ) sin( 2π x m ) cos( 4π x m ) )

Despejamos el vector a columna de las incógnitas utilizando el operador \

x=0:0.125:0.875;
y=[2.2,2.8,6.1,3.9,0,-1.1,0.6,1.1];
V=zeros(length(x),4);
V(:,1)=1;
V(:,2)=cos(2*pi*x);
V(:,3)=sin(2*pi*x);
V(:,4)=cos(4*pi*x);
a=(V'*V)\(V'*y');
hold on
plot(x,y,'ro','markersize',4,'markerfacecolor','r')
z=@(x) a(1)+a(2)*cos(2*pi*x)+a(3)*sin(2*pi*x)+a(4)*cos(4*pi*x);
fplot(z,[x(1),x(end)])
xlabel('x')
ylabel('y')
grid on
title('Ajuste no lineal')
hold off
>> (V'*V)
ans =
    8.0000   -0.0000   -0.0000   -0.0000
   -0.0000    4.0000    0.0000   -0.0000
   -0.0000    0.0000    4.0000    0.0000
   -0.0000   -0.0000    0.0000    4.0000
>> a
a =
    1.9500
    0.7445
    2.5594
   -1.1250

En la siguiente página, estudiaremos cómo se ajustan los datos a funciones armónicas, cuando las abscisas xj están igualmente espaciadas, lo que nos permitirá prescindir de la matriz V calculando las incógnitas a,b,c,d de una forma más sencilla ya que como podemos observar la matriz V'*V es diagonal

Ajuste mediante nlinfit

Ejemplo 1

En un experimento hemos obtenido los siguientes datos:

x 0.008 0.009 0.015 0.027 0.0427 0.068 0.139 0.197
y 7 8.5 11 14.5 16.5 19 21 21.5

calculamos los parámetros a y b a partir del ajuste de datos por el procedimiento de mínimos cuadrados implementado en la función MATLAB nlinfit, tomando como modelo la función

y= ax b+x

x=[0.008,0.009,0.015,0.027,0.0427,0.068,0.139,0.197];
y=[7,8.5,11,14.5,16.5,19,21,21.5];
hold on
plot(x,y,'ro','markersize',4,'markerfacecolor','r')
f=@(a,x) a(1)*x./(a(2)+x);      
a0=[23,0.02];  %valor inicial
af=nlinfit(x,y,f,a0)
g=@(x) f(af,x);
fplot(g,[x(1),x(end)])
title('Ajuste no lineal, y=ax/(b+x)')
xlabel('x')
ylabel('y')
grid on
hold off 

Corremos el script en la ventana de comandos

af =   23.5104    0.0171

Ejemplo 2

En un experimento hemos obtenido los siguientes datos:

x 1 2 3 4 5 6 7 8 9 10 11 12
y 0.07 0.21 0.35 0.485 0.605 0.71 0.79 0.847 0.89 0.924 0.946 0.963

Calculamos los parámetros a y b a partir del ajuste de datos por el procedimiento de mínimos cuadrados implementado en la función MATLAB nlinfit, tomando como modelo la función

y=1 aexp(x/a)bexp(x/b) ab

x=1:12;
y=[0.07 0.21 0.35 0.485 0.605 0.71 0.79 0.847 0.89 0.924 0.946 0.963];
hold on
plot(x,y,'ro','markersize',4,'markerfacecolor','r')
f=@(a,x) 1-(a(1)*exp(-x/a(1))-a(2)*exp(-x/a(2)))/(a(1)-a(2));      
a0=[2 1];  %valor inicial
af=nlinfit(x,y,f,a0)
g=@(x) f(af,x);
fplot(g,[x(1),x(end)])
title('Ajuste no lineal')
xlabel('x')
ylabel('y')
grid on
hold off 

Corremos el script en la ventana de comandos

af =    2.4120    2.4119

Ejemplos en el curso de Física. Ajuste no lineal: nlinfit

Péndulo compuesto.

El prisma de vidrio

El arco iris

Carga de un condensador

Fuerza magnética sobre conductor rectilíneo. La balanza de corriente

Datos de las olas del mar

Análisis de las alturas y periodos de las olas

Función de distribución de Weibull

Función de distribución de Rayleigh

Energía producida por un aerogenerador