Reacciones con autocatálisis. Grupo D12

Trabajo realizado por estudiantes
Título	Reacciones con autocatálisis. Grupo D12
Asignatura	Ecuaciones Diferenciales
Curso	Curso 2014-15
Autores	Javier Ruiz de Galarreta López, Argimiro Martínez López, Eduardo Moyano, Alberto Rodríguez Ruiz.
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 Introducción: objetivos y metodología

El objetivo de este trabajo es estudiar las concentraciones de los reactivos y productos de una reacción química a lo largo de la evolución en el tiempo de dicha reacción. Para ello modelizaremos el proceso empleando ecuaciones diferenciales y teniendo en cuenta la Ley de Acción de Masas y el Principio de Conservación de la materia.

Al ser las ecuaciones y sistemas que se nos presentan de difícil resolución analítica, utilizaremos el programa informático Matlab para resolverlos mediante métodos numéricos como el de Euler, el del Trapecio, el de Runge-Kutta y el de Heun.

Plasmaremos las soluciones en gráficos que nos permitirán dar una interpretación físico química a los resultados.

2 Reacción I

En primer lugar consideramos una reacción química irreversible en una solución bien mezclada, que se desarrolla para un volumen y una temperatura constantes, y en la que uno de los componentes (B, en este caso), hace de efecto catalítico.

[math] A + B →_{k1} 2B [/math]

Para desarrollar nuestro análisis nos apoyaremos en dos leyes físico químicas:

Ley de acción de masas: Enunciada por Guldberg y Waage en 1864, establece que dada una reacción química reversible en equilibrio, a temperatura constante, la relación de concentraciones de los reactivos y productos tiene un valor constante.

Principio de conservación de la masa: Establece que en una reacción química ordinaria la masa permanece constante, siendo la masa consumida de los reactivos igual a la masa obtenida de los productos. Fue elaborada independientemente por Míjail Lomonósov en 1745 y Antoine Lavoisier en 1785.

2.1 Interpretación del problema y deducción de las ecuaciones diferenciales

Deducir a partir de la ley de accion de masas y del principio de conservacion de la masa que las concentraciones de A y B deben satisfacer las ecuaciones

[math]x'(t)+y'(t)=0[/math]
[math]y'(t)=k_1x(t)y(t)[/math]

.....................

2.2 Obtención del problema de valor inicial

Una vez demostrado el origen de las ecuaciones (1) y (2) en el apartado anterior, obtendremos el problema de valor inicial (PVI) mediante las manipulaciones detalladas a continuación:

[math]x'(t)+y'(t)=0\quad\quad(1)[/math]

[math]y'(t)=k_1x(t)y(t)\quad\quad(2)[/math]

Integrando (1) obtenemos (3):

[math]\int_\mathbb{D} (\frac{dx(t)}{dt} + \frac{dy(t)}{dt})\,dt =\int_\mathbb{D} 0\,dt \quad \Rightarrow \quad x(t)+y(t)=k_2\quad \Rightarrow \quad x(t)=k_2-y(t)\quad\quad(3)[/math]

Sustituyendo (3) en (2), obtenemos (4):

[math]y'(t)=k_1x(t)y(t)\quad \Rightarrow \quad y'(t)=k_1(k_2-y(t))y(t)\quad \Rightarrow \quad y'(t)=k_1k_2y(t)-k_1y^2(t)\quad\quad(4)[/math]

Conocemos las concentraciones iniciales de A y de B, por lo que podemos obtener el valor de la constante k₂:

[math][A]:\quad x(0)=1 mol/l;\quad\quad [B]:\quad y(0)=0,01 mol/l[/math]

[math]x(t)+y(t)=k_2\quad \Rightarrow \quad x(0)+y(0)=k_2\quad \Rightarrow \quad 1+0,01=k_2\quad \Rightarrow \quad k_2=1,01[/math]

Por último, conocemos también el valor de k₁, por lo que podremos plantear a continuación el correspondiente PVI:

[math]k_1=1[/math]

[math](PVI)\begin{cases}y'(t)=1,01y(t)-y^2(t) \\ \\ y(0) = 0,01 \end{cases}\quad (t,y)\in[0,10]\times\mathbb{R}[/math]

2.3 Unicidad de la solución

Intersección B con dominio

Aplicamos a continuación el Teorema de Existencia y Unicidad:

Sea [math]f(t,y)=y'(t) \Rightarrow f(t,y)=1.01y(t)-y^2(t)\quad[/math], y [math]\quad\frac{\partial }{\partial y} f(t,y) = 1.01-2y(t)[/math].

Y sea [math]B[(t_0,y_0),r],\quad r\gt0,\quad en (t_0,y_0)=(0,0.01)[/math].

[math]f(t,y)[/math] es continua en [math] \mathbb{B}\cap \mathbb{D}\quad \Rightarrow \quad \exists y(t)\quad /\quad y'(t)=1,01y(t)-y^2(t)[/math]

[math]\frac{\partial }{\partial y} f(t,y)[/math] es continua en [math] \mathbb{B}\cap \mathbb{D}\quad \Rightarrow \quad \exists ! y(t)\quad /\quad y'(t)=1,01y(t)-y^2(t)[/math]

Con lo que queda demostrado que el (PVI) tiene solución, y además esta es única.

2.4 Resolución numérica del problema de valor inicial

2.4.3 Método de Runge-Kutta

Representación gráfica y(t). Método Runge-Kutta

%Método de Runge-Kutta
clc

%Función f(t)=y'(t)
fyt=inline('1.01*y-y^2','y');

%Datos del problema
y0=0.01;
t0=0;
tN=10;
h=0.1;

%Variable independiente (t)
t=t0:h:tN;

%Variable dependiente (y)
y=zeros(1,length(t));
y(1)=y0;

for i=1:length(t)-1
    k1=fyt(y(i));
    k2=fyt(y(i)+k1/2*h);
    k3=fyt(y(i)+k2/2*h);
    k4=fyt(y(i)+k3*h);
    y(i+1)=y(i)+h/6*(k1+2*k2+2*k3+k4);
end

%Representación gráfica
hold on
plot(t,y)
legend('Sustancia B: y(t)')
xlabel('Tiempo (s)')
ylabel('Concentración (mol/l)')
grid on

2.5 Resolución en forma de sistema de ecuaciones

A continuación vamos a plantear el problema de valor inicial con un sistema de dos ecuaciones equivalente al problema resuelto anteriormente.

Estas ecuaciones serán:

[math](PVI)\begin{cases}y'(t)=k_1 x(t)y(t) \\ \\ x'(t)=-k_1 x(t)y(t) \\ \\ x(0)=1 \\ \\ y(0) = 0.01\end{cases}\quad t\in[0,10][/math]

Resolveremos el problema por el método de Euler y por el método de Runge-Kutta, plasmando las soluciones en gráficos.

2.5.1 Método de Euler

Representación gráfica de las concentraciones de A y B. Resolución del sistema por Euler

%Resolución del sistema por el método de Euler

clear all

%Definimos el vector de la variable independiente, el número de
%subintervalos y la longitud de paso

t0 = 0;
tN = 10;
h = 0.1;
t = t0:h:tN
N = (tN-t0)/h;

%Definimos los vectores de las variables dependientes

x = zeros(1, N+1);
y = zeros(1, N+1);

x(1) = 1;
y(1) = 0.01;

for i = 1:N
    
    x(i+1)=x(i)+h*(-x(i).*y(i));
    y(i+1)=y(i)+h*(y(i).*x(i));
    
end

%Representamos gráficamente las soluciones

hold on

plot(t,x,'b');
plot(t,y,'r');
grid on

legend('Sustancia A: x(t)','Sustancia B: y(t)')
xlabel('Tiempo (s)')
ylabel('Concentración (mol/l)')


hold off

2.5.2 Método de Runge-Kutta

El método de resolución numérica de Runge-Kutta, involucra un gran número de parámetros. Para simplificar la implementación de este método en Matlab, el sistema se tratará como una matriz [math]z[/math], y para la obtención de los parámetros [math]k_i[/math] utilizaremos la función inline, con la que declararemos previamente las expresiones necesarias para el cálculo de los parámetros [math]k_i[/math], ahorrándonos tener que escribir innumerables veces dichas expresiones.

