Reaccion Autocatalisis Grupo 2B

De MateWiki
Saltar a: navegación, buscar
Trabajo realizado por estudiantes
Título Reacción con autocatálisis. Grupo 2-B
Asignatura Ecuaciones Diferenciales
Curso Curso 2016-17
Autores Gonzalo Royo Navajas, Marta Caracuel Mateos, Carlota Sánchez Martínez, Abid Al-Akioui Sanz, Alejandro Prieto Martínez, Pablo Retamar Leboutet
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 Introducción

Una reacción de autocatálisis es un proceso mediante el cual un compuesto químico induce y controla una reacción química sobre sí mismo. En este caso, se estudiará una reacción regida por la constante k, bimolecular e irreversible, en la cual se suponen unos reactivos A y B que presentan un volumen y una temperatura constantes:

                           				A+B → C

Dichas condiciones satisfacen la Ley de acción de masas, justifica que "La velocidad de la reacción en cuestión será proporcional al producto de las concentraciones de A y B". Para desarrollar el análisis, es necesario apoyarse en el Principio de conservación de la materia, por el cual "la suma de la masa de los reactivos y la de los productos se mantendrá constante a lo largo del tiempo". Así pues, se obtendrá un tercer compuesto que actuará como producto:

                                                      A+B → k * 2B

