Reacciones de autocatalisis Grupo 9A

De MateWiki
Saltar a: navegación, buscar
Trabajo realizado por estudiantes
Título Reacciones con autocatálisis. Grupo A2
Asignatura Ecuaciones Diferenciales
Curso Curso 2014-15
Autores David Carmona Rodriguez,Alejandro Muñoz Cotter, Daniel Alonso Palop, Luis Bermeosolo Echeverria
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 Introducción, comentarios generales y planteamiento de la primera reacción

La autocatálisis es el proceso mediante el cual un compuesto químico induce y controla una reacción química sobre sí mismo. Los compuestos autocatalíticos no son catalizadores en sentido estricto ya que su estructura química resulta alterada durante el proceso.

Consideramos una solución bien mezclada a temperatura y volumen constantes. En esta solución tiene lugar una reacción química en la que en el momento inicial se encuentran dos reactivos A y B. A medida que avanza el tiempo se forma el producto 2B, teniendo en cuenta que la presencia de B hace de efecto catalítico en la reacción y satisfaciendo la ley de acción de masas que establece que la velocidad de reacción es proporcional al producto de las concentraciones de los reactivos.

Reacción bimolecular: [math] A + B \rightarrow _{k1} 2B [/math]


Para tiempo [math](t=0)[/math] nombramos e identificamos las variables:


A: concentracion inicial [math]x_0[/math] (mol/l)


B:concentracion inicial [math]y_0[/math] (mol/l)


Sucede la reaccion[math] \rightarrow [/math] El tiempo comienza [math](t\gt0)[/math]


A: concentracion en funcion de t. [math]x(t)[/math] (mol /l)

B:concentracion en función de t. [math]y(t)[/math] (mol/l)

Como el volumen se mantiene constante [math](V=cte)[/math]


[math]V=volumen[/math]


[math]M_A(t)[/math] = masa del compuesto A


[math]M_B(t)[/math] = masa del compuesto B


Según la ley de concentración de la masa: [math]M_A(t)+M_B(t)=k[/math] siendo [math]k=''cte''[/math] en el proceso, si dividimos por V:

[math]\frac{M_A(t)}{V} + \frac{M_B(t)}{V} =\frac{k}{V}[/math]          (renombramos [math]k^*=\frac{K}{V}[/math])


Sustituimos por nuestros términos:

[math]x(t)+y(t) = k^*[/math]


Derivando [math]x'(t)=+y'(t)=0[/math] ec.(1)


Segun la ley de acción de masas:


[math]\mbox{Velocidad de reacción = (cte)·(Cantidad de reactivo A)·(cantidad de reactivo de B)}[/math]


[math]y'(t)=k1*x(t)*y(t)[/math] ec.(2)


si integramos la ec.(1) [math]x(t)+y(t)=k^*[/math] (con [math]k^*=\frac{k}{V}[/math]) despejamos [math]x(t)=k^*-y(t)[/math]

y sustituyendo en ec.(2) ya tenemos planteado el P.V.I con las condiciones iniciales dadas en el enunciado :


