Revisión del 15:02 6 mar 2015

1 Introducción y explicación del PVI

El problema que vamos a tratar es un problema que aborda el tema de Reacciones Complejas. Para ello consideramos, según el enunciado dado, una reacción química irreversible en una solución bien mezclada. Tal y como se nos indica inicialmente en los datos del problema, suponemos también que dicha reacción dada ocurre a un volumen y temperatura constantes. La reacción química con la que trabajaremos responde a la estructura que a continuación se presenta:

                                                      A + B →[math]k_{1}[/math] C 

Ésta nos indica que tenemos inicialmente 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, tal y como expresamos en el esquema de la reacción.

Otro dato también a suponer es que dicha reacción química satisface la Ley de Acción de Masas, ley que establece que la velocidad de reacción es proporcional al producto de las concentraciones de los reactivos.

Se nos plantea la siguiente ecuación diferencial en el enunciado que, supuestamente, corresponde con la variación de la concentración de C a lo largo del tiempo:

                                                            [math]
   \left . 
 \begin{matrix} 
y'(t) = k_1 · (a_0 - y(t)) · (b_0 - y(t)),\\
  t\gt0
 \end{matrix}
   \right \}
[/math]

Tal y como podemos comprobar, si que podemos obtener la variación de la concentración del producto C a lo largo del tiempo a partir de dicha ecuación diferencial. Esto es así debido a que nuestra reacción química se trata de una reacción bimolecular y satisface a su vez la ley de acción de masas, conceptos que ya se han explicado.

Una vez conocidos todos estos datos, planteamos ahora nuestro Problema de Valor Inicial con el que trabajaremos en los siguientes apartados, que es el siguiente:

                                                            [math]
   \left . 
      \begin{matrix} 
                                                                y'(t) = k_1 · (a_0 - y(t)) · (b_0 - y(t)) \\
       A(t) = a_0 - y(t)\\B(t) = b_0 - y(t)
 \end{matrix}
   \right \}
[/math]

donde [math]a_0[/math]=concentración inicial del reactivo A, [math]b_0[/math]=concentración inicial del reactivo B y [math]k_1[/math]= constante de proporcionalidad de la velocidad de variación de la concentración del producto C.

Como hemos visto, la concentración de A [math](A(t) = a_0 - y(t))[/math] y la concentración de B [math](B(t) = b_0 - y(t))[/math] son inversamente proporcional a la concentración de C [math](C(t)=y(t))[/math]: las concentraciones de los dos reactivos A y B parten de su concentración inicial [math]a_0[/math] y [math] b_0[/math] van disminuyendo a medida que la concentración del producto C aumenta con el tiempo.

1.1 ¿Qué sucedería si el proceso fuese reversible?

Hasta ahora hemos planteado el PVI suponiendo el proceso irreversible. En cambio, si lo suponemos reversible, únicamente tendríamos que añadir un componente negativo en la velocidad de reacción del compuesto C, proporcional a la concentración de C, destinada a la creación de A y B.

Concentración reactivos k1=1 k2=2

Concentración reactivos k1=1 k2=3

En las gráficas vemos cómo las reactivos acaban estabilizándose.

Si [math]k_2[/math]=∞ ó [math]k_1[/math]=0, no habría reacción. En el primer caso se crearía A y B más rápido que su destrucción. En el segundo nunca se llegaría a destruir nada de A y B.

Si [math]k_2[/math]=0 ó [math]k_1[/math]=∞, la reacción no sería reversible. En el primer caso no se crearía nada de A y B por lo que sólo se crearía C. En el segundo caso la creación de C es mucho mayor que la de A y B.

2 Resolución Numérica del problema asociado

Ahora, nos disponemos a resolver numéricamente el problema de valor inicial descrito en el apartado anterior, suponiendo:

           [math]
          \left . 
      \begin{matrix} 
a_0 = 3 mol/l \\  
b_0 = 1 mol/l \\
k_1 = 1 mol/s 
     \end{matrix}
        \right \}
           [/math]

