Relaciones potenciales y cambio de variable

Técnicas Experimentales 2021-1

Facultad de Ciencias, UNAM

En este notebooks aprenderemos:

  • Cómo graficar logaritmos.
  • Reforzaremos los conceptos aprendidos en el capitulo 4 de Introducción al análisis gráfico de datos experimentales de B. Oda Noda.
  • Aplicaremos los conceptos aprendidos en un ejemplo acerca de la Tercera Ley de Kepler.
  • Practicaremos el método de cambio de variable para obtener las constantes de una ecuación potencial a partir de un conjunto de datos.

Por hacer:

  • Lee este notebook, corre y comprende el código y teoría de la sección 1.
  • Resuelve la sección 2.

Por entregar:

  • Este mismo notebook contestado individualmente.
In [1]:
# Aqui pondermos las bibliotecas que necesitaremos en este notebook. No es necesario importarlas más de una vez. 
# Por convención se ponen en orden alfabético

import matplotlib.pyplot as plt
import numpy as np

1. La tercera ley de Kepler

Johannes Kepler publicó sus leyes de movimiento planetario entre 1609 y 1619. Estas leyes describen las órbtas de los planetas alrededor del Sol. Las tres leyes dicen que:

  1. La orbita de un planeta es un aelipse con el Sol en uno de los focos.
  2. Un segnemto de línea uqe une al planeta con el Sol barre areas iguales en tiempos iguales.
  3. El cuadrado del periodo orbital de un planeta es proporcional al cubo del radio del eje mayor de la órbita. Es decir:
\begin{equation*} \frac{R^3}{T^2}=C, \label{eq:Kepler} \tag{1} \end{equation*}

en donde $C$ es una constante, $T$ es el periodo orbital y $R$ es el radio del eje mayor de la órbita.

Utilizando mediciones del perído y el radio de la órbita de los satélites de Urano, verificaremos que se cumpla la tercera ley de Kepler para estos cuerpos celestes (Ejercicio sugerido en el capítulo 4 de Oda Noda, B. Introducción al análisis gráfico de datos experimetales).

Para ello seguiremos el siguinete procedimiento:

  1. Graficaremos las observaciones de los periodos y radios de las órbitas de las lunas de Urano para corroborar que se trata de una relación potencial.
  2. Haremos un cambio de variable para obtener una recta auxiliar.
  3. Encontraremos la pendiente y ordenada al origen de dicha recta (las constantes de la ecuación).
  4. Finalmente regresaremos a las variables originales para encontrar la ecuación que desribe la relación empírica entre radio de las órbitas y su periodo.

Los datos son:

Satélite Período (días) Radio de la órbita ($1\times10^4$ km)
I 1986U7 0.33 4.93
II 1986U8 0.37 5.33
III 1986U9 0.43 5.91
IV 1986U3 0.46 6.18
V 1986U6 0.48 6.27
VI 1986U2 0.49 6.44
VII 1986U1 0.51 6.61
VIII 1986U4 0.56 6.99
IX 1986U5 0.62 7.51
X 1985U1 0.76 8.59
XI Miranda 1.41 12.94
XII Ariel 2.52 19.10
XIII Umbriel 4.14 26.60
XIV Titania 8.71 43.59
XV Oberon 13.46 58.35

Datos tomados de Astronomy, Mayo 1986, 14, 6-22

In [2]:
import numpy as np

R = np.array([4.93, 5.33, 5.91, 6.18, 6.27, 6.44, 6.61, 6.99, 
              7.51, 8.59, 12.94, 19.10, 26.60, 43.59, 58.35])  # Radios, esta es nuestra variable dependiente, R.

T = np.array([0.33,0.37,0.43,0.46,0.48,0.49,0.51,0.56,0.62,
              0.76,1.41,2.52,4.14,8.71,13.46]) # Períodos,  esta es la variable independiente, T.
In [3]:
# No se especifica la incertidumbre , así que tomaremos la mitad de la mínima escala

dR = 0.005 # 10^4 km, incertidumbre en R
dT = 0.005 # días, incertidumbre en T

Grafiquemos los datos para inspeccionar qué tipo de relacion existe entre las variables:

