Revisión actual del 13:53 6 mar 2015

1 Introducción

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 [math]A[/math] y [math]B[/math], que van formando un producto [math]C[/math] en lo que se conoce como una reacción bimolecular, es decir, una molécula de [math]A[/math] y una de [math]B[/math] producen una de [math]C[/math].:

                                                                    [math]A + B → C[/math]

Supondremos también que se satisface 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.

2 Concentración del reactivo [math]C[/math] a lo largo del tiempo

En primer lugar, suponiendo que la reacción es irreversible, es decir, que [math]C[/math] no puede volver a formar [math]A[/math] y [math]B[/math], se puede expresar la situación mediante el siguiente problema de valor inicial, cuando t>0.:

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

Donde [math] {a_{0}}[/math] y [math]{b_{0}}[/math] representan las concentraciones iniciales de A y B respectivamente y [math]{k_{1}}[/math] es la constante de proporcionalidad que indica la ley de acción de masas. [math]y’(t)[/math] representa la velocidad de reacción, es decir, la velocidad con que se forma el producto [math]C[/math] ; mientras que la solución de la ecuación diferencial, [math]y(t)[/math] , representa la concentración de [math]C[/math], la cantidad de compuesto que se ha formado. Para comprobar que la solución de este problema de Cauchy es única, comenzamos definiendo.:

[math]f(t,y)={k_{1}}({a_{0}}-y(t))({b_{0}}-y(t))=0,   t\gt0 [/math]

Observamos que [math]f(t,y)[/math] es continua en el intervalo [math] I=(0,∞)∩B((0,0),r\gt0)[/math], por lo que ya aseguramos la existencia de al menos una solución. Examinando la derivada parcial.:

[math]\frac{∂f}{∂y}=k_1(-a_0-b_0+2y)[/math]

Se deduce que al ser un polinomio no va a haber problemas de continuidad para ninguna bola. De esta forma, podemos asegurar que existe una única solución.

2.1 Concentración del reactivo C a lo largo del tiempo en un proceso reversible

Ahora bien, ¿qué pasaría si la reacción fuese reversible? Con esto lo que queremos indicar es que además de formarse [math]C [/math], este también se descompone de nuevo en [math]A [/math] y [math]B [/math]. Por lo tanto, nos quedaría una ecuación diferencial con la siguiente forma, cuando [math]t\gt0 [/math]. :

[math] y'(t)={k_{1}}({a_{0}}-y(t))({b_{0}}-y(t))-Cy(t)[/math]

Si nos fijamos, la modificación es sencillamente restar [math]y(t) [/math] multiplicada por una constante, la cual refleja qué porcentaje de [math]C [/math] se va descomponiendo a medida que se da, simultáneamente, la reacción [math]A + B → C [/math]. Sin embargo, esta constante no puede tener cualquier valor, ya que en una reacción se cumple la conservación de masa. Por ello, a [math]C [/math] le asociaremos el siguiente dominio::

[math]C∈(0,1)[/math]

Obsérvese que [math]C [/math] tiene que ser mayor que 0 para que sea reversible, pero al mismo tiempo inferior que 1, ya que si no se estaría descomponiendo más compuesto [math]C [/math] del que realmente se ha formado, en otras palabras, estaría apareciendo masa de la nada.

3 Método de Euler

Una vez llegados a este punto, se procede a ver una interpretación gráfica de la reacción. Para ello, se refleja en la gráfica (en un intervalo de 2 segundos), cómo varían las concentraciones tanto de [math]A[/math] ([math]a_0[/math]=3 mol/L) y [math]B[/math] ([math]b_0[/math]=1 mol/L) como de [math] C[/math], para una [math]k_1[/math]=1 mol/s. Se comienza, por el momento, aplicando el método de Euler. Posteriormente se utilizarán métodos más precisos.