[math] \left \{ \begin{matrix} y'(t) = k_{1}*(K_{2} - y(t))*y(t) \\ y(0)=0.01 \end{matrix} \right . [/math]

Ahora procedemos a comprobar si tiene solucion y ademas es unica mediante la aplicacion del teorema de la existencia y unicidad:

Existe solución para el PVI planteado si existe una "Bola" de radio [math] r [/math] alrededor del punto de estudio [math] (t_{0},y_{0}) [/math] en la que la función sea continua. Esta solución será además única cuando también se cumpla que [math] \displaystyle{ \delta f \over \delta y} [/math] sea también contínua.

Observamos que nuestra funcion: [math] f(t,y(t))= k_{1}*(K_{2} - y(t))*y(t) [/math] Y: [math] \displaystyle{ \delta f \over \delta y} = k_{1}(K_{2} - 2y) [/math]

La función [math] f(t,y(t))= k_{1}*(K_{2} - y(t))*y(t) [/math] es de clase [math]C^1[/math] y su derivada es continua siempre, entonces podemos afirmar que existe solución única.

2 Primera reacción propuesta

Procedemos a resolver el PVI mediante metodos numericos estudiados en la asignatura interpretando los resultados relacionados con el proceso quimico:

2.1 Ecuación

2.1.1 Método de Euler

Resolvemos la ecuacion mediante el programa de EULER que consiste en un algoritmo basado en la formula: [math]y_{n+1} = y_n + h f (t_n,y_n)[/math]

que permite dar en un numero finito de pasos [math](N= \frac{t_n-t_0}{V})[/math] un aproximacion numerica a la solucion del problema.

%METODO DE EULER
%DATOS DEL PROBLEMA
t0=0; tN=10; y0=0.01; K1=1; K2=1.01;h=0.1;
%ELECCION DEL PASO
%GENERACION DEL VECTOR TIEMPO t EN FUNCION DE h
t=t0:h:tN; 
%PREPARACION DEL VECTOR SOLUCION APROXIMADA
N=(tN-t0)/h;
y=zeros(1,N+1); 
y(1)=y0;
for i=1:N
y(i+1)=y(i)+h*(K1*(K2-y(i))*y(i)); %METODO DE EULER
x(i+1)=1.01-y(i+1);
end 
%FINALIZO EL PROGRAMA
%GRAFICAS
hold on
plot(t,y,'b');
plot(t,x,'r');
hold off
xlabel('segundos (s)');
ylabel('Concentracion (mol/L)');
legend('Concentracion de B','Concentracion de A','Location','Best');


Método Numérico de Euler

Obtenemos esta grafica que será comentada en la conclusion de este apartado junto con las demas graficas obtenidas.

2.1.2 Método del Trapecio

Resolvemos la ecuacion mediante el metodo del trapecio. Otro metodo numerico para aproximar la solucion de la ecuacion con menor error que el metodo de EULER. Es un metodo implicito, por tanto habra que despejar el termino [math]y_{n-1}[/math] y aplicar la formula: [math] \begin{array}{c}\\ y_{n+1}=y_n+{h \over 2}*[f(t_n,y_n)+f(t_{n+1},y_{n+1})]\end{array} [/math]

Despejando analiticamente:

[math] \begin{array}{c}y_{n+1}=y_n+{h \over 2}*[y_n*(K_2-y_n)+y_{n+1}*(K_2-y_{n+1})]\\y_{n+1}=y_n+{h \over 2}*y_n*(K_2-y_n)+{h \over 2}*y_{n+1}*(K_2-y_{n+1})\\y_{n+1}-{h \over 2}*y_{n+1}*(K_2-y_{n+1})=y_n*[1+{h \over 2}*(K_2-y_n)]\\{h \over 2}*(y_{n+1})^2+y_{n+1}*(1-{K_2*h \over 2}+[-y_n*[1+{h \over 2}*(K_2-y_n)]])=0\\y_{n+1}={-(1-{K_2*h \over 2})+\sqrt{(1-{K_2*h \over 2})^2-4*{h \over 2}*[-y_n*[1+{h \over 2}*(K_2-y_n)]]} \over 2*{h \over 2}}\\y_{n+1}={-1+{K_2*h \over 2}+\sqrt{(1-{K_2*h \over 2})^2-2*h*[-y_n*[1+{h \over 2}*(K_2-y_n)]]} \over h}\end{array} [/math]


Codigo Matlab:


%METODO DEL TRAPECIO
%DATOS DEL PROBLEMA
t0=0; tN=10;h=0.1; y0=0.01; K1=1; K2=1.01;
%generacion del vector tiempo, conocido h
t=t0:h:tN;
%calculo del numero de intervalos
N=(tN-t0)/h; %calculo del numero de intervalos
%preparacion del vector solucion aproximada
y=zeros(1,length(t));
%inicializacion
y(1)=y0;
%aplicacion del metodo numerico del trapecio
for i=1:N 
y(i+1)=(1/(h*K1))*((0.5*h*K1*K2-1)+sqrt((1-0.5*h*K1*K2)^2-2*h*K1*(-y(i)-(h/2)*y(i)*(K2-y(i)))));
end
%calculamos ahora la concentracion de A mediante la ley de conservacion de masa
x=1.01-y;
%Graficas.
hold on
plot(t,y,'b');
plot(t,x,'r');
hold off
xlabel('segundos (s)');
ylabel('Concentracion (mol/L)');
legend('Concentracion de B','Concentracion de A','Location','Best');


Método Numérico del Trapecio

Obtenemos la misma gráfica.

2.1.3 Método de Runge-Kutta

Resolvemos mediante Rounge Kutta de 4º orden

%METODO DE RUNGE-KUTTA
%DATOS DEL PROBLEMA
t0=0; tN=10;h=0.1; y0=0.01; x0=1; K1=1; K2=1.01;
%generacion del vector tiempo, conocido h
t=t0:h:tN;
%calculo del numero de intervalos
N=(tN-t0)/h; %calculo del numero de intervalos
%preparacion del vector solucion aproximada
y=zeros(1,length(t));
%inicializacion
y(1)=y0;
x(1)=x0;
x=y;
%inicializacion
y(1)=y0;
x(1)=x0;
U=inline('x*y','t','y','x');
V=inline('-x*y','t','y','x');
%aplicacion del metodo numerico del trapecio
for k=1:N
  %Runge-Kutta de orden 4 (RK4):
  K1_y=U(t(k),y(k),x(k));
  K1_x=V(t(k),y(k),x(k));
  K2_y=U(t(k)+(h/2),y(k)+h*K1_y/2,x(k)+h*K1_x/2);
  K2_x=V(t(k)+(h/2),y(k)+h*K1_y/2,x(k)+h*K1_x/2);
  K3_y=U(t(k)+(h/2),y(k)+h*K2_y/2,x(k)+h*K2_x/2);
  K3_x=V(t(k)+(h/2),y(k)+h*K2_y/2,x(k)+h*K2_x/2);
  K4_y=U(t(k)+h,y(k)+K3_y*h,x(k)+K3_x*h);
  K4_x=V(t(k)+h,y(k)+K3_y*h,x(k)+K3_x*h);
 
  y(k+1)=y(k)+(h/6)*(K1_y+2*K2_y+2*K3_y+K4_y);
  x(k+1)=x(k)+(h/6)*(K1_x+2*K2_x+2*K3_x+K4_x);
end
%Graficas.
hold on
plot(t,y,'b');
plot(t,x,'r');
hold off
xlabel('segundos (s)');
ylabel('Concentracion (mol/L)');
legend('Concentracion de B','Concentracion de A','Location','Best');


Método Numérico de Runge-Kutta

Obtenemos también la misma gráfica.

2.1.4 Interpretacion de las graficas y soluciones

En una reacción autocatalítica si comenzamos con una cantidad pequeña de B, la velocidad de reacción aumentará a medida que se vaya formando más B. En el otro extremo, cuando haya desaparecido prácticamente todo el componente A, la velocidad ha de tender a cero. Este comportamiento se puede apreciar en la gráfica anterior, en la que la velocidad varía a lo largo de una función cuyo máximo corresponde a concentraciones iguales de A y de B.


Las funciones por las que se rige la reaccion son simetricas, a la vez que una disminuye, la otra crece, lo que era de esperar ya que a medida que se forma el compuesto B disminuye la cantidad de compuesto A que tenemos.

2.2 Sistema de ecuaciones

Ahora escribimos las ecuaciones en forma de sistema:

[math] \left \{ \begin{matrix} y'(t)=k_{1}x(t)y(t) \\ x'(t)=-k_{1}x(t)y(t) \end{matrix} \right . [/math]

Con las condiciones iniciales: La concentracion de A para tiempo inicial cero es [math] 1\frac{mol}{l} [/math], mientras que la concentracion de B es mucho menor [math] 0.01\frac{mol}{l} [/math]. EL P.V.I nos queda:

[math] \left \{ \begin{matrix} y'(t)=k_{1}x(t)y(t) \\ x'(t)=-k_{1}x(t)y(t) \\ x(0)=1 \\ y(0)=0.01 \end{matrix} \right . [/math] Con [math] k=1\frac{mol}{s} [/math]

Resolvemos el PVI mediante el metodo de EULER y Runge Kutta 4º orden teniendo en cuenta que son ecuaciones no lineales:

%METODO DE EULER PARA SISTEMAS
k=1;t0=0; tN=10;h=0.1; y0=0.01; x0=1;
%generacion del vector tiempo, conocido h
t=t0:h:tN;
%calculo del numero de intervalos
N=(tN-t0)/h; 
%preparacion del vector solucion aproximada
y=zeros(1,length(t));
x=y;
%inicializacion
y(1)=y0;
x(1)=x0; 
for i=1:N  
    y(i+1)=y(i)+h*(x(i)*y(i));
    x(i+1)=x(i)+h*(-x(i)*y(i));
end
%Graficas
hold on
plot(t,y,'b');
plot(t,x,'r');
hold off
xlabel('segundos (s)');
ylabel('Concentracion (mol/L)');
legend('Concentracion de B','Concentracion de A','Location','Best');


Método Numérico de Euler en Sistema

Como podemos observar, el resultado es el mismo que en los casos anteriores, sólo que esta vez se ha usado un sistema de ecuaiones para resolverlo.

%METODO DE RUNGE-KUTTA PARA SISTEMA
t0=0; tN=10; y0=0.01; x0=1; h=0.1;
%generacion del vector tiempo, conocido h
t=t0:h:tN;
%calculo del numero de intervalos
N=(tN-t0)/h; %calculo del numero de intervalos
%preparacion del vector solucion aproximada
y=zeros(1,length(t));
x=y;
%inicializacion
x(1)=x0;
y(1)=y0;
%preparamos las funciones a usar en el bucle
U=inline('x*y','t','y','x');
V=inline('-x*y','t','y','x');
%Recorremos, con un bucle, los vectores, y en ellos almacenamos los
%resultados obtenidos con RK4.
for k=1:N
 %Runge-Kutta de orden 4 (RK4):
 k1_y=U(t(k),y(k),x(k));
 k1_x=V(t(k),y(k),x(k));
 k2_y=U(t(k)+(h/2),y(k)+h*k1_y/2,x(k)+h*k1_x/2);
 k2_x=V(t(k)+(h/2),y(k)+h*k1_y/2,x(k)+h*k1_x/2);
 k3_y=U(t(k)+(h/2),y(k)+h*k2_y/2,x(k)+h*k2_x/2);
 k3_x=V(t(k)+(h/2),y(k)+h*k2_y/2,x(k)+h*k2_x/2);
 k4_y=U(t(k)+h,y(k)+k3_y*h,x(k)+k3_x*h);
 k4_x=V(t(k)+h,y(k)+k3_y*h,x(k)+k3_x*h);
 
 y(k+1)=y(k)+(h/6)*(k1_y+2*k2_y+2*k3_y+k4_y);
 x(k+1)=x(k)+(h/6)*(k1_x+2*k2_x+2*k3_x+k4_x);
end
%Graficas
hold on
plot(t,y,'b');
plot(t,x,'r');
hold off
xlabel('segundos (s)');
ylabel('Concentracion (mol/L)');
legend('Concentracion de B','Concentracion de A','Location','Best');


Método Numérico de Runge-Kutta en Sistema

Como podemos observar de nuevo, el resultado es el mismo que en los casos anteriores. Esta vez se ha usado un sistema de ecuaciones, como en el caso de Euler con sistema, pero esta vez resuelto por el método numérico para sistemas de Runge-Kutta.

3 Segunda reacción propuesta: Reacción consecutiva propuesta por Lotka

3.1 Deducción e interpretación de las ecuaciones diferenciales (Apartado 5)

Consideramos ahora la siguiente reacción consecutiva propuesta por Lotka (1920):

A+X→2X (Con cte k1) ; X+Y→2Y (Con cte k2)  ; Y→B (Con cte k3)

Donde A, X, B e Y son sustacias distintas. Observamos que las etapas 1 y 2 son autocatalíticas ya que vemos autocatálisis en los sustancias X e Y respectivamente. La reacción transcurre consumiendo A para producir el producto final B, de acuerdo con la reacción global: A→B.

Los intermedios X e Y dominan la velocidad y la composición de la mezcla reactiva en las fases intermedias, pero acaban por desaparecer como se observa.

Denotamos x=x(t), y=y(t), A=A(t) y B=B(t) las concentraciones de las sustancias X, Y, A y B respectivamente. Según el principio de conservación de la masa sabemos que:

                                 A+x+y+B=constante

Si derivamos la expresión anterior respecto del tiempo se llega a la primera ecuación:

                                   A'+x'+y'+B'=0

Por otro lado si observamos la primera y la segunda etapa y nos fijamos en que le ocurre a la sustancia X se llega a la siguiente conclusión:

                                     x'=k1*A*x-k2*x*y

La cual nos dice que la variación de la concentración respecto del tiempo de la sustancia X es proporcional a la concentración de A y de X (con cte k1 (etapa 1)) y a la concentracion de X e Y (con cte k2 (etapa 2)). Teniendo en cuenta además que el primer sumando es positivo pues se está creando sustancia X y el segundo negativo pues se está elimando sustancia X (para crear la sustancia Y).

Análogo razonamiento haríamos con la sustancia Y y la sustancia B mirando repectivamente las etapas 2,3 y 3; llegando así a las siguientes ecuaciones:

                         y'=k2*x*y-k3*y       ;       B'=k3*y

3.2 Planteamiento, resolución por Euler e interpretación del PVI (Apartado 6)

Aplicamos el metodo de Euler para resolver el siguiente sistema:


\begin{matrix} x'=k1Ax-k2xy\\ y'=k2xy-k3y\\ A'=-k1Ax\\ B'=k3y\\ \end{matrix} Teniendo en cuenta las siguientes condiciones:

\begin{matrix} k1=k2=2k3=0.5\\ A=5\\ B=0\\ x=5·10^{-4}\\ y=10^{-5}\\ \end{matrix}

Codigo matlab:

% Introducimos las constantes iniciales dadas en el enunciado
t0=0; tf=200;
k1=0.1; k2=0.1; k3=0.05;
A0=5;
x0=5*10^(-4);
y0=10^(-5);
B0=0;

%Metodo de Euler para h=0.01
h1=0.01;
t1=[t0:h1:tf]; % Introducimos el vector independiente
N1=(tf-t0)/h1; % Numero de subintervalos
A1=linspace(0,0,N1+1); % Creamos los vectores dependientes para h=0.01
x1=linspace(0,0,N1+1);
y1=linspace(0,0,N1+1);
B1=linspace(0,0,N1+1);
A1(1)=A0; x1(1)=x0; y1(1)=y0; B1(1)=B0; % Introducimos las constantes iniciales en el primer termino de los vectores dependientes
for i=1:N1 % Metodo de Euler
        A1(i+1)=A1(i)+h1*((-k1)*x1(i)*A1(i));
        x1(i+1)=x1(i)+h1*(k1*A1(i)*x1(i)-k2*y1(i)*x1(i));
        y1(i+1)=y1(i)+h1*(k2*x1(i)*y1(i)-k3*y1(i));
        B1(i+1)=B1(i)+h1*(k3*y1(i));
end
%Metodo de Euler para h=0.001
h2=0.001;
t2=[t0:h2:tf]; % Introducimos el vector independiente
N2=(tf-t0)/h2; % Numero de subintervalos
A2=linspace(0,0,N2+1); % Creamos los vectores dependientes para h=0.001
x2=linspace(0,0,N2+1);
y2=linspace(0,0,N2+1);
B2=linspace(0,0,N2+1);
A2(1)=A0; x2(1)=x0; y2(1)=y0; B2(1)=B0;
for j=1:N2 % Metodo de Euler
         A2(j+1)=A2(j)+h2*((-k1)*x2(j)*A2(j));
        x2(j+1)=x2(j)+h2*(k1*A2(j)*x2(j)-k2*y2(j)*x2(j));
        y2(j+1)=y2(j)+h2*(k2*x2(j)*y2(j)-k3*y2(j));
        B2(j+1)=B2(j)+h2*(k3*y2(j));
end
subplot(1,2,1) % Comandos para la visualizacion
hold on
plot(t1,A1,'r','linewidth',2)
plot(t1,B1,'k','linewidth',2)
plot(t1,x1,'b','linewidth',2)
plot(t1,y1,'g','linewidth',2)
xlabel('tiempo')
ylabel('concentracion')
title('Grafico del metodo de Euler con h=0.01')
legend('Concentración de A','Concentración de B','Concentración de X','Concentración de Y')
hold off
subplot(1,2,2)
hold on
plot(t2,A2,'r','linewidth',2)
plot(t2,B2,'k','linewidth',2)
plot(t2,x2,'b','linewidth',2)
plot(t2,y2,'g','linewidth',2)
xlabel('tiempo')
ylabel('concentracion')
title('Grafico del metodo de Euler con h=0.001')
legend('Concentración de A','Concentración de B','Concentración de X','Concentración de Y')
hold off

Mediante este procedimiento se obtienen las siguientes graficas: centro

3.3 Planteamiento, resolución por Heun e interpretación del PVI (Apartado 7)

El método de Heun es un método explícito, debemos primeramente definir una serie de constantes K1 y K2.

[math] \left \{ \begin{matrix} y0,t0\\ y_{(n+1)}=y_n+\frac{h}{2}(K1+K2)\\ K1=f(t_n,y_n)\\ K2=f(t_n+h,y_n +K1h) \end{matrix} \right . [/math]

Aplicando el metodo mediante el programa Matlab:

%Datos del problema de valor inicial
t0=0;
tN=200;
A0=5;
X0=5*(10^-4);
Y0=10^(-5);
B0=0;
k1=0.1;
k2=0.1;
k3=0.05;
h=0.01;
%Vector de tiempo t de t0 a tN con paso h
t=t0:h:tN;
%N=número de intervalos
N=(tN-t0)/h;
%Vectores vacíos donde se almacena solución
A=zeros(1,N+1);
X=zeros(1,N+1);
Y=zeros(1,N+1);
B=zeros(1,N+1);
%Valores de arranque
A(1)=A0;
X(1)=X0;
Y(1)=Y0;
B(1)=B0;
for i=1:N
%Definimos las constantes
%Valor de K1 y K2 para la ecuacion de X
K1X=k1.*A(i).*X(i)-k2.*X(i).*Y(i);
K2X=k1.*(A(i)+K1X.*h).*(X(i)+K1X.*h)-k2.*(X(i)+K1X.*h).*(Y(i)+K1X.*h);
%Valor de K1 y K2 para la ecuacion de Y
K1Y=k2.*X(i).*Y(i)-k3.*Y(i);
K2Y=k2.*(X(i)+K1Y.*h).*(Y(i)+K1Y.*h)-k3.*(Y(i)+K1Y.*h);
%Valor de K1 y K2 para la ecuacion de B
K1B=k3.*Y(i);
K2B=k3.*(Y(i)+K1B.*h);
%Valor de K1 y K2 para la ecuacion de A
K1A=-k1.*X(i).*A(i);
K2A=-k1.*(X(i)+K1A.*h).*(A(i)+K1A.*h);
%Resolucion X Y B A
X(i+1)=X(i)+0.5*h.*(K1X+K2X);
Y(i+1)=Y(i)+0.5*h.*(K1Y+K2Y);
B(i+1)=B(i)+0.5*h.*(K1B+K2B);
A(i+1)=A(i)+0.5*h.*(K1A+K2A);
end

%Gráficas separadas para interpretar.
%Gráficas separadas para interpretar.
subplot(2,2,1)
plot(t,A,'r')
title('Concentracion de A en funcion de t');
subplot(2,2,2)
plot(t,X)
title('Concentracion de X en funcion de t');
subplot(2,2,3)
plot(t,Y,'g')
title('Concentracion de Y en funcion de t');
subplot(2,2,4)
plot(t,B,'b');
title('Concentracion de B en funcion de t');


centro

3.3.1 INTERPRETACIÓN DE LAS GRÁFICAS

La concentracion de A tiende a 0 a medida que avanza el tiempo hasta llegar a un tiempo de aproximadamente 50 segundos.


[math] A + X \rightarrow _{k1} 2X [/math] (1)


Segun A disminuye, por tanto reacciona, se va formando el compuesto X en la primera reaccion autocatalitica, cuando el total de la concentracion de A ha reaccionado y se ha formado X, hay un periodo de tiempo en el que la concentracion de X se mantiene (podemos verlo en la grafica mediante la pendiente que es practicamente nula en el intervalo de tiempo [35,45]. En el instante en el que X tiene mayor cantidad se da con mayor velocidad la reaccion de autocatalisis 2, (ley de accion de masas) en la que el compuesto Y interviene para su propia formacion.


[math] X+ Y \rightarrow _{k2} 2Y [/math] (2)


En el instante t=48 (aproximadamente) empieza a haber una concentracion considerable de Y y se empieza a formar el compuesto B cuya concentracion final coincide con la concentracion inicial de A.


[math] Y \rightarrow _{k3} B [/math]