In [1]:
"""
IPython Notebook v4.0 para python 2.7
Librerías adicionales: numpy, matplotlib
Contenido bajo licencia CC-BY 4.0. Código bajo licencia MIT. (c) Sebastian Flores.
"""
# Configuracion para recargar módulos y librerías
%reload_ext autoreload
%autoreload 2
from IPython.core.display import HTML
HTML(open("style/mat281.css", "r").read())
Out[1]:
Porque área y volumen tienen potencias enteras, y las leyes son conservativas. Al integrar sobre un área o volumen, se deben obtener constantes.
Existen 7 dimensiones fundamentales, cuyas unidades básicas son definidas por el SI (Sistema Internacional de Unidades):
Existen 7 dimensiones fundamentales, cuyas unidades b asicas son definidas por el SI (Sistema Internacional de Unidades):
Para construir estos parámetros no-dimensionales:
Adimensionalizando tenemos: $$\Pi_1 = h g^x t^y$$ Luego se tiene la siguiente relación entre las dimensiones $$\begin{align}[\Pi_1] &= [h g^x t^y] = [h] [g]^x [t]^y \\ &= L \ L^x \ T^{-2x} \ T^y = L^0 T^0 \end{align}$$ Es decir, debemos resolver el sistema: $$\begin{align} 1 + x & = 0 \\ −2x + y &= 0 \end{align}$$
Se obtiene: $x = −1$ e $y = −2$, con lo cual $$ \Pi_1 = h g^x t^y = \frac{h}{gt^2}$$
A partir de eso, podemos establecer las siguientes relaciones: $$ \begin{align} h &= c g t^2 \\ g &= \frac{1}{c} \frac{h}{t^2} \\ t &= \sqrt{\frac{1}{c}\frac{h}{g}} \end{align}$$ Y todo eso ¡sin ningún conocimiento físico excepto las unidades!
Adimensionalizando la energía $E$ tenemos: $$\Pi_1 = E r^x t^y \rho^z$$ Luego se tiene la siguiente relación entre las dimensiones $$\begin{align} [\Pi_1] &= [E r^x t^y \rho^z] = [E] [r]^x [t]^y [\rho]^y \\ &= \Big(M L^2/ T^2 \Big) \ L^x \ T^y \Big(M/ L^3 \Big)^z= M^0 L^0 T^0 \end{align}$$ Es decir, debemos resolver el sistema: $$\begin{align} 1 + z & = 0 \\ 2 + x -3z &= 0 \\ -2 + y &= 0 \end{align}$$
Se obtiene: $x = −5$, $y = 2$ y $z = −1$, con lo cual $$ \Pi_1 = E r^x t^y \rho^z = \frac{E t^2}{r^5 \rho}$$
A partir de eso, podemos establecer las siguientes relaciones: $$ \begin{align} E &= c_E \frac{r^5\rho}{t^2} \\ r &= c_r \sqrt[5]{\frac{E t^2}{\rho}} \end{align}$$ Y todo eso ¡sin ningún conocimiento físico excepto las unidades!
La anécdota cuenta que en 1945 Estados Unidos hizo explotar una bomba (nombre clave Trinity) en el desierto de Nuevo México. Luego en 1947, se liberó una secuencia de fotos de la explosión. Utilizando las fotos y el análisis dimensional, el científico británico Goeffrey Taylor estimó la energía liberarada con sólo un 10% de error.
Conclusión: Coeficiente de fricción ya es un coeficiente adimensional.
Adimensionalizando la distancia de detención $d$ tenemos: $$\Pi_2 = d m^x t^y v^z$$ Luego se tiene la siguiente relación entre las dimensiones $$\begin{align} [\Pi_2] &= [d m^x t^y v^z] = [d] [m]^x [t]^y [v]^y \\ &= L M^x \ T^y \Big(L/T\Big)^z= M^0 L^0 T^0 \end{align}$$ Es decir, debemos resolver el sistema: $$\begin{align} x & = 0 \\ 1 + z &= 0 \\ y -z &= 0 \end{align}$$
Se obtiene: $x = 0$, $y = z$ y $z = −1$, con lo cual $$ \Pi_2 = d m^x t^y v^z = \frac{d}{v t}$$
Adimensionalizando la fuerza de frenado $f$ tenemos: $$\Pi_3 = f m^x t^y v^z$$ Luego se tiene la siguiente relación entre las dimensiones $$[\Pi_3] = [f m^x t^y v^z] = [f] [m]^x [t]^y [v]^y = M L/T^2 M^x \ T^y \Big(L/T\Big)^z= M^0 L^0 T^0$$ Es decir, debemos resolver el sistema: $$\begin{align} 1+x & = 0 \\ 1 + z &= 0 \\ y -z -2 &= 0 \end{align}$$
Se obtiene: $x = -1$, $y = 1$ y $z = −1$, con lo cual $$ \Pi_3 = f m^x t^y v^z = \frac{f t}{mv}$$
El teorema de Buckingham nos dice que existe por tanto una relación del tipo: $$\Phi(\Pi_1,\Pi_2,\Pi_3) = 0$$ es decir, podemos encontrar una relación para $\Pi_2$ del tipo : $$\Pi_2 =\phi(\Pi_1, \Pi_3)$$ Es decir $$\frac{d}{vt} =\phi(\mu, \frac{ft}{mv})$$ por lo tanto $$ d= vt \ \phi(\mu, \frac{ft}{mv})$$
Análisis dimensional no nos dice nada más. ¡¡Necesitamos datos!!
Considere el clásico Cherry Tree Dataset, descargable desde http://www.statsci.org/data/general/cherry.html.
Atkinson, A. C. (1982) Regression diagnostics, transformations and constructed variables (with discussion). Journal of the Royal Statistical Society, Series B, 44, 1-36.
Ryan, T. A. Jr., Joiner, B. L. and Ryan, B. F. (1985) The Minitab Student Handbook, Boston: Duxbury Press, 328-329.
Hand D. J., Daly F., Lunn A. D., McConway K. J., Ostrowski E. (1994) A Handbook of Small Data Sets. Chapman and Hall, London. Data set 210.
Los datos son:
Los datos corresponden a una muestra de 31 arboles de black cherry, en Allegheny National Forest, Pennsylvania. Los datos fueron recolectados para estimar el volumen de un árbol en función de su altura y diámetro.
In [7]:
%%bash
head data/cherry.txt
In [12]:
from matplotlib import pyplot as plt
import numpy as np
data = np.loadtxt("data/cherry.txt", skiprows=1)
D, H, V = data.T
fig = plt.figure(figsize=(16,8))
plt.plot(D,H,'o')
plt.ylim(ymin=0)
plt.show()
In [13]:
from matplotlib import pyplot as plt
import numpy as np
data = np.loadtxt("data/cherry.txt", skiprows=1)
D, H, V = data.T
fig = plt.figure(figsize=(16,8))
plt.subplot(3,1,1)
plt.plot(D,V,'r>')
plt.xlabel("Diametro")
plt.ylabel("Volumen")
plt.ylim(ymin=0)
plt.subplot(3,1,2)
plt.plot(H,V,'bo')
plt.xlabel("Altura")
plt.ylabel("Volumen")
plt.ylim(ymin=0)
plt.subplot(3,1,3)
plt.plot(D,H,'bs')
plt.xlabel("Diametro")
plt.ylabel("Altura")
plt.ylim(ymin=0)
fig.tight_layout()
plt.show()
In [14]:
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
data = np.loadtxt("data/cherry.txt", skiprows=1)
D, H, V = data.T
fig = plt.figure(figsize=(16,8))
ax = fig.add_subplot(111, projection='3d')
ax.scatter(D, H, V, s=60)
plt.xlabel("Diametro [pulgadas]")
plt.ylabel("Altura [pies]")
plt.title("Volumen [pies cubicos]")
plt.show()
¿Qué tipo de modelamiento es posible realizar con los datos?
$$¿v = c_0 + c_1 h^1 + c_2 h^2 + c_3 h^3 + c_4 d + c_5 d^2 + c_6 d^3 ?$$$$¿v = c_0 h d^2 ?$$$$¿v = c_0 h^2 d?$$$$¿v = c_0 h^3 d^3 ?$$$$¿v = c_0 + c_1 h^3 + c_2 d^3 ?$$$$¿v = c_0 d_3 h^0 + c_1 d^{5/2} h^{1/2} + ... + c_2 d^{1/2} h^{5/2} + c_3 d^0 h^3 ?$$Utilizando el teorema $\Pi$:
El teorema de Buckingham permite por tanto establecer la siguiente relación: $$ \Pi_1 = f(\Pi_2)$$ es decir $$ \frac{v}{h^3} = f\Big( \frac{d}{h}\Big)$$
El teorema de Buckingham establece: $$ \Pi_1 = f(\Pi_2)$$ PERO NADA MÁS.
In [15]:
from matplotlib import pyplot as plt
import numpy as np
data = np.loadtxt("data/cherry.txt", skiprows=1)
D, H, V = data.T
Pi_1 = V/H**3
Pi_2 = D/H
fig = plt.figure(figsize=(16,8))
plt.plot(Pi_2,Pi_1,'ob', alpha=0.5, ms=12)
plt.ylim(ymin=0)
plt.xlabel('$\Pi_2$', fontsize=20)
plt.ylabel('$\Pi_1$', fontsize=20)
plt.show()
In [18]:
from matplotlib import pyplot as plt
import numpy as np
data = np.loadtxt("data/cherry.txt", skiprows=1)
D, H, V = data.T
Pi_1 = V/H**3
Pi_2 = D/H
# Regresion lineal
m, b = np.polyfit(Pi_2, Pi_1, 1)
# Plotting
fig = plt.figure(figsize=(16,8))
plt.plot(Pi_2,Pi_1,'ob', alpha=0.5, ms=12)
x = np.sort(Pi_2)
label = "y = {0:.6f} + {1:.4} x".format(b,m)
plt.plot(x,m*x+b,'r', lw=2, label=label)
plt.xlabel('$\Pi_2$', fontsize=20)
plt.ylabel('$\Pi_1$', fontsize=20)
plt.ylim(ymin=0)
plt.legend()
plt.show()
Para obtener los coeficientes de la relación no lineal de tipo "potencia", utilizamos el clásico truco de laboratorio de física: tomar logaritmos y luego obtener los coeficientes de una relación lineal. Tomando logaritmo a: $$ \Pi_1 = k (\Pi_2)^n$$ Obtenemos $$ \log \Pi_1 = \log k + n \log \Pi_2$$
Es decir, podemos realizar una regresión lineal a $\log \Pi_1$ y $\log \Pi_2$, de modo de encontrar $\log \Pi_1 = b + m \log \Pi_2$ y luego calcular $n=m$ y $k=e^b$.
In [22]:
from matplotlib import pyplot as plt
import numpy as np
data = np.loadtxt("data/cherry.txt", skiprows=1)
D, H, V = data.T
Pi_1 = V/H**3
Pi_2 = D/H
# Regresion de potencia
x = np.log(Pi_2)
cm, cb = np.polyfit(np.log(Pi_2), np.log(Pi_1), 1)
n, k = cm, np.exp(cb)
# Plotting
fig = plt.figure(figsize=(16,8))
plt.plot(Pi_2,Pi_1,'ob', alpha=0.5, ms=12)
x = np.sort(Pi_2)
label = "y = {0:.6f} x^{1:.4}".format(k,n)
plt.plot(x,k*x**n,'r', lw=2, label=label)
#plt.plot(x,b + m*x,'g', lw=2, label=label)
plt.xlabel('$\Pi_2$', fontsize=20)
plt.ylabel('$\Pi_1$', fontsize=20)
plt.ylim(ymin=0)
plt.legend()
plt.show()
Sólo podemos saberlo analizando el error de nuestro modelo.
In [23]:
# Error modelo lineal
V_pred_lineal = (b + m*Pi_2) * H**3
error_pred_lineal = V - V_pred_lineal
print "Predicción Lineal de Volumen"
print "\tError promedio:", np.abs(error_pred_lineal).mean()
print "\tError cuadrático medio:", (error_pred_lineal**2).sum()**(0.5)/len(error_pred_lineal)
print "\tError máximo:", np.abs(error_pred_lineal).max()
In [24]:
# Error modelo no lineal
V_pred_potencia = (k*Pi_2**n) * H**3
error_pred_potencia = V - V_pred_potencia
print "Predicción No Lineal de Volumen"
print "\tError promedio:", np.abs(error_pred_potencia).mean()
print "\tError cuadrático medio:", (error_pred_potencia**2).sum()**(0.5)/len(error_pred_potencia)
print "\tError máximo:", np.abs(error_pred_potencia).max()