La resolución de la ecuación diferencial se apoya en las dos leyes mencionadas con anterioridad. De esta forma, es posible determinar las siguientes ecuaciones: [math] \left\{\begin{matrix} x=A+B\\ y=2*B \end{matrix}\right. [/math]

Por el método de sustitución se deducen las variables dependientes A y B: [math] \left\{\begin{matrix} A=x-\frac{y}{2}\\ B=\frac{y}{2} \end{matrix}\right. [/math]

En base al Principio de conservación de la masa, se conoce la igualdad:

                                                      x+y=cte                                                                       (1)

A su vez, también se sabe que la velocidad es proporcional al producto de las concentraciones de los reactivos: v=k*A*B, por lo tanto, v=y’=-x’. De esta forma se concluye que y’=k*A*B en la que se sustituyen los valores de los reactivos previamente determinados:

    y’=k*[x-y/2]*[y/2] → y'=k*[x-y/2]*y → y'=k*[cte-y-y/2]*y → y'=k*[cte-3/2*y]*y → y'=3/2*k*[2/3*cte-y]*y → y'=k*[cte-y]*y → y'=k*x*y            

De (1) se puede afirmar que:

                                                     x=cte-y

Y sustituyendo en la fórmula anteriormente demostrada, se obtiene que:

                                                 y'=-k*y^2+k*y*cte

A partir de la anterior ecuación diferencial se puede formar el PVI asociado: [math] \left\{\begin{matrix} y'=-k*y^2+k*y*cte\\ y(0)=y_0 \end{matrix}\right. [/math]

Para poder resolver el PVI se necesita calcular el valor de la cte:

                                                 x=A+B → A=x-y/2
                                                   y=2B → B=y/2

A partir de los datos numéricos, A=2 y B=0.05 mol/L, y despejando en las ecuaciones anteriores, se obtienen los valores de x e y, siendo 1.95 y 0.1 respectivamente. Por lo tanto:

                                              cte=x+y=1.95+0.1=2.05

Para comprobar que el PVI tiene solución única se aplica el Teorema de Picard, siendo:

                                              f=-k*y^2+k*y*cte

Al ser f una función polinómica está definida y es continua en RxR. A continuación, se calcula la derivada parcial de f con respecto a y, dando:

                                              f'=-2*k*y+k*cte

Análogamente, al ser la derivada una función polinómica, está definida y es continua en RxR, concluyendo, así, que el PVI tiene solución y es única.

2 Resolución del PVI

La resolución del PVI se ha llevado acabo por diferentes métodos. Para ello, se ha escogido un paso h=0.1 y una k=1.4 mol/s. El estudio se realiza durante los primeros 10 segundos en los que se produce la reacción.

Método de Euler, cuyo esquema numérico se indica a continuación:

\[ PVI = \begin{cases} y_0 \\ y_{n+1} = y_n + h*f(t_n,y_n) \\ \end{cases}\quad \]

caption
 1 % Datos
 2 A=2;
 3 B=0.05;
 4 k=1.4;
 5 h=0.1;
 6 t0=0;
 7 tN=10;
 8 c=2.05;
 9 f=@(t,y)(-k*y^2+k*c*y);
10 % Discretizamos
11 N=round((tN-t0)/h);
12 % Creamos el vector de tiempo t 
13 t=linspace(t0,tN,N+1);
14 % Vector con la solucion por aproximacion
15 y=zeros(size(t));
16 y(1)=B;
17 % Utilizamos el metodo de Euler 
18 for i=1:N
19     y(i+1)=y(i)+h*f(t(i),y(i));
20 end
21 x=c-y;
22 % Dibujamos la solucion aproximada
23 hold on
24 plot(t,y,'b')
25 plot(t,x,'g')
26 hold off
27 legend('Cantidad de productos','Cantidad de reactivos')
28 xlabel('Tiempo (s)');
29 ylabel('Concentración (mol/L)');
30 % Calculamos el momento en el que coinciden las dos concentraciones:
31 for i=1:N
32     if abs(y(i)-x(i))<0.01
33         fprintf('Corte en punto (%.2f,%.2f)\n',t(i),y(i))
34     end
35 end


Métodos del Trapecio y Runge-Kutta:

El método del Trapecio sigue siguiente método numérico: \[ 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 \]

Mientras que el de Runge-Kutta: \[ PVI = \begin{cases} y_0 \\ y_{n+1} = y_n + \frac{h}{6}*(K_1+2*K_2+2*K_3+k_4) \\ K_1=f(t_n,y_n) \\ K_2=f(t_n+\frac{1}{2}*h,y_n+\frac{1}{2}*K_1*h) \\ K_3=f(t_n+\frac{1}{2}*h,y_n+\frac{1}{2}*K_2*h) \\ K_4=f(t_n+h,y_n+K_3*h)\\ \end{cases}\quad \]

caption
 1 % Datos
 2 A=2;
 3 B=0.05;
 4 k=1.4;
 5 h=0.1;
 6 t0=0;
 7 tN=10;
 8 c=2.05;
 9 
10 f=@(t,y)(-k*y^2+k*c*y);
11 
12 % Discretizamos
13 N=round((tN-t0)/h);
14 
15 % Creamos el vector de tiempo t 
16 t=linspace(t0,tN,N+1);
17 
18 % Metodo del trapecio
19 yt=zeros(1,N+1);
20 yt(1)=B;
21 for i=1:N
22     yt(i+1)=(-1+0.5*h*k*c+sqrt((1-0.5*h*k*c)^2+2*h*k*(yt(i)+0.5*h*f(t(i),yt(i)))))/(h*k);
23 end
24 
25 xt=c-yt;
26 
27 hold on
28 plot(t,yt,'k')
29 plot(t,xt,'b')
30 hold off
31 legend('Cantidad de productos (trapecio)','Cantidad de reactivos (trapecio)')
32 
33 % Metodo de Runge-Kutta
34 yr=zeros(1,N+1);
35 yr(1)=B;
36 for j=1:N
37     K1=f(t(j),yr(j));
38     K2=f(t(j)+0.5*h,yr(j)+0.5*h*K1);
39     K3=f(t(j)+0.5*h,yr(j)+0.5*h*K2);
40     K4=f(t(j)+h,yr(j)+K3*h);
41     yr(j+1)=yr(j)+h*(K1+2*K2+2*K3+K4)/6;
42 end
43 xr=c-yr;
44 
45 % Dibujamos la solucion aproximada
46 figure
47 hold on
48 plot(t,yr,'r')
49 plot(t,xr,'g')
50 hold off
51 legend('Cantidad de productos (Runge-kutta)','Cantidad de reactivos (Runge-kutta)')


Una vez realizados los programas y obtenidas las gráficas correspondientes,se puede concluir que la velocidad de la reacción aumentará hasta llegar a su máximo valor en la intersección de ambas curvas, disminuyendo dicha velocidad hasta agotar la concentración de los reactivos.

caption

Tras realizar el sistema se procede al cálculo correspondiente, representando gráficamente ambos métodos, donde se puede observar una variación de ambos. El punto de cruce de las curvas es diferente en respectivos métodos, siendo en Euler (1.40,1.02) y en Runge Kutta 4 (1.30,1.02), como se puede observar en la gráfica. Se argumenta que el método Euler aporta un desarrollo del sistema más preciso debido a que este método es de orden 1, dando menos error que el método de Runge Kutta 4, de orden 4. Sin embargo, esta diferencia entre ambas gráficas queda reducida a medida que se aumenta el paso h.

3 Reacciones consecutivas de Lotka

Lotka propone un modelo de reacciones consecutivas en las que se consume A, y tras determinadas reacciones se produce B. En las reacciones intermedias se crean dos nuevos compuestos, X e Y que marcarán la velocidad de reacción.

La reacción global queda definida por tres reacciones intermedias que suceden de forma consecutiva. Cabe destacar, por un lado, el hecho de que dos de ellas son autocatalíticas y, por otro lado, que todas ellas están determinadas por diferentes constantes ki que determinarán la velocidad de cada reacción. Las tres reacciones son: [math] \left\{\begin{matrix} A+X→k_1*2*X\\ X+Y→k_2*2*Y\\ Y→k_3*B \end{matrix}\right. [/math]

Para determinar las ecuaciones diferenciales correspondientes a la velocidad de formación (o descomposición) de cada sustancia, se hace un breve estudio de cada una de ellas. En nuestro caso, se ha seguido el orden contrario, y se empieza por la última de ellas para acabar con la primera. Además, es necesario utilizar el Principio de conservación de masas para poder determinar las cuatro ecuaciones correspondientes a los cuatro compuestos presentes. El análisis seguido es el siguiente:

  • Y→k3*B. A partir de Y se formará B, regido por una constante k3. Mediante dicha reacción se podrá definir una velocidad de transformación de B: B'=k3*Y.
  • X+Y→k2*2*Y. En la formación de Y hay que tener en cuenta la parte proporcional a la descomposición de X (k2*X*Y) y descontar la cantidad de Y que se descompone para la formación de B (k3*B). Por tanto, la velocidad de transformación de Y queda determinada por: Y'=k2*Y*X-k3*Y.
  • A+X→k1*2*X. La formación de X tiene una parte proporcional a la descomposición de A (k1*A*X), al que se descuenta aquella parte de X que reacciona para formar el compuesto Y (k2*X*Y). Así pues, se obtendrá una velocidad de transformación regida por X'=k1*A*X-k2*X*Y.
  • Derivado la ecuación correspondiente al Principio de conservación de masas (A+X+Y+B=cte) se obtiene la siguiente ecuación: A'+X'+Y'+B'=0. Sustituyendo los valores obtenidos anteriormente de X', Y' y B' se obtiene la ecuación diferencial que rige la velocidad de A: A'+x'+y'+B'=0 → A'=-x'-y'-B'=-k1*A*x+k2*x*y-k2*x*y+k3*y-k3*y=-k1*A*x → A'=-k1*A*X

Por lo tanto, a partir de este sistema de ecuaciones se puede formar el siguiente PVI asociado, al que hay que añadir las condiciones iniciales, que para este caso son las concentraciones iniciales de cada compuesto: \[ PVI = \begin{cases} X'=k_1*A*X-k_2*X*Y\\Y'=k_2*X*Y-k_3*Y\\B'=k_3*Y\\A'=-k_1*A*X\\X(0)=X_0\\Y(0)=Y_0\\B(0)=B_0\\A(0)=A_0\\ \end{cases}\quad \]


Este PVI queda a continuación resuelto numéricamente mediante el método de Euler. Para ello, las condiciones iniciales elegidas son: A0=4.7, B0=2, X0=2 e Y0=1. El estudio se realiza en los primeros 200 segundos en los que se producen las reacciones. Además, está resuelto usando dos pasos, primero h=0.01 y después h=0.005.


caption
 1 % Datos
 2 A0=4.7;
 3 B0=2;
 4 X0=2;
 5 Y0=1;
 6 k=1.4;
 7 t0=0;
 8 tN=200;
 9 c=2.05;
10 k1=0.1;
11 k2=0.15;
12 k3=0.06;
13 
14 y0=[A0;B0;X0;Y0];
15 ne=4;
16 
17 % Funcion:
18 f=@(t,y)[-k1*y(1)*y(3);k3*y(4);k1*y(1)*y(3)-k2*y(3)*y(4);k2*y(3)*y(4)-k3*y(4)];
19 
20 % Usamos h=0.01:
21 h=0.01;
22 % Discretizamos
23 N=round((tN-t0)/h);
24 
25 % Creamos el vector de tiempo t 
26 t=linspace(t0,tN,N+1);
27 
28 % Matriz de soluciones
29 y=zeros(ne,N+1);
30 y(:,1)=y0;
31 
32 % Aplicamos el método de Euler:
33 for k=1:N
34     y(:,k+1)=y(:,k)+h*f(t(k),y(:,k));
35 end
36 
37 % Dibujamos la solucion aproximada
38 subplot(2,1,1)
39 hold on
40 plot(t,y(1,:),'r')
41 plot(t,y(2,:),'b')
42 plot(t,y(3,:),'g')
43 plot(t,y(4,:),'k')
44 legend('A','B','X','Y');
45 hold off
46 
47 % Usamos h=0.005:
48 h=0.005;
49 % Discretizamos
50 N=round((tN-t0)/h);
51 
52 % Creamos el vector de tiempo t 
53 t=linspace(t0,tN,N+1);
54 
55 % Matriz de soluciones
56 y=zeros(ne,N+1);
57 y(:,1)=y0;
58 
59 % Aplicamos el método de Euler:
60 for k=1:N
61     y(:,k+1)=y(:,k)+h*f(t(k),y(:,k));
62 end
63 
64 % Dibujamos la solucion aproximada
65 subplot(2,1,2)
66 hold on
67 plot(t,y(1,:),'r')
68 plot(t,y(2,:),'b')
69 plot(t,y(3,:),'g')
70 plot(t,y(4,:),'k')
71 legend('A','B','X','Y');
72 hold off


caption
caption










Como se puede observar en los gráficos de detalle, existe una pequeña perturbación entre ambas gráficas, las cuales han sido tomadas arbitrariamente como ejemplo. Esta perturbación es producto de las condiciones iniciales que se introducen en el programa realizado. A pesar de la existencia de dichas discordancias y a la inexistencia de un incremento del error, se podría decir, incluso concluir, que las reacciones estudiadas son estables. Esto se debe a la necesidad de aplicar un gran zoom para poder llegar a apreciar la diferencia entre las curvas de las gráficas una vez superpuestas.

A continuación se procede a hacer un breve estudio comparativo de la reacción que tiene lugar cuando, manteniendo las ecuaciones que rigen las reacciones, se cambia el valor que toma la constante k1 y se mantienen iguales el resto de valores. Para este estudio se han tomado cuatro valores de k1 en el intervalo (0.02,2), y se ha resuelto el PVI numéricamente mediante el método de Heun.

caption
 1 % Datos
 2 A0=4.7;
 3 B0=2;
 4 X0=2;
 5 Y0=1;
 6 k=1.4;
 7 t0=0;
 8 tN=200;
 9 c=2.05;
10 k1=0.1;
11 k2=0.15;
12 k3=0.06;
13 
14 y0=[A0;B0;X0;Y0];
15 ne=4;
16 
17 % Usamos h=0.01:
18 h=0.01;
19 % Discretizamos
20 N=round((tN-t0)/h);
21 
22 % Creamos el vector de tiempo t 
23 t=linspace(t0,tN,N+1);
24 
25 % Matriz de soluciones
26 y=zeros(ne,N+1);
27 y(:,1)=y0;
28 
29 % Aplicamos el método de Heun
30 a=1;
31 b=linspace(0.02,2,4);
32 for j=1:4
33     k1=b(j);
34     f=@(t,y)[-k1*y(1)*y(3);k3*y(4);k1*y(1)*y(3)-k2*y(3)*y(4);k2*y(3)*y(4)-k3*y(4)];
35     for i=1:N
36         K1=f(t(i),y(:,i));
37         K2=f(t(i)+h,y(:,i)+K1*h);
38         y(:,i+1)=y(:,i)+0.5*h*(K1+K2);
39     end
40     subplot(2,2,a)
41     hold on
42     plot(t,y(1,:),'r')
43     plot(t,y(2,:),'b')
44     plot(t,y(3,:),'g')
45     plot(t,y(4,:),'k')
46     legend('A','B','X','Y');
47     title(k1);
48     hold off
49     a=a+1;
50 end


La resolución de este problema se caracteriza por los distintos valores que toma k1 en el intervalo (0.002-2). A raíz de esto se aprecian en las gráficas, también adjuntas, los cambios que se producen en las reacciones. Además, en una de las ecuaciones definidas por Lotka se puede observar la presencia de k1: A+X=k1*2*X. Debido a esta relación, y a partir del estudio de las gráficas, se puede concluir que a medida que el valor de la constante k1 aumenta:

  • El valor de X experimenta un crecimiento instantáneo, que es reducido inmediatamente y cuya curva va transformándose en un asíntota.
  • El valor de Y experimenta una reacción muy parecida a la de X pero presenta una forma menos brusca.
  • El valor de A experimenta una desaparición prácticamente instantánea al inicio la reacción .
  • El valor de B experimenta un crecimiento bastante rápido.

Estas diferencias se pueden encontrar claramente entre cada una de las cuatro gráficas, de forma que se podría afirmar que a partir de k1= 0.68 la reacción mantiene tanto las velocidades como los cambios prácticamente iguales. Para valores pequeños de k1, como en el caso en el que su valor es 0.02 o incluso 0.1 (estudiado anteriormente), se puede apreciar una inestabilidad que comienza en torno a los 120 segundos de los dos compuestos intermedios, X e Y. Esta inestabilidad se reduce conforme el valor de k1 aumenta, hasta hacerse inapreciable cuando la constante alcanza valores de 0.68 y superiores. Además, se observa que el tiempo en el que ocurre la totalidad de la reacción se ve reducidoa medida que el valor de k1 aumenta. Por lo tanto, cuanto mayor sea dicho valor, menos tiempo necesita la reacción global para llevarse a cabo en su totalidad, es decir, que todo el compuesto A se transforme en compuesto B, y los compuestos intermedios X e Y desaparezcan por completo.