Revisión del 23:37 4 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 A y B, que van formando un producto C en lo que se conoce como una reacción bimolecular, es decir, una molécula de A y una de B producen una de C,

                                                                    A + B → C

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 C a lo largo del tiempo

En primer lugar, suponiendo que la reacción es irreversible, es decir, que C no puede volver a formar A y B, se puede expresar la situación mediante el siguiente problema de valor inicial: {█(y^' (t)=k_1 (a_0-y(t) )(b_0-y(t) )@y(0)=0)┤ , t>0

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

f(t,y)=k_1 (a_0-y(t) )(b_0-y(t) ),t>0

∂f/∂y=k_1 (〖-a〗_0 〖-b〗_0+2y)

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

∂f/∂y=k_1 (〖-a〗_0 〖-b〗_0+2y)

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 C, este también se descompone de nuevo en A y B. Por lo tanto, nos quedaría una ecuación diferencial con la siguiente forma:

y^' (t)=k_1 (a_0-y(t) )(b_0-y(t) )-Cy(t),t>0

Si nos fijamos, la modificación es sencillamente restar y(t) multiplicada por una constante, la cual refleja qué porcentaje de C se va descomponiendo a medida que se da, simultáneamente, la reacción A + B → C. 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 C le asociaremos el siguiente dominio:

C∈(0,1)

Obsérvese que C 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 C 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 A (a_0=3 mol/L) y B (b_0=1 mol/L) como de C, para una k_1=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, a y b
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 A y B y del producto C, se ha diseñado un programa que, mediante diferencia de valores sucesivos de la variable y(t), 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 t→∞:

%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
        limiteC=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
        limiteA=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
        limiteB=b_(i+1)
        break
    end
end

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

limiteA=2.0036≈2 mol/L limiteB=0.0036≈0mol/L limiteC=0.9964≈1mol/L

Lo cual es lógico, ya que A y B se combinan en una proporción de 1:1, luego B, que tiene una concentración inicial de b_0=1, será el reactivo limitante (el que desaparecerá por completo), mientras que A, con a_0=3, acabará con una concentración molar final igual a 2.

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 Trapecios (segundo orden) y de Runge-Kutta (cuarto orden). Desgraciadamente, no es posible comparar los errores cometidos por cada método con la expresión real de la concentración de C, ya que la ecuación diferencial no es de fácil resolución analítica, luego resulta imposible alcanzar una expresión explícita de la concentración en términos de nuestros datos conocidos: a_0 〖,b〗_0 〖 ,k〗_1,t (tiempo). Sin embargo, podemos confiar con bastante seguridad en el método sobre todo de Runge-Kutta, debido a su cuarto orden.

Para los trapecios, al ser un método implícito, es necesario despejar la componente y_(n+1):

y_(n+1)=((1+h/2∙k_1∙(a_0+b_0 ) )-√((1+h/2∙k_1∙(a_0+b_0 ) )^2-2∙h∙k_(1∙) (y_n+h/2∙k_1∙(a_0-y_n )∙(b_0-y_n )+〖h/2∙k〗_1∙a_0∙b_0 ) ))/(k_1*h)

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:

%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 dos, en azul las soluciones de los Trapecios, y en negro las
%de Runge-Kutta. 
figure(2)
hold on
plot(t,y2)
plot(t,a2)
plot(t,b2)
plot(t,y3,'k')
plot(t,a3,'k')
plot(t,b3,'k')
hold off

Que proporciona el resultado que se puede observar en la primera gráfica adjunta. Como podemos ver, es prácticamente indistinguible una gráfica de la otra. Sin embargo, si la juntamos con las gráficas proporcionadas por el método de Euler, que se representa ahora en rojo, la diferencia es palpable (ver gráfica 2), lo cual evidencia la peor precisión del método de Euler.

Reacción compleja

Gráfico 2

5 Reacción consecutiva

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

                                                                  A + B → k1C → k2D

De tal modo que los reactivos A y B se combinan en la proporción 1:1 para formar compuesto C, el cual, a su vez, se va descomponiendo para formar compuesto D. 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 C es proporcional a las concentraciones existentes de A y de B y del mismo modo la de formación de la sustancia D lo es a la concentración de C que se ha ido formando.

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

[math] \begin{cases} x=cosh(u)cos(v) \\ y=sinh(u)sin(v) \end{cases} [/math]

[math] \begin{cases} \ltmath\gty'\ltsub\gt1\lt/sub\gt(t)=k\ltsub\gt1\lt/sub\gt(a\ltsub\gt0\lt/sub\gt − y\ltsub\gt1\lt/sub\gt(t))(b\ltsub\gt0\lt/sub\gt − y\ltsub\gt1\lt/sub\gt(t))[/math] \\ [math]y'\ltsub\gt2\lt/sub\gt(t)=k\ltsub\gt2\lt/sub\gty\ltsub\gt1\lt/sub\gt(t)[/math] \\ [math]y\ltsub\gt1\lt/sub\gt(0)=0, y\ltsub\gt2\lt/sub\gt(0)=0[/math] \end{cases} </math>

Donde a₀ y b₀ representan las concentraciones iniciales de A y B respectivamente, k₁ es la constante de proporcionalidad que rige la primera reacción y k₂ la que rige la segunda reacción o reacción consecutiva. y₁’(t) e y₂’(t) representan las velocidades de las reacciones, es decir, las velocidades con que se forman los productos C y D. y₂(t) representa la concentración del producto D pero es importante apreciar que en este caso y₁(t) no nos da el valor de la concentración de C directamente, sino que ésta será y₁(t)- y₂’(t), la cantidad que se forma en la primera reacción menos la que se emplea para formar el producto D 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 k₂ mayor que k₁. Y otro en el cual se considera que es más lenta, lo que implica que k₂ es menor que k₁. Concretamente vamos a tomar los siguientes valores para resolver el problema numéricamente: a₀=3, b₀=1, k₁=1, k₂= . 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 k₁=5 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 a medida que avanza el tiempo, las reacciones C y D aumentan sus concentraciones. Antes de que pase 1 segundo, la concentración C alcanza un máximo (de unos 0.3 molar) y va disminuyendo progresivamente hasta alcanzar un valor de 0 molar. La concentración D sigue creciendo hasta llegar a 1 molar. Una vez que tanto C como D alcanzan estos valores, se estabiliza la reacción.

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

CASO k₁=5 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 k₁=1, k₂=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 C es más rápida que D, 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. C 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. D 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 C 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.