Reacciones con autocatálisis. Grupo D12

Trabajo realizado por estudiantes
Título	Reacciones con autocatálisis. Grupo D12
Asignatura	Ecuaciones Diferenciales
Curso	Curso 2014-15
Autores	Javier Ruiz de Galarreta López, Argimiro Martínez López, Eduardo Moyano, Alberto Rodríguez Ruiz.
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 Introducción

1.1 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.

1.2 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.

1.3 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 Reacción I

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

Deducir a partir de la ley de accion de masas y del principio de conservacion de la masa que las concentraciones de A y B deben satisfacer las ecuaciones

[math]x'(t)+y'(t)=0[/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.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 Resolución por 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 Resolución por Runge-Kutta

El método de resolución numérica de Runge-Kutta, involucra un gran número de parámetros.

Definiendo un problema de valor inicial como:

[math] y' = f(x, y), \quad y(x_0) = y_0[/math]

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

[math]y_{i+1} = y_i + {1 \over 6}h\left ( k_1 + 2k_2 + 2k_3 + k_4 \right ) [/math]

Donde

[math] \begin{cases} k_1 & = f \left( x_i, y_i \right) \\ k_2 & = f \left( x_i + {1 \over 2}h, y_i + {1 \over 2}k_1 h \right) \\ k_3 & = f \left( x_i + {1 \over 2}h, y_i + {1 \over 2}k_2 h \right) \\ k_4 & = f \left( x_i + h, y_i + k_3h \right) \\ \end{cases} [/math]

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. [math] :\ltmath\gt(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] \quad \Rightarrow \quad (PVI)\begin{cases}z'(t)= \begin{bmatrix} x'(t) \\ y'(t) \\ \end{bmatrix} = \begin{bmatrix} -k_1x(t)y(t) \\ k_1x(t)y(t) \\ \end{bmatrix} \\ \\ 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]

Función F:

[math] F(x,y) = \begin{bmatrix} x'(t) \\ y'(t) \\ \end{bmatrix} = \begin{bmatrix} -k_1x(t)y(t) \\ k_1x(t)y(t) \\ \end{bmatrix} \quad \Rightarrow \quad k_1=1 \quad \Rightarrow \quad F(x,y) = \begin{bmatrix} -x(t)y(t) \\ x(t)y(t) \\ \end{bmatrix} [/math]

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('[-x*y;x*y]','x','y');

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

%Matriz 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(x(i),y(i));
    k2=F(x(i)+k1(1)/2*h,y(i)+k1(2)/2*h);
    k3=F(x(i)+k2(1)/2*h,y(i)+k2(2)/2*h);
    k4=F(x(i)+k3(1)*h,y(i)+k3(2)*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