% Euler
clear all
% datos del problema
t0=0;
tN=2;
y0=0;
h=0.1;
k1=1;
a0=3;
b0=1;
% Definimos la variable independiente
t=t0:h:tN;
% Guardamos los valores de la solución aproximada en el vector "y1" y
% establecemos los vectores de concentraciones A y B
y1=zeros(1,length(t)); 
y1(1)=y0;
a1=zeros(1,length(t));
a1(1)=a0;
b1=zeros(1,length(t));
b1(1)=b0;
% A continuación, calculamos las componentes de los vectores y1, a1 y b1
for i=1:length(t)-1
    y1(i+1)=y1(i)+h*(k1*(a0-y1(i))*(b0-y1(i)));    %Euler
    a1(i+1)=a0-y1(i+1); % evolución de la concentración de A
    b1(i+1)=b0-y1(i+1); % evolución de la concentración de B
end
%grafico 
figure(1)
hold on
plot(t,y1,'g')
plot(t,a1,'r')
plot(t,b1,'b')
legend('concentración de C','concentración de A','concentración de B','location','best')
hold off 
%en verde se muestra la concentración de C, en rojo la de A, y en azul la de B

Reacción compleja

Aunque con observar los resultados ya podemos decir con casi absoluta seguridad cuáles son los límites de las concentraciones de cada reactivo [math]A[/math] y [math]B[/math] y del producto [math]C[/math], se ha diseñado un programa que, mediante diferencia de valores sucesivos de la variable [math]y(t)[/math], nos devuelve el valor de la función cuando dicha diferencia es muy pequeña (se ha puesto una diferencia de 0.001), es decir, cuando tenemos una asíntota al tender [math]t→∞[/math]:

%Estudio de las asíntotas
t1=t0:h:10;%se aumenta la longitud de t, y se estima que se alcanzarán las asíntotas antes de t=10
y_=zeros(1,length(t1)); %se define un vector y_ idéntico a y1 pero con t hasta 10
a_=y_;
b_=y_;
y_(1)=y0;
a_(1)=a0;
b_(1)=b0;
for i=1:length(t1)-1
    y_(i+1)=y_(i)+h*(k1*(a0-y_(i))*(b0-y_(i)));
    diferencia=y_(i+1)-y_(i);
    if abs(diferencia)<=0.001
        limC=y_(i+1)
        break
    end
end
for i=1:length(t1)-1
    y_(i+1)=y_(i)+h*(k1*(a0-y_(i))*(b0-y_(i)));    
    a_(i+1)=a0-y_(i+1);
    diferencia2=a_(i+1)-a_(i);
    if abs(diferencia2)<=0.001
        limA=a_(i+1)
        break
    end
end
for i=1:length(t1)-1
    y_(i+1)=y_(i)+h*(k1*(a0-y_(i))*(b0-y_(i)));    
    b_(i+1)=b0-y_(i+1);
    diferencia3=b_(i+1)-b_(i);
    if abs(diferencia3)<=0.001
        limB=b_(i+1)
        break
    end
end

Que nos devuelve unos valores de las concentraciones para el límite de:

[math]\lim A=2,0036≈2 mol/L[/math]: [math]\lim B=0,0036≈0mol/L[/math]: [math]\lim C=0,9964≈1mol/L[/math]:

Lo cual es lógico, ya que [math]A[/math] y [math]B[/math] se combinan en una proporción de 1:1, luego [math]B[/math], que tiene una concentración inicial de [math]b_0[/math]=1, será el reactivo limitante (el que desaparecerá por completo), mientras que [math]A[/math], con [math]a_0[/math]=3, acabará con una concentración molar final igual a 2, y [math]C[/math], del que no habia nada,alcanza una concentraccion de [math]1mol/L[/math].

4 Método del Trapecio y Método de Runge-Kutta

