Solución Numérica de Ecuaciones Diferenciales con ode45 en Matlab, Ejemplos.

[Recuerda que en este Blog los enlaces para la descarga del código se encuentran al final del artículo.]

En muchas ocasiones, nos enfrentamos a ecuaciones diferenciales de las cuales es muy difícil o simplemente no puede obtener una solución analítica, para esto existen los métodos numéricos, y en esta ocasión explicaremos un método numérico (en realidad su función) extremadamente conocido entre la comunidad usuaria de Matlab y en general en el ámbito científico llamado ode45() el cual es una variación un optimización del método de Runge-Kutta modificado.

LA SINTAXIS.

LO QUE RETORNA.

En esta entrada nos centraremos en la sintaxis básica ya que no necesitamos de las estadísticas de rendimiento, etc… este método retorna un vector columna con los valores de la variable independiente (de ahora en adelante para nuestros ejemplos será ‘t‘) y un vector columna (que llamaremos ‘y‘) y que contiene cada uno de los valores para los que se va dando la solución en cada instante ‘t‘, esto lo ilustraremos en un momento con un ejemplo.

LO QUE RECIBE.

Esta función recibe 3 parámetros esenciales, el nombre de un archivo .m en el que se define la ecuación diferencial acompañado por un símbolo ‘@‘, el siguiente argumento, corresponde a un vector fila con dos valores entre los cuales queremos conocer la solución de la ED y como último parámetro recibe un vector de valores iniciales que contiene los valores iniciales necesarios para una correcta solución de la ED, En otras palabras el prototipo básico para usar ode45() es el siguiente:

y en este caso la solución se numérica se almacenará en el vector ‘y‘ para cada uno de los instantes de tiempo presentes en el vector ‘t‘.

¡UN EJEMPLO!

Todo es más fácil de entender con ejemplos, y para este primer ejemplo encontraremos la solución numérica para la siguientrama que la masa se encuentra en la posición X(0)=0 es decir, en su posición de equilibrio, y que lleva una velocidad dx(0)=2 e ecuación diferencial de primer orden, con condición inicial y(0)=0.

Bueno, para dar solución a esta ecuación siempre será necesario dos scripts, uno donde definimos los parámetros de ‘tiempo_de_solución‘ y ‘condiciones_iniciales‘, y pues, posteriormente hacemos el llamado a ode45( ), como el de la figura siguiente:

En este caso decidimos que queremos visualizar la solución en el intervalo de 0 a 3 para la variable independiente (en este caso x), en la linea siguiente dijimos que la condición inicial sería y(0)=0, en la linea siguiente se hace uso de ode45() sabiendo que los argumentos que retorna es la variable independiente, seguida de la variable dependiente que contiene la solución y finalmente hacemos un plot() de ambos vectores para visualizar la solución.

Cabe resaltar que en la imagen anterior está presente el símbolo @ que acompaña el nombre del archivo donde se define (escribe) nuestra ecuación diferencial, y esto se debe a que con el @ sucede como si le pasáramos un puntero o dirección de donde se encuentra el archivo que define a la ED y de está forma ode45() hará llamadas recursivas cuantas veces sea necesario, en nuestro caso el nombre del archivo es EDO1 y su contenido es el siguiente:

Como se puede ver el nombre de la función (y por ende del archivo) es EDO1 y dentro de este se define la ecuación diferencial que planteamos al principio, luego de guardados guardados ambos archivos sin errores, procedemos a ejecutar el script encargado de llamar a ode45() con sus parámetros, obteniendo la solución como se ve en la figura 1.

Fig 1. Solución de la Primera Ecuación Diferencial Propuesta.

 

¡UN EJEMPLO DE SEGUNDO ORDEN!

Para este siguiente ejemplo usaremos un conocido ejercicio de sistemas físicos que genera una ED de segundo orden a la cual le hallaremos la solución con ode45() para ver el comportamiento de su movimiento a través del tiempo, este sistema es el conocido masa-resorte de la figura 2.