Utilizaremos, en este caso, el Método de Euler como método numérico para la resolución del problema, eligiendo un paso [math]h[/math]=0.1 en los primeros 2 segundos. El código Matlab resultante tras la aplicación numérica de este método a nuestro problema de valor inicial es el siguiente:

%PVI por Euler con una h = 0.1
clear all

%condiciones dadas
y0 = 0;
t0 = 0;
tN = 2;
h = 0.1;
k = 1;
%calculamos N
N = (tN - t0)/h;
%definimos la variable t
t = linspace(t0,tN,N+1);
%funcion y
y = zeros(1,N+1);

y(1)= y0;

%solucion obtenida
for i=1:N
    y(i+1)=y(i)+ h*(k*(3 - y(i))*(1 - y(i)));
end

%Concentracion de los reactivos
ya = zeros(1,N+1);
yb = zeros(1,N+1);

ya = 3.- y;
yb = 1.- y; 


%solucion grafica
hold on
plot(t,y,'linewidth',2)
plot(t,ya,'g','linewidth',2)
plot(t,yb,'r','linewidth',2)
legend('Concentracion de C','Concentracion de A','Concentracion de B','Location','best');
hold off

Una vez establecido dicho código Matlab, dibujamos en una misma gráfica tanto las concentraciones de los dos reactivos como la del producto, en función de su evolución en el tiempo, quedando las mismas de la siguiente manera:

Gráfica de las dos concentraciones y el producto. Utilizando el método de Euler y con una h = 0.1

Observamos ahora en la gráfica la variación de las tres concentraciones según el paso del tiempo. Como era de esperar, las concentraciones de A y B disminuirán exponencialmente, mientras que la concentración de C, por el contrario, irá aumentando de manera logarítmica. Esto es debido a que los reactivos necesitan irse consumiendo para poder dar lugar a C, que va creciendo y se va formando más cantidad del mismo. El punto inicial de las tres concentraciones para un tiempo [math]t[/math]=0 representado en la gráfica nos indica la concentración de cada elemento que tenemos en el instante inicial, antes de que comience a darse lugar la reacción. Este valor coincide con los datos que se nos planteaban en el problema de valor inicial:

      [math]
          \left . 
      \begin{matrix} 
a_0 = 3 mol/l \\  
b_0 = 1 mol/l \\
k_1 = 1 mol/s 
     \end{matrix}
        \right \}
           [/math]

Las gráficas de las funciones de las concentraciones tienden a valores diferentes a medida que se avanza en el tiempo. Observamos que la gráfica que representa la variación de la concentración de A con el tiempo (color verde) tiende a 2, la de B (color rojo) tiende a 0 y la de C (color azul) tiende a 1. Esto es así debido a tratarse de una reacción bimolecular, que, como ya se ha explicado anteriormente, significa que una molécula de A y una de B producen una de C. Por tanto, esa es la razón por la que la concentración de A pasa de [math]a_0[/math]=3 mol/l a [math]A(t)[/math]=2 mol/l en función del tiempo t transcurrido, y así ocurre también con B y C, disminuyendo B de [math]b_0[/math]=1 mol/l a [math]B(t)[/math]=0 mol/l y, por el contrario, aumentando C de [math]c_0[/math]=0 mol/l a [math]C(t)[/math]=1 mol/l.

Cuando t tiende a infinito:

Una vez analizada cómo es la variación de las concentraciones en función del transcurso del tiempo, resolveremos numéricamente el mismo sistema para un tiempo lo suficientemente grande como para ser capaces de establecer una aproximación de los límites de las concentraciones cuando t → ∞. Para ello, lo único a modificar en el programa anterior es el dato del tiempo, que lo aumentaremos razonablemente para lograr ver hacia dónde tienden las gráficas con el paso del tiempo, mientras que el resto de datos se mantienen.

Ahora, la gráfica del Método de Euler para un tiempo t → ∞ queda de la siguiente manera:

Gráfica con un tiempo mayor para observar los límites de las concentraciones. Utilizando el método de Euler