A continuación, para obtener con más precisión las representaciones gráficas de las respectivas concentraciones, se opta por utilizar los métodos de los [math]Trapecios[/math] (segundo orden) y de [math]Runge-Kutta[/math] (cuarto orden). Además, se ha obtenido la expresión explícita de la solución de nuestro problema de valor inicial de manera analítica, mediante variables separables y fracciones simples, resultando la función real:: [math]y_e(t)=\frac{3e^{2t}-3}{3e^{2t}-1}[/math] De manera que las concentraciones de A y de B en cada valor de t quedarían como:: [math]a_e(t)=a_0-y_e(t)[/math]: [math]b_e (t)=b_0-y_e (t)[/math]: Para los trapecios, al ser un método implícito, es necesario despejar la componente [math]y_{n+1}[/math]:

[math]y_{n+1}=\frac{(1+\frac{h}{2}∙k_1∙(a_0+b_0))-\sqrt{(1+\frac{h}{2}∙k_1∙(a_0+b_0))^2 -2∙h∙k_1(y_n+\frac{h}{2}∙k_1∙(a_0-y_n )∙(b_0-y_n)+\frac{h}{2}∙k_1∙a_0∙b_0) }}{k_1∙h}[/math]

Cabe destacar que este programa parte de los datos ya establecidos en el programa de Euler para t0, tN, y0, a0, b0, k1 y h:

%Función real
ye=(3*exp(2.*t)-3)./(3*exp(2.*t)-1);
ae=a0-ye;
be=b0-ye;

%Trapecios
y2=zeros(1,length(t)); 
a2=y2;
b2=y2;
y2(1)=y0;
a2(1)=a0;
b2(1)=b0;
for i= 1:length(t)-1
    y2(i+1)=((1+0.5*h*k1*(a0+b0))-sqrt((1+0.5*h*k1*(a0+b0))^2-2*h*k1*(y2(i)+0.5*h*k1*(a0-y2(i))*(b0-  y2(i))+0.5*h*k1*a0*b0)))/(k1*h);
    a2(i+1)=a0-y2(i+1);
    b2(i+1)=b0-y2(i+1);
end

%Runge-Kutta
y3=zeros(1,length(t)); 
a3=y3;
b3=y3;
y3(1)=y0;
a3(1)=a0;
b3(1)=b0;
for i=1:length(t)-1 
    r1=k1*(a0-y3(i))*(b0-y3(i));
    r2=k1*(a0-y3(i)-0.5*r1*h)*(b0-y3(i)-0.5*r1*h);
    r3=k1*(a0-y3(i)-0.5*r2*h)*(b0-y3(i)-0.5*r2*h);
    r4=k1*(a0-y3(i)-r3*h)*(b0-y3(i)-r3*h);
    y3(i+1)=y3(i)+h/6*(r1+2*r2+2*r3+r4);
    a3(i+1)=a0-y3(i+1);
    b3(i+1)=b0-y3(i+1);
end

%se pintan las cuatro gráficas, con la función real en verde y los puntos con asteriscos, 
%Euler en rojo, los Trapecios en azul, y Runge-Kutta en negro
figure(2)
hold on

plot(t,ye,'g*-')
plot(t,y2)
plot(t,y3,'k')
plot(t,y1,'r')
legend('Función real','Trapecios', 'Runge-Kutta', 'Euler','Location','best');
plot(t,ae,'g*')
plot(t,a2)
plot(t,a3,'k')
plot(t,a1,'r')
plot(t,be,'g*')
plot(t,b2)
plot(t,b3,'k')
plot(t,b1,'r')
hold off

Los resultados del programa se muestran en las dos siguientes figuras. Como podemos observar, el método de Euler, en Rojo, rápidamente se separa de los demás y de la función real, lo que es debido a que es un método de orden 1 y su precisión es bastante reducida. Para analizar el comportamiento de los otros dos métodos, se ha hecho un zoom en una zona aleatoria de la gráfica, y se muestra en la figura 2. Lo lógico, debido a que es de orden 2, es que el método de los trapecios sea el siguiente que se aleje de la solución real, y esto es precisamente lo que se evidencia en la gráfica, pues como vemos, la línea Azul se separa de la negra y de la Verde. Es por tanto que el método de Runge-Kutta es el más preciso de todos, haciendo honor a su cuarto orden, y como podemos observar, las líneas Verde y negra son casi indistinguibles.

