Reacción Autocatálisis A2

Trabajo realizado por estudiantes
Título	Reacción Atocatálisis. Grupo 2-A
Asignatura	Ecuaciones Diferenciales
Curso	Curso 2015-16
Autores	Ana Arcos González Cristina Sáez del Moral Nerea Sauto VegaPablo Merayo Vidal María Arias Casado Borja Rodríguez Monereo
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 Introducción

A continuación se presenta un estudio sobre las reacciones de autocatálisis, basadas en el principio de conservación de la materia y en la Ley de Acción de masas.

• Principio de conservación de la materia: En una reacción química ordinaria, la masa permanece constante, es decir, la masa consumida de los reactivos es igual a la masa obtenida de los productos.

• Ley de Acción de masas: establece que para una reacción química reversible en equilibrio a una temperatura constante, una relación determinada de concentraciones de reactivos y productos, tienen un valor constante.

El enunciado del estudio a realizar será: Se considera una reacción química irreversible en una solución bien mezclada. Supondremos que la reacción ocurre para un volumen y temperatura constantes. Al inicio se encuentran dos reactivos A y B que van formando un producto C en lo que se conoce como una reacción bimolecular, es decir, una molécula de A y una de B producen una de C. En este ejercicio analizaremos el caso particular en el que A se transforma en B pero suponiendo que la presencia de B hace de efecto catalítico en la reacción. Escribiremos este proceso como una reacción bimolecular

A+B->k12B

2 Deducción Problema de valor inicial

Deducir a partir de la ley de acción de masas y del principio de conservación de la masa que las concentraciones de A y B deben satisfacer las ecuaciones

y’(t) + x’(t) = 0

y’(t) = k1 x(t) y(t)

Integrando la primera ecuación y sustituyendo el valor de x(t) en la segunda, calcular la ecuación diferencial para y(t). Escribir el PVI asociado y decidir si tiene solución ´única.

Lo primero es denominar x (t) e y (t) como las concentraciones de los reactivos que aparecen en la ecuación de arriba. Para encontrar la ecuación en función de una única variable y (t), se parte de la ley de conservación de la masa, obteniendo como resultado

X (t)+y (t)=cte.=c

Si derivamos, obtenemos la primera ecuación que nos piden demostrar que se cumple:

x’(t)+y’(t)=0.

La segunda ecuación la podemos obtener, directamente, a partir de la Ley de Acción de masas que no dice que:

y’(t)=k1*x(t)*y(t)

Integrando y sustituyendo la variable x (t) de la primera ecuación en la segunda obtenemos la ecuación diferencial deseada:

y′ (t)=k1∗ (c−y(t))∗y(t)

Si se parte de que la concentración inicial de B en t=0 es y0, obtendremos el problema de valor inicial pedido.

3 Resolución Numérica PVI

Suponiendo que las concentraciones iniciales de A y B son 1 mol/l y 0.01 mol/l respectivamente, y k = 1 mol/s, resolver el PVI por el método de Euler, eligiendo un paso h = 0.1, en los primeros 10 segundos. Dibujar en una misma gráfica las concentraciones de las dos sustancias. Interpretar las gráficas.

• Método Euler

El método de Euler, llamado así en honor al matemático Leonhard Euler, es un procedimiento de integración numérica. Su expresión matemática es\[ PVI = \begin{cases} y_0 \\ y_{n+1} = y_n + h*f(t_n,y_n) \\ \end{cases}\quad \]

A continuación emplearemos esta expresión en un programa de Matlab para obtener una aproximación de las soluciones. El código de dicho programa es el siguiente

%DATOS DEL PROBLEMA

texto alternativo

 %DATOS DEL PROBLEMA
%x(t)= Concentración de A a lo largo del tiempo en mol/L = 1
%y(t)= Concentración de B a lo largo del tiempo en mol/L = 0.01
%k = 1 mol/s
%h = 0.1 
%En los primeros diez segundos t=[0,10]s
k=1;
h=0.1;
t0=0;
tN=10;
y0=0.01;
%Definimos el vector tiempo con h=0.1
t=t0:h:tN;
%La constante c del P.V.I. es la suma de las concentraciones como resultado del Principio de conservación de la masa.c=0.1+y0=1.01;
y=zeros(1, length(t));
%La variable y será la solución de la concentración B.
y(1)=y0;
%Definimos la función para aplicar, a continuación, el Método de Euler.
f=inline('y*(1.01-y)','t','y');
for i=1:length(t)-1
   y(i+1)=y(i)+h*f(t(i),y(i)); %Aplicamos Euler
