Reacciones de autocatalisis Grupo 9A

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 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}*(cte - 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}*(cte - y(t))*y(t) [/math] Y: [math] \displaystyle{ \delta f \over \delta y} = k_{1}(cte - 2y) [/math]

Ambas funciones son de clase [math]C^\infty[/math] podemos afirmar entonces, 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

%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); %TAMBIEN PUEDO USAR zeros(1,length(t))
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
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');

2.2 Método del Trapecio

Resolvemos la ecuacion mediante el metodo del trapecio.

%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');

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');