Es importante también atender a la escala que se muestra en la figura 2, que nos da una idea de los rangos en el eje de ordenadas de los que estamos hablando. Cabe decir que estamos trabajando con un paso de h=0.1, lo que significa que con un paso menor, estas diferencias serian mucho menos perceptibles.

Fifura 1

Figura 2

5 Reacción consecutiva

Consideremos ahora que se produce una reacción consecutiva de la forma.:

                                                                  [math]A + B → k_1C → k_2D[/math]

De tal modo que los reactivos [math]A[/math] y [math]B[/math] se combinan en la proporción 1:1 para formar compuesto [math]C[/math], el cual, a su vez, se va descomponiendo para formar compuesto [math]D[/math]. Al igual que se ha hecho antes al analizar las otras reacciones, se considerará que se cumple la ley de acción de masas. Es decir, la velocidad con la que se forma compuesto [math]C[/math] es proporcional a las concentraciones existentes de [math]A[/math] y de [math]B[/math] y del mismo modo la de formación de la sustancia [math]D[/math] lo es a la concentración de [math]C[/math] que se ha ido formando.

Esta reacción consecutiva está representada por el siguiente problema de valor inicial::

[math] \left . \begin{matrix} y'_1(t)={k_{1}}({a_{0}}-y_1(t))({b_{0}}-y_1(t))\\ y'_2(t)={k_{2}}(y_1(t)-y_2(t))\\ y_1(0)=0\\y_2(0)=0 \end{matrix} \right \} [/math] Donde [math]a_0[/math] y [math]b_0[/math] representan las concentraciones iniciales de [math]A[/math] y [math]B[/math] respectivamente, [math]k_1[/math] es la constante de proporcionalidad que rige la primera reacción y [math]k_2[/math] la que rige la segunda reacción o reacción consecutiva. [math]y_1’(t)[/math] e [math]y_2’(t)[/math] representan las velocidades de las reacciones, es decir, las velocidades con que se forman los productos [math]C[/math] y [math]D[/math]. [math]y_2(t)[/math] representa la concentración del producto [math]D[/math] pero es importante apreciar que en este caso [math]y_1(t)[/math] no nos da el valor de la concentración de [math]C[/math] directamente, sino que ésta será [math]y_1(t)[/math]- [math]y_2(t)[/math], la cantidad que se forma en la primera reacción menos la que se emplea para formar el producto [math]D[/math] en la segunda.

Este problema se va a resolver para dos casos diferentes. Uno en el cual se considera que la segunda reacción es más rápida que la primera, lo que supone tomar un valor de la constante [math]k_2[/math] mayor que [math]k_1[/math]. Y otro en el cual se considera que es más lenta, lo que implica que [math]k_2[/math] es menor que [math]k_1[/math]. Concretamente vamos a tomar los siguientes valores para resolver el problema numéricamente; [math]a_0[/math]=3 , [math]b_0[/math]=1, [math]k_1[/math]=1, [math]k_2[/math]=5, [math]k_2[/math]=1/5, .En ambos casos se resuelve para los primeros 10 segundos.

5.1 La reacción D es más rápida que la reacción C

CASO [math]k_1=5 [/math] h=0.1

%datos
y0=[0;0];
z0=[0;0];
t0=0;
tN=10;
h=0.1;
 
%variable independiente t
t=t0:h:tN;
 
%matriz de las soluciones
y=zeros(2,length(t)); %euler
y(:,1)=y0;
 
z=zeros(2,length(t)); %runge-kutta
z(:,1)=y0;
 
