Aproximación por minimos cuadrados grupo 23

Se está testeando un nuevo metamaterial sintético en el laboratorio frente a la fatiga. Para ello se golpea periódicamente el material con un dispositivo automático, manteniendo la fuerza pero cambiando la frecuencia. Se registran los datos de tiempo de rotura tr en horas para cada frecuencia de golpeo ω (número de golpes por minuto). Los datos obtenidos son:

Frecuencia/tiempo de rotura

1.33	1040
2.67	1036
4	999
5.33	1017
6.67	1049
8.67	1000
10.67	1007
12	NaN
14.67	NaN
16	NaN
17.3	NaN
18.67	943
20	916
21.33	944
22	929
22.67	955
24	940
24,67	911
25.33	785
25.67	703
26	634
26.33	326
26.67	7

Algunos ensayos no tienen datos por un error en el laboratorio.

1 Dibujar los puntos y generar la gráfica que los representa

Comenzamos por introducir en Matlab los datos y a continuación generar la gráfica que los representa para su fácil visualización:

w=[1.33 2.67 4 5.33 6.67 8.67 10.67 18.67 20 21.33 22 22.67 23.33 24 24.67 25.33 25.67 26 26.33 26.67]  %vector de la frecuencia
tr=[1040 1036 999 1017 1049 1000 1007 943 916 944 929 955 931 940 911 785 703 624 326 7]  %vector del tiempo de rotura
plot(w.',tr.','k--o')
xlabel('frecuencia')
ylabel('tiempo de rotura')
title('ensayo fatiga')
grid  %dibujamos la gráfica tiempo de rotura/frecuencia y le damos formato

2 Cálculo y representación de una recta por mínimos cuadrados

2.1 Cálculo de la recta

Procedemos a ajustar los datos a una recta, para ello introducimos en Matlab:

% añadir al programa anterior
m1r=[ones(20,1),w.'] %matriz de los coeficientes, tamaño 20x2, siendo la primera fila una columna de unos y la segunda el vector de las frecuencias traspuesto
m2r=[m1r,tr.']  %matriz ampliada de la matriz anterior junto con el vector de los tiempos de rotura traspuesto
rank(m1r)
rank(m2r) %usamos la función que nos calcula el rango de las matrices

Si ejecutásemos el programa en este momento obtenemos que el rango de m1 es 2 y el de m3 es 3, por lo que sabemos que el sistema es incompatible, pasamos a calcular la resolución del sistema, que Matlab calcula automáticamente por mínimos cuadrados y a continuación definimos la recta y la añadimos en el gráfico creado anteriormente escribiendo:

% añadir al programa anterior
xr=m1r\tr.'  %definimos xr como el vector con las soluciones del sistema de la recta
yr=xr(1)+xr(2)*w 
hold on  %le pedimos que mantenga la gráfica anterior
plot(w,yr,'b-')  %añade la recta calculada por mínimos cuadrados
legend('datos','recta de minimos') %añade una leyenda

2.2 Cálculo del error cuadrático medio de la recta

Queremos tener una referencia del error cometido por la aproximación, por lo que calculamos el error cuadrático medio de la siguiente forma

% añadir al programa anterior
errorr=(sum((tr-yr).^2))/20

2.3 Aproximación del valor 16 golpes/min

Para saber qué valor predice la recta para una frecuencia de 16 golpes por minuto ejecutamos:

% añadir al programa anterior
zr=xr(1)+xr(2)*16

2.4 Resultado

una vez escrito el código lo ejecutamos y Matlab nos da los siguientes resultados:

xr=[1.146876510867884e+03;-16.980810431368106], por lo que la ecuación de la recta es: tr=1146.9-17*ω

errorr=4.310542091992706e+04, el error cuadrático medio es de 43105

zr=8.751835439659939e+02, es decir, para una frecuencia de 16 golpes por minuto se estima un tiempo de rotura de 875 horas 11 minutos y además representa la siguiente gráfica:

gráfica de los datos de partida junto con la ecuación lineal de mínimos cuadrados

3 Cálculo y representación de una parábola por mínimos cuadrados

3.1 Cálculo de la parábola

Procedemos a ajustar los datos a una parábola, para ello introducimos en Matlab:

%añadir al programa anterior
m1p=[ones(20,1),w.',w.'.*w.']  %matriz de los coeficientes, tamaño 20x3, siendo la primera fila una columna de unos, la segunda el vector de las frecuencias traspuestas y la tercera es el vector de las frecuencias al cuadrado
m2p=[m1p,tr.']  %matriz ampliada con los datos del tiempo de rotura
rank(m1p)
rank(m2p)  %usamos la función que nos calcula el rango de las matrices

Si ejecutásemos el programa en este momento nos calcularía el rango de cada una de las matrices, siendo 3 para m1 y 4 para m2, por lo que sabemos que el sistema es incompatible, pasamos a calcular la resolución del sistema, que Matlab calcula automáticamente por mínimos cuadrados y a continuación definimos la parábola y la añadimos en el gráfico creado anteriormente escribiendo:

%añadir al programa anterior
xp=m1p\tr.'
y1p=xp(1)+xp(2)*w+xp(3)*w.^2   %ecuación de la parábola calculada por mínimos cuadrados
h=1:2:27
y2p=xp(1)+xp(2)*h+xp(3)*h.^2  %representamos una parábola equivalente usando más datos
hold on
plot(h,y2p,'b-')
legend('datos','parabola de minimos')  %añadimos a la gráfica anterior la parábola junto con una leyenda

3.2 Cálculo del error cuadrático medio de la parábola

Queremos tener una referencia del error cometido por la aproximación, por lo que calculamos el error cuadrático medio de la siguiente forma

%añadir al programa anterior
errorp=(sum((tr-y1p).^2))/20

3.3 Aproximación del valor 16 golpes/min

Para saber qué valor predice la recta para una frecuencia de 16 golpes por minuto ejecutamos:

%añadir al programa anterior
Zp=xp(1)+xp(2)*16+xp(3)*16.^2

3.4 Resultado

una vez escrito el código lo ejecutamos y Matlab nos da los siguientes resultados:

xp=[8.301635504889812e+02;57.398707088415634;-2.553441460170276] por lo que la ecuación de la parábola es: tr=830.16-57.39*ω-2.55*ω^2

errorp=2.947828440943336e+04 por lo que el error cuadrático medio es de 29478

zp=1.094861850100041e+03, es decir, para una frecuencia de 16 golpes por minuto se estima un tiempo de rotura de 1094 horas, y además representa la siguiente gráfica:

gráfica de los datos de partida junto con la ecuación parabólica de mínimos cuadrados

4 Cálculo y representación de una ecuación exponencial por mínimos cuadrados

4.1 Cálculo de la ecuación exponencial

Procedemos a ajustar los datos a una ecuación exponencial, para ello introducimos en Matlab:

%añadir al programa anterior
m1e=[ones(20,1),w.',exp(w.')] %matriz de los coeficientes, tamaño 20x3, siendo la primera fila una columna de unos, la segunda el vector de las frecuencias traspuestas y la tercera es el vector de los exponentes de las frecuencias
m2e=[m1e,tr.'] %matriz ampliada con los datos del tiempo de rotura
rank(m1e)
rank(m2e)  %calculamos los rangos de las matrices

Si ejecutásemos el programa en este momento nos calcularía el rango de cada una de las matrices y obtenemos que es 3 para m1 y 4 para m2 por lo que sabemos que el sistema es incompatible, pasamos a calcular la resolución del sistema, que Matlab calcula automáticamente por mínimos cuadrados y a continuación definimos la parábola y la añadimos en el gráfico creado anteriormente escribiendo:

%añadir al programa anterior
xe=m1e\tr.'  %definimos x como el vector con las soluciones del sistema
ye=xe(1)+xe(2)*w+xe(3)*exp(w)  %ecuación exponencial calculada por mínimos cuadrados
hold on
plot(w,ye,'b-')
legend('datos','exponencial de minimos')  %añadimos a la gráfica anterior la parábola junto con una leyenda

4.2 Cálculo del error cuadrático medio de la parábola

Queremos tener una referencia del error cometido por la aproximación, por lo que calculamos el error cuadrático medio de la siguiente forma

%añadir al programa anterior
errore=sum((tr-ye).^2)/20

4.3 Aproximación del valor 16 golpes/min

Para saber qué valor predice la recta para una frecuencia de 16 golpes por minuto ejecutamos:

%añadir al programa anterior
ze=xe(1)+xe(2)*16+xe(3)*exp(16)

4.4 Resultado

una vez escrito el código lo ejecutamos y Matlab nos da los siguientes resultados:

xe=[1.025976427703857e+03;-1.974268082492748;-2.319344866513919e-09] por lo que la ecuación exponencial es tr=1025.97-1.97*w-2.319e-9*exp(w)

errore=1.908149369400242e+03 el error cuadrático medio vale 1908.1

ze=9.943675284291536e+02 para una frecuencia de 16 golpes por minuto se estima un tiempo de rotura de 994.36 horas y además genera la siguiente gráfica:

gráfica de los datos de partida junto con la ecuación exponencial de mínimos cuadrados

5 Conclusiones

La distribución con el menor error cuadrático medio es la función exponencial, además es la que más coincide, visualmente, con los datos de partida por lo que nos parece que 1098.1 es la predicción más fiable de las 3