¿Cómo se puede conocer el comportamiento de ecuaciones diferenciales que representan un sistema dinámico? Una forma de hacerlo es simulando numéricamente (obteniendo la solución numérica).
Para quienes no están familiarizados con la programación en MATLAB o en GNU Octave, les presento un minicurso de simulación de ecuaciones diferenciales que preparé en 2006.
Scripts
Los scripts o guiones son archivos que se ejecutan línea a línea y que pueden a su vez llamar a funciones de MATLAB/Octave1 o funciones definidas por el usuario. Considere por ejemplo, hacer el gráfico de la función $\sin(x)$, sin(x)en el intervalo [0 6]. Podemos escribir el código necesario en el archivo primer.m:
|
1 2 3 4 5 6 7 8 9 |
clear all close all x = 0:0.1:6; plot(x,sin(x)) title('Funcion seno de x') xlabel('x [seg]') ylabel('sin(x)') grid |
Si desde la ventana de comandos de MATLAB se escribe:
>> primer
O desde GNU Octave:
octave:1> primer
se ejecuta el contenido del script primer. Al ejecutar el código anterior se genera la siguiente gráfica:
![Gráfica de la función $sin(x)$ en el intervalo $[0,2\pi]$](http://www.controldesistemasnolineales.com/blog/wp-content/uploads/2015/01/primer.png)
El siguiente es un script para MATLAB el cual permite simular la ecuación diferencial $dx/dt = -x$, con condición inicial $x(0) = -1.5$. Recuerde que ésta es una ecuación diferencial de primer orden.
MATLAB usa la función ode23 para calcular numéricamente la solución $x(t)$ de la ecuación diferencial. El tiempo de simulación es entre [0 5] s y la condición inicial es x0.
|
1 2 3 4 5 6 7 8 9 10 |
clear all close all x0 = -1.5; [t,x] = ode23('priorden',[0 5],x0); plot(t,x) title('Solucion de la ED de primer orden') xlabel('t [seg]') ylabel('x(t)') |
La ecuación diferencial es simulada mediante el comando ode23. ¿Cómo? Este comando requiere que le pasemos la ecuación diferencial $\dot x = -x$. Esto se hace en el archivo priorden.m que contiene solamente una función, la función function derx = priorden(t,x) donde está escrita la ecuación diferencial:
|
1 2 3 |
function derx = priorden(t,x) a = -1; derx = a*x; |
De nuevo, priorden es el nombre de la función o function donde se define la ecuación diferencial $\dot x = -x$ y asimismo es el nombre del archivo priorden.m. La ecuación diferencial depende de (t,x). La ecuación derx = a*x; representa la ecuación diferencial con coeficiente constante a = -1.
Al ejecutar simu se obtiene la gráfica:
La solución analítica, dada por $x(t) = -1.5e^{-t}$, escrita en MATLAB/Octave x(t) = -1.5*exp(-t) la podemos sobreescribir en la gráfica anterior usando el comando hold:
|
11 12 13 14 15 |
% se sobrescribe la solucion x(t) = -1.5*exp(-t) hold plot(t,-1.5*exp(-t),'xr') |
Las opciones 'xr' significan que va a graficar cada punto con una 'x' y en color rojo. El resultado:
En GNU Octave
El código usado en GNU es equivalente, con las siguientes salvedades
- la función
priordenposee dos particularidades: se puede incluir dentro del script, el orden de los parámetros de entrada se invierte(x,t), - ahora se usa la función
lsodeen lugar deode23para el cálculo de la solución numérica de la ecuación diferencial, - la función
lsoderequiere definir de forma precisa el vector de tiempo $t\in [0,5]$, esto se hace con la funciónlinspace, en la cual además se indica el número de puntos que se calcularán2.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 |
clear all close all function derx = priorden(x,t) a = -1; derx = a*x; endfunction x0 = -1.5; t = linspace(0,5,25); x = lsode('priorden',x0,t); plot(t,x) title('Solucion de la ED de primer orden') xlabel('t [seg]') ylabel('x(t)') % se sobrescribe la solucion x(t) = -1.5*exp(-t) hold plot(t,-1.5*exp(-t),'xr') |
O en forma breve (sin hold):
|
1 2 3 4 5 6 7 8 9 10 11 12 |
clear all close all f=@(x,t)[-x]; # ecuacion de primer orden x0 = -1.5; t = linspace(0,5,25); x = lsode(f,x0,t); plot(t,x,t,-1.5*exp(-t),'xr') title('Solucion de la ED de primer orden') xlabel('t [seg]') ylabel('x(t)') |
Recuerde que:
- MATLAB y GNU Octave están basados en vectores.
- Las funciones definidas por el usuario son estructuras básicas de programación en MATLAB/Octave.
- Para usar adecuadamente MATLAB/Octave reconozca las funciones básicas que necesita y estúdielas.
Simule en línea en GNU Octave ecuaciones diferenciales:
- Solución de la ecuación diferencial $\dot x = -x$, con $x_0 = 1$: enlace
- Otro ejemplo avanzado
- Un ejemplo de simulación de ecuaciones diferenciales no lineales: el modelo predador-presa (predator-prey model)
Gracias
Saludos