Como ya se podía intuir en el anterior gráfico a este, las gráficas de las funciones de las concentraciones tienen sus límites para un tiempo t → ∞ determinados según la reacción bimolecular en la que están participando. Por tanto, tal y como ya se ha comentado, la concentración de A tiende a 2 mol/l, la de B tiende a 0 mol/l y la de C tiende a 1 mol/l, cuando la reacción tiende a un tiempo t → ∞.

2.1 Método del Trapecio

Resolvemos ahora nuestro mismo problema de valor inicial por otros dos métodos numéricos de aproximación, que son el Método del Trapecio y el Método de Runge Kutta. De esta manera, podemos comparar la solución al problema obtenida por los tres métodos numéricos empleados, debido a que, al ser métodos aproximados, las soluciones serán más precisas en función del método empleado. Mantenemos todos los datos iguales, utilizando también el mismo paso de tiempo.

Tendremos, por un lado, el problema de valor inicial resuelto por el Método del Trapecio.

Al ser un método de aproximación implícito, tenemos que despejar primero la variable [math]y_{n+1}[/math] que queremos calcular, que es la que necesitamos para operar el método numéricamente. Aplicando el método:

Esta variable, una vez despejada, quedaría de la siguiente forma:

[math]y_{n+1} = \frac{-Q+sqrt{Q^2-4·W·E}}{2·W}[/math]

%PVI por Trapecio con una h = 0.1
clear all

%condiciones dadas
y0 = 0;
t0 = 0;
tN = 2;
h = 0.1;
k = 1;
a = 3;
b = 1;
%calculamos N
N = (tN - t0)/h;
%definimos la variable t
t = linspace(t0,tN,N+1);
%funcion y
y = zeros(1,N+1);

y(1)= y0;

%solucion obtenida
for i=1:N
    %Definimos tres variables para reducir la ecuación
    Q = 1 + (h*k*(a+b))/2;
    W = -(h*k)/2;
    E = -y(i)- h*k*(1/2)*((a - y(i))*(b - y(i)) + a*b);
    y(i+1)=(-Q + sqrt(Q^2- (4*W*E)))/(2*W);
end

%Concentracion de los reactivos
ya = zeros(1,N+1);
yb = zeros(1,N+1);

ya = 3.- y;
yb = 1.- y; 


%solucion grafica
hold on
plot(t,y,'linewidth',2)
plot(t,ya,'g','linewidth',2)
plot(t,yb,'r','linewidth',2)
legend('Concentracion de C','Concentracion de A','Concentracion de B','Location','best');
hold off

Acompañado de su gráfica:

Gráfica de las dos concentraciones y el producto. Utilizando el método del Trapecio y con una h = 0.1

2.2 Método de Runge-Kutta de Orden 4

Y, por otro lado, el problema de valor inicial resuelto empleando el Método de Runge Kutta, quedando así:

%PVI por Runge Kutta orden 4 con una h = 0.1
clear all

%condiciones dadas
y0 = 0;
t0 = 0;
tN = 2;
h = 0.1;
k = 1;
a = 3;
b = 1;
%calculamos N
N = (tN - t0)/h;
%definimos la variable t
t = linspace(t0,tN,N+1);
%funcion y
y = zeros(1,N+1);

y(1)= y0;

%solucion obtenida
for i=1:N
    K1 = k*(a - y(i))*(b - y(i));
    K2 = k*(a-(y(i)+K1*h/2))*(b - (y(i)+K1*h/2));
    K3 = k*(a-(y(i)+K2*h/2))*(b - (y(i)+K2*h/2));
    K4 = k*(a-(y(i)+(K3*h)))*(b - (y(i)+K3*h));
    y(i+1)= y(i) + h/6*(K1+(2*K2)+(2*K3)+K4);
end

%Concentracion de los reactivos
ya = zeros(1,N+1);
yb = zeros(1,N+1);

ya = 3.- y;
yb = 1.- y; 


%solucion grafica
hold on
plot(t,y,'linewidth',2)
plot(t,ya,'g','linewidth',2)
plot(t,yb,'r','linewidth',2)
legend('Concentracion de C','Concentracion de A','Concentracion de B','Location','best');
hold off

