Diferencia entre revisiones de «Reacciones de autocatalisis Grupo 9A»

Revisión del 15:03 5 mar 2015

1 Introducción y planteamiento

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 en [math] D \cap B((t_{0},y_{0}),r) [/math].

Observamos que nuestra funcion: [math] f(t,y(t))=f= 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))=f= 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 Ecuación diferencial

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

2.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

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 parábola cuyo máximo corresponde a concentraciones iguales de A y de B?.

2.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.3 Método de Runge-Kutta

Resolvemos mediante Rounge Kutta

%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.

3 Reacción consecutiva propuesta por Lotka

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

3.2 Método de Runge-Kutta

3.2.1 Interpretación

3.3 Método de Euler

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',1)
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|thumb

3.4 Método de Heun

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

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

Aplicando el metodo mediante el programa Matlab:

4 Conclusión