Muestra las diferencias entre dos versiones de la página.
Ambos lados, revisión anterior Revisión previa Próxima revisión | Revisión previa | ||
acemu:ensayos:jornadas:03_2010:runge-kutta [2012/06/18 16:44] luis |
acemu:ensayos:jornadas:03_2010:runge-kutta [2012/06/19 16:11] (actual) luis [Inclusión de la fricción del aire en la simulación] |
||
---|---|---|---|
Línea 1: | Línea 1: | ||
[[ACEMU: | [[ACEMU: | ||
- | ====== Simulación por el Método Runge-Kutta ====== | + | ====== Simulación por el Método Runge-Kutta |
==== Introducción ==== | ==== Introducción ==== | ||
- | En este artículo, intentamos trasmitir | + | En esta sección, intentamos trasmitir |
- | Para no dejar ideas fuera de contexto, dejamos el posteado tal cual fué realizado, por eso el diálogo en primera persona. | + | |
- | ==== La clase de Kenneth | + | ==== Análisis |
- | Estimados: aprovechando el feriado, me puse a jugar con la curva de empuje del motor (la " | + | Trabajando |
\\ | \\ | ||
- | Para empezar, aclaremos cuales fueron las aproximaciones | + | Para comenzar, aclaremos cuales fueron las aproximaciones |
\\ | \\ | ||
- | - usé la curva de empuje obtenida, pasando los datos a Newton. | + | - 1° se tomó la curva de empuje obtenida, pasando los datos a Newton. |
- | - asumí | + | - 2° se asumieron |
* Etapa 1: con el motor haciendo empuje sobre el vector. | * Etapa 1: con el motor haciendo empuje sobre el vector. | ||
* Etapa 2: vuelo libre (sólo actúa la fuerza de la gravedad) | * Etapa 2: vuelo libre (sólo actúa la fuerza de la gravedad) | ||
- | También | + | También |
< | < | ||
Línea 33: | Línea 32: | ||
m es la masa del vector | m es la masa del vector | ||
\\ | \\ | ||
- | Asumí | + | Se asumió |
\\ | \\ | ||
- | Como E(t) lo tenemos | + | Como E(t) se encuentra |
+ | Todo fue programado empleando el lenguaje AWK [[http:// | ||
\\ | \\ | ||
- | Asumí | + | También se asumió |
\\ | \\ | ||
- | == Resultados: == | + | ==== Resultados: |
//Tiempo de vuelo total//: aprox 20 s\\ | //Tiempo de vuelo total//: aprox 20 s\\ | ||
- | //Velocidad máxima alcanzada//: | + | //Velocidad máxima alcanzada//: |
//Apogeo//: 485 m a los 10 s\\ | //Apogeo//: 485 m a los 10 s\\ | ||
- | //Velocidad final del vector//: -97 m/s (dejó un lindo agujerito en la tierra)\\ | + | //Velocidad final del vector//: -97 m/s (recordemos que no se considera apertura de paracaídas)\\ |
\\ | \\ | ||
- | == Mejoras | + | ==== Como mejoras |
\\ | \\ | ||
- | 1) Usar RK de orden 4\\ | + | 1) Usar Runge-Kutta |
2) Modificar la ecuación para introducir fricción con el aire\\ | 2) Modificar la ecuación para introducir fricción con el aire\\ | ||
3) Incluir la apertura del paracaídas en la simulación\\ | 3) Incluir la apertura del paracaídas en la simulación\\ | ||
- | 4) Arrancar la simulación cuando el **Empuje** supera el **Peso** del cohete, | + | 4) Arrancar la simulación cuando el **Empuje** supera el **Peso** del cohete, |
\\ | \\ | ||
- | == El fuente del programa AWK: == | + | ==== El fuente del programa AWK ==== |
+ | < | ||
Código: [Seleccionar] | Código: [Seleccionar] | ||
# Programa AWK para calcular altura y velocidad de un cohete a partir de | # Programa AWK para calcular altura y velocidad de un cohete a partir de | ||
Línea 99: | Línea 99: | ||
} | } | ||
} | } | ||
+ | </ | ||
- | Adjunto los resultados obtenidos: tabla de texto con los datos de tiempo (s), altura (m) y velocidad (m/s) | + | ==== La Ecuación ==== |
- | En cuanto a la ecuación | + | En cuanto a la ecuación la misma se basa en la [[http:// |
- | En este caso F(t)=E(t)-P(t) es la fuerza resultante que actúa sobre el cohete, que sería la resta del empuje del motor menos el peso del cohete (suponemos además que el cohete despega y vuela de manera completamente vertical, eso me olvidé de aclararlo). | + | En este caso F(t)=E(t)-P(t), es la fuerza resultante que actúa sobre el cohete, que sería la resta del empuje del motor, menos el peso del cohete (suponemos además que el cohete despega y vuela de manera completamente vertical). |
- | Todas son funciones del tiempo, porque el empuje del motor varía a medida que transcurre el tiempo y lo mism ocurre con la masa del cohete, de ahi que el peso P también varíe con el tiempo. Una aproximación | + | Todas son funciones del tiempo, porque el empuje del motor varía a medida que transcurre el tiempo y lo mismo ocurre con la masa del cohete, de ahí que el peso P también varíe con el tiempo.\\ Una aproximación fue suponer, que P(t)=constante, |
- | Del otro lado de la ecuación tenés la masa del cohete multimplicada por la aceleración que sufre el cohete. | + | < |
+ | F(t)=E(t)-P | ||
+ | </ | ||
- | La aceleración no es otra cosa que la derivada segunda | + | Del otro lado de la ecuación tenemos |
- | Por eso es que podés despejar | + | La aceleración, |
- | y(t)'' | + | < |
+ | (y'' | ||
+ | </code> | ||
- | o sea | + | para indicar que es la derivada segunda. |
+ | Por eso es que podemos despejar la aceleración que sufre el cohete en la ecuación de la siguiente manera: | ||
+ | |||
+ | < | ||
+ | y(t)'' | ||
+ | </ | ||
+ | o sea | ||
+ | < | ||
y(t)'' | y(t)'' | ||
+ | </ | ||
- | donde los valores de E(t) son los que obtuvimos en la toma de datos del ensayo del motor. | + | donde los valores de E(t) son los que obtuvimos en la toma de datos del ensayo del motor.\\ |
- | Teniendo en cuenta que el peso P=mg, en realidad podríamos escribir: | + | Teniendo en cuenta que el peso P=mg, en realidad podríamos escribir:\\ |
+ | < | ||
y'' | y'' | ||
+ | </ | ||
- | Esa es la ecuación diferencial que estoy resolviendo numéricamente. OK, cómo hago eso no lo voy a explicar acá, porque el fundamento es más complicado, pero es una simple receta, una vez que entendiste de donde sale la receta, la aplicás a ojos cerrados y no tiene ningún misterio. | + | Esa es la ecuación diferencial que estamos |
- | Lo único que determina la solución serán las condiciones iniciales (altura y velocidad inicial) que asumí que ambas valían cero. | + | Lo único que determina la solución serán las condiciones iniciales (altura y velocidad inicial) que asumimos |
Una mejora que puede agregarse es tener en consideración la fricción del aire. Generalmente la fricción es una función de la velocidad: a mayor velocidad, mayor fricción. | Una mejora que puede agregarse es tener en consideración la fricción del aire. Generalmente la fricción es una función de la velocidad: a mayor velocidad, mayor fricción. | ||
Línea 133: | Línea 148: | ||
Aparecería entonces una nueva fuerza: | Aparecería entonces una nueva fuerza: | ||
+ | < | ||
F(t)=E(t) - P - f(y') | F(t)=E(t) - P - f(y') | ||
+ | </ | ||
- | donde f(y') es la fuerza de fricción, que depende de la velocidad y' (derivada primera de la distancia) | + | donde **f(y')** es la fuerza de fricción, que depende de la velocidad |
Así la ecuación diferencial quedaría: | Así la ecuación diferencial quedaría: | ||
+ | < | ||
y'' | y'' | ||
+ | </ | ||
- | y se resuelve de manera similar a la anterior. Si la altura que alcanza el cohete es MUY alta, entonces la fricción también dependerá de la altura (ya que varía la densidad del aire y ese es uno de los factores que influye en la fricción, pero para bajas alturas se puede suponer constante y por lo tanto independiente de la altura que alcance el cohete. | + | y se resuelve de manera similar a la anterior.\\ Si la altura que alcanza el cohete es **MUY** alta, entonces la fricción también dependerá de la altura (ya que varía la densidad del aire y ese es uno de los factores que influye en la fricción), pero para bajas alturas se puede suponer constante y por lo tanto independiente de la altura que alcance el cohete. |
- | Hasta ahi es realmente sencillo de simular | + | Hasta este punto la simulación resulta bastante sencilla, |
+ | Se complica un poco la estimación | ||
- | Igual no creo que la influencia de la variación de masa sea demasiado importante. | + | Igual no creemos |
- | Hasta ahi la simulación es 1D... suponemos que el cohete sube derechito | + | Hasta ahí la simulación es 1D... suponemos que el cohete sube derecho |
- | Con el método de Runge-Kutta, | + | Con el método de [[http:// |
Por supuesto que las condiciones iniciales (altura y velocidad inicial) determinan la solución. | Por supuesto que las condiciones iniciales (altura y velocidad inicial) determinan la solución. | ||
- | Lo que no aclaré | + | Lo que no aclaramos, |
En nuestro caso el cambio de variable lo podríamos escribir así: | En nuestro caso el cambio de variable lo podríamos escribir así: | ||
+ | < | ||
y'=v | y'=v | ||
+ | </ | ||
- | (lo cual es obvio, pues uno lo que lee es: la derivada primera de la altura es la velocidad) | + | (la derivada primera de la altura es la velocidad) |
- | Por lo tanto y'' | + | Por lo tanto **y'' |
- | y eso da lugar al sistema de 2 EDOs de 1er orden: | + | y eso da lugar al sistema de 2 **EDOs** de 1er orden: |
+ | < | ||
y'=v | y'=v | ||
v' | v' | ||
+ | </ | ||
La solución de la primera ecuación es la altura y(t) y le corresponde en nuestro caso la condición inicial y(0)=0 | La solución de la primera ecuación es la altura y(t) y le corresponde en nuestro caso la condición inicial y(0)=0 | ||
Línea 174: | Línea 198: | ||
Podríamos haber puesto otras condiciones iniciales, pero lo normal en nuestro caso sería que el cohete arranque desde el piso con velocidad cero. | Podríamos haber puesto otras condiciones iniciales, pero lo normal en nuestro caso sería que el cohete arranque desde el piso con velocidad cero. | ||
- | Ahora que tenés | + | Ahora que tenemos |
- | Supongamos que tenés | + | Supongamos que tenemos |
- | Esto es que el término de la derecha, | + | < |
+ | y'=f(t,y) | ||
+ | </ | ||
- | Entonces, si conocés un valor de la solución a tiempo | + | Esto es que el término de la derecha, f(t,y) es una función que depende exclusivamente |
- | Para eso, hacés el siguiente cálculo (esto tomalo como viene, otro día te explico | + | Entonces, si conocemos un valor de la solución a tiempo t, aplicando la receta |
+ | Para eso, hacemos el siguiente cálculo : | ||
+ | |||
+ | < | ||
K1=h*f(t0, | K1=h*f(t0, | ||
K2=h*f(t0+h, | K2=h*f(t0+h, | ||
y1=y0+(K1+K2)/ | y1=y0+(K1+K2)/ | ||
t1=t0+h | t1=t0+h | ||
+ | </ | ||
- | ¡Listo! | + | **t0** e **y0** son las condiciones iniciales, con eso encontramos |
- | + | ||
- | OK, ¿qué hago ahora? | + | |
- | Pues muy simple, ahora conozco la solución a tiempo t1, puedo volver a aplicar la receta RK para encontrar la solución a tiempo t2=t1+h. | + | OK, ¿qué hacemos |
- | Llamale | + | Pues muy simple, ahora conocemos la solución |
- | Ahora con t2 e y2, vuelvo a aplicar la receta y calculo y3 para el tiempo t3 | + | A esto lo denominamos **y2**. |
+ | Ahora con **t2** e **y2**, volvemos a aplicar la receta y calculo **y3** para el tiempo **t3**. | ||
- | ... y así sigo, todo lo que quiera. Es realmente una papa 8) | + | ... y así seguimos, todo lo necesario. |
- | Ahora bien, la receta que te pasé es para aplicar a 1 EDO de 1er orden, pero nosotros tenemos un sistema de 2 EDOs de 1er orden. | + | Ahora bien, la receta que vimos, |
- | No pasa nada, supongamos que el sistema de EDOs es: | + | No pasa nada, supongamos que el sistema de [[http:// |
+ | < | ||
y' | y' | ||
v' | v' | ||
+ | </ | ||
- | con condciones | + | con condiciones |
- | Le aplicás | + | Le aplicamos |
+ | Para no confundir, identificaremos con **K** las Kes de la 1er ecuación y con **L** las Kes de la segunda ecuación. Así que al aplicar la receta queda así: | ||
+ | < | ||
K1=h*f(t0, | K1=h*f(t0, | ||
L1=h*g(t0, | L1=h*g(t0, | ||
Línea 220: | Línea 253: | ||
v1=v0+(L1+L2)/ | v1=v0+(L1+L2)/ | ||
t1=t0+h | t1=t0+h | ||
+ | </ | ||
- | ¡y listo! Si tenés | + | Tenemos |
+ | Claro que si tenemos** | ||
- | A esta altura, si mirás | + | A esta altura, si miramos |
+ | No pretendamos | ||
+ | Por ejemplo la parte que se lee la tabla de datos debe parecer muy oscura (es una sóla linea), pero eso es porque | ||
- | Luego que termina la parte de lectura de datos, el programa está dividido en 2 secciones: mientras | + | Luego que termina la parte de lectura de datos, el programa está dividido en 2 secciones: mientras |
+ | Vemos que en el cálculo | ||
+ | Por esto consideramos al programa | ||
- | OK, volviendo | + | Volviendo |
- | Te paso la receta de orden 4: | + | Veamos ahora la receta de orden 4: |
+ | < | ||
K1=h*f(t,y) | K1=h*f(t,y) | ||
K2=h*f(t+h/ | K2=h*f(t+h/ | ||
Línea 237: | Línea 277: | ||
t=t+h | t=t+h | ||
y=y+(K1+2*K2+2*K3+K4)/ | y=y+(K1+2*K2+2*K3+K4)/ | ||
+ | </ | ||
- | al final soy un nabo, demoro más en escribir toda esta explicación que en reescribir el programa para hacerlo con RK4 e incluir | + | ==== Inclusión de la fricción del aire en la simulación ==== |
- | Bueno, acabo de hacer la simulación incluyendo la fricción del aire. | + | La fórmula de fuerza de fricción que estamos |
- | + | ||
- | La fórmula de fuerza de fricción que estoy usando es: | + | |
+ | < | ||
f(v)=0.5*Cd*da*A*v² | f(v)=0.5*Cd*da*A*v² | ||
+ | </ | ||
- | donde | + | donde: |
- | Cd: es el coeficiente de fricción del cuerpo, que en el caso del cohete | + | **Cd:** es el coeficiente de fricción del cuerpo, que en el caso del cohete |
+ | Este coeficiente es muy variable, depende de varios factores, no sólo de la forma del cuerpo, sino también de la velocidad (a través del [[http:// | ||
+ | Acá es donde toma importancia la terminación del cohete: que tan áspera o suave sea la superficie, que elementos extraños perturben el flujo del aire, etc. | ||
- | da: es la densidad del aire, que tomé como constante y de magnitud 1.29 kg/m³ | + | **da**: es la densidad del aire, tomada |
- | A: es el área de la sección transversal del cohete, que es Pi*r², medido en m² , por lo tanto tomé r=0.025 m | + | **A**: es el área de la sección transversal del cohete, que es Pi*r², medido en m² , por lo tanto r=0.025 m |
- | v: es la velocidad del vector. | + | **v**: es la velocidad del vector. |
Es evidente que esta fuerza varía también en el tiempo, a medida que varía la velocidad del vector. | Es evidente que esta fuerza varía también en el tiempo, a medida que varía la velocidad del vector. | ||
Línea 260: | Línea 303: | ||
Así que ahora las ecuaciones quedan así: | Así que ahora las ecuaciones quedan así: | ||
+ | < | ||
y'=v | y'=v | ||
v' | v' | ||
+ | </ | ||
- | + | A continuación el link del fuente AWK para la simulación en cuestión, utilizando Runge-Kutta de orden 4 | |
[[ACEMU: | [[ACEMU: |