end
%Calculamos la concentración de A a partir de las concentraciones de B y de C, como la diferencia entre ambas (siendo la concentración de B "y").
x=1.01-y;
%Dibujamos ambas gráficas
plot(t,y,'k x');
hold on
plot (t,x,'g x'); 
hold off
title('Gráfica Concentración-Tiempo');
xlabel('tiempo(s)');
ylabel('Concentración (mol/l)');
legend('Concentración de B', 'Concentración de A', 'Location', 'Best');

• Método Trapecio

Vamos a aplicar un método de integración numérica conocido como el Método del Trapecio por aproximar el valor de la integral al del área de un trapecio. Su expresión matemática es\[ PVI = \begin{cases} y_0 \\ y_{n+1} = y_n + \frac{h}{2}*(f(t_n,y_n)+f(t_{n+1},y_{n+1}))\\ \end{cases}\quad \]

Introducimos esta expresión en el código y obtenemos la gráfica de la solución.

texto alternativo

 dy=k1*(c-y)*y   x0=1   y0=0.01   k1=1   h=0.1   [0,10]   c=1.01
 % Datos del problema
 k=1;
 c=1.01;
 t0=0;
 tN=10; 
 y0=0.01;
 h=0.1; % tamaño de paso
 t=t0:h:tN; 
 N=(tN-t0)/h; % número de intervalos
 y=zeros(1,length(t)); % matriz de 1 por N+1
 y(1)=y0;
 % Trapecio
for i=1:N
y(i+1)=(1/(h*k))*((0.5*h*k*c-1)+sqrt((1-0.5*h*k*c)^2-2*h*k*(-y(i)-(h/2)*y(i)*(c-y(i)))));
end
x=1.01-y;
% Gráfico
hold on
plot(t,y,'k x')
plot(t,x,'g x')
hold off
legend('Concentración de B','Concentración de A', 'location',   'best')

• Método Runge Kutta

Resolver el PVI numéricamente usando el método del trapecio y el método de Runge Kutta de cuarto orden con el mismo paso de tiempo. 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, utilizaremos la función inline, con la que declararemos previamente la función algebraica necesaria para el cálculo de los parámetros \(k_i\), ahorrándonos tener que escribir innumerables veces la correspondiente expresión matemática.

texto alternativo

 % dy=k1*(c-y)*y   x0=1   y0=0.01   k1=1   h=0.1   [0,10]   c=1.01
 % Datos del problema
 k=1;
 c=1.01;
 t0=0;
 tN=10; 
 y0=0.01;
 h=0.1; % tamaño de paso
 t=t0:h:tN; 
 N=(tN-t0)/h; % número de intervalos
 y=zeros(1,length(t)); % matriz de 1 por N+1
 y(1)=y0;
 % RK4
 for i=1:N
    K1=k*(c-y(i))*y(i);
    K2=k*(c-y(i)+0.5*K1*h)*(y(i)+0.5*K1*h);
    K3=k*(c-y(i)+0.5*K2*h)*(y(i)+0.5*K2*h);
    K4=k*(c-y(i)+K3*h)*(y(i)+K3*h);
    y(i+1)=y(i)+h/6*(K1+2*K2+2*K3+K4); 
 end
 x=1.01-y;
 % Gráfica
 hold on
 plot(t,y,'k x')
 plot(t,x,'g x')
 hold off
 legend('Concentración de B','Concentración de A’, ‘location’, ’best')

4 Sistema de Ecuaciones

Escribiremos ahora las ecuaciones (1) en forma de sistema

y’(t) = k1 x(t) y(t)

x’(t) = −k1 x(t) y(t)

Plantear y resolver el PVI asociado al sistema anterior por el método de Euler y el de Runge Kutta de orden 4 para tiempo t = 10 segundos. Dibujar las concentraciones en una sola figura e interpretar el resultado. Otra forma de plantear la ecuación anteriormente estudiada es mediante un sistema de ecuaciones, definiendo el problema de valor inicial asociado a este sistema, que es x(0)=1 e y(0)=0.01 )

