acemu:ensayos:jornadas:03_2010:runge-kutta
Diferencias
Muestra las diferencias entre dos versiones de la página.
| Ambos lados, revisión anteriorRevisión previaPróxima revisión | Revisión previa | ||
| acemu:ensayos:jornadas:03_2010:runge-kutta [2012/06/18 16:57] – luis | acemu:ensayos:jornadas:03_2010:runge-kutta [2012/06/19 16:11] (actual) – [Inclusión de la fricción del aire en la simulación] luis | ||
|---|---|---|---|
| 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] | ||
| Línea 101: | Línea 101: | ||
| </ | </ | ||
| - | En cuanto a la ecuación no hay misterio alguno, es física de liceo, la vieja ley de Newton: F=ma | + | ==== La Ecuación ==== |
| - | En este caso F(t)=E(t)-P(t), | + | En cuanto a la ecuació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 que hice fue suponer que P(t)=constante, | + | 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). |
| - | Del otro lado de la ecuación tenés | + | 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, |
| - | La aceleración, | + | < |
| + | F(t)=E(t)-P | ||
| + | </ | ||
| - | Por eso es que podés despejar | + | Del otro lado de la ecuación tenemos la masa del cohete multiplicada por la aceleración que sufre el mismo. |
| - | y(t)'' | + | La aceleración, |
| - | o sea | + | < |
| + | (y'' | ||
| + | </ | ||
| + | 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 132: | 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 173: | 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é | + | OK, ¿qué |
| - | Pues muy simple, ahora conozco | + | Pues muy simple, ahora conocemos |
| - | Llamale a eso y2. | + | A esto lo denominamos **y2**. |
| - | Ahora con t2 e y2, vuelvo | + | Ahora con **t2** e **y2**, volvemos |
| + | ... y así seguimos, todo lo necesario. | ||
| - | ... y así sigo, todo lo que quiera. Es realmente una papa 8) | + | Ahora bien, la receta que vimos, es para aplicar a un [[http://es.wikipedia.org/ |
| - | 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. | + | No pasa nada, supongamos que el sistema de [[http:// |
| - | + | ||
| - | No pasa nada, supongamos que el sistema de EDOs es: | + | |
| + | < | ||
| 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 219: | 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 236: | 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 259: | 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: | ||
acemu/ensayos/jornadas/03_2010/runge-kutta.1340063839.txt.gz · Última modificación: por luis