In [4]:
fig = plt.figure(figsize=(6,4))
ax1= fig.add_subplot(1,1,1)

ax1.errorbar(T, R, yerr=dR, xerr=dT, 
             elinewidth=2, linewidth=0,
             marker='o', markersize=10, 
             markerfacecolor='y',markeredgecolor='0.5',ecolor='0.5')
ax1.set_xlabel('Período (días)')
ax1.set_ylabel('Radio de la órbita ($10^4$ km)')
ax1.set_title('Radio de la órbita como función del período')
plt.show()

Nota: La incertidumbre es demasiado pequeña comparada con el rango de los datos, por eso no se ven las barras de error.

Al graficar los datos vemos que la relación entre T y R es potencial; definitivamente no es una recta. Esto quiere decir que sigue una ecuación de la forma:

\begin{equation*} R=aT^n, \label{eq:potencial} \tag{2} \end{equation*}

en donde $a$ y $n$ son constantes. Si sacamos el logaritmo de la ecuación anterior obtenemos:

\begin{eqnarray*} \log{R} & = & \log{aT^n} \\ & = & \log{a}+n \log{T} \label{eq:recta} \tag{3} \end{eqnarray*}

Esta es la ecuación de una recta en donde la variable independiente es $\log{T}$, la dependiente es $\log{R}$, la pendiente es la constante $n$ y la ordenada al origen es la constante $\log{a}$. Así que si queremos encontrar las constantes de la ecuación potencial $a$ y $n$, basta con encontrar la pendiente y ordenada al origen de la ecuación de la recta. Renombrando a las variables (esto se llama cambio de variable) obtenemos:

\begin{equation*} R' = a'+nT' \label{eq:cambio_var} \tag{4} \end{equation*}

en donde $R'=\log{R}$, $T=\log{T}$, $a'=\log{a}$ y $n=n$ se queda igual. Nota que solo renombramos a las varibles originales para que sea más claro que tenemos una recta.

Para obtener la incertidumbre asociada a las variables $R'$ y $T'$ recordemos que la medición se reporta como $R\pm\delta R$, por lo cual el valor de $R'$ estará dentro del rango $[\log{(R - \delta R)}, \log{(R + \delta R)}]$. Similarmente, $T'$ estará dentro del intervalo $[\log(T - \delta T),\log(T + \delta T)]$.

Grafiquemos las variables $T'$ y $R'$:

In [5]:
# Nuevas variables
T_p = np.log(T) # usamos la función log de numpy para calular el logaritmo natural de las variables T y R
R_p = np.log(R)

# Incertidumbres - Necesitamos un valor inferior y uno superior. 
# Pondremos ambos en un arreglo de 2x15 (2 renglones y 15 columnas). El primer renglón contiene 
# los valores inferiores de la incertidumbre y el segundo los valores superiores de la incertibumbre

dR_superior = abs(np.log(R) - np.log(R+dR)) # Diferencia entre el valor de log(R) y la cota superior del intervalo
dR_inferior = abs(np.log(R) - np.log(R-dR)) # Diferencia entre el valor de log(R) y la cota inferior del intervalo
dR_p = np.array([dR_superior, dR_inferior]) 
print(np.shape(dR_p)) # Para checar el tamaño del arreglo

dT_superior = abs(np.log(T) - np.log(T+dT)) # Diferencia entre el valor de log(R) y la cota superior del intervalo
dT_inferior = abs(np.log(T) - np.log(T-dT)) # Diferencia entre el valor de log(R) y la cota inferior del intervalo
dT_p = np.array([dT_superior, dT_inferior]) 
print(np.shape(dT_p)) # Para checar el tamaño del arreglo

(2, 15)
(2, 15)
In [6]:
# Graficamos
fig = plt.figure(figsize=(6,4))
ax1= fig.add_subplot(1,1,1)

ax1.errorbar(T_p, R_p, yerr=dR_p, xerr=dT_p, 
             capsize=10, elinewidth=2, linewidth=0,
             marker='o', markersize=10, 
             markerfacecolor='g',markeredgecolor='0.5',ecolor='0.5')