%Partimos del sistema de ecuaciones
%y'(t)=k*x(t)*y(t)
%x'(t)=-k*x(t)*y(t)
%Donde:
%x(t)= concentración de A a lo largo del tiempo en mol/L
%y(t)= concentración de B a lo largo del tiempo en mol/L
%x'(t)= velocidad con la que decrece la concentración de A en mol/(L*s)
%y'(t)= velocidad con la que aumenta la concentración de B en mol/(L*s)
clear, close all
%Introducimos los datos del problema
t0=0;
tN=10;
y0=0.01; %Concentración de B en t=0
x0=1;    %Concentración de A en t=0
h=0.1;   %Paso de discretización
k=1;
%Calculamos N=número de intervalos
N=(tN-t0)/h;
%Discretización temporal
t=t0:h:tN;
%Preparación de un vector solución
y=zeros(1,length(t));
x=zeros(1,length(t));
%Tanto x como y son vectores de dimensiones 1xN+1
%Inicialización
y(1)=y0;
x(1)=x0;

• Aplicamos el método de Euler, previamente:

%Creamos el bucle aplicando el método de Euler
for n=1:N
   y(n+1)=y(n)+h*(k*y(n)*x(n));
   x(n+1)=x(n)+h*(-k*y(n)*x(n));
end
%Guardamos los resultados en un matriz A
A=[t;x;y]';
%Dibujamos la solución
figure(1);
hold on
plot(t,y,'gx-')
plot(t,x,'kx-')
hold off
legend('Concetración de B','Concentración de A','location','best')
title('Método de Euler')
xlabel('Tiempo (s)');
ylabel('Concentración (mol/L)');
axis([0 10 0 1.2])
%Representamos los puntos con cruces que luego hemos unido con una línea

• Aplicamos el método de Runge-Kutta, que también hemos usado con anterioridad:

%Creamos el bucle aplicando el método de Runge-Kutta
for n=1:N
   ky1=k*x(n)*y(n);
   kx1=-k*x(n)*y(n);
   ky2=k*(x(n)+h/2*kx1)*(y(n)+h/2*ky1);
   kx2=-k*(x(n)+h/2*kx1)*(y(n)+h/2*ky1);
   ky3=k*(x(n)+h/2*kx2)*(y(n)+h/2*ky2);
   kx3=-k*(x(n)+h/2*kx2)*(y(n)+h/2*ky2);
   ky4=k*(x(n)+h*kx3)*(y(n)+h*ky3);
   kx4=-k*(x(n)+h*kx3)*(y(n)+h*ky3);
   y(n+1)=y(n)+h/6*(ky1+2*ky2+2*ky3+ky4);
   x(n+1)=x(n)+h/6*(kx1+2*kx2+2*kx3+kx4);
end 
%Guardamos los resultados en un matriz A
A=[t;x;y]';
%Dibujamos la solución
figure(1);
hold on
plot(t,y,'gx-')
plot(t,x,'kx-')
hold off
legend('Concetración de B','Concentración de A','location','best')
title('Método de Euler')
xlabel('Tiempo (s)');
ylabel('Concentración (mol/L)');
axis([0 10 0 1.2])
%Representamos los puntos con cruces que luego hemos unido con una línea

Podemos observar como la representación gráfica coincide con la obtenida en el apartado anterior, ya que se trata del mismo problema resuelto por métodos diferentes. Un aspecto importante a la hora de aplicar Runge-Kutta a sistemas de ecuaciones diferenciales no lineales es que, al depender K2 de K1, K3 de K2, y K4 de K3, es necesario definir primero las Ki de ambas incógnitas antes que las Kj sindo i<j, de lo contrario, a la hora de aplicar el bucle, los valores de las Ki utilizados serían los de la anterior iteración del bucle, obteniendo, por lo tanto, resultados erróneos.

5 Conclusiones

