Integración de ecuaciones cinéticas con Octave

Solapas principales

Octave permite resolver las ecuaciones de un sistema cinético, de forma muy sencilla y eficaz. Aplicaremos Octave a dos sistemas cinéticos complejos: reacciones consecutivas y reacciones oscilantes (mecanismo de Lotka)

a) Reacciones consecutivas:

$A\stackrel{k_1}{\rightarrow} B \stackrel{k_2}{\rightarrow} C$, resolvemos con Octave las siguientes ecuaciones diferenciales: \begin{eqnarray} \frac{d[A]}{dt}& = & -k_1[A]\\ \frac{d[B]}{dt}& = & k_1[A]-k_2[B]\\ \frac{d[C]}{dt}& = & k_2[B] \end{eqnarray}

Creamos una función que proporcione la velocidad de los tres componentes en función de las concentraciones:

octave> function xdot = f(x,t)  # x representa las concentraciones de A, B y C, t es el tiempo en un instante dado y xdot son las derivadas de [A], [B] y [C] respecto al tiempo.

> global k  # k es un vector que representa las constantes cinéticas $k_1$ y $k_2$ se define de forma global para que sea visible en la función anterior.

> xdot(1) = -k(1)*x(1);  # define la primera ecuación

> xdot(2) = k(1)*x(1) - k(2)*x(2);  # define la segunda ecuación

> xdot(3) = k(2)*x(2);  # define la tercera ecuación.

> endfunction

Ahora pasamos a indicar los valores de constantes de velocidad y concentraciones inicales:

octave> global k; k(1) = 1.0; k(2) = 1.0;   # define k(1) y k(2)

octave> conc0(1) = 5.0;  # define la concentración inicial de A

octave> conc0(2) = 0.0;  # define la concentración inicial de B

octave> conc0(3) = 0.0;  # define la concentración inicial de C

A continuación, indicamos los tiempos entre los que vamos a realizar la integración y el número de etapas que debe realizar octave durante la integración.

octave> t_inicial = 0.0;

octave> t_final = 10.0;

octave> n_pasos = 200;

octave> t = linspace(t_inicial, t_final, n_pasos);  #La integración de las ecuaciones diferenciales se realiza llamando a la rutina lsode incluida en octave.

octave> conc = lsode("f", conc0, t);  #La salida es la matriz conc, cuyas columnas son concetraciones de los tres componentes y el tiempo. Sin embargo, resulta más visual una grafica de concentraciones frente a tiempos. Escribimos el código que hace la gráfica:

octave> gset xlabel "tiempo";

octave> gset ylabel "concentraciones";

octave> plot (t, conc);

Mecanismo de Lotka

$A+X\stackrel{k_1}{\rightarrow} 2X; \;\;\;\; X+Y\stackrel{k_2}{\rightarrow} 2Y; \;\;\;\; Y\stackrel{k_3}{\rightarrow} B$

Las ecuaciones que nos dan las variaciones en el tiempo de los distintos reactantes son: \begin{eqnarray} \frac{d[A]}{dt} & = & 0.0\\ \frac{d[X]}{dt} & = & k_1[A][X]-k_2[X][Y]\\ \frac{d[Y]}{dt} & = & k_2[X][Y]-k_3[Y]\\ \frac{d[B]}{dt} & = & k_3[Y] \end{eqnarray}

Resolvemos estas ecuaciones numéricamente utilizando octave.

octave> function xdot = f(x,t)  #  x representa las concentraciones de A, X , Y y B, t es el tiempo en un instante dado y xdot son las derivadas de [A], [X], [Y] y [B] respecto al tiempo.

> global k

> xdot(1) = 0.0;   # define la primera función

> xdot(2) = k(1)*x(1) *x(2) - k(2)*x(2)*x(3); # define la segunda función.

> xdot(3) = k(2)*x(2)*x(3) - k(3)*x(3);  # define la tercera ecuación.

> xdot(4) = k(3)*x(3);  # define la cuarta ecuación.

> endfunction

octave> global k; k(1) = 0.05; k(2) = 0.1; k(3) =0.05;   #$ define k(1) y k(2)

octave> conc0(1) = 5.0;   # define la concentración inicial de A

octave> conc0(2) = 5.0e-4;   # define la concentración inicial de X

octave> conc0(3) = 1.0e-5;   #define la concentración inicial de Y

octave> conc0(3) = 0.0;   # define la concentración inicial de B

octave> t_inicial = 0.0;

octave> t_final = 650.0;

octave> n_pasos = 10000;

octave> t = linspace(t_inicial, t_final, n_pasos);

octave> conc = lsode("f", conc0, t);

octave> plot (t, conc);