ax1.set_xlabel('$\log{T}$' )
ax1.set_ylabel('$\log{R}$')
ax1.set_title('Cambio de variable - ecuación lineal asociada a T y R')
plt.show()

OJO: En esta figura exageré el tamaño de las barras de error solo para que veas que funcionan.

Efectivamente, cuando graficamos las variables T' y R' obtenemos una recta. Ahora, utilizando el método de mínimos cuadrados que aprendimos hace algunas semanas podemos estimar el valor de la pendiente $n$ y la ordenada al origen $a'$. Voy a utilizar la función min_cuadrados que hice para el notebook pasado:

In [7]:
def min_cuadrados(X,Y):
    '''Esta función calcula la pendiente y la ordenada al origen de la recta y=mx+b por mínimos 
    cuadrados a partir de los vectores de mediciones X y Y.
    Input: 
    X - arreglo de numpy 1D
    Y - arreglo de numpy 1D del mismo tamaño que X.
    Output:
    m, b : Escalares, la pendiente y ordenada al origen.
    '''
    N = len(X) # numero de mediciones
    sum_xy = np.sum(X*Y) # suma de todos los Xi*Yi
    sum_x = np.sum(X) # suma de todas las X
    sum_y = np.sum(Y) # suma de todas las Y
    sum_x2 = np.sum(X**2) # suma de todas las Xˆ2

    m = ((N*sum_xy) - (sum_x*sum_y)) / ((N*sum_x2) - (sum_x**2))
    b = (sum_y - (m*sum_x)) / N
    return(m, b)
In [8]:
n, a_p = min_cuadrados(T_p, R_p)  # llamo a la función min_cuadrados con las variables T_p y R_p como argumento
print('la pendiente es: %1.2f' %n)
print('la ordenada es: %1.2f' %a_p)

la pendiente es: 0.67
la ordenada es: 2.33
In [9]:
fig = plt.figure(figsize=(6,4))
ax1= fig.add_subplot(1,1,1)

ax1.errorbar(T_p, R_p, yerr=dR_p, xerr=dT_p, 
             capsize=10, elinewidth=2, linewidth=0,
             marker='o', markersize=10, 
             markerfacecolor='g',markeredgecolor='0.5',ecolor='0.5', label='observaciones')
ax1.plot(T_p, n*T_p+a_p, '0.5', label='ajuste min. cuadrados')
ax1.set_xlabel('$\log{T}$' )
ax1.set_ylabel('$\log{R}$')
ax1.set_title('Cambio de variable - mínimos cuadrados')
ax1.legend()
plt.show()

Ahora que conocemos $a'$ y $n$ podemos regresar a la ecuación original (Eq. 2). Recordemos que $a'=\log(a)$, por lo que $a$ es el logaritmo inverso de $a'$ (antilog o $e^{a'}$ para base e como en nuestro caso) $a=e^(a')$; y $n=n$. Sustituyendo en la ecuación potencial (Eq. 2):

\begin{equation*} R = aT^n \end{equation*}

Para comparar nuestra ecuaci'on con la 3ra Ley de Kepler, debemos escribirla en los mismos términos que la ecuaci on 1, por lo que reacomodamos la ecuación anterior:

$$\frac{R}{T^n}=a$$

Elevando al cubo:

$$\frac{R^3}{T^{3n}}=a^3$$

.

Sustituyendo el valor de n=0.67:

\begin{equation*} \frac{R^3}{T^{2.01}}= C, \label{eq:Kepler_exp} \tag{5} \end{equation*}

que es la tercera ley de Kepler (Ecuación 1). La constante C es $a^3$.

2. Ejercicio de aplicación

Ahora repite el procedimiento anterior para encontrar la ecuación que describe la relación entre los siguientes datos:

X Y
0.3 13.0
0.6 9.5
0.9 8.5
1.2 7.5
1.5 6.5
3.0 5.0

Considera las incertidumbres. No tienen unidades.

Por hacer:

  1. Grafica los datos.
In [10]:
X = np.array([0.3, 0.6, 0.9, 1.2, 1.5, 3.0])
Y = np.array([13.0, 9.5, 8.5, 7.5, 6.5, 5.0])
dX = 0.05
dY = 0.05
In [11]:
fig = plt.figure(figsize=(6,4))
ax1= fig.add_subplot(1,1,1)

ax1.errorbar(X, Y, yerr=dY, xerr=dX, 
             elinewidth=2, linewidth=0,
             marker='o', markersize=10, 
             markerfacecolor='purple',markeredgecolor='0.5',ecolor='0.5')
ax1.set_xlabel('X' )
ax1.set_ylabel('Y')
plt.show()
  1. Si la relación es potencial, haz un cambio de variables.
    • Grafica las nuevas variables.
    • Calcula la pendiente y ordenada al origen de la recta auxiliar.
    • Grafica la recta auxiliar junto con las nuevas variables.
In [12]:
X_p = np.log(X)
Y_p = np.log(Y)

dX_superior = abs(np.log(X) - np.log(X+dX)) # Diferencia entre el valor de log(X) y la cota superior del intervalo
dX_inferior = abs(np.log(X) - np.log(X-dX)) # Diferencia entre el valor de log(X) y la cota inferior del intervalo
dX_p = np.array([dX_superior, dX_inferior]) 
print(np.shape(dX_p)) # Para checar el tamaño del arreglo

dY_superior = abs(np.log(Y) - np.log(Y+dY)) # Diferencia entre el valor de log(Y) y la cota superior del intervalo
dY_inferior = abs(np.log(Y) - np.log(Y-dY)) # Diferencia entre el valor de log(Y) y la cota inferior del intervalo
dY_p = np.array([dY_superior, dY_inferior]) 
print(np.shape(dY_p)) # Para checar el tamaño del arreglo

(2, 6)
(2, 6)
In [13]:
fig = plt.figure(figsize=(6,4))
ax1= fig.add_subplot(1,1,1)

ax1.errorbar(X_p, Y_p, yerr=dY_p, xerr=dX_p, 
             elinewidth=2, linewidth=0,
             marker='o', markersize=10, 
             markerfacecolor='orchid',markeredgecolor='0.5',ecolor='0.5')
ax1.set_xlabel('Xp' )
ax1.set_ylabel('Yp')
plt.show()
In [14]:
m2, b2 = min_cuadrados(X_p, Y_p)
In [15]:
fig = plt.figure(figsize=(6,4))
ax1= fig.add_subplot(1,1,1)

ax1.errorbar(X_p, Y_p, yerr=dY_p, xerr=dX_p, 
             elinewidth=2, linewidth=0,
             marker='o', markersize=10, 
             markerfacecolor='orchid',markeredgecolor='0.5',ecolor='0.5',label='datos')
ax1.plot(X_p,m2 * X_p + b2,'-', color='0.3', label='min cuadrados')
ax1.set_xlabel('Xp')
ax1.set_ylabel('Yp')
ax1.legend()
plt.show()
  1. Regresa a las variables originales
    • Escribe la ecuación que relaciona a X y Y.
    • Grafica la curva ajustada junto con los datos originales.
In [16]:
print('La ecuación que relaciona a X y Y es: Y=%1.1fX^{%1.1f}' %(np.exp(b2), m2 ))

La ecuación que relaciona a X y Y es: Y=7.9X^{-0.4}

La ecuación que relaciona a $X$ y $Y$ es:

$$Y=7.9X^{-0.4}$$
In [19]:
fig = plt.figure(figsize=(6,4))
ax1= fig.add_subplot(1,1,1)

X_ajuste = np.linspace(0.3,3.0,30)

ax1.errorbar(X, Y, yerr=dY, xerr=dX, 
             elinewidth=2, linewidth=0,
             marker='o', markersize=10, 
             markerfacecolor='purple',markeredgecolor='0.5',ecolor='0.5',label='datos')
ax1.plot(X_ajuste,np.exp(b2)*(X_ajuste**m2),'-', color='0.3', label='ajuste cambio de variable')
ax1.set_xlabel('X')
ax1.set_ylabel('Y')
ax1.legend()
plt.show()
In [ ]:
In [ ]: