Diferencias

Muestra las diferencias entre dos versiones de la página.

--- acemu:ensayos:jornadas:03_2010:runge-kutta [2012/06/18 17:39] – 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:Ensayos:Jornadas:03_2010|Retorno a página anterior]]
-====== Simulación por el Método Runge-Kutta ======
+====== Simulación por el Método Runge-Kutta - EN REVISION ======
 ==== Introducción ====
-En este artículo, intentamos trasmitir en forma de artículo, un posteado en el Foro de A.C.E.M.U., realizado por el maestro Kenneth Irving, sobre las simulaciones antes vistas, pero resueltas a través del Método [[http://es.wikipedia.org/wiki/Runge-Kutta|Runge-Kutta]] de resolución de ecuaciones diferenciales.\\
+En esta sección, intentamos trasmitir algunas ideas volcadas en el foro de A.C.E.M.U., sobre simulaciones de vuelo resueltas a través del Método [[http://es.wikipedia.org/wiki/Runge-Kutta|Runge-Kutta]] de resolución de ecuaciones diferenciales.\\
-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 "normal"), y llegué a unos número interesantes en cuanto a predicción de velocidad y altura del vector.\\
+Trabajando con las curvas de empuje del motor (ACEMU MX#001), se obtuvieron unos datos interesantes en cuanto  a predicción de las posibles velocidades y alturas del vector.\\
 \\
-Para empezar, aclaremos cuales fueron las aproximaciones que usé:\\
+Para comenzar, aclaremos cuales fueron las aproximaciones utilizadas: \\
 \\
-  - 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í dos etapas de vuelo:
+  - 2° se asumieron dos etapas de vuelo:
      * Etapa 1: con el motor haciendo empuje sobre el vector.
      * Etapa 2: vuelo libre (sólo actúa la fuerza de la gravedad)
-También asumí que no hay rozamiento con el aire, o sea que la ecuación que usé es:
+También se asumió que no hay rozamiento con el aire, o sea que la ecuación utilizada fué:
 <file>
@@ Línea 33: / Línea 32: @@
   m es la masa del vector
 \\
-Asumí que la masa del cohete no varía con el tiempo y la fijé en 1.8 kg, por lo tanto  el peso es P=m*9.8\\
+Se asumió que la masa del cohete no varía con el tiempo, la cual se fijó en 1.8 kg, por lo tanto  el peso es P=m*9.8\\
 \\
-Como E(t) lo tenemos bajo la forma de una tabla de valores, no tenemos otra opción que resolver numéricamente el problema, para lo cual usé [[http://es.wikipedia.org/wiki/Runge-Kutta|Runge-Kutta]] de orden 2. Todo fue programado empleando el lenguaje AWK [[http://es.wikipedia.org/wiki/Awk|Referencia al lenguaje AWK en wikipedia]].\\
+Como E(t) se encuentra bajo la forma de una tabla de valores, no cabría otra opción que resolver numéricamente el problema, para lo cual se utilizó [[http://es.wikipedia.org/wiki/Runge-Kutta|Runge-Kutta]] de orden 2.\\
+Todo fue programado empleando el lenguaje AWK [[http://es.wikipedia.org/wiki/Awk|Referencia al lenguaje AWK en wikipedia]].\\
 \\
-Asumí que el cohete lleva paracaídas... pero que no abre.\\
+También se asumió que lleva paracaídas, pero no se consideró su apertura.\\
 \\
-== Resultados: ==
+==== Resultados: ====
 //Tiempo de vuelo total//: aprox 20 s\\
-//Velocidad máxima alcanzada//: 90 m/s a los 1.3 segundos de vuelo (en ese momento se encuentra a unos 65 m de altura)\\
+//Velocidad máxima alcanzada//: 90 m/s a los 1.3 segundos de vuelo (en ese momento se encuentra a unos 65 mts. de altura)\\
 //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 a introducir: ==
+==== Como mejoras a introducir a este primer análisis: ====
 \\
-) Usar RK de orden 4\\
+) Usar Runge-Kutta de orden 4\\
 ) Modificar la ecuación para introducir fricción con el aire\\
 ) Incluir la apertura del paracaídas en la simulación\\
-) Arrancar la simulación cuando el **Empuje** supera el **Peso** del cohete, así no hay que andar recortando a mano la tabla de datos de empuje, y evitamos tonterías como que al principio el cohete "cae" (obviamente la simulación no tiene en cuenta que el cohete está apoyado)\\
+) Arrancar la simulación cuando el **Empuje** supera el **Peso** del cohete, a los efectos de no recortar a mano la tabla de datos de empuje, y evitar consideraciones tales como que al principio el cohete "cae" (obviamente la simulación no tiene en cuenta que el cohete está apoyado)\\
 \\
-== El fuente del programa AWK: ==
+==== El fuente del programa AWK ====
 <code>
 Código: [Seleccionar]
@@ Línea 101: / Línea 101: @@
 </code>
-== La Ecuación ==
+==== La Ecuación ====
-En cuanto a la ecuación no hay misterio alguno, es física de liceo, la vieja [[http://es.wikipedia.org/wiki/Leyes_de_newton|Ley de Newton]]: F=ma
+En cuanto a la ecuación la misma se basa en la [[http://es.wikipedia.org/wiki/Leyes_de_newton|Ley de Newton]]: **F=ma**
-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 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, por eso
+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, por eso
 <code>
@@ Línea 113: / Línea 113: @@
 </code>
-Del otro lado de la ecuación tenés la masa del cohete multiplicada por la aceleración que sufre el cohete.
+Del otro lado de la ecuación tenemos la masa del cohete multiplicada por la aceleración que sufre el mismo.
 La aceleración, no es otra cosa que la derivada segunda de la distancia, (en nuestro caso la altura 'y' del cohete), y que suele representarse con dos apóstrofes
@@ Línea 123: / Línea 123: @@
 para indicar que es la derivada segunda.
-Por eso es que podés despejar la aceleración que sufre el cohete en la ecuación de la siguiente manera:
+Por eso es que podemos despejar la aceleración que sufre el cohete en la ecuación de la siguiente manera:
 <code>
@@ Línea 140: / Línea 140: @@
 </code>
-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 resolviendo numéricamente.
-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 valían cero.
 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 162: / Línea 162: @@
 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 aqui es realmente sencillo de simular, porque los datos o los tenés, o te los va dando la simulación (como la altura y la velocidad en cada instante).\\ Se me complica un poco la estimación del ritmo de variación de masa del motor, para tenerlo en cuenta de manera de mejorar la simulación, pero ya le voy a encontrar una vuelta.\\ Seguramente se podrá estimar la pérdida de masa en cada instante, mediante el cociente de la fuerza de empuje en el instante t, dividido la integral de la fuerza de empuje, multiplicado por la masa de combustible.\\ Eso me daría la fracción de masa perdido en ese delta t..., suponiendo que la fuerza de empuje es proporcional a la masa eyectada, lo cual no es estrictamente cierto, porque también es proporcional a la velocidad de eyección de esa masa, pero más o menos sirve).
+Hasta este punto la simulación resulta bastante sencilla, porque los datos o los tenemos, o los va dando la simulación (como la altura y la velocidad en cada instante).\\
+Se complica un poco la estimación del ritmo de variación de masa del motor, para tenerlo en cuenta de manera de mejorar la simulación.\\ Seguramente se podrá estimar la pérdida de masa en cada instante, mediante el cociente de la fuerza de empuje en el instante t, dividido la integral de la fuerza de empuje, multiplicado por la masa de combustible.\\ Eso daría la fracción de masa perdida en ese delta t..., suponiendo que la fuerza de empuje es proporcional a la masa eyectada, lo cual no es estrictamente cierto, porque también es proporcional a la velocidad de eyección de esa masa, pero a los efectos de este análisis lo tomaremos como un dato válido.
-Igual no creo que la influencia de la variación de masa sea demasiado importante. Calculo que el término de fricción debe ser más importante y es realmente una pavada de incluir.
+Igual no creemos que la influencia de la variación de masa sea demasiado importante. Creemos que el término de de fricción debe ser más importante y sería muy sencillo incluir.
-Hasta ahí la simulación es 1D... suponemos que el cohete sube derechito y baja derechito (y de culete)... habría que, al menos llegar a una simulación 2D, que debe incluir el efecto del viento, o por lo menos un ángulo de despegue.
+Hasta ahí la simulación es 1D... suponemos que el cohete sube derecho y baja derecho (y de cola),habría que, al menos llegar a una simulación 2D, que debe incluir el efecto del viento, o por lo menos un ángulo de despegue.
 Con el método de [[http://es.wikipedia.org/wiki/Runge-Kutta|Runge-Kutta]], en lugar de resolver la ecuación uno lo que obtiene, son aproximaciones numéricas a la solución de la ecuación para ciertos valores de la variable independiente, en este caso el tiempo.
@@ Línea 172: / Línea 173: @@
 Por supuesto que las condiciones iniciales (altura y velocidad inicial) determinan la solución.
-Lo que no aclaré es que al resolver la ecuación diferencial de segundo orden, la solución es la altura y(t).\\ Como el método de [[http://es.wikipedia.org/wiki/Runge-Kutta|Runge-Kutta]] se puede aplicar solamente a EDOs ([[http://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial_ordinaria|Ecuaciones Diferenciales Ordinarias]]) de primer orden, lo que tenés que hacer para aplicarla a una de segundo orden, es aplicar un cambio de variable y transformar la EDO de 2do orden en un sistema de 2 EDOs de 1er orden.
+Lo que no aclaramos, es que al resolver la ecuación diferencial de segundo orden, la solución es la altura y(t).\\ Como el método de [[http://es.wikipedia.org/wiki/Runge-Kutta|Runge-Kutta]] se puede aplicar solamente a EDOs ([[http://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial_ordinaria|Ecuaciones Diferenciales Ordinarias]]) de primer orden, lo que tenemos que hacer para aplicarla a una de segundo orden, es aplicar un cambio de variable y transformar la EDO de 2do orden en un sistema de 2 EDOs de 1er orden.
-----
 En nuestro caso el cambio de variable lo podríamos escribir así:
+<code>
 y'=v
+</code>
-(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''=v' (la aceleración es la derivada primera de la velocidad)
+Por lo tanto **y''=v'** (la aceleración es la derivada primera de la velocidad)
-y eso da lugar al sistema de 2 EDOs de 1er orden:
+y eso da lugar al sistema de 2 **EDOs** de 1er orden:
+<code>
 y'=v
 v'=E(t)-P
+</code>
 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 193: / 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.
-Ahora que tenés eso, le aplicás la "receta Runge-Kutta" que es MUY simple.
+Ahora que tenemos eso, le aplicamos la "receta [[http://es.wikipedia.org/wiki/Runge-Kutta|Runge-Kutta]]" que es muy simple.
-Supongamos que tenés una Ecuacion Ddiferencial Ordinaria genérica,  expresada de la siguiente manera: y'=f(t,y)
+Supongamos que tenemos una [[http://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial_ordinaria|Ecuación Diferencial Ordinaria]] genérica,  expresada de la siguiente manera:
-Esto es que el término de la derecha, f(t,y) es una función que depende exclusivamente de la variable independiente (el tiempo) y la solución de la ecuación diferencial, y(t) que es lo que queremos encontrar.
+<code>
+y'=f(t,y)
+</code>
-Entonces, si conocés un valor de la solución a tiempo t, aplicando la receta de RK podés encontrar una aproximación a la solución a tiempo t+h, donde h es el incremente, en nuestro caso 1/60 s.
+Esto es que el término de la derecha, f(t,y) es una función que depende exclusivamente de la variable independiente (el tiempo), y la solución de la ecuación diferencial, y(t) que es lo que queremos encontrar.
-Para eso, hacés el siguiente cálculo (esto tomalo como viene, otro día te explico de donde sale):
+Entonces, si conocemos un valor de la solución a tiempo t, aplicando la receta de [[http://es.wikipedia.org/wiki/Runge-Kutta|Runge-Kutta]], podemos encontrar una aproximación a la solución a tiempo t+h, donde h es el incremente, en nuestro caso 1/60 s.
+Para eso, hacemos el siguiente cálculo :
+<code>
 K1=h*f(t0,y0)
 K2=h*f(t0+h,y0+K1)
 y1=y0+(K1+K2)/2
 t1=t0+h
+</code>
-¡Listo! t0 e y0 son las condiciones iniciales, con eso encontraste un valor aproximado a la solución para t1=t0+h, y a ese valor le llamás y1.
+**t0** e **y0** son las condiciones iniciales, con eso encontramos un valor aproximado a la solución para **t1=t0+h**, y a ese valor le llamas **y1**.
-OK, ¿qué hago ahora?
+OK, ¿qué hacemos 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.
+Pues muy simple, ahora conocemos la solución a tiempo **t1**, podemos volver a aplicar la receta [[http://es.wikipedia.org/wiki/Runge-Kutta|Runge-Kutta]] para encontrar la solución a tiempo **t2=t1+h**.
-Llamale a eso y2.
+A esto lo denominamos **y2**.
-Ahora con t2 e y2, vuelvo a aplicar la receta y calculo y3 para el tiempo t3
+Ahora con **t2** e **y2**, volvemos a aplicar la receta y calculo **y3** para el tiempo **t3**.
+... 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/wiki/Ecuaci%C3%B3n_diferencial_ordinaria|EDO]] de 1er orden, pero nosotros tenemos un sistema de 2 [[http://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial_ordinaria|EDO]]s de 1er orden.
-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://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial_ordinaria|EDO]]s es:
-No pasa nada, supongamos que el sistema de EDOs es:
+<code>
 y'=f(t,y,v)
 v'=g(t,y,v)
+</code>
-con condciones iniciales t0, y0 y v0 (fijate que a cada ecuación debe corresponderle una condición inicial de la correspondiente solución).
+con condiciones iniciales **t0**, **y0** y **v0** (veamos que a cada ecuación debe corresponderle una condición inicial de la correspondiente solución).
-Le aplicás la receta a CADA ecuación en ORDEN. 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í:
+Le aplicamos la receta a **CADA** ecuación en **ORDEN**.\\
+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í:
+<code>
 K1=h*f(t0,y0,v0)
 L1=h*g(t0,y0,v0)
@@ Línea 239: / Línea 253: @@
 v1=v0+(L1+L2)/2
 t1=t0+h
+</code>
-¡y listo! Si tenés un sistema de N EDOs... aplicás la receta N veces, una a cada EDO... una verdarera chotada. Claro que si tenés N EDOs lo más probable es que te convenga trabajar matricialmente, porque te podés llegar a hacer un entrevero lindo  :-X
+Tenemos un sistema de ** N** [[http://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial_ordinaria|EDO]]s; aplicamos la receta **N** veces, una a cada [[http://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial_ordinaria|EDO]], realmente sencillo.\\
+ Claro que si tenemos** N** [[http://es.wikipedia.org/wiki/Ecuaci%C3%B3n_diferencial_ordinaria|EDO]]s lo más probable es que convenga trabajar matricialmente.
-A esta altura, si mirás el programa de AWK que escribí, vas a empezar a reconocer la "receta de RK" por ahi. No pretendas entender al detalle el funcionamiento del programa, porque si sabés algo de programación te va a aprecer que le faltan cosas. Por ejemplo la parte que se lee la tabla de datos debe parecer muy oscura (es una sóla linea), pero eso es porque AWK ya tiene implícito la lectura de datos: los tipos que lo diseñaron sabían bien lo que hacían. AWK es un programa para procesar datos de texto, así que para qué molestarse en que el usuario tenga que programar eso: se da por defecto.
+A esta altura, si miramos el programa de [[http://es.wikipedia.org/wiki/Awk|AWK]] que vimos más arriba, reconoceremos la "**receta de Runge-Kutta**" por ahí.\\
+No pretendamos entender al detalle el funcionamiento del programa, porque para quien sabe algo de programación le va a parecer que le faltan cosas.\\
+Por ejemplo la parte que se lee la tabla de datos debe parecer muy oscura (es una sóla linea), pero eso es porque **AWK** ya tiene implícito la lectura de datos: quienes lo diseñaron sabían bien lo que hacían. **AWK** es un programa para procesar datos de texto, así que para qué molestarse en que el usuario tenga que programar eso: se da por defecto.
-Luego que termina la parte de lectura de datos, el programa está dividido en 2 secciones: mientras tengo datos de la tabla de empuje, calculo las ecuaciones considerando el empuje, cuando se me terminó la tabla de valores de empuje (el motor se apagó) arranca la segunda sección, de vuelo libre. Viste que en el cálculo desapareció el empuje E, y sólo queda el peso P (en realidad como va dividido por m, alcanzaría con poner g es lugar de P/m).  Pero viste que el programa es rechoto, igual que la receta.
+Luego que termina la parte de lectura de datos, el programa está dividido en 2 secciones: mientras tenemos datos de la tabla de empuje, se calculan las ecuaciones considerando el empuje, cuando se termina la tabla de valores de empuje (el motor se apagó), arranca la segunda sección de vuelo libre.\\
+Vemos que en el cálculo **desapareció el empuje E**, y **sólo queda el peso P** (en realidad como va dividido por m, alcanzaría con poner g es lugar de P/m).\\
+Por esto consideramos al programa tan sencillo, al igual que la receta.
-OK, volviendo al método numérico, RK de orden 2 no es tan bueno como las recetas de RK de orden superior, pero para empezar está bien.
+Volviendo al método numérico, **Runge-Kutta** de orden 2 no es tan bueno como las recetas de Runge-Kutta de orden superior, pero para comenzar está bien.
-Te paso la receta de orden 4:
+Veamos ahora la receta de orden 4:
+<code>
 K1=h*f(t,y)
 K2=h*f(t+h/2,y+K1/2)
@@ Línea 256: / Línea 277: @@
 t=t+h
 y=y+(K1+2*K2+2*K3+K4)/6
+</code>
-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 la fricción del aire
+==== 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 usando es:
-La fórmula de fuerza de fricción que estoy usando es:
+<code>
 f(v)=0.5*Cd*da*A*v²
+</code>
-donde
+donde:
-Cd: es el coeficiente de fricción del cuerpo, que en el caso del cohete consideré que valía 0.7, según una publicación de NAR. 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 número de Mach y el número de Reynolds que mide la turbulencia del fluido que atraviesa), pero para bajas velocidades (subsónicas) se puede considerar independiente de Re y Ma. 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.
+**Cd:** es el coeficiente de fricción del cuerpo, que en el caso del cohete consideramos valía 0.7, según una publicación de [[http://www.nar.org/|NAR]].\\
+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://es.wikipedia.org/wiki/N%C3%BAmero_Mach|número de Mach]] y el [[http://es.wikipedia.org/wiki/Numero_de_Reynolds|número de Reynolds]], que mide la turbulencia del fluido que atraviesa), pero para bajas velocidades (subsónicas) se puede considerar independiente de **Re** y **Ma**.\\
+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 como constante y de magnitud 1.29 kg/m³
-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.
@@ Línea 279: / Línea 303: @@
 Así que ahora las ecuaciones quedan así:
+<code>
 y'=v
 v'=(E(t)-P-f(v))/m
+</code>
+A continuación el link del fuente AWK para la simulación en cuestión, utilizando Runge-Kutta de orden 4   [[ACEMU:Ensayos:Jornadas:03_2010:software|Simulación - Programa escrito en AWK para Runge-Kutta 4]].
 [[ACEMU:Ensayos:Jornadas:03_2010|Retorno a página anterior]]