Solución Numérica de Ecuaciones Diferenciales con Ode en Scilab con Ejemplos.

En muchas ocasiones, a lo largo de nuestra carrera, nos podemos enfrentar a derivadas o sistemas de ecuaciones diferenciales, para las cuales es muy difícil o simplemente no es posible hallar una solución analítica, o en otra palabras, no es posible expresar su solución mediante funciones elementales, en este tipo de situaciones debemos recurrir a los métodos numéricos que nos permiten aproximar las soluciones, el comportamiento de dichos sistemas de ecuaciones. En esta entrada dedicaremos algo de tiempo para explicar uno de ellos (en realidad su implementación en Scilab).

Para los que anteriormente hayan trabajado con Matlab seguramente sabrán que hay diferentes funciones para este tipo de usos y que la más conocida es la famosa ode45( ), pues muchos vendrán aquí buscando su equivalente en Scilab, pues somos todos afortunados ya que Scilab tiene su equivalente a ode45( ) y es la conocida función ode( ) cuya sintaxis no dista demasiado de su competencia en Matlab.

¿La Sintaxis de Ode( )?

En esta entrada nos centraremos en la sintaxis básica de la función, es decir, dejaremos de lado las opciones de estadísticas de rendimiento, etc…y mencionaremos lo que es necesario para solucionar una ecuación diferencial de orden N (en otras entradas hablamos de como resolver sistemas de EE.DD).

La función ode( ) tiene el siguiente prototipo:

ode_scilab_prototipo

Donde:

  •  Y es el vector que contendrá las soluciones de la ecuación que deseamos solucionar, para cada instante de tiempo t.
  • Yo es un vector columna que debe contener los valores iniciales respectivos de la ecuación diferencial.
  • to es un escalar que indica en donde fueron evaluados lo valores iniciales presentes en Yo, es decir, que si por ejemplo tengo una ecuación diferencial de segundo orden, necesitaré dos valores iniciales (a saber, Y(to) y Y'(to)) y con este parámetro, Scilab nos dice que ambas condiciones iniciales deben ser tomadas en el mismo punto to.
  • t corresponde  a un vector que al igual que los parámetros anteriores debe ser definido por le usuario, y es sobre este vector t que se calculará la solución de la ecuación diferencial, es decir que si quiero observar la solución de una ecuación diferencial de 2do orden entre 0 y 4, entonces solo bastará con generar un vector t (o con cualquier otro nombre) que vaya desde 0 a 4, con el número de puntos que queramos.
  • f es el nombre de una función (archivo con extensión .sci) que contiene la definición de nuestra ecuación diferencial, y será ejemplificada con algunos ejemplos mas adelante.

Hay que tener algo muy en cuenta para evitar problemas posteriores, y es que cuando se definen los valores iniciales, estos pueden ser tomados antes de que inicie el rango que queremos visualizar, o incluso sobre este, pero no después…en otras palabras si como en el ejemplo anterior, quiero visualizar la solucion de una ecuación diferencial entre t=3 y t=6 puedo pasar como parámetros Yo con valores iniciales tomados en To=2 o To=3 por ejemplo, pero no con To mayores o iguales a 3, y tiene bastante lógica ya que si quiero hacer la visualización entre t=3 y t=6, si paso valores iniciales por ejemplo en To=4, pues ode() no podrá generar soluciones para tiempos anteriores a To=4 y Scilab nos indicará esto con un mensaje de error.

 ¡Un Ejemplo de Primer Orden!

La verdad creo que uno entiende mejor las cosas cuando ve ejemplos, y por esta razón daremos solución al siguiente ejemplo de primer orden. Resolver la siguiente ecuación diferencial y visualizar la solución para un 0=<t<=4;

edo1_scilab

Con condición inicial:

edo1_condicion_inicial

Lo que indica que:

edo1_tiempo_inicial

Teniendo lo anterior claro, entonces abrimos un nuevo script en Scilab y ponemos lo siguiente, la explicaciones se encuentra en los propios comentarios:

ode_ejemplo_1

Como se puede ver, al momento de hacer el llamado a la función ode(), uno de los parámetros es (como explicábamos antes) la función F donde se define la forma de la ecuación diferencial y su estructura es la siguiente:

funcion_ejemplo_1

Como pueden ver el nombre de la función es por simplicidad ‘f’ aunque claramente se puede escoger otro nombre.

Después de desarrollar y guardar lo anterior, solo basta con dar clic en la opción ‘ejecutar’ (el símbolo con forma de play) para ejecutar el script que llamará a la función ode().

ejecutary el resultado obtenido es el siguiente:

sol_edo_1

Fig 1. Solución con ode( ) a la ecuación diferencial dy=9t^2-4t-3sin(t) con y(0)=3.

Para comprobar la efectividad de la función ode( ), la solución analítica para este ejemplo es:

original_1

si graficamos la función original con puntos rojos, sobre la linea azul continua (que es la obtenida con ode( ) ) podremos apreciar la precisión de la función y además que no cometimos ningún error.

Fig 2. Comparación soluciones analítica y numérica.

Fig 2. Comparación soluciones analítica y numérica.

 

¡Mejor un Ejemplo de Segundo Orden!

La función ode(), al igual que la función ode45() en Matlab permite resolver ecuaciones de orden N, pero para esto hay que escribirla de forma que la máxima derivada sea de orden 1, entonces en estos casos debemos realizar una reducción de orden, en otras palabras, si queremos resolver una ecuación de orden N, debemos convertirla en un sistema de N ecuaciones de orden 1. (pero no os preocupéis, con el siguiente ejemplo ilustraremos lo anterior).

Sea el sistema Masa-Resorte que se ve a continuación:

Fig 3. Sistema Masa-Resorte.

Fig 3. Sistema Masa-Resorte.

Y cuyo movimiento a través del tiempo está determinado por la ecuación diferencial de segundo orden que se presenta a continuación, esta ecuación puede ser fácilmente obtenida mediante un análisis de fuerzas de Newton, (recordemos que la derivada se puede indicar con un ”.” y se lee ‘prima’).

 edo1

[Recordemos que en este caso, la variable dependiente es la X, entonces también la podemos indicar como X(t)].

Como mencionamos anteriormente debemos llevar esta ecuación de segundo orden, a dos ecuaciones de orden 1, en otras palabras debemos seguir un procedimiento para realizar la reducción de orden, y lo indicamos a continuación.

reduccion

Después de haber hecho este cambio de variables, si reescribimos la ecuación diferencial original (y lo debemos hacer para un mejor entendimiento) quedaría de la siguiente forma.

nueva_edo

Hay que recordar que a la hora de trabajar con ode() debemos tener una ED de la forma y’=f(t,y) por lo tanto debemos tener despejada a la primera derivada, quedando la ecuación anterior (y a su vez el sistema) de la siguiente forma:

reescrito

Teniendo en cuenta que:

prima

¿Notan como logramos reescribir una ecuación diferencial de segundo orden en un par de ecuaciones de primer orden?

No escogimos al azar las nuevas variables X1 y X2 pudiendo haber escogido U1 y U2, etc. La elección se hizo con el fin de lograr un mayor entendimiento a la hora de pasar estos parámetros a ode() ya que lo que recibe esta función es un vector columna, entonces le pasaremos un vector de nombre X que contendrá a X1 y X2 dentro de él, en otras palabras, los subíndices 1 y 2 harán referencias a las columnas del vector principal donde se almacenan los valores de X1 y X2 respectivamente, lo ilustramos de la siguiente forma:

vec_ode_2

Y para las derivadas es básicamente lo mismo, aunque en este caso se trata de un vector columna:

vector_dev

y será este vector de derivadas el que usaremos en la función que describe el sistema.

Después de haber comprendido los procesos anteriores se procede a implementar los scripts en Scilab (¿no es eso lo que queremos?). El primer paso será hacer el script general, donde definiremos las condiciones iniciales (de velocidad y posición), el tiempo inicial (en este caso To=0), el tiempo en el que queremos la solución al que llamaremos ‘t’ como en los casos anteriores, cargaremos la función xpp.sci al workspace de Scilab como siempre, y despues de hacer el respectivo llamado a la función ode() para obtener nuestras soluciones, las graficaremos en un plot() y escribiremos sus características, todo esto como se muestra a continuación.

script

El siguiente archivo es el que define el sistema de ecuaciones de primer orden que quedó luego de la reducción que se hizo en pasos anteriores, este archivo lleva el mismo nombre con el que se debe llamar en la función ode() y su extensión es .sci

archivo_funcion

con lo anterior, solo hace falta ejecutar el script en el que hacemos el llamado de la función ode(), si queremos entender como es el retorno de esta función, entonces podemos ilustrar como se almacenan las soluciones en la variable X como sigue:

Como se almacenan los resultados en la variable X.

Como se almacenan los resultados en la variable X.

en este caso vemos que al contrario de Matlab y su ode45(), en Scilab, la solución es almacenada en las filas, decir que si queremos ‘plotear’ la posición de la masa, basta entonces con tomar todos los valores de la primera fila, ya que a esta corresponde la X primitiva (la de posición) y si queremos graficar la velocidad, solo necesitamos hacer uso de todos los valores presentes en la fila 2 ya que como vimos cuando hicimos la reducción X2 era la primera derivada de X, es decir la velocidad.

Graficando la posición y la velocidad se obtiene lo siguiente:

Grafica de la posición y la velocidad, accediendo a los respectivos valores de la solución X.

Grafica de la posición y la velocidad, accediendo a los respectivos valores de la solución X.

De la misma forma en que se resolvió esta ecuación diferencial de segundo orden para el sistema Masa-Resorte, se puede resolver cualquier ecuación diferencial de orden N, obviamente llevándola por medio de “reducción de orden” a un sistema de N ecuaciones diferenciales de orden 1 y pasándola a ode() para que nos retorne las soluciones correspondientes.

Hasta aquí este tutorial sobre el uso de la función ode() en Scilab, como hacer uso de él y como interpretar los resultados que nos entrega. los scripts que desarrollamos en este articulo se encuentran aquí para su descarga, análisis y modificación.

Espero que todo esto haya sido de ayuda para ustedes, cualquier sugerencia o pregunta es bienvenida, saludos y éxitos..no olvides comentar.

Autor: Julio E. Marulanda.

Responder

Introduce tus datos o haz clic en un icono para iniciar sesión:

Logo de WordPress.com

Estás comentando usando tu cuenta de WordPress.com. Cerrar sesión / Cambiar )

Imagen de Twitter

Estás comentando usando tu cuenta de Twitter. Cerrar sesión / Cambiar )

Foto de Facebook

Estás comentando usando tu cuenta de Facebook. Cerrar sesión / Cambiar )

Google+ photo

Estás comentando usando tu cuenta de Google+. Cerrar sesión / Cambiar )

Conectando a %s