Utilización de MATLAB en la solución de ecuaciones diferenciales con una introducción a la modelizacón matematica Muchas de las ecuaciones diferenciales con las cuales nos encontramos no son manejables en forma exacta como se aprende en los cursos sobre ODE. Una herramienta en el análisis de estas ecuaciones es el uso de paquetes de solución numérica de ecuaciones diferenciales. Otras metodologías proporcionan un analisis cualitativo de las soluciones y asi obtenemos alguna información acerca de las soluciones aunque no exista una forma explicita para la misma..Veremos como Matlab puede ser utilizado como una herramienta para encontrar una solución aproximada y así poder dibujarla.

Ecuación diferencial escalar. Primeramente estudiaremos el caso mas simple con una sola variable, la cual tiene la forma general: dz = F (t , z ) dt z (t0 ) = z0 Pensamos a la variable t como el tiempo y z como la variable que describe el estado del sistema. En general uno trata con sistemas ¨autónomos¨ (independientes del tiempo), pero aquí es donde matlab se hace fuerte pudiendo tratar con sistemas mas generales en los cuales el miembro derecho dependa del tiempo y de las variables de estado del sistema. Veamos Matlab puede ser usado para resolver el modelo de Verhulst dN N = rN (1 − ) dt K

(1.1)

N (t0 ) = N 0 La variable de estado es escrita como N en vez de z. Los parámetros r y K son contantes que representan el crecimiento sin obtaculos (por ejemplo sin limitaciones en los recursos o poblaciones pequeñnas o en estado insipiente) y la capacidad de transporte. Supongamos por un momento de cómo se resuelve la ecuación de Verhults por separación de variables. He aquí como Matlab puede evaluar aproximadamente y dibujar sin conocer en forma explícita la formula para la solución exacta. Primero escribimos una Mfile que compute F, el lado derecho de la ecuación diferencial. Para el caso del modelo de Verhults

F = F (t , N ) = rN (1 −

N ) K

A continuación se reproduce el código de la función que resuelve el modelo de Verhulst