En esta reacción biomolecular, tenemos inicialmente A en concentración \(1\frac{mol}{l}\) y B en concentración \(0.1\frac{mol}{l}\), en ella B actúa como catalizador, lo que implica que a medida que su concentración va aumentando, la velocidad de la reacción aumenta de manera exponencial, hasta llegar a un punto en que la cantidad de sustancia que queda por reaccionar es demasiado pequeña y la velocidad de reacción empieza a disminuir. La máxima velocidad de la reacción se alcanzará a los 4,5610 segundos del inicio de la reacción, instante en el que las concentraciones de A y B son iguales (0.5050 mol/L) a partir de ese momento, como ya hemos mencionado, la velocidad de la reacción empieza a disminuir, es decir, la pendiente de las curvas de las concentraciones se va acercando a la horizontal, haciéndose nula en el instante en que la concentración de A ha desaparecido prácticamente por completo (a los diez segundo de iniciarse la reacción), siendo éste el final de la reacción autocatalítica. Podemos observar que nuestras funciones tienen asíntotas horizontales en 0 y en 1.01. Estas asíntotas nos muestran los valores máximos y mínimos que pueden presentar ambas concentraciones. Por otra parte, la función y(t) es siempre creciente, mientras que la función x(t) es siempre decreciente. El crecimiento o decrecimiento de las funciones nos muestra como la concentración de A al principio es muy alta, mientras que la de B es mínima y, según transcurre el tiempo, las moléculas de A se van transformando en moléculas de B, por lo que la concentración de esta última aumenta a costa de la disminución de la otra (cumple la ley de concentración de masa). También podemos observar una simetría entre ambas funciones. La recta x=4.5610 es una recta de simetría, al igual que la y=0.5050. La simetría de las funciones nos muestra que B se produce al mismo ritmo que disminuye A, sumando ambos valores 1.01 en cualquier instante, lo que verifica el principio de la conservación de masa.

6 Deducción e interpretación de las ecuaciones diferenciales a partir de la reacción consecutiva de Lotka.

Partimos de la reacción consecutiva propuesta por Lotka:

A + X → k₁*2X

X + Y → k₂*2Y

Y → k₃*B

Con contantes: k₁, k₂ y k₃.

Donde A, X, Y, D son sustancias distintas. Observamos que las dos primeras reacciones son autocatalíticas. La reacción consume A para producir B mientras que X y Y dominan la velocidad y mezcla en los estados intermedios. Siguiendo la estrategia de los apartados 1 y 2 deducir las siguientes ecuaciones diferenciales para las concentraciones interpretando adecuadamente los términos:

x’(t) = k1 A x − k2 x y y’= k2 x y − k3 y B’= k3 y A’ + x’ + y’ + B’ = 0

Observamos que las dos primeras ecuaciones son autocatalíticas, la reacción transcurre consumiendo A y originándose B, como se ve en la reacción global

A + B → k₁*2B.

Para deducir las ecuaciones diferenciales nos basamos en los dos principios físicos ya nombrados, ley de acción de masas y principio de conservación de la masa.

Las concentración de las cuatro sustancias depende del tiempo, B(t), A(t), X(t), Y(t) pero las simplificamos a A, X, B, Y para un mejor manejo. Para deducir las tres primeras ecuaciones nos apoyamos en la ley de acción de masas, que dice: la velocidad de reacción es proporcional al producto de la concentración de los reactivos.

De la primera ecuación de Lotka sacamos que x’(t) = k1 A x, [1] de la segunda sacamos que y’(t)= k2 x y, [2] y, de la tercera, B’= k3 y , [3].

Además, sabemos, en la segunda reacción, que la velocidad con la que se forma y es igual desaparece x, por lo que podemos escribir x’(t)= -k2 x y. [4]. De la misma manera, a partir de la tercera reacción sacamos y’= -k3 y [5].

A partir de [1] y de [4] obtenemos la primera ecuación a demostrar:

x’ = k₁*A*x - k₂*x*y

A partir de [2] y [5] obtenemos la segunda ecuación a demostrar:

y’ = k₂*x*y - k₃*y

La tercera ecuación es más sencilla a deducir puesto que es, directamente, [3]:

B’ = k₃*y

Para la última ecuación diferencial a deducir partimos del principio de conservación de la masa, que establece que la suma de las concentraciones de los compuestos es constante (las pérdidas de masa son despreciables): A+B+X+Y= cte.

Lo que derivando nos da la cuarta ecuación a deducir:

x’+y’+A’+B’= 0

Quedando así deducidas las cuatro ecuaciones diferenciales pedidas, que se emplearan en los siguientes apartados.

7 Resolución del sistema de la reacción de Lotka por el método de Euler

Plantear el PVI asociado al sistema anterior y resolverlo en el tiempo t (0, 200) para los siguientes valores k1 = k2 = 2k3 = 0.1, y concentraciones iniciales de A, X, Y, B dadas por 5, 5*10−4, 10−5, 0 respectivamente. Usar el método de Euler con h = 0.01 y h = 0.001. ¿Es estable?.