Fig 2. Sistema Masa-Resorte.

Y elegimos este sistema porque tenemos la ecuación diferencial de segundo orden que modela su movimiento (se puede obtener fácilmente por análisis de fuerzas de Newton) y es la siguiente.

[Cabe recordar que el operador . (punto) también se usa para indicar derivadas se lee (prima)].

 

* Antes de continuar debemos aclarar algunas cosas, primero la variable dependiente en este caso es X  [X(t)], en otras palabras X es función de ‘t‘.

* ode45() solo resuelve ecuaciones del tipo y’=f(t,y) por lo tanto si vamos a resolver ecuaciones de orden superior (orden>1) entonces si tenemos una ecuación de orden N, deberemos llevar esta ecuación a un sistema de N ecuaciones diferenciales de primer orden, y esto se logra mediante algo conocido como reducción de orden de ecuaciones diferenciales.

Continuando con la solución a  nuestro problema propuesto, debemos entonces reducir el orden de nuestra ecuación diferencial a un sistema de dos ecuaciones diferenciales de primer orden, ¿cómo lo hacemos?.

Definimos dos variables X1 y X2, y se sigue el procedimiento a continuación:

Entonces habiendo 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.

Hay que recordar que a la hora de trabajar con ode45() 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 de la siguiente forma, y a su vez el sistema quedaría descrito como:

Teniendo en cuenta que:

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

Ahora, no escogimos 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 ode45() ya que lo que recibe esta función es un vector, 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 X donde se almacenan los valores de X1 y X2 respectivamente, lo ilustramos de la siguiente forma:

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

y será este el vector que se usará tanto en el archivo que define el sistema de ecuaciones como en el retorno de ode45() el cual recordemos que retorna [t,X].

Teniendo lo anterior en cuenta y habiendo comprendido comprendido el proceso, procedemos entonces a escribir los scripts, el primero como siempre, define el tiempo_de_solución y las condiciones_iniciales, y de la misma forma que antes, hace el llamado a ode45() y recibe su solución en los vectores [t,X] quedando dicho script de la siguiente forma:

En este caso en las condiciones iniciales le decimos al programa que la masa se encuentra en la posición X(0)=0 es decir, en su posición de equilibrio, y que lleva una velocidad dx(0)=2 (hay que aclarar que las condiciones iniciales se dan con el formato CI=[x(0) dx(0)] ).

En el siguiente script se define el sistema equivalente de ecuaciones diferenciales de primer orden y cuyo archivo lleva el mismo nombre que se usó en el llamado de ode45() que es EDO2.m en este caso, quedando como se muestra a continuación.

Para entender el manejo de índices que se hizo en esta función volvemos a visualizar cada variable:

Cuando se realiza la ejecución del código que se encarga de llamar a la función ode45() se logra el siguiente resultado, donde se aprovechó el hecho de que la función ode45() retorna la solución para X1 y X2, recuerde que X1=x lo que quiere decir que la solución de X1 (en este caso la primera columna de la matriz principal X) representará la solución de x -> posición, y X2 corresponde a x’ lo que quiere decir que la columna 2 de la matriz principal X contendrá a x’ o que quiere decir que X(:,2) contiene la velocidad.

Para ilustrar mejor lo retornado por ode45() entonces mostramos en la siguiente tabla el contenido de la matriz X:

Tabla 1. Contenido de la matriz X retornada por ode45()

Es por esta razón que al graficar la matriz X de la la siguiente forma, se ve en la figura 3 la posición y la velocidad, y que además pueden ser identificadas fácilmente ya que cumplen con las condiciones iniciales dadas.

Fig 3. Gráfica de x(t) y dx/dt, posición y velocidad retornada por ode45().

 

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 ode45() para que el nos retorne las soluciones correspondientes.

 

Hasta aquí este tutorial sobre el uso de ode45(), 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.

Autor: Julio Marulanda.