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 García, Alberto Rodríguez Ruiz.
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

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 de 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, entre otras cosas, que la velocidad de reacción es proporcional al producto de las concentraciones de los reactivos.

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]

La otra ecuación que necesitamos la obtenemos directamente de aplicar la Ley de acción de masas, que establece que la velocidad de la reacción (la velocidad a la que aparece el producto o desaparece el reactivo) es igual al producto de las concentraciones por una constante.

[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 operaciones 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].

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

Y siendo [math](t,y) \in \mathbb{D} = [0,10]\times \mathbb{R}^+ [/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

En este apartado resolveremos el problema de valor inicial planteado anteriormente con tres métodos de integración numérica: el método de Euler, el método del Trapecio y el método de Runge-Kutta de orden 4.

2.4.1 Método de Euler

El método de Euler, llamado así en honor al matemático Leonhard Euler, es un procedimiento de integración numérica. Su expresión matemática es: [math] PVI = \begin{cases} y_0 \\ y_{n+1} = y_n + h*f(t_n,y_n) \\ \end{cases}\quad [/math]

A continuación emplearemos esta expresión en un programa de Matlab para obtener una aproximación de las soluciones. El código de dicho programa es el siguiente.

Representación gráfica y(t). Método de Euler

%Método de Euler
clc

%Datos iniciales
t0=0;
tN=4;
y0=0.01;
h=0.1;
t0=0;
tN=10;

%Cálculo de subintervalos
N=round((tN-t0)/h);
t=t0:h:tN;

%Creamos el vector donde vamos a guardar las soluciones
y=zeros(1,N+1);
y(1)=y0;

%Euler
for i=1:N;
    y(i+1)=y(i)+h*(1.01*y(i)-(y(i)).^2);
end
%Tabla de Resultados
[t',y']
x=1.01-y;
hold on 
plot(t,y,'r')
plot(t,x,'b')
hold off
legend('Sustancia B: y(t)','Sustancia A: x(t)','Location','best');
xlabel('Tiempo (s)')
ylabel('Concentración (mol/l)')
grid on

2.4.2 Método del trapecio

Vamos a aplicar un método de integración numérica conocido como el Método del Trapecio por aproximar el valor de la integral al del área de un trapecio. Su expresión matemática es:

[math] PVI = \begin{cases} y_0 \\ y_{n+1} = y_n + \frac{h}{2}*(f(t_n,y_n)+f(t_{n+1},y_{n+1}))\\ \end{cases}\quad [/math]

Introducimos esta expresión en el código y obtenemos la gráfica de la solución.

Representación gráfica y(t). Método del Trapecio

%Método del Trapecio
clc

%Datos iniciales
t0=0;
tN=4;
y0=0.01;
h=0.1;
t0=0;
tN=10;

%Cálculo de subintervalos
N=round((tN-t0)/h);
t=t0:h:tN;

%Creamos el vector donde vamos a guardar las soluciones
y=zeros(1,N+1);
y(1)=y0;

for i=1:N
y(i+1)=(((1-(h/2)*1.01)-sqrt(((1-(h/2)*1.01)^2)+4*(h/2)*((1+(h/2)*1.01)*y(i)-(h/2)*(y(i))^2)))/-h);
end
%Sacamos la tabla de resultados
[t',y']
x=1.01-y;
%Gráfico
hold on
plot(t,y,'r')
plot(t,x,'b')
legend('Sustancia B: y(t)','Sustancia A: x(t)','Location','best');
xlabel('Tiempo (s)')
ylabel('Concentración (mol/l)')
grid on

2.4.3 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, utilizaremos la función inline, con la que declararemos previamente la función algebraica necesaria para el cálculo de los parámetros [math]k_i[/math], ahorrándonos tener que escribir innumerables veces la correspondiente expresión matemática.

A continuación, detallamos la función [math]f(y)[/math] de la que dependerán los valores de los parámetros [math]k_i[/math]:

[math]f(y)=y'(t)=1.01y(t)-y(t)^2 \quad → \quad \begin{cases} k_1 & = f(y_i) \\ k_2 & = f(y_i + {1 \over 2}k_1h) \\ k_3 & = f(y_i + \frac{1}{2}k_2h) \\ k_4 & = f(y_i + k_3h) \\ \end{cases}\quad [h=0.1] [/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(k_1+2k_2+2k_3+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 y(t). Método Runge-Kutta

%Método de Runge-Kutta
clc

%Función f(t)=y'(t)
fy=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=fy(y(i));
    k2=fy(y(i)+k1/2*h);
    k3=fy(y(i)+k2/2*h);
    k4=fy(y(i)+k3*h);
    y(i+1)=y(i)+h/6*(k1+2*k2+2*k3+k4);
end

%Variable dependiente (x)
x=1.01-y;

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

2.4.4 Comparación de métodos

Podemos observar fácilmente que las gráficas obtenidas con los diferentes métodos son iguales a primera vista. Sin embargo, superponiéndolas queda patente que no son exactamente iguales debido a que los métodos no poseen la misma precisión y , al operar todos con la misma longitud de paso, unos aproximan la solución exacta mejor que otros. Así se observa en estas imágenes:

Comparación de métodos

Comparación de métodos

Como podemos comprobar en las imágenes, el método de Euler es menos preciso que el del trapecio o el de Runge-Kutta.

2.5 Resolución en forma de sistema de ecuaciones

A continuación vamos a plantear el problema de valor inicial con un sistema no lineal 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

En este caso aplicaremos en el siguiente código la expresión para EDO a cada una de las ecuaciones que conforma el sistema, obteniendo la correspondiente gráfica.

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]. Además, utilizaremos la función inline, con la que declararemos previamente la función algebraica necesaria para el cálculo de los parámetros [math]k_i[/math], ahorrándonos tener que escribir innumerables veces la correspondiente expresión matemática.

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 \→ \quad k_1=1\quad \Rightarrow \quad (PVI)\begin{cases}\vec{z}'(t)= \begin{bmatrix} x'(t) \\ y'(t) \\ \end{bmatrix} = \begin{bmatrix} -x(t)y(t) \\ x(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)=\vec{F}(\vec{z})= \begin{bmatrix} f(\vec{z}) \\ g(\vec{z}) \\ \end{bmatrix}= \begin{bmatrix} x'(t) \\ y'(t) \\ \end{bmatrix}= \vec{z}'(t)= \begin{bmatrix} -x(t)y(t) \\ x(t)y(t) \\ \end{bmatrix} \quad → \quad \begin{cases} [\vec{k_1}] & = F(\vec{z}_i) \\ [\vec{k_2}] & = F(\vec{z}_i + {1 \over 2}\vec{k_1}h) \\ [\vec{k_3}] & = F(\vec{z}_i + \frac{1}{2}\vec{k_2}h) \\ [\vec{k_4}] & = F(\vec{z}_i + \vec{k_3}h) \\ \end{cases}\quad [h=0.1] [/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

Como cabía esperar, obtenemos la misma solución resolviendo el problema como sistema o como EDO. 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. Ambos compuestos alcanzan los extremos de sus concentraciones en el mismo instante de tiempo

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

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

[math] A + X →_{k1} 2X;[/math]  :[math] X + Y →_{k2} 2Y;[/math]  :[math]Y →_{k3} 2B [/math]

Como se puede observar, se trata de una reacción compleja integrada por 3 reacciones simples (las dos primeras autocatalíticas) en las que intervienen un total de 4 compuestos, A, B, X e Y. Siguiendo el desarrollo de la reacción, es lógico suponer que A se agotará al comienzo de la reacción, que los máximos de concentración de X e Y se darán en los estados intermedios y que al final el compuesto B será el único restante, siendo el producto final de la reacción.

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

Para averiguar cómo evolucionan las concentraciones de los cuatro compuestos durante el desarrollo de la reacción necesitaremos obtener cuatro ecuaciones diferenciales. Para ello nos volveremos a basar en los principios físicos que ya empleamos en la anterior reacción. Definimos A(t), x(t), y(t), y B(t) (A, x, y, B; para ahorrar notación) como las concentraciones de A, X, Y, y B, respectivamente.

En primer lugar aplicamos el Principio de Conservación de la masa, que establece que la suma de las concentraciones de los compuestos es constante (no desaparece masa en el proceso), lo que derivando nos indica que:

[math] A' + x' + y' + B' = 0 \quad\quad(I)[/math]

A continuación aplicaremos la Ley de Acción de Masas (velocidad de reacción proporcional al producto de la concentración de los reactivos).

Aplicando esta ley a la última reacción simple obtenemos que:

[math]B'=k_3y \quad\quad(II)[/math]

En el caso de las concentraciones de X e Y las expresiones resultantes son algo más complicadas ya que estos dos compuestos actúan como reactivos y productos en más de una reacción simple, por lo que la variación de sus concentraciones será la suma de las variaciones de las reacciones en las que aparezcan:

[math] x' = x'_{1} + x'_{2} → x' = k_{1}Ax - k_{2}xy \quad\quad(III)[/math]

[math] y' = y'_{2} + y'_{3} → y' = k_{2}xy − k_3{y} \quad\quad(IV)[/math]

Por último sustituimos (II), (III) y (IV) en (I) obteniendo:

[math] A' + k_{1}Ax - k{2}xy + k_{2}xy − k_3{y}+ k_3y= 0 : [/math]

De donde obtenemos que:

[math] A' = -k_{1}Ax [/math]

Ahora que ya tenemos las variaciones de concentración de los compuestos en función de las concentraciones podemos aplicar métodos numéricos que nos lleven a una solución aproximada.

3.2 Resolución numérica del sistema de ecuaciones

Una vez obtenidas las ecuaciones y conociendo las condiciones iniciales y los valores de las constantes podemos plantear el problema de valor inicial para luego resolverlo.

[math](PVI) \begin{cases} A' = -k_{1}Ax\\ x' = k_{1}Ax - k{2}xy\\ y' = k_{2}xy − k_3{y}\\ B'=k_3y \\ \\ A(0)=5 \\ x(0)=5*10^{-4} \\ y(0)=10^{-5}\\ B(0)=0 \\ \end{cases} \quad t\in[0,200][/math]

Siendo [math] k_{1}=K_{2}=2k_{3}=0.1 [/math] y las unidades de la concentración mol/l.

3.2.1 Método de Euler

Aplicamos el método de Euler de forma análoga a como lo hicimos en apartado 2.5.1

Representación gráfica de las concentraciones de A, B, X, e Y a lo largo del tiempo. Resoluciónd delsistema por el método de Euler con longitud de paso h=0.01

Representación gráfica de las concentraciones de A, B, X, e Y a lo largo del tiempo. Resoluciónd delsistema por el método de Euler con longitud de paso h=0.001

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

clear all

%Definimos la variable independiente y establecemos los datos del problema
%de valor inicial. Para la longitud de paso h introducimos en primer lugar
%0.01 y luego 0.001

t0=0;
tN=200;
h= input('Introduzca una longitud de paso ');
t = t0:h:tN;
N=(tN-t0)/h;
k1 = 0.1;
k2 = 0.1;
k3 = 0.05;
A0 = 5;
B0 = 0;
X0 = 5*10^(-4);
Y0 = 10^(-5);


%Definimos las variables dependientes.

A = zeros(1, N-1);
B = A;
X = A;
Y = A;

A(1) = A0;
B(1) = B0;
X(1) = X0;
Y(1) = Y0;

for i = 1:N
    
    X(i+1) = X(i) + h*(k1*A(i)*X(i)-k2*Y(i)*X(i));
    Y(i+1)=Y(i)+h.*(k2.*X(i).*Y(i)-k3.*Y(i));
    B(i+1) = B(i) + h*(k3*Y(i));
    A(i+1) = A(i) + h*(-k1*A(i)*X(i));
   
end

%Representamos gráficamente la solución

hold on

plot (t,A,'b')
plot(t, X, 'r')
plot(t,Y, 'g')
plot(t, B, '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

A continuación tendremos que comprobar si el método es estable.

Un método es estable cuando el error local es pequeño y además este error no se amplifica a lo largo del tiempo. Para calificar la estabilidad del método de Euler resolveremos el sistema anterior con este método, con pasos "h=0.01" y "h=0.001", y lo compararemos con la resolución del sistema con el método de Heun, que se asemeja más a la solución exacta.
En la comparación representaremos en azul el método de Heun y en rojo el método de Euler.
Para "h=0.01" tenemos lo siguiente:

Comparación entre Método de Heun (Azul) con el Método de Euler con paso h=0.01 (Rojo)

Comparación entre Método de Heun (Azul) con el Método de Euler con paso h=0.01 (Rojo)

Para "h=0.001" tenemos lo siguiente:

Comparación entre Método de Heun (Azul) con el Método de Euler con paso h=0.001 (Rojo)

Comparación entre Método de Heun (Azul) con el Método de Euler con paso h=0.001 (Rojo)

Como podemos ver en las imágenes de los dos procesos las gráficas se superponen y no se aprecian diferencias hasta que no ampliamos la imagen. Y esto sucede en todo momento con lo que podemos garantizar que en este caso el método de Euler es estable para los pasos que hemos usado.

3.2.2 Método de Heun

Si bien el método de resolución numérica de Heun no involucra un número de parámetros tan notable como el método de Runge-Kutta, en este caso tenemos un sistema de 4 ecuaciones diferenciales, por lo que, en aras de simplificar, utilizaremos el mismo criterio que en el apartado 2.5.2.: 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].

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,200][/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}(A,x,y,B)=\vec{F}(\vec{z})= \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} -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} \quad → \quad \begin{cases} [\vec{k_1}] & = F(\vec{z}_i) \\ [\vec{k_2}] & = F(\vec{z}_i + {1 \over 2}\vec{k_1}h) \\ \end{cases}\quad [h=0.01] [/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 Interpretación de las gráficas

El proceso comienza con el predominio absoluto del compuesto A en la concentración total, siendo la presencia de los otros tres compuestos despreciable. Durante los primeros segundos la situación apenas cambia, pero tras los diez primeros segundos la concentración de A cae bruscamente mientras que la concentración de X aumenta a la misma velocidad; antes de llegar a los veinte segundos transcurridos las concentraciones de ambos compuestos se igualan y otros quince segundos más tarde el compuesto A prácticamente ha desaparecido y el X ha alcanzado su concentración máxima (5 mol/litro).

Es entonces cuando concluye la primera reacción y comienzan las otras dos, al no haber más compuesto A que cree compuesto X. La velocidad de aparición de Y es inicialmente mayor que la de B, superando al poco tiempo en concentración a la sustancia X. Sin embargo, llega un momento en el que la concentración de X es tan baja que se ralentiza la creación de Y y pronto empieza a disminuir su concentración tras alcanzar un máximo en torno a los 56 segundos, perdiendo peso la segunda reacción en beneficio de la tercera.

Tras los 60 segundos de reacción el compuesto X pasa a tener una concentración despreciable, por lo que la única reacción simple que se sigue desarrollando es la tercera, destruyéndose reactivo Y y creándose producto B. La concentración de Y va reduciéndose paulatinamente hasta que en torno a los 160 segundos de reacción alcanza niveles insignificantes mientras que la concentración de B alcanza su máximo de 5 mol/litro.

Es interesante observar que en todo momento se cumple el Principio de Conservación de la masa: la suma en cualquier momento de las concentraciones de las sustancias es 5 mol /litro. Esto es fácilmente comprobable en los momentos en que las concentraciones de A, X y B alcanzan su máximo de 5 mol/litro y el resto de concentraciones son nulas o despreciables.