Primero, expresamos el sistema en forma matricial:

[math](PVI)\begin{cases}x'(t)=-k_1 x(t)y(t) \\ \\ y'(t)=k_1 x(t)y(t) \\ \\ x(0)=1 \\ \\ y(0) = 0.01\end{cases}\quad t\in[0,10] \quad \Rightarrow \quad (PVI)\begin{cases}\vec{z}'(t)= \begin{bmatrix} x'(t) \\ y'(t) \\ \end{bmatrix} = \begin{bmatrix} -k_1x(t)y(t) \\ k_1x(t)y(t) \\ \end{bmatrix} \\ \\ \vec{z}(0)= \begin{bmatrix} x(0) \\ y(0) \\ \end{bmatrix} = \begin{bmatrix} 1 \\ 0.01 \\ \end{bmatrix} \end{cases}\quad t\in[0,10][/math]

A continuación, detallamos la función [math]\vec{F}[/math] de la que dependerán los valores de los parámetros [math]k_i[/math]:

[math]\vec{F}(x,y)=\begin{bmatrix}f(x) \\ f(y) \\ \end{bmatrix}=\begin{bmatrix}x'(t) \\ y'(t) \\ \end{bmatrix}=\vec{z}'(t)=\begin{bmatrix}-k_1x(t)y(t) \\ k_1x(t)y(t) \\ \end{bmatrix}\quad \Rightarrow \quad k_1=1\quad \Rightarrow \quad \vec{F}(x,y)=\begin{bmatrix}-x(t)y(t) \\ x(t)y(t) \\ \end{bmatrix}\quad \Rightarrow \quad \begin{cases} [\vec{k_1}] & = F(\vec{z}_i) \\ [\vec{k_2}] & = F(\vec{z}_i + {1 \over 2}k_1h) \\ [\vec{k_3}] & = F(\vec{z}_i + \frac{1}{2}k_2h) \\ [\vec{k_4}] & = F(\vec{z}_i + k_3h) \\ \end{cases} [/math]

Entonces el método RK4 para este problema está dado por la siguiente ecuación:

[math]\vec{z}_{i+1} = \vec{z}_{i} + {1 \over 6}h(\vec{k_1}+2\vec{k_2}+2\vec{k_3}+\vec{k_4})[/math]

Mediante el siguiente programa en Matlab calculamos numéricamente la solución y la representamos en un gráfico:

Representación gráfica de las concentraciones de A y B. Resolución del sistema por Runge-Kutta

%Método de Runge-Kutta
clc

%Función F=[x'(t); y'(t)]
F=inline('[-z(1)*z(2);z(1)*z(2)]','z');

%Datos del problema 
x0=1;
y0=0.01;
t0=0;
tN=10;
h=0.1;

%Matriz de variable independiente (t)
t=t0:h:tN;

%Matriz de variables dependientes z=[x;y]
z=zeros(2,length(t));
z(:,1)=[x0;y0];

for i = 1:length(t)-1
    k1=F(z(:,i))
    k2=F(z(:,i)+k1.*h/2)
    k3=F(z(:,i)+k2.*h/2)
    k4=F(z(:,i)+k3.*h)
    z(:,i+1)=z(:,i)+h/6.*(k1+2.*k2+2.*k3+k4);
end
 
%Representación gráfica
hold on
plot(t,z(1,:))
plot(t,z(2,:),'r')
legend('Sustancia A: x(t)','Sustancia B: y(t)')
xlabel('Tiempo (s)')
ylabel('Concentración (mol/l)')
grid on
%hold off

2.6 Interpretación de las gráficas

Podemos observar por los resultados obtenidos que la velocidad de aparición de B es igual a la velocidad de desaparición de A (como dicta el principio de conservación de la masa), que A acaba desapareciendo prácticamente por completo y que la concentración de B termina por estabilizarse alrededor de 1 mol/l.

También cabe destacar que cuando el tiempo transcurrido alcanza 4,5 segundos las concentraciones de ambos compuestos son las mismas.

3 Reacción II

3.1 Interpretación del problema y deducción de las ecuaciones diferenciales

En este apartado vamos a estudiar la reacción propuesta por Alfred J. Lotka en 1920, que es la siguiente: