Hola, gente.
Hoy traigo un post de un tema que mezcla Física, Matemáticas, Métodos numéricos y algo de Programación en Python; espero que les guste.
¡Por cierto!, dedico este post a @Hazail que es una empedernida mujer de ciencia y la dueña de mis suspiros.
Para entender este post necesitamos tener algunas nociones sobre calculo infinitesimal, mecánica elemental y lo básico de programación. Sin embargo, no mencionare algo sin antes dar una breve explicación de lo que estoy usando, esperando ser lo suficientemente claro.
Nuestro propósito es calcular de alguna manera la posición del extremo de un resorte dispuesto en forma horizontal y que esta sujeto del otro extremo a un pared, usando ciertas condiciones iniciales. Luego, haremos una sencilla animación usando Python.
Pero para lograr todo ello, primero estudiaremos lo que hay detrás.
Todos hemos visto un resorte en acción. Su comportamiento es hipnótico. El va y ven de su danza interminable puede convertir a cualquier persona en un gato curioso, deseoso de jugar con el.
Pero ademas de hermoso, el movimiento de un resorte guarda comportamientos que a primera vista parecen imposibles.
Pero tienen explicación, solo es cuestión de tener un verdadero deseo de conocerla.
Para nuestro objetivo, primero debemos aclarar como se comporta, a grandes rasgos, un resorte.
Imaginemos nuestro resorte atado a la pared. Este resorte se encuentra en un estado llamado "estado de equilibrio" o también se dice que el resorte esta en su "posición de equilibrio". Esto significa que no esta ni estirado, ni comprimido, o sea, esta en su forma "natural".
No es difícil pensar en lo que pasa cuando lo estiramos. Cuando nosotros hacemos eso, el resorte nos "jala" en dirección contraria.
y lo mismo sucede si en vez de estirar, lo comprimimos.
También es casi seguro, que conozcas el hecho de que entre mas estiras (o comprimes) el resorte, mayor es la fuerza con la que este se opone. Esto significa que la fuerza cambia con respecto a la distancia en la que se encuentra estirado (o comprimido). En física se conoce cuantitativamente esta propiedad y esta dada por la siguiente ecuación:
En la formula, "x" significa la posición en el marco de referencia usado. ( Aquí, nuestro marco tendrá su origen en la posición de equilibro del resorte.) y k es una constante que tiene que ver con el resorte (material del que esta hecho, calibre del alambre, etc).
Por ejemplo, supongamos que comprimimos un resorte de acero hecho de alambre muy delgado. En este caso, la contante k tiene un valor de 1.5 Newtons por cada metro. Digamos que lo comprimimos 0.3 metros. Para calcular la fuerza con la que solo debemos sustituir los valores en la formula:
El valor negativo de x es porque, en nuestro sistema de coordenadas, el cero esta en la posición de equilibro, entonces, si el resorte se mueve hacia la izquierda, las posiciones son negativas (pueden pasar a ver mi post de vectores si esto no queda muy claro).
Con el calculo anterior obtenemos un resultado coherente, pues la fuerza positiva resultante significa que el resorte "empuja" hacia la derecha.
Ahora, pensemos un una configuración diferente de nuestro sistema de estudio. Ahora pondremos un cuerpo con masa en el extremo libre del resorte.
El movimiento de ese cuerpo estará directamente ligado al comportamiento del resorte y gracias Newton, sabemos calcularlo.
La segunda ley de Newton nos dice que la fuerza total que siente un cuerpo, determina su movimiento de acuerdo a la siguiente relación:
donde "m" es la masa del cuerpo que "siente" la fuerza (sería ideal un chiste de star wars aquí) y a con flechita arriba es el vector aceleración (ver mi post de vectores si hay dudas de eso).
Ahora, en nuestro sistema vamos a despreciar la fricción (¬¬ típico de físicos novatos como yo). El vector fuerza total que "siente" el cuerpo tiene dos componentes: una vertical una horizontal. Pero la fuerza total vertical es nula, pues el peso del cuerpo es contrarrestada con la fuerza que ejerce la superficie sobre la que esta el sistema.
Así, al final de cuentas, podemos trabajar solo con la componente horizontal. Y, como no hay fricción, la única fuerza que siente el cuerpo, es la que el resorte ejerce, la cual es la que determinara el movimiento.
Ya vimos como se calcula esa fuerza, asi que podemos reescribir la segunda ley de Newton, aplicado a nuestro problema de la siguiente forma:
Una vez entendida la parte física, la parte matemática es muy simple. Nos apoyaremos en algunas definiciones de mecánica para describir el movimiento de nuestro objeto rosa en forma de una ecuación diferencial (signifique lo que signifique).
De las clases de física básica, sabemos que la tasa de cambio de la posición mientras pasa el tiempo, le llamamos velocidad. Una forma elegante de decir lo mismo es así:
Podemos decir algo similar para definir la aceleración. Esta es cómo va cambiando la velocidad mientras pasa el tiempo:
No son mas que simples derivadas. Ahora, como se pudo ver, la aceleración la obtuvimos derivando dos veces seguidas la posición. Esto lo podemos escribir así:
Y con esto podemos rescribir la segunda ley de Newton aplicado a nuestro problema ("otra vez" Forrest Gump 1994).
Ahora solo aplicamos un poquito de álgebra para darle una presentación mas pomposa:
A nuestra ultima ecuación la podemos llamar "ecuación diferencial" (no importa ahora porque), y describe exactamente el movimiento que tendrá nuestro objeto rosa. En este caso, la ecuación tiene solución analítica. Pero para entenderla necesitamos teoría de ecuaciones diferenciales y eso es del diablo, así que mejor usaremos otra estrategia.
La idea básica los "métodos numéricos" es calcular "paso a paso". Recomiendo ampliamente que investiguen un poco sobre este tema. Por ahora tendrán que conformarse con una probadita: El método de Euler.
Este método surge de una idea bastante simple. De geometría analítica sabemos que la pendiente de un recta la podemos calcular conociendo dos de sus puntos:
donde los puntos tienes coordenadas
Luego, una de las primeras cosas que aprendemos de calculo diferencial es que la derivada de y con respecto a x, evaluada en determinado punto, nos da la pendiente de la recta tangente a ese punto, es decir:
podemos combinar estas dos ecuaciones interpretando (a,b) como (x1,y1):
Ahora despejemos y2:
Esta ultima ecuación nos dice que, si conocemos (X1,Y1) y tenemos la expresión de la derivada, podemos obtener el valor de y2.
Con un ligero cambio en la notación, para posteriores usos, al final nos queda:
Con la anterior formula, si conocemos un primer punto, entonces podemos calcular el que sigue. Luego podríamos usar ese punto calculado y repetir el proceso para calcular un tercero. Luego podemos hacer lo mismo y calcular un cuarto punto y así sucesivamente. Si graficáramos esos puntos ser vería un comportamiento similar al que de la siguiente imagen.
Y podemos ver también que entre mas pequeño es el intevalo X2 - X1 mas nos acercamos a dibujar una linea continua. Una linea continua representaría las soluciones exactas de la ecuación diferencial.
Si ustedes un día estudian un poco mas a fondo el tema de métodos numéricos, verán que uno de los problemas con estos es que son susceptibles a errores. El método de Euler es el más sencillo de ellos pero también es el que más acumula errores.
Una ligera modificación reduce considerablemente este hecho: en vez de usar la pendiente del punto anterior, hay que usar la pendiente del punto de destino. La ecuación queda así:
Usando definiciones mencionadas más arriba llegamos a lo siguiente:
Nota: Al aplicar el método, nos quedan dos ecuaciones porque para calcular la posición necesitamos su derivada, lo cual es la velocidad. Pero ademas necesitamos calcular la velocidad, que es la segunda derivada de la posición, de la cual tenemos una expresión dada por la ecuación diferencial del movimiento del resorte. Se puede ver también que la pendiente usada para la posición es la pendiente del punto de destino y la que se usa para calcular la velocidad es la pendiente del punto anterior. Esto quiere decir que usamos Euler para calcular la velocidad y Euler mejorado para calcular la posición (sí, es posible).
El codigo completo lo pondré al final
Las siguientes imágenes png serán las utilizadas para el programa:
https://k60.kn3.net/taringa/4/C/E/5/A/C/mladiciones/099.png
https://k60.kn3.net/taringa/B/E/8/5/4/7/mladiciones/BF2.png
https://k60.kn3.net/taringa/D/2/C/C/3/B/mladiciones/F5D.png
https://k60.kn3.net/taringa/2/F/B/4/0/1/mladiciones/8E3.png
Deben guardarlas con el nombre correcto (piso, masarosa,resorte, linea)
Iniciamos un script de python en algún editor de texto plano. Primero importamos dos librerias que utilizaremos. (previa instalación)
La imagen de la masarosa.png mide 100x100, resorte.png 260x100, el piso tiene 600px de largo y la linea tiene 400 de alto. Con base en estas dimensiones y teniendo en cuenta como usa las coordenadas pygame, creamos las siguientes "constantes":
Para la simulación necesitamos conocer cuatro valores iniciales y el tamaño de los intervalos de tiempo. En este caso, los intervalos de tiempo duraran el tiempo que un fotograma esté en pantalla.
Ahora, con pygame creamos una ventana con titulo y color de fondo, un reloj y la es estructura basica de un loop principal de un juego:
Ahora creamos unas superficies a partir de unas imágenes. También creamos una variable "x" y una "v" y les asignamos los valores de posición y velocidad inicial. Colocamos el código justo antes del bloque "while":
Ahora debemos "pintar" todo en la ventana que creamos. Esto lo hacemos dentro del loop principal antes de "pygame.display.flip()":
Si corremos nuestro código nos queda así:
Donde se ve que el resorte no esta "estirado" como debería. Para simular esto, vamos a escalar la imagen original aumentando (o disminuyendo) su tamaño horizontal lo mismo que la masa rosa se desplazó. Para esto creamos una nueva variable "resorte" dentro del loop la cual guardara la superficie escalada. Esta nueva superficie sera la que "pintaremos" en vez de "resorte original, asi que habra que cambiar la linea "ventana.blit(resorte_original, (coordenada_equilibro_resorte, 150))" por "ventana.blit(resorte, (coordenada_equilibro_resorte, 150))"
Si corremos el programa ahora nos aparece el sistema en su configuración inicial:
Para animar, solo hay que actualizar la posición. Esta la calcularemos con el método que aprendimos mas arriba. Primero necesitamos calcular la velocidad, luego usamos esa velocidad para calcular la nueva posición. Colocamos el código después de la ultima sentencia en el bucle:
Y el resultado al correr el programa es:
CODIGO COMPLETO:
Eso fue todo. Comenten.
Hoy traigo un post de un tema que mezcla Física, Matemáticas, Métodos numéricos y algo de Programación en Python; espero que les guste.
¡Por cierto!, dedico este post a @Hazail que es una empedernida mujer de ciencia y la dueña de mis suspiros.
Para entender este post necesitamos tener algunas nociones sobre calculo infinitesimal, mecánica elemental y lo básico de programación. Sin embargo, no mencionare algo sin antes dar una breve explicación de lo que estoy usando, esperando ser lo suficientemente claro.
Nuestro propósito es calcular de alguna manera la posición del extremo de un resorte dispuesto en forma horizontal y que esta sujeto del otro extremo a un pared, usando ciertas condiciones iniciales. Luego, haremos una sencilla animación usando Python.
Pero para lograr todo ello, primero estudiaremos lo que hay detrás.
Todos hemos visto un resorte en acción. Su comportamiento es hipnótico. El va y ven de su danza interminable puede convertir a cualquier persona en un gato curioso, deseoso de jugar con el.
Pero ademas de hermoso, el movimiento de un resorte guarda comportamientos que a primera vista parecen imposibles.
Pero tienen explicación, solo es cuestión de tener un verdadero deseo de conocerla.
Para nuestro objetivo, primero debemos aclarar como se comporta, a grandes rasgos, un resorte.
Imaginemos nuestro resorte atado a la pared. Este resorte se encuentra en un estado llamado "estado de equilibrio" o también se dice que el resorte esta en su "posición de equilibrio". Esto significa que no esta ni estirado, ni comprimido, o sea, esta en su forma "natural".
No es difícil pensar en lo que pasa cuando lo estiramos. Cuando nosotros hacemos eso, el resorte nos "jala" en dirección contraria.
y lo mismo sucede si en vez de estirar, lo comprimimos.
También es casi seguro, que conozcas el hecho de que entre mas estiras (o comprimes) el resorte, mayor es la fuerza con la que este se opone. Esto significa que la fuerza cambia con respecto a la distancia en la que se encuentra estirado (o comprimido). En física se conoce cuantitativamente esta propiedad y esta dada por la siguiente ecuación:
En la formula, "x" significa la posición en el marco de referencia usado. ( Aquí, nuestro marco tendrá su origen en la posición de equilibro del resorte.) y k es una constante que tiene que ver con el resorte (material del que esta hecho, calibre del alambre, etc).
Por ejemplo, supongamos que comprimimos un resorte de acero hecho de alambre muy delgado. En este caso, la contante k tiene un valor de 1.5 Newtons por cada metro. Digamos que lo comprimimos 0.3 metros. Para calcular la fuerza con la que solo debemos sustituir los valores en la formula:
El valor negativo de x es porque, en nuestro sistema de coordenadas, el cero esta en la posición de equilibro, entonces, si el resorte se mueve hacia la izquierda, las posiciones son negativas (pueden pasar a ver mi post de vectores si esto no queda muy claro).
Con el calculo anterior obtenemos un resultado coherente, pues la fuerza positiva resultante significa que el resorte "empuja" hacia la derecha.
Ahora, pensemos un una configuración diferente de nuestro sistema de estudio. Ahora pondremos un cuerpo con masa en el extremo libre del resorte.
El movimiento de ese cuerpo estará directamente ligado al comportamiento del resorte y gracias Newton, sabemos calcularlo.
La segunda ley de Newton nos dice que la fuerza total que siente un cuerpo, determina su movimiento de acuerdo a la siguiente relación:
donde "m" es la masa del cuerpo que "siente" la fuerza (sería ideal un chiste de star wars aquí) y a con flechita arriba es el vector aceleración (ver mi post de vectores si hay dudas de eso).
Ahora, en nuestro sistema vamos a despreciar la fricción (¬¬ típico de físicos novatos como yo). El vector fuerza total que "siente" el cuerpo tiene dos componentes: una vertical una horizontal. Pero la fuerza total vertical es nula, pues el peso del cuerpo es contrarrestada con la fuerza que ejerce la superficie sobre la que esta el sistema.
Así, al final de cuentas, podemos trabajar solo con la componente horizontal. Y, como no hay fricción, la única fuerza que siente el cuerpo, es la que el resorte ejerce, la cual es la que determinara el movimiento.
Ya vimos como se calcula esa fuerza, asi que podemos reescribir la segunda ley de Newton, aplicado a nuestro problema de la siguiente forma:
Una vez entendida la parte física, la parte matemática es muy simple. Nos apoyaremos en algunas definiciones de mecánica para describir el movimiento de nuestro objeto rosa en forma de una ecuación diferencial (signifique lo que signifique).
De las clases de física básica, sabemos que la tasa de cambio de la posición mientras pasa el tiempo, le llamamos velocidad. Una forma elegante de decir lo mismo es así:
Podemos decir algo similar para definir la aceleración. Esta es cómo va cambiando la velocidad mientras pasa el tiempo:
No son mas que simples derivadas. Ahora, como se pudo ver, la aceleración la obtuvimos derivando dos veces seguidas la posición. Esto lo podemos escribir así:
Y con esto podemos rescribir la segunda ley de Newton aplicado a nuestro problema ("otra vez" Forrest Gump 1994).
Ahora solo aplicamos un poquito de álgebra para darle una presentación mas pomposa:
A nuestra ultima ecuación la podemos llamar "ecuación diferencial" (no importa ahora porque), y describe exactamente el movimiento que tendrá nuestro objeto rosa. En este caso, la ecuación tiene solución analítica. Pero para entenderla necesitamos teoría de ecuaciones diferenciales y eso es del diablo, así que mejor usaremos otra estrategia.
La idea básica los "métodos numéricos" es calcular "paso a paso". Recomiendo ampliamente que investiguen un poco sobre este tema. Por ahora tendrán que conformarse con una probadita: El método de Euler.
Este método surge de una idea bastante simple. De geometría analítica sabemos que la pendiente de un recta la podemos calcular conociendo dos de sus puntos:
donde los puntos tienes coordenadas
Luego, una de las primeras cosas que aprendemos de calculo diferencial es que la derivada de y con respecto a x, evaluada en determinado punto, nos da la pendiente de la recta tangente a ese punto, es decir:
podemos combinar estas dos ecuaciones interpretando (a,b) como (x1,y1):
Ahora despejemos y2:
Esta ultima ecuación nos dice que, si conocemos (X1,Y1) y tenemos la expresión de la derivada, podemos obtener el valor de y2.
Con un ligero cambio en la notación, para posteriores usos, al final nos queda:
Con la anterior formula, si conocemos un primer punto, entonces podemos calcular el que sigue. Luego podríamos usar ese punto calculado y repetir el proceso para calcular un tercero. Luego podemos hacer lo mismo y calcular un cuarto punto y así sucesivamente. Si graficáramos esos puntos ser vería un comportamiento similar al que de la siguiente imagen.
Y podemos ver también que entre mas pequeño es el intevalo X2 - X1 mas nos acercamos a dibujar una linea continua. Una linea continua representaría las soluciones exactas de la ecuación diferencial.
Si ustedes un día estudian un poco mas a fondo el tema de métodos numéricos, verán que uno de los problemas con estos es que son susceptibles a errores. El método de Euler es el más sencillo de ellos pero también es el que más acumula errores.
Una ligera modificación reduce considerablemente este hecho: en vez de usar la pendiente del punto anterior, hay que usar la pendiente del punto de destino. La ecuación queda así:
Usando definiciones mencionadas más arriba llegamos a lo siguiente:
Nota: Al aplicar el método, nos quedan dos ecuaciones porque para calcular la posición necesitamos su derivada, lo cual es la velocidad. Pero ademas necesitamos calcular la velocidad, que es la segunda derivada de la posición, de la cual tenemos una expresión dada por la ecuación diferencial del movimiento del resorte. Se puede ver también que la pendiente usada para la posición es la pendiente del punto de destino y la que se usa para calcular la velocidad es la pendiente del punto anterior. Esto quiere decir que usamos Euler para calcular la velocidad y Euler mejorado para calcular la posición (sí, es posible).
El codigo completo lo pondré al final
Las siguientes imágenes png serán las utilizadas para el programa:
https://k60.kn3.net/taringa/4/C/E/5/A/C/mladiciones/099.png
https://k60.kn3.net/taringa/B/E/8/5/4/7/mladiciones/BF2.png
https://k60.kn3.net/taringa/D/2/C/C/3/B/mladiciones/F5D.png
https://k60.kn3.net/taringa/2/F/B/4/0/1/mladiciones/8E3.png
Deben guardarlas con el nombre correcto (piso, masarosa,resorte, linea)
Iniciamos un script de python en algún editor de texto plano. Primero importamos dos librerias que utilizaremos. (previa instalación)
[color=#000000]
#!/usr/bin/env python
# -*- coding: utf-8 -*-
import numpy
import pygame
import sys
[/color]
La imagen de la masarosa.png mide 100x100, resorte.png 260x100, el piso tiene 600px de largo y la linea tiene 400 de alto. Con base en estas dimensiones y teniendo en cuenta como usa las coordenadas pygame, creamos las siguientes "constantes":
[color=#000000]
dimension_ventana = (600,400)
coordenada_equilibro_masa = 250
coordenada_equilibro_resorte = -10
coordenada_linea = (300,0)
coordenada_piso = (0,250)
ancho_resorte = 260
color_fondo = (255,255,255)
fps = 60
[/color]
Para la simulación necesitamos conocer cuatro valores iniciales y el tamaño de los intervalos de tiempo. En este caso, los intervalos de tiempo duraran el tiempo que un fotograma esté en pantalla.
[color=#000000]
dt = 1/float(fps)
posicion_incial = 100
velocidad_inicial = 0
k = 9
masa = 1
[/color]
Ahora, con pygame creamos una ventana con titulo y color de fondo, un reloj y la es estructura basica de un loop principal de un juego:
[color=#000000]
ventana = pygame.display.set_mode(dimension_ventana)
pygame.display.set_caption("Simulación")
reloj = pygame.time.Clock()
while True:
for evento in pygame.event.get():
if evento.type == pygame.QUIT:
sys.exit(0)
pygame.display.flip()
ventana.fill(color_fondo)
reloj.tick(fps)
[/color]
Ahora creamos unas superficies a partir de unas imágenes. También creamos una variable "x" y una "v" y les asignamos los valores de posición y velocidad inicial. Colocamos el código justo antes del bloque "while":
[color=#000000]
piso = pygame.image.load('piso.png')
linea = pygame.image.load('linea.png').convert_alpha()
masarosa = pygame.image.load('masarosa.png')
resorte_original = pygame.image.load('resorte.png').convert_alpha()
x = posicion_incial
v = velocidad_inicial
[/color]
Ahora debemos "pintar" todo en la ventana que creamos. Esto lo hacemos dentro del loop principal antes de "pygame.display.flip()":
[color=#000000]
ventana.blit(resorte_original, (coordenada_equilibro_resorte, 150))
ventana.blit(masarosa, (coordenada_equilibro_masa + x, 150))
ventana.blit(linea, coordenada_linea)
ventana.blit(piso, coordenada_piso)
[/color]
Si corremos nuestro código nos queda así:
Donde se ve que el resorte no esta "estirado" como debería. Para simular esto, vamos a escalar la imagen original aumentando (o disminuyendo) su tamaño horizontal lo mismo que la masa rosa se desplazó. Para esto creamos una nueva variable "resorte" dentro del loop la cual guardara la superficie escalada. Esta nueva superficie sera la que "pintaremos" en vez de "resorte original, asi que habra que cambiar la linea "ventana.blit(resorte_original, (coordenada_equilibro_resorte, 150))" por "ventana.blit(resorte, (coordenada_equilibro_resorte, 150))"
[color=#000000]
resorte = pygame.transform.scale(resorte_original, (int(numpy.abs(ancho_resorte+x)),100))
ventana.blit(resorte, (coordenada_equilibro_resorte, 150))
[/color]
Si corremos el programa ahora nos aparece el sistema en su configuración inicial:
Para animar, solo hay que actualizar la posición. Esta la calcularemos con el método que aprendimos mas arriba. Primero necesitamos calcular la velocidad, luego usamos esa velocidad para calcular la nueva posición. Colocamos el código después de la ultima sentencia en el bucle:
[color=#000000]
v -= (k/masa)*x*dt
x += v*dt
[/color]
Y el resultado al correr el programa es:
CODIGO COMPLETO:
#!/usr/bin/env python
# -*- coding: utf-8 -*-
import numpy
import pygame
import sys
dimension_ventana = (600,400)
coordenada_equilibro_masa = 250
coordenada_equilibro_resorte = -10
coordenada_linea = (300,0)
coordenada_piso = (0,250)
ancho_resorte = 260
color_fondo = (255,255,255)
fps = 60
dt = 1/float(fps)
posicion_incial = 100
velocidad_inicial = 0
k = 9
masa = 1
ventana = pygame.display.set_mode(dimension_ventana)
pygame.display.set_caption("Simulación")
ventana.fill(color_fondo)
reloj = pygame.time.Clock()
piso = pygame.image.load('piso.png')
linea = pygame.image.load('linea.png').convert_alpha()
masarosa = pygame.image.load('masarosa.png')
resorte_original = pygame.image.load('resorte.png').convert_alpha()
x = posicion_incial
v = velocidad_inicial
while True:
for evento in pygame.event.get():
if evento.type == pygame.QUIT:
sys.exit(0)
resorte = pygame.transform.scale(resorte_original, (int(numpy.abs(ancho_resorte+x)),100))
ventana.blit(resorte, (coordenada_equilibro_resorte, 150))
ventana.blit(masarosa, (coordenada_equilibro_masa + x, 150))
ventana.blit(linea, coordenada_linea)
ventana.blit(piso, coordenada_piso)
pygame.display.flip()
ventana.fill(color_fondo)
reloj.tick(fps)
v -= (k/masa)*x*dt
x += v*dt
Eso fue todo. Comenten.