Revisión del 01:30 3 mar 2015

1 Introducción: objetivos y metodología

El objetivo de este trabajo es estudiar las concentraciones de los reactivos y productos de una reacción química a lo largo de la evolución en el tiempo de dicha reacción. Para ello modelizaremos el proceso empleando ecuaciones diferenciales y teniendo en cuenta la Ley de Acción de Masas y el Principio de Conservación de la materia.

Al ser las ecuaciones y sistemas que se nos presentan de difícil resolución analítica, utilizaremos el programa informático Matlab para resolverlos mediante métodos numéricos como el de Euler, el del Trapecio, el de Runge-Kutta y el de Heun.

Plasmaremos las soluciones en gráficos que nos permitirán dar una interpretación físico química a los resultados.

2 Reacción I

En primer lugar consideramos una reacción química irreversible en una solución bien mezclada, que se desarrolla para un volumen y una temperatura constantes, y en la que uno de los componentes (B, en este caso), hace de efecto catalítico.

[math] A + B →_{k1} 2B [/math]

Para desarrollar nuestro análisis nos apoyaremos en dos leyes físico químicas:

Ley de acción de masas: Enunciada por Guldberg y Waage en 1864, establece que dada una reacción química reversible en equilibrio, a temperatura constante, la relación de concentraciones de los reactivos y productos tiene un valor constante.

Principio de conservación de la masa: Establece que en una reacción química ordinaria la masa permanece constante, siendo la masa consumida de los reactivos igual a la masa obtenida de los productos. Fue elaborada independientemente por Míjail Lomonósov en 1745 y Antoine Lavoisier en 1785.

2.1 Interpretación del problema y deducción de las ecuaciones diferenciales

Para hallar las ecuaciones diferenciales que nos describan cómo se desarrolla la reacción en el tiempo nos basaremos en las dos leyes enunciadas anteriormente. En primer lugar, definimos nuestra variable independiente como el tiempo [math]t[/math] en segundos. Las dos variables dependientes serán las concentraciones de los compuestos A y B, que denominaremos [math]x(t)[/math] e [math]y(t)[/math], ambas expresadas en mol/l.

El principio de conservación de la masa expone que la suma de las concentraciones de los reactivos y los productos es constante. Esto equivale, derivando la expresión [math]x(t)+y(t)=cte[/math], a que desaparece tanto reactivo como producto aparece, o lo que es lo mismo, que la velocidad de aparición del producto es igual a la velocidad de desaparición (velocidad de aparición negativa) del reactivo:

[math]x'(t)=-y'(t)[/math]

[math]y'(t)=k_1x(t)y(t)[/math]

.....................

2.2 Obtención del problema de valor inicial

Una vez demostrado el origen de las ecuaciones (1) y (2) en el apartado anterior, obtendremos el problema de valor inicial (PVI) mediante las manipulaciones detalladas a continuación:

[math]x'(t)+y'(t)=0\quad\quad(1)[/math]

[math]y'(t)=k_1x(t)y(t)\quad\quad(2)[/math]

Integrando (1) obtenemos (3):

[math]\int_\mathbb{D} (\frac{dx(t)}{dt} + \frac{dy(t)}{dt})\,dt =\int_\mathbb{D} 0\,dt \quad \Rightarrow \quad x(t)+y(t)=k_2\quad \Rightarrow \quad x(t)=k_2-y(t)\quad\quad(3)[/math]

Sustituyendo (3) en (2), obtenemos (4):

[math]y'(t)=k_1x(t)y(t)\quad \Rightarrow \quad y'(t)=k_1(k_2-y(t))y(t)\quad \Rightarrow \quad y'(t)=k_1k_2y(t)-k_1y^2(t)\quad\quad(4)[/math]

Conocemos las concentraciones iniciales de A y de B, por lo que podemos obtener el valor de la constante k₂:

[math][A]:\quad x(0)=1 mol/l;\quad\quad [B]:\quad y(0)=0,01 mol/l[/math]

[math]x(t)+y(t)=k_2\quad \Rightarrow \quad x(0)+y(0)=k_2\quad \Rightarrow \quad 1+0,01=k_2\quad \Rightarrow \quad k_2=1,01[/math]

Por último, conocemos también el valor de k₁, por lo que podremos plantear a continuación el correspondiente PVI:

[math]k_1=1[/math]

[math](PVI)\begin{cases}y'(t)=1,01y(t)-y^2(t) \\ \\ y(0) = 0,01 \end{cases}\quad (t,y)\in[0,10]\times\mathbb{R}[/math]

2.3 Unicidad de la solución

Intersección B con dominio

Aplicamos a continuación el Teorema de Existencia y Unicidad:

Sea [math]f(t,y)=y'(t) \Rightarrow f(t,y)=1.01y(t)-y^2(t)\quad[/math], y [math]\quad\frac{\partial }{\partial y} f(t,y) = 1.01-2y(t)[/math].

Y sea [math]B[(t_0,y_0),r],\quad r\gt0,\quad en (t_0,y_0)=(0,0.01)[/math].

[math]f(t,y)[/math] es continua en [math] \mathbb{B}\cap \mathbb{D}\quad \Rightarrow \quad \exists y(t)\quad /\quad y'(t)=1,01y(t)-y^2(t)[/math]

[math]\frac{\partial }{\partial y} f(t,y)[/math] es continua en [math] \mathbb{B}\cap \mathbb{D}\quad \Rightarrow \quad \exists ! y(t)\quad /\quad y'(t)=1,01y(t)-y^2(t)[/math]

Con lo que queda demostrado que el (PVI) tiene solución, y además esta es única.

2.4 Resolución numérica del problema de valor inicial

2.4.1 Método de Euler

2.4.2 Método del trapecio

2.4.3 Método de Runge-Kutta

Representación gráfica y(t). Método Runge-Kutta

%Método de Runge-Kutta
clc

%Función f(t)=y'(t)
fyt=inline('1.01*y-y^2','y');

%Datos del problema
y0=0.01;
t0=0;
tN=10;
h=0.1;

%Variable independiente (t)
t=t0:h:tN;

%Variable dependiente (y)
y=zeros(1,length(t));
y(1)=y0;

for i=1:length(t)-1
    k1=fyt(y(i));
    k2=fyt(y(i)+k1/2*h);
    k3=fyt(y(i)+k2/2*h);
    k4=fyt(y(i)+k3*h);
    y(i+1)=y(i)+h/6*(k1+2*k2+2*k3+k4);
end

%Representación gráfica
hold on
plot(t,y)
legend('Sustancia B: y(t)')
xlabel('Tiempo (s)')
ylabel('Concentración (mol/l)')
grid on

2.5 Resolución en forma de sistema de ecuaciones

A continuación vamos a plantear el problema de valor inicial con un sistema de dos ecuaciones equivalente al problema resuelto anteriormente.

Estas ecuaciones serán:

[math](PVI)\begin{cases}y'(t)=k_1 x(t)y(t) \\ \\ x'(t)=-k_1 x(t)y(t) \\ \\ x(0)=1 \\ \\ y(0) = 0.01\end{cases}\quad t\in[0,10][/math]

Resolveremos el problema por el método de Euler y por el método de Runge-Kutta, plasmando las soluciones en gráficos.

2.5.1 Método de Euler

Representación gráfica de las concentraciones de A y B. Resolución del sistema por Euler

%Resolución del sistema por el método de Euler

clear all

%Definimos el vector de la variable independiente, el número de
%subintervalos y la longitud de paso

t0 = 0;
tN = 10;
h = 0.1;
t = t0:h:tN
N = (tN-t0)/h;

%Definimos los vectores de las variables dependientes

x = zeros(1, N+1);
y = zeros(1, N+1);

x(1) = 1;
y(1) = 0.01;

for i = 1:N
    
    x(i+1)=x(i)+h*(-x(i).*y(i));
    y(i+1)=y(i)+h*(y(i).*x(i));
    
end

%Representamos gráficamente las soluciones

hold on

plot(t,x,'b');
plot(t,y,'r');
grid on

legend('Sustancia A: x(t)','Sustancia B: y(t)')
xlabel('Tiempo (s)')
ylabel('Concentración (mol/l)')


hold off

2.5.2 Método de Runge-Kutta

El método de resolución numérica de Runge-Kutta, involucra un gran número de parámetros. Para simplificar la implementación de este método en Matlab, el sistema se tratará como una matriz [math]z[/math], y para la obtención de los parámetros [math]k_i[/math] utilizaremos la función inline, con la que declararemos previamente las expresiones necesarias para el cálculo de los parámetros [math]k_i[/math], ahorrándonos tener que escribir innumerables veces dichas expresiones.

Primero, expresamos el sistema en forma matricial:

[math](PVI)\begin{cases}x'(t)=-k_1 x(t)y(t) \\ \\ y'(t)=k_1 x(t)y(t) \\ \\ x(0)=1 \\ \\ y(0) = 0.01\end{cases}\quad t\in[0,10] \quad \Rightarrow \quad (PVI)\begin{cases}\vec{z}'(t)= \begin{bmatrix} x'(t) \\ y'(t) \\ \end{bmatrix} = \begin{bmatrix} -k_1x(t)y(t) \\ k_1x(t)y(t) \\ \end{bmatrix} \\ \\ \vec{z}(0)= \begin{bmatrix} x(0) \\ y(0) \\ \end{bmatrix} = \begin{bmatrix} 1 \\ 0.01 \\ \end{bmatrix} \end{cases}\quad t\in[0,10][/math]

A continuación, detallamos la función [math]\vec{F}[/math] de la que dependerán los valores de los parámetros [math]k_i[/math]:

[math]\vec{F}(x,y)=\begin{bmatrix}f(x,y) \\ g(x,y) \\ \end{bmatrix}=\begin{bmatrix}x'(t) \\ y'(t) \\ \end{bmatrix}=\vec{z}'(t)=\begin{bmatrix}-k_1x(t)y(t) \\ k_1x(t)y(t) \\ \end{bmatrix}\quad \Rightarrow \quad k_1=1\quad \Rightarrow \quad \vec{F}(x,y)=\begin{bmatrix}-x(t)y(t) \\ x(t)y(t) \\ \end{bmatrix}\quad \Rightarrow \quad \begin{cases} [\vec{k_1}] & = F(\vec{z}_i) \\ [\vec{k_2}] & = F(\vec{z}_i + {1 \over 2}k_1h) \\ [\vec{k_3}] & = F(\vec{z}_i + \frac{1}{2}k_2h) \\ [\vec{k_4}] & = F(\vec{z}_i + k_3h) \\ \end{cases} [/math]

Entonces el método RK4 para este problema está dado por la siguiente ecuación:

[math]\vec{z}_{i+1} = \vec{z}_{i} + {1 \over 6}h(\vec{k_1}+2\vec{k_2}+2\vec{k_3}+\vec{k_4})[/math]

Mediante el siguiente programa en Matlab calculamos numéricamente la solución y la representamos en un gráfico:

Representación gráfica de las concentraciones de A y B. Resolución del sistema por Runge-Kutta

%Método de Runge-Kutta
clc

%Función F=[x'(t); y'(t)]
F=inline('[-z(1)*z(2);z(1)*z(2)]','z');

%Datos del problema 
x0=1;
y0=0.01;
t0=0;
tN=10;
h=0.1;

%Vector de variable independiente (t)
t=t0:h:tN;

%Matriz de variables dependientes z=[x;y]
z=zeros(2,length(t));
z(:,1)=[x0;y0];

for i = 1:length(t)-1
    k1=F(z(:,i));
    k2=F(z(:,i)+k1.*h/2);
    k3=F(z(:,i)+k2.*h/2);
    k4=F(z(:,i)+k3.*h);
    z(:,i+1)=z(:,i)+h/6.*(k1+2.*k2+2.*k3+k4);
end
 
%Representación gráfica
hold on
plot(t,z(1,:))
plot(t,z(2,:),'r')
legend('Sustancia A: x(t)','Sustancia B: y(t)')
xlabel('Tiempo (s)')
ylabel('Concentración (mol/l)')
grid on
%hold off

2.6 Interpretación de las gráficas

Podemos observar por los resultados obtenidos que la velocidad de aparición de B es igual a la velocidad de desaparición de A (como dicta el principio de conservación de la masa), que A acaba desapareciendo prácticamente por completo y que la concentración de B termina por estabilizarse alrededor de 1 mol/l.

También cabe destacar que cuando el tiempo transcurrido alcanza 4,5 segundos las concentraciones de ambos compuestos son las mismas.

3 Reacción II

3.1 Interpretación del problema y deducción de las ecuaciones diferenciales

En este apartado vamos a estudiar la reacción propuesta por Alfred J. Lotka en 1920, que es la siguiente:

3.2 Resolución numérica del sistema de ecuaciones

3.2.1 Método de Euler

3.2.2 Método de Heun

El método de resolución numérica de Runge-Kutta, involucra un gran número de parámetros. Para simplificar la implementación de este método en Matlab, el sistema se tratará como una matriz [math]z[/math], y para la obtención de los parámetros [math]k_i[/math] utilizaremos la función inline, con la que declararemos previamente las expresiones necesarias para el cálculo de los parámetros [math]k_i[/math], ahorrándonos tener que escribir innumerables veces dichas expresiones.

Primero, expresamos el sistema en forma matricial:

[math](PVI)\begin{cases} A'(t)=-k_1 A(t)x(t) \\ x'(t)=k_1 A(t)x(t)-k_2 x(t)y(t) \\ y'(t)=k_2 x(t)y(t)-k_3 y(t) \\ B'(t)=k_3 x(t)y(t) \\ \\ A(0)=5 \\ x(0)=5·10^{-4} \\ y(0)=10^{-5} \\ B(0) = 0 \end{cases}\quad t\in[0,200] \quad \Rightarrow \quad \begin{cases} k_1=0.1 \\ k_2=0.1 \\ k_3=0.05 \\ \end{cases} \quad \Rightarrow \quad (PVI)\begin{cases}\vec{z}'(t)= \begin{bmatrix} A'(t) \\ x'(t) \\ y'(t) \\ B'(t) \\ \end{bmatrix} = \begin{bmatrix} -0.1 A(t)x(t) \\ 0.1 A(t)x(t)-0.1 x(t)y(t) \\ 0.1 x(t)y(t)-0.05 y(t) \\ 0.05 y(t)x(t)\\ \end{bmatrix} \\ \\ \vec{z}(0)= \begin{bmatrix} A(0) \\ x(0) \\ y(0) \\ B(0) \\ \end{bmatrix} = \begin{bmatrix} 5 \\ 5·10^{-4} \\ 10^{-5} \\ 0 \end{bmatrix} \end{cases}\quad t\in[0,10][/math]

A continuación, detallamos la función [math]\vec{F}[/math] de la que dependerán los valores de los parámetros [math]k_i[/math]:

[math]\vec{F}(x,y)=\begin{bmatrix}f(\vec{z}) \\ g(\vec{z}) \\ h(\vec{z}) \\ i(\vec{z}) \\ \end{bmatrix}=\begin{bmatrix}A'(t) \\ x'(t) \\ y'(t) \\ B'(t) \\ \end{bmatrix}=\vec{z}'(t)=\begin{bmatrix}-k_1x(t)y(t) \\ k_1x(t)y(t) \\ \end{bmatrix}\quad \Rightarrow \quad k_1=1\quad \Rightarrow \quad \vec{F}(x,y)=\begin{bmatrix}-x(t)y(t) \\ x(t)y(t) \\ \end{bmatrix}\quad \Rightarrow \quad \begin{cases} [\vec{k_1}] & = F(\vec{z}_i) \\ [\vec{k_2}] & = F(\vec{z}_i + {1 \over 2}k_1h) \\ \end{cases} [/math]

Entonces el método de Heun para este problema está dado por la siguiente ecuación:

[math]\vec{z}_{i+1} = \vec{z}_{i} + {1 \over 2}h(\vec{k_1}+\vec{k_2})[/math]

Mediante el siguiente programa en Matlab calculamos numéricamente la solución y la representamos en un gráfico:

Representación gráfica de las concentraciones de A,X,Y y B. Resolución del sistema por Heun

%Método de Heun
clc

%Función F=[A'(t);x'(t); y'(t);B'(t)]
F=inline('0.1.*[-z(1)*z(2);z(1)*z(2)-z(2)*z(3);z(2)*z(3)-z(3)/2;z(3)/2]','z');

%Datos del problema
A0=5;
x0=5*10^-4;
y0=10^(-5);
B0=0;
t0=0;
tN=200;
h=0.1;

%Vector de variable independiente (t)
t=t0:h:tN;

%Matriz de variables dependientes z=[A;x;y;B]
z=zeros(4,length(t));
z(:,1)=[A0;x0;y0;B0];

for i = 1:length(t)-1
    k1=F(z(:,i));
    k2=F(z(:,i)+k1.*h/2);
    z(:,i+1)=z(:,i)+h/2.*(k1+k2);
end
 
%Representación gráfica
hold on
plot(t,z(1,:),'b')
plot(t,z(2,:),'r')
plot(t,z(3,:),'g')
plot(t,z(4,:),'m')
legend('Sustancia A: A(t)','Sustancia X: x(t)','Sustancia Y: y(t)','Sustancia B: B(t)','Location','Best')
xlabel('Tiempo (s)')
ylabel('Concentración (mol/l)')
grid on
hold off

3.3 Comentario de las gráficas