function F=verhulst(t,N) % esta funcion devuelve la velocidad de cambio de la poblacion para el % modelo de verhulst % el usuario debe ajustar parametros en la seccion de USER PARAMETERS % % dN/dt=F(N)=r N (1 - N/K) % % Input variables % t=time % N=current population size % % Input parameters (externos): % r=unhitered per capita growth rate of population % N=current population size % output variables % F=velocidad de cambio de la poblacion (dN/dt) % parametros del usuario r=0.55; K=665; % computo principal F=r*(1-N/K)*N; return Hay que hacer notar que aunque F no depende explícitamente del tiempo se ha incluido como argumento en la entrada de la primera linea. Esto es necesario para que el ODE solver pueda trabajar apropiadamente El primer argumento de la función debe ser la variable tiempo, y el segundo argumento de la funcion debe ser la variable de estado, la cual en este caso es llamada N. Matlab por supuesto requiere números específicos para todos los parámetros de entrada, de modo que se deben proporcionar valores para r y K . Es conveniente recolectar todas las constantes de entrada en un lugar como por ejemplo se hizo en la seccion ¨%User parameters¨ de forma que el usuario pueda cambiar los valores con facilidad. Ahora supondremos que queremos resolver el modelo de verhults (1.1) sobre el intervalo de tiempo 0< t < 18 horas (se debe estar alerta hacerca de la consistencia de las unidades tal como horas para t, pero por supuesto que a Matlab no le interesa el sistema de unidades que uno use. Pero si es importante que los cálculos sean consistentes y asi poder interpretar los resultados en termino de resultados del mundo real. Supongamos que en t=0 la población de tamaño N0=9.6 (medida con repecto unidad biomasica del sistema). A continuación se lista el programa que puede llevar a cabo tal tarea en modo comando es >>tspan=[0 18]; >>N0=9.6 >>[t,N]=ode45('verhulst',tspan,N0); >>plot(t,N) >>title('matlab solution of Verhulst model') >>xlabel('t (hours')

>>ylabel('N (biomass)') Las variables t y N son realmente vectores columna que guardan los valores de t y N en cada paso de tiempo que Matlab utiliza para resolver la ecuación diferencial. La grafica de la solucion de la ecuación se muestra en la figura 1. Notemos que la solucion numérica aproximada tiene la misma forma que la solución exacta

Figura 1

Sistemas de ecuaciones diferenciales En muchas ocasiones necesitados describir la evolución de un estado el cual requiere de 2 o mas variables para poderlo describir. Si la velocidad de cambio del sistema puede ser expresada como una funcion instantanea del tiempoy del estado del sistema , entonces se puede modelar matemáticamente como una ecuación diferencial para dicho vector:

dz = F (t , z ) dt

(1.2)

z (t0 ) = z (0)

Donde z es un vector que describe el estado y F es una funcion vectorial del tiempo En componentes esto puede ser escrito como sigue

dz1 = F1 (t , z1 , z2 ,......, zn ) dt dz2 = F2 (t , z1 , z2 ,......, zn ) dt .................................... .................................... dzn = Fn (t , z1 , z2 ,......, zn ) dt Con condiciones iniciales z1 (t0 ) = z1(0) z2 (t0 ) = z2 (0) ............................. zn (t0 ) = zn (0)

Suponemos aquí que el estado requiere de n variables z1 , z2 ,….., zn para su descripción. Como ejemplo se considerará el caso de dos especies que compiten . El estado del sistema está descrito por dos variables x e y las cuales representas las poblaciones de cada una de las especies. El sistema de ecuaciones diferenciales que describen el modelo es: dx = (a − by ) x, dt (1.3) dy = (m − nx) y, dt Donde a,b,m y n son parámetros constantes que necesitan ser especificados externamente al modelo y debemos aportar condiciones iniciales de la forma x(t0 ) = x0 y (t0 ) = y0

(1.4)

No tenemos ninguna solución explicita disponible para este sistema de ecuaciones diferenciales. Pero podemos preguntarle a MATLAB si nos puede proporcionar una solución aproximada con el paquete ODE paquete ode45. Conceptuallmente es mas de lo mismo que en el caso de una sola variable. Pero debemos ser cuidadosos con la notación. El ODE solver de MATLAB espera que se describa el estado del sistema con un vector. Por lo tanto para resolver 2.3 con MATLAB debemos primero describir el estado del sistema con un vector de estado con componentes igual a las variables de estado

⎡z ⎤ z = ⎢ 1⎥, ⎣ z2 ⎦ z1 = x, z2 = y ,

Reescribiendo el sistema de ecuaciones diferenciales 2.3 en término de este vector de estado tenemos: ⎡ (a − bz2 ) z1 ⎤ dz = F (t , z ) = ⎢ ⎥, − ( ) m nz z dt ⎣ 1 2⎦ ⎡x ⎤ z (t0 ) = ⎢ 0 ⎥ , ⎣ y0 ⎦ La cual es de la forma 2-2. Ahora podemos escribir una M-file que calcule el lado derecho del sistema de ecuaciones diferenciales en la forma de un vector F(t,z).

function F=compspec(t,z) % evalua la velocidad de cambio de una poblacion de un modelo simple de % especies en competencia % El usuario debe ajustar los parámetros en la seccion user parameters % dx/dt = a x +b x y % dy/dt = m y -m x y % % Variables de entrada % t= tiempo % z= vector de estado con componentes z(1)=x, z(2)=y, % con x e y represent5ando la poblacion de las especies % Parametros de entrada (externos) % Parametros de entrada % a,m= velocidad de crecimiento aistado per capita de cada especie % parametros de interaccion % variables de salida % F=vector columna con F(1)=dx/dt y F(2)=dy/dt %USER PARAMETERS a=3; b=5/2; m=2; n=1; % traduccion del vector de estado a variables mas convenientes x=z(1); y=z(2); % ajusta de F como un vector columna con dos entradas F=zeros(2,1); % Computo principal F(1)=a*x-b*x*y; f(2)=m*y-n*x*y; return Ahora supondremos que queremos resolver la evolución de la población de especies en competencia sobre el intervalo de tiempo 0=< t =< 2.5, comenzando con poblaciones iniciales x0=0.2 e y0=0.18 (imaginando que el tamaño de la población digamos que está en unidades de millar digamos)

Entonces aquí es como MATLAB puede ser instruido para realizar esta tarea tspan=[0 2.5]; z0=[0.2 0.18]; [t z]=ode45('compspec',tspan,z0); La variable t es nuevamente un vector columna que guarda los valores de t que utilize MATLAB para discretizar el tiempo. La variable z es ahora una matriz. La primera coluna se refiere a los valores de z1 en varios tiempos, y la segunda columna se refiere a los valores de z2 en los distintos tiempos. Hay un número de graficos que pueden utilizar esta información. Quizas estemos interesados en dibujar x versus t fig.2 plot(t,z(:,1)) title('poblacion de la primera especie vs tiempo') xlabel('t') ylabel('x')

Fig.2

O quizas deseemos graficar y versus t (fig.3). Para lo cual corremos el script siguiente tspan=[0 2.5]; z0=[0.2 0.18]; [t z]=ode45('compspec',tspan,z0); plot(t,z(:,2)) title('poblacion de la segunda especie vs tiempo') xlabel('t') ylabel('y')

FIG.3

Y tambien podemos plotear todas las variables de estado simultáneamente como una funcion del tiempo(fig.4) mediante el script. tspan=[0 2.5]; z0=[0.2 0.18]; [t z]=ode45('compspec',tspan,z0); plot(t,z) title('poblacion de la segunda especie vs tiempo') xlabel('t') ylabel('y')

FIG.4

Finalmente para el sistema de dos variables , es frecuentemente interesante haccer un dibujo del espacio de faces de una variable en funcion de la otra(fig.5) cuyo strip es: tspan=[0 2.5]; z0=[0.2 0.18]; [t z]=ode45('compspec',tspan,z0); plot(z(:,1),z(:,2)) title('Ploteo del espacio de estados') xlabel('x') ylabel('y')

Fig.5 : Solucion del modelo de especies en competencia (2.3). Parámetros a=3,b=5/2,m=2,n=1. Condiciones iniciales x0=0.2,y0=0.18