Acompañado también de su gráfica:

Gráfica de las dos concentraciones y el producto. Utilizando el método de Runge Kutta de orden 4

3 Segunda reacción. A partir de C obtenemos otro compuesto D

Suponemos ahora, según el enunciado dado, que está ocurriendo una reacción consecutiva que tiene la siguiente estructura:

                                                    A + B → C → D

La forma en la que esta reacción química funciona es muy similar a la reacción anterior con la que hemos trabajado: es igualmente una reacción bimolecular, luego sigue ocurriendo que una molécula de A y una de B van a producir una de C, pero ahora hay que considerar también el hecho de que este compuesto C se irá descomponiendo en la misma proporción en la que lo estaba haciendo anteriormente, dando lugar, por tanto, al compuesto D.

Para poder calcular la velocidad con la que se forma el compuesto D, tenemos también en cuenta la conservación de la Ley de Acción de masas en esta reacción consecutiva que estamos analizando.

La [math]y(t)[/math] calculada anteriormente, es decir, la concentración del compuesto C en función del tiempo, ahora pasaría a ser la suma de las concentraciones y(t)= yreal(t)+z(t) es decir: y(t)= [C] + [D]

                                  yreal(t)=y(t)- z(t)
 

Por tanto, después de analizar nuestra nueva reacción, definimos el nuevo problema de valor inicial asociado a esta reacción, que es:

[math] \left . \begin{matrix} yreal'(t)=y'(t)- z'(t)\\ z'(t)=k·yreal(t)\\ yreal'(t)=k_1·(a_0-y(t))·(b_0-y(t))-k_2·(y(t)- z(t))\\ a_0=3\\b_0=1\\y(0)=0\\z(0)=0 \end{matrix} \right \} [/math]

Será con este PVI con el que trabajaremos a partir de ahora.

3.1 Resolución numérica

Para la resolución numérica de este PVI, vamos a proponer dos casos diferentes con los que trabajar. Suponemos los mismos datos iniciales que hemos utilizado en la primera reacción, que son, a saber [math]a_0[/math]=3 , [math]b_0[/math]=1 y [math]k_1[/math]=1, para los primeros 10 segundos de duración de la reacción.

1) Utilizaremos, en el primer caso un valor de [math]k_2[/math]=5, lo que supone que [math]k_2[/math] es mayor que [math]k_1[/math], luego la segunda reacción es mucho más rápida que la primera.

2) En el segundo caso, utilizaremos un valor de [math]k_2[/math]=1/5, por lo que [math]k_2[/math] es menor que [math]k_1[/math], luego la segunda reacción es mucho más lenta que la primera.

Emplearemos dos Métodos Numéricos diferentes para resolver ambos casos, siendo estos el Método de Euler y el Método de Runge-Kutta. Esto nos va a permitir comparar los resultados obtenidos en los distintos procedimientos y así establecer una solución lo más precisa posible.

3.1.1 Euler

A la hora de crear el código Matlab con el método de Euler, ponemos una entrada de datos (input) tanto en el valor de la constante [math]k_2[/math], para poder realizar con el mismo código los dos casos descritos anteriormente, como en el tamaño de longitud de paso [math]h[/math], ya que vamos a utilizar varios valores también de la misma para comparar las distintas soluciones obtenidas.

El código numérico que resuelve nuestro PVI con el método de Euler es el siguiente:

%Establecemos los valores iniciales
t0=0;
tN=10;
y0=0;
k1=1;
a=3;
b=1;
k2=input('introduce k2: ');
%Discretizamos
h=input('introduce longitud de paso: ');
N=round((tN-t0)/h);
t=t0:h:tN;
%Creamos el vector y, donde almacenaremos las soluciones
y1=zeros(1,N+1);
yd=zeros(1,N+1);

y1(1)=y0;
yd(1)=y0;

yy1=y0;
yyd=y0;

%Utilizamos el 'for' para crear los vectores solución
for n=1:N
  yy1=yy1+h*(k1*(a-yy1)*(b-yy1)-k2*(yy1-yyd));
  yyd=yyd+h*(k2*(yy1-yyd));
  y1(n+1)=yy1;
  yd(n+1)=yyd;
end

%Definimos las variables de las concentraciones para dibujarlas
ya=3.-y1;
yb=1.-y1;
yc=y1-yd;

hold on

plot(t,yc,'b','linewidth',2)
plot(t,ya,'g','linewidth',2)
plot(t,yb,'r','linewidth',2)
plot(t,yd,'k','linewidth',2)
legend('Concentracion de C','Concentracion de A','Concentracion de B','Concentracion de D','Location','best');
hold off

Una vez hemos descrito el código, nos disponemos a pintar en gráficas las distintas evoluciones de las concentraciones de los elementos a lo largo del tiempo según nuestra reacción química consecutiva.

La primera gráfica que representamos es la correspondiente al primer caso que hemos descrito anteriormente, es decir, a un valor de [math]k_2[/math]=5 y a una longitud de paso [math]h[/math]=0.1, quedando de la siguiente manera:

Gráfica de las dos concentraciones A y B, el producto y la concentración D originada a partir de C. Utilizando el método de Euler.h=0.1

Como es lógico, C y D empiezan a crecer. D siempre crecerá tendiendo a 1, y C seguirá un recorrido decreciente, puesto que tiene que acabar tendiendo a 0.

C empieza a decrecer antes del primer segundo, y su máximo está tan solo a poco más de 0.3 molar. Este hecho es debido a que la velocidad de reacción de D es mayor que la velocidad de reacción de C, que era lo que suponía el utilizar una [math]k_2[/math] superior a la [math]k_1[/math] de la primera reacción, tal y como se ha explicado previamente.

A continuación, vamos a estudiar qué ocurriría si, utilizando la misma [math]k_2[/math] anterior, decidimos utilizar una longitud de paso [math]h[/math]=0.3. Introduciendo nuestros nuevos datos en el input y dibujando la gráfica solución del problema, observamos lo que ocurre con la misma:

Gráfica de las dos concentraciones A y B, el producto y la concentración D originada a partir de C. Utilizando el método de Euler.h=0.3

Tal y como se puede observar en la gráfica, es un tamaño de paso demasiado grande para resolver este problema, ya que sitúa a D por encima de 1 y a B y C con valores negativos. Cabe destacar que, a pesar de que transcurridos los 10 segundos de duración de la reacción se puede observar que los límites de las gráficas de las concentraciones tienden a los mismos valores con los que lo hacía para el tamaño de longitud de paso anterior, no es un tamaño de h que permita ver claramente cuál es la solución del problema.

Tras este pequeño análisis, vamos a resolver ahora el segundo caso inicialmente descrito, es decir, emplearemos ahora una [math]k_2[/math]=1/5 y mantendremos la misma longitud de paso que hemos utilizado en la primera gráfica, [math]h[/math]=0.1.

La gráfica nos quedaría así:

Gráfica de las dos concentraciones A y B, el producto y la concentración D originada a partir de C. Utilizando el método de Euler.h=0.1

Al ser la velocidad de la reacción C 5 veces mas rápida que la de D, observamos que C tiene un máximo bastante mayor que D, en este caso de valor 0.9 ( siendo 1 el límite al que tiende la concentración de C), y también se observa, en comparación con el primer caso, que la gráfica de la concentración de C decrece con menos pendiente hacia 0.

3.1.2 Runge-Kutta de Orden 4

Ahora nos disponemos a desarrollar el problema con el segundo método numérico anteriormente propuesto, el Método de Runge-Kutta. Igual que ocurría para el método de Euler, en este caso también ponemos una entrada de datos (input) tanto en el valor de la constante [math]k_2[/math] como en el tamaño de longitud de paso [math]h[/math], ya que también vamos a utilizar varios valores diferentes de las mismas.

El código numérico que resuelve nuestro PVI con el método de Runge-Kutta es el siguiente:

%C como reactivo y D como producto
t0=0;
tn=10;
h=input('introduce longitud de paso: ');
N=round((tn-t0)/h);
t=t0:h:tn;
y1=zeros(1,N+1);
yd=zeros(1,N+1);
k1=1;
k2=input('introduce k2: ');
a=3;
b=1;
for n=1:N
    
    %runge kutta 
    kc1=k1*(a-y1(n))*(b-y1(n));
    kc2=k1*(a-y1(n)-h/2*kc1)*(b-y1(n)-h/2*kc1);
    kc3=k1*(a-y1(n)-h/2*kc2)*(b-y1(n)-h/2*kc2);
    kc4=k1*(a-y1(n)-h*kc3)*(b-y1(n)-h*kc3);
    y1(n+1)=y1(n)+(h/6)*(kc1+2*kc2+2*kc3+kc4);
    kd1=k2*(y1(n)-yd(n));
    kd2=k2*(y1(n)-yd(n)-h/2*kd1);
    kd3=k2*(y1(n)-yd(n)-h/2*kd2);
    kd4=k2*(y1(n)-yd(n)-h*kd3);
    yd(n+1)=yd(n)+(h/6)*(kd1+2*kd2+2*kd3+kd4);
 
end
%concentraciones
ya=3.-y1;
yb=1.-y1;
yc=y1-yd;

hold on

plot(t,yc,'b','linewidth',2)
plot(t,ya,'g','linewidth',2)
plot(t,yb,'r','linewidth',2)
plot(t,yd,'k','linewidth',2)
legend('Concentracion de C','Concentracion de A','Concentracion de B','Concentracion de D','Location','best');
hold off

El Método de Runge-Kutta es una forma numérica de resolver problemas de valor inicial que trabaja de la misma manera que lo hace el Método de Euler, dando soluciones aproximadas al problema de valor inicial propuesto. Por tanto, las gráficas que nos van a salir con este método van a ser muy parecidas a las obtenidas con el Método de Euler, pero sabiendo de antemano que el Método de Runge-Kutta consigue hacer una mejor aproximación de las soluciones de la que puede hacer Euler.

Una vez conocido esto, nos disponemos a pintar las mismas gráficas representadas con el método anterior pero utilizando ahora Runge-Kutta:

Primer caso: valor de [math]k_2[/math]=5 y longitud de paso [math]h[/math]=0.1, quedando de la siguiente manera:

Gráfica de las dos concentraciones A y B, el producto y la concentración D originada a partir de C. Utilizando el método de Runge Kutta h=0.1 k2=5

Las gráficas evolucionan de la misma manera que lo hacían antes, dando sin embargo una solución más aproximada que la primera obtenida.

Estudiamos también qué ocurriría si decidimos utilizar una longitud de paso [math]h[/math]=0.3, con las mismas condiciones con las que lo hemos hecho con el anterior método pero usando ahora Runge-Kutta. Representamos la gráfica obtenida:

Gráfica de las dos concentraciones A y B, el producto y la concentración D originada a partir de C. Utilizando el método de Runge Kutta h=0.3 k2=5

El error para este método es notable, ya que C llega a valer 0.6 cuando no debería sobrepasar 0.5. Por lo tanto, este tamaño de longitud de paso [math]h[/math]=0.3 no nos sirve tampoco para el Método de Runge-Kutta. Aunque no nos vale para ninguno de los dos métodos, ya que obtenemos errores en ambos dos, si que es preciso comentar que error obtenido con el Método de Runge-Kutta para [math]h[/math]=0.3 es un error menor al obtenido con el Método de Euler para [math]h[/math]=0.3.

Ahora el segundo caso: una [math]k_2[/math]=1/5 y mantendremos la misma longitud de paso que hemos utilizado en la primera gráfica, [math]h[/math]=0.1. La gráfica sería:

Gráfica de las dos concentraciones A y B, el producto y la concentración D originada a partir de C. Utilizando el método de Runge Kutta h=0.1 k2=1/5

Las gráficas evolucionan igual, siendo los valores mucho mejor aproximados.

3.2 Discusión y análisis