for i=1:length(t)-1
 %euler
 y(:,i+1)=y(:,i)+h*[1*(3-y(1,i))*(1-y(1,i)); 5*(y(1,i)-y(2,i))];
 %runge-kutta
 K1=[1*(3-z(1,i))*(1-z(1,i)); 5*(z(1,i)-z(2,i))];
 K2=[1*(3-(z(1,i)+h/2*K1(1)))*(1-(z(1,i)+h/2*K1(1))); 5*((y(1,i)+h/2*K1(2))-(y(2,i)+h/2*K1(2)))];
 K3=[1*(3-(z(1,i)+h/2*K2(1)))*(1-(z(1,i)+h/2*K2(1))); 5*((y(1,i)+h/2*K2(2))-(y(2,i)+h/2*K2(2)))];
 K4=[1*(3-(z(1,i)+h*K3(1)))*(1-(z(1,i)+h*K3(1))); 5*((y(1,i)+h*K3(2))-(y(2,i)+h*K3(2)))];
 z(:,i+1)=z(:,i)+h/6*(K1+2*K2+2*K3+K4);
end

% se calcula la concentración de C=y1t)-y2(t)
Ay=y(1,:)-y(2,:); %euler
Az=z(1,:)-z(2,:); %runge-kutta
 
figure(1)
hold on
plot(t,Ay) %concentración de C por euler
plot(t,y(2,:),'r') %concentración de D por euler
plot(t,Az,'k') %concentración de C por runge-kutta
plot(t,z(2,:),'g') %concentración de D por runge-kutta
legend('concentración de C por euler','concentración de D por euler','concentración de C por runge-kutta','concentración de D por runge-kutta','Location','best')
hold off

Reacción compleja

Se observa que en en los primeros instantes aumentan las concentraciones de [math]C[/math] y [math]D[/math] aumentan sus concentraciones. Antes de que pase 1 segundo, la concentración [math]C[/math] alcanza un máximo (de unos 0.3 molar) y va disminuyendo progresivamente hasta alcanzar un valor de 0 molar. La concentración [math]D[/math] sigue creciendo hasta llegar a 1 molar. Una vez que tanto [math]C[/math] como [math]D[/math] alcanzan estos valores, se estabiliza la reacción.

Como hemos impuesto que la reacción [math]D[/math] es más rápida que [math]C[/math], se observa que [math]C[/math] empieza creciendo más rápidamente, pero como reacciona más despacio que [math]D[/math] su concentración disminuye progresivamente, mientras que la de [math]D[/math] aumenta.

CASO [math]k_2=5 [/math] h=0.3

%datos
y0=[0;0];
z0=[0;0];
t0=0;
tN=10;
h=0.3;
 
%variable independiente t
t=t0:h:tN;
 
%matriz de las soluciones
y=zeros(2,length(t)); %euler
y(:,1)=y0;
 
z=zeros(2,length(t)); %runge-kutta
z(:,1)=y0;
 
for i=1:length(t)-1
 %euler
 y(:,i+1)=y(:,i)+h*[1*(3-y(1,i))*(1-y(1,i)); 5*(y(1,i)-y(2,i))];
 %runge-kutta
 K1=[1*(3-z(1,i))*(1-z(1,i)); 5*(z(1,i)-z(2,i))];
 K2=[1*(3-(z(1,i)+h/2*K1(1)))*(1-(z(1,i)+h/2*K1(1))); 5*((y(1,i)+h/2*K1(2))-(y(2,i)+h/2*K1(2)))];
 K3=[1*(3-(z(1,i)+h/2*K2(1)))*(1-(z(1,i)+h/2*K2(1))); 5*((y(1,i)+h/2*K2(2))-(y(2,i)+h/2*K2(2)))];
 K4=[1*(3-(z(1,i)+h*K3(1)))*(1-(z(1,i)+h*K3(1))); 5*((y(1,i)+h*K3(2))-(y(2,i)+h*K3(2)))];
 z(:,i+1)=z(:,i)+h/6*(K1+2*K2+2*K3+K4);