clear all
%Definimos la variable independiente y establecemos los datos del problema
%de valor inicial. Para la longitud de paso h introducimos en primer lugar
%0.01 y luego 0.001
t0=0;
tN=200;
h= input('Introduzca una longitud de paso ');
t = t0:h:tN;
N=(tN-t0)/h;
k1 = 0.1;
k2 = 0.1;
k3 = 0.05;
A0 = 5;
B0 = 0;
X0 = 5*10^(-4);
Y0 = 10^(-5);
%Definimos las variables dependientes.
A = zeros(1, N+1);
B = A;
X = A;
Y = A;
A(1) = A0;
B(1) = B0;
X(1) = X0;
Y(1) = Y0;
for i = 1:N
    X(i+1) = X(i) + h*(k1*A(i)*X(i)-k2*Y(i)*X(i));
    Y(i+1)=Y(i)+h.*(k2.*X(i).*Y(i)-k3.*Y(i));
    B(i+1) = B(i) + h*(k3*Y(i));
    A(i+1) = A(i) + h*(-k1*A(i)*X(i));
End
%Representamos gráficamente la solución
hold on
plot (t,A,'c')
plot(t, X, 'r')
plot(t,Y, 'g')
plot(t, B, 'k')
legend('Sustancia A: A(t)','Sustancia X: x(t)','Sustancia Y: y(t)','Sustancia B: B(t)','Location','Best')
xlabel('Tiempo (s)')
ylabel('Concentración (mol/l)')
grid on
hold off

Un método es estable cuando el error local es pequeño y además este error no se amplifica a lo largo del tiempo. Para calificar la estabilidad del método de Euler resolveremos el sistema anterior con este método, con pasos "h=0.01" y "h=0.001", y lo compararemos con la resolución del sistema con el método de Heun, que se asemeja más a la solución exacta. En la comparación representaremos en azul el método de Heun y en rojo el método de Euler. Para "h=0.01" tenemos lo siguiente: Como no se separan ambos métodos podemos afirmar que es estable.

8 Resolución por Heun

%Intorducimos datos del tiempo y calculamos su vector
t0=0; tf=200;
h=0.01; N=(tf-t0)/h;
t=t0:h:tf;
% introducimos datos aportados por el enunciado
A0=5; X0=5*(10^-4); Y0=10^(-5); B0=0;
%creamos las matrices A B X Y que tendrán guardados los puntos solución
%para cada instante de tiempo
A=zeros(1,N+1);
X=zeros(1,N+1);
Y=zeros(1,N+1);
B=zeros(1,N+1);
% introducimos los valores iniciales dentro de las matrices
A(1)=A0; X(1)=X0; Y(1)=Y0; B(1)=B0;
aa=A0; xx=X0; yy=Y0; bb=B0;
% introducimos los valores de las constantes dadas en el enunciado
K1=0.1; K2=K1; K3=(K1/2);

% creamos un bucle for que nos va a ir solucionando el PVI por el método de
% Heun para cada valor de t y metiendo este valor en cada uno de los
% lugares correspondientes a cada matriz
for n=1:N
   K1A=-K1*aa*xx;
   K2A=K1A+K1A*h;
   K1X=K1*aa*xx-K2*xx*yy;
   K2X=K1X+K1X*h;
   K1Y=K2*xx*yy-K3*yy;
   K2Y=K1Y+K1Y*h;
   K1B=K3*yy;
   K2B=K1B+K1B*h;
   aa=aa+0.5*h*(K1A+K2A);
   xx=xx+0.5*h*(K1X+K2X);
   yy=yy+0.5*h*(K1Y+K2Y);
   bb=bb+0.5*h*(K1B+K2B);
   A(n+1)=aa;
   X(n+1)=xx;
   Y(n+1)=yy;
   B(n+1)=bb;
   end
% una vez tenemos todas las matrices completas dibujamos, para ello usa bos
% la función plot. Para poder dibujar varias gráficas dentro de una misma
% pestaña usamos el comando hold on // hold off
plot(t,A,'*r');
hold on
plot(t,X,'*g');
plot(t,Y,'*y');
plot(t,B,'*');
hold off
xlabel('Tiempo(s)','FontSize',11);
ylabel('concentración (mol/L)','Fontsize',11);
title('CONCENTRACIÓN - TIEMPO (Heun)','FontName','Berlin sans FB','FontSize',11);
legend('concentración de A','concentración de X','contenctración de Y','concentración de B','location','best');
%xlabel, ylabel, title y legend son comandos que nos permiten escribir en
%la gráfica qué es cada eje, un titulo y la legenda para la mejor
%comprensión de la gráfica.