end

% se calcula la concentración de C=y1t)-y2(t)
Ay=y(1,:)-y(2,:); %euler
Az=z(1,:)-z(2,:); %runge-kutta
 
figure(1)
hold on
plot(t,Ay) %concentración de C por euler
plot(t,y(2,:),'r') %concentración de D por euler
plot(t,Az,'k') %concentración de C por runge-kutta
plot(t,z(2,:),'g') %concentración de D por runge-kutta
legend('concentración de C por euler','concentración de D por euler','concentración de C por runge-kutta','concentración de D por runge-kutta','Location','best')
hold off

Reacción compleja

Repetimos el programa pero para un tamaño de paso mayor, de h=0.3, y observamos que no es suficientemente preciso, sobre todo al principio, no se puede analizar bien lo que está pasando en la reacción.

5.2 La reacción D es más lenta que la reacción C

CASO [math]k_1 [/math]=1, [math]k_2 [/math]=1/5

%datos
y0=[0;0];
z0=[0;0];
t0=0;
tN=10;
h=0.1;
 
%variable independiente t
t=t0:h:tN;
 
%matriz de las soluciones
y=zeros(2,length(t)); %euler
y(:,1)=y0;
 
z=zeros(2,length(t)); %runge-kutta
z(:,1)=y0;
 
for i=1:length(t)-1
 %euler
 y(:,i+1)=y(:,i)+h*[1*(3-y(1,i))*(1-y(1,i)); 1/5*(y(1,i)-y(2,i))];
 %runge-kutta
 K1=[1*(3-z(1,i))*(1-z(1,i)); 1/5*(z(1,i)-z(2,i))];
 K2=[1*(3-(z(1,i)+h/2*K1(1)))*(1-(z(1,i)+h/2*K1(1))); 1/5*((y(1,i)+h/2*K1(2))-(y(2,i)+h/2*K1(2)))];
 K3=[1*(3-(z(1,i)+h/2*K2(1)))*(1-(z(1,i)+h/2*K2(1))); 1/5*((y(1,i)+h/2*K2(2))-(y(2,i)+h/2*K2(2)))];
 K4=[1*(3-(z(1,i)+h*K3(1)))*(1-(z(1,i)+h*K3(1))); 1/5*((y(1,i)+h*K3(2))-(y(2,i)+h*K3(2)))];
 z(:,i+1)=z(:,i)+h/6*(K1+2*K2+2*K3+K4);
end

% se calcula la concentración de C=y1t)-y2(t)
Ay=y(1,:)-y(2,:); %euler
Az=z(1,:)-z(2,:); %runge-kutta
 
figure(1)
hold on
plot(t,Ay) %concentración de C por euler
plot(t,y(2,:),'r') %concentración de D por euler
plot(t,Az,'k') %concentración de C por runge-kutta
plot(t,z(2,:),'g') %concentración de D por runge-kutta
legend('concentración de C por euler','concentración de D por euler','concentración de C por runge-kutta','concentración de D por runge-kutta','Location','best')
hold off

Reacción compleja

Si ahora imponemos que la reacción [math]C [/math] es más rápida que [math]D [/math], obtenemos lo siguiente. Se parece a la gráfica que obtuvimos en el programa anterior, pero más extendida en el tiempo, más lenta. [math]C [/math] alcanza una concentración máxima de unos 0.8 molar antes de empezar a decrecer a medida que avanza el tiempo y acercándose a un valor de 0 molar. [math]D [/math] crece más lentamente que en la gráfica anterior, pero se aproxima al mismo valor, 1 molar. Se observa que aunque hagamos la reacción [math]C [/math] más rápida, tenemos las mismas concentraciones finales, pero tardan más tiempo en obtenerse, ya que la reacción tarda más en estabilizarse.