Diferencia entre revisiones de «Reacciones con autocatálisis. Grupo C2»

Revisión del 12:38 3 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,: [math] A + B \rightarrow C [/math] 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. En este ejercicio analizaremos el caso particular en el que A se transforma en B pero suponiendo que la presencia de B hace de efecto catalítico en la reacción. Escribiremos este proceso como una reacción bimolecular: [math] A + B \rightarrow _{k1} 2B [/math]

2 Primera reacción

2.1 Interpretación

Teniendo que x e y son las concentraciones de los reactivos que intervienen en la reacción y tomando la reacción que describe la transformación de A y B para producir 2B a una velocidad k1. Basándonos en la ley de la conservación de la masa (suma de concentraciones es siempre constante). [math] x(t) + y(t) = cte [/math] Y derivando en ambas partes con respecto al tiempo: [math] x'(t) + y'(t) = 0 [/math] Por otro lado, usando la ley de conservación de masas, que establece que la velocidad de reacción es proporcional al producto de las concentraciones de los reactivos, en este caso el factor de proporcionalidad es k1. [math] y'(t) = k_{1}*x(t)*y(t) [/math] De la primera igualdad obtenemos que: [math] x(t) = cte - y(t) [/math] y sustituyendo en la igualdad de la ley de conservación de masas resulta: [math] y'(t) = k_{1}*(cte - y(t))*y(t) [/math]

2.2 Ecuación diferencial

2.2.1 Método de Euler

Vamos a resolver el problema por el método de Euler: [math] \begin{array}{c} y_n\\ y_{n+1}=y_n+h*f(t_n,y_n)\end{array} [/math] Aplicado a nuestra ecuación y sustituyendo por los datos del problema, tendríamos que: [math] \begin{array}{c} y_n\\ y_{n+1}=y_n+h*[y_n*(1.01-y_n)]\end{array} [/math] Y resolvemos:

Gráfica Método de Euler h=0.1

%Datos:

%Concentración A: 1mol/L,
%Concentración B: 0.01mol/L,
%k=1mol/s, h=0.1, t=[0,10] (s)

%-----------------------------------------------------------------------

k=1;
t0=0; tN=10;
h=0.1; y0=0.01;
%Definimos el vector del tiempo, con un paso de 0.1.
t=t0:h:tN;
%La constante C, del P.V.I., es la suma de las concentraciones (Principio
%de Conservación de la Masa).
% C=0.1+y0; ----> C=1.01;
%Preasignamos el vector y, que será nuestra solución de la concentración B.
y=zeros(1,length(t));
%Introducimos el valor incial en la primera posición del vector.
y(1)=y0;
%Definimos la función, para aplicarla a continuación en el Método de Euler.
F=inline('y*(1.01-y)','t','y');
%Recorremos, con un bucle, nuestro vector y, y en él almacenamos los
%resultados obtenidos con Euler.
for i=1:length(t)-1   
    y(i+1)=y(i)+h*F(t(i),y(i)); %Euler.
end
%De nuevo, por el Principio de Conservación de la Masa,como la suma de las
%concetraciones será C (cte.), podemos calcular la concentración de A como 
%la diferencia entre C y la concentración de B (el vector y).
x=1.01-y;
%Dibujamos ambas gráficas.
plot(t,y,'r');
hold on
plot(t,x,'g');
hold off
title('Gráfica Concentración-Tiempo');
xlabel('tiempo (segundos)');
ylabel('concentración (mol/L)');
legend('Concentración de B','Concentración de A','Location','Best');

En la gráfica se puede observar que la concentración inicial de la sustancia A es de 1 [math]\frac{mol}{L}[/math] , mientras que la de B es de 0.01 [math] \frac{mol}{L} [/math]. Debido a que B hace efecto catalítico en la reacción, las curvas representativas de las concentraciones, creciente para la sustancia B (la concentración es mayor que la de A con el paso del tiempo) y decreciente para la sustancia A, son exponenciales. Una vez pasados los primeros 6-7 segundos, la velocidad de la reacción disminuye pues queda poca concentración de A para reaccionar y, a los 10 segundos, prácticamente todo la cantidad de sustancia es de B. Podemos calcular las cantidades finales concretas:

%Concentraciones transcurridos 10 segundos de la reacción:
CF_A=x(length(x)); CF_B=y(length(y));
fprintf('La concentración final de A es de %.4f mol/L,y la de B de %.4f mol/L\n',CF_A,CF_B);

Cuyo resultado en pantalla será:

La concentración final de A es de 0.0039 mol/L,y la de B de 1.0061 mol/L

Sin embargo, no podemos concretar exactamente cuándo ambas concentraciones son iguales. Esto es debido a la discretización, pues las gráficas no se pueden pintar como curvas continuas, si no como puntos muy próximos y por tanto, no se puede decir exactamente el valor donde [math] \left [ A \right ] [/math] = [math] \left [ B \right ] [/math] , ya que lo más probable es que en ninguno de los vectores que representan el valor de las concentraciones coincidan. Sin embargo, podemos hacer una estimación del intervalo de tiempo donde ocurra, mirando la gráfica: el tiempo transcurrido será de entre 4.5 y 5 segundos. Para más precisión basta con ver los vectores de las concentraciones, y vemos que deben cortarse pasados entre 4.8 y 4.9 segundos, y la concentración estará entre 0.4999 y 0.5254 .

Sustancia B:

Columns 45 through 55

    0.4244    0.4493    0.4745    0.4999    0.5254    0.5508    0.5761    0.6011    0.6257    0.6497    0.6731

Sustancia A:

Columns 45 through 55

    0.5856    0.5607    0.5355    0.5101    0.4846    0.4592    0.4339    0.4089    0.3843    0.3603    0.3369

2.2.2 Método del Trapecio

Ahora vamos a resolver el mismo problema de valor inicial, pero esta vez por el método del trapecio: [math] \begin{array}{c} y_n \\ y_{n+1}=y_n+{h \over 2}*[f(t_n,y_n)+f(t_{n+1},y_{n+1})]\end{array} [/math] Como podemos observar, se trata de un método implícito. Esto quiere decir que nuestra incógnita depende de una función en la que aparece también. Para solucionarlo, aplicaremos el método a nuestra ecuación y despejaremos la incógnita ([math] y_{n+1} [/math] ) en función de lo demás: [math] \begin{array}{c}y_{n+1}=y_n+{h \over 2}*[y_n*(C-y_n)+y_{n+1}*(C-y_{n+1})]\\y_{n+1}=y_n+{h \over 2}*y_n*(C-y_n)+{h \over 2}*y_{n+1}*(C-y_{n+1})\\y_{n+1}-{h \over 2}*y_{n+1}*(C-y_{n+1})=y_n*[1+{h \over 2}*(C-y_n)]\\{h \over 2}*(y_{n+1})^2+y_{n+1}*(1-{C*h \over 2}+[-y_n*[1+{h \over 2}*(C-y_n)]])=0\\y_{n+1}={-(1-{C*h \over 2})+\sqrt{(1-{C*h \over 2})^2-4*{h \over 2}*[-y_n*[1+{h \over 2}*(C-y_n)]]} \over 2*{h \over 2}}\\y_{n+1}={-1+{C*h \over 2}+\sqrt{(1-{C*h \over 2})^2-2*h*[-y_n*[1+{h \over 2}*(C-y_n)]]} \over h}\end{array} [/math]

Gráfica Método del Trapecio h=0.1

%Datos:

%Concentración A: 0.1mol/L,
%Concentración B: 0.01mol/L,
%k=1mol/s, h=0.1, t=[0,10] (s)

%-----------------------------------------------------------------------

k=1;
t0=0; tN=10;
h=0.1; y0=0.01; c=1.01;
%Definimos el vector del tiempo, con un paso de 0.1.
t=t0:h:tN;
%La constante C, del P.V.I., es la suma de las concentraciones (Principio
%de Conservación de la Masa).
% C=0.1+y0; ----> C=1.01;
%Preasignamos el vector y, que será nuestra solución de la concentración B.
y=zeros(1,length(t));
%Introducimos el valor incial en la primera posición del vector.
y(1)=y0;
%Recorremos, con un bucle, nuestro vector y, y en él almacenamos los
%resultados obtenidos con el Método del Trapecio.
for i=1:length(t)-1   
    %Trapecio:
    y(i+1)=(1/(h*k))*((0.5*h*k*c-1)+sqrt((1-0.5*h*k*c)^2-2*h*k*(-y(i)-(h/2)*y(i)*(c-y(i)))));
end
%De nuevo, por el Principio de Conservación de la Masa,como la suma de las
%concetraciones será C (cte.), podemos calcular la concentración de A como 
%la diferencia entre C y la concentración de B (el vector y).
x=1.01-y;
%Dibujamos ambas gráficas.
plot(t,y,'r');
hold on
plot(t,x,'g');
hold off
title('Gráfica Concentración-Tiempo');
xlabel('tiempo (segundos)');
ylabel('concentración (mol/L)');
legend('Concentración de B','Concentración de A','Location','Best');

Se observa que la solución obtenida es la misma que por el método de Euler.

2.2.3 Método de Runge-Kutta

Otra vez, vamos a resolver el mismo problema de valor inicial, ayudándonos en este caso del método de Runge-Kutta de orden 4.

Gráfica Método de Runge-Kutta de cuarto orden h=0.1

%Datos:

%Concentración A: 1mol/L,
%Concentración B: 0.01mol/L,
%k=1mol/s, h=0.1, t=[0,10] (s)

%-----------------------------------------------------------------------

t0=0; tN=10;
h=0.1; y0=0.01;
%Definimos el vector del tiempo, con un paso de 0.1.
t=t0:h:tN;
%La constante C, del P.V.I., es la suma de las concentraciones (Principio
%de Conservación de la Masa).
% C=0.1+y0; ----> C=1.01;
%Preasignamos el vector y, que será nuestra solución de la concentración B.
y=zeros(1,length(t));
%Introducimos el valor incial en la primera posición del vector.
y(1)=y0;
%Definimos la función, para aplicarla a continuación en el Método de 
%Runge-Kutta.
F=inline('1*y*(1.01-y)','t','y');
%Recorremos, con un bucle, nuestro vector y, y en él almacenamos los
%resultados obtenidos con RK4.
for k=1:length(t)-1
    %Runge-Kutta de orden 4 (RK4):
    K1=F(t(k),y(k));
    K2=F(t(k)+(h/2),y(k)+h*K1/2);
    K3=F(t(k)+(h/2),y(k)+h*K2/2);
    K4=F(t(k)+h,y(k)+K3*h);
    y(k+1)=y(k)+(h/6)*(K1+2*K2+2*K3+K4);
end
%De nuevo, por el Principio de Conservación de la Masa,como la suma de las
%concetraciones será C (cte.), podemos calcular la concentración de A como 
%la diferencia entre C y la concentración de B (el vector y).
x=1.01-y;
%Dibujamos ambas gráficas.
plot(t,y,'r');
hold on
plot(t,x,'g');
hold off
title('Gráfica Concentración-Tiempo');
xlabel('tiempo (segundos)');
ylabel('concentración (mol/L)');
legend('Concentración de B','Concentración de A','Location','Best');

De nuevo, ambas soluciones coinciden.

2.3 Sistema de ecuaciones

Otra forma de plantear la resolución de la reacción bimolecular de autocatálisis anterior es plantear tanto la concentración de A como la de B como las variables de un sistema de ecuaciones diferenciales. Este sistema es: [math] \left \{ \begin{matrix} y'(t)=k_{1}x(t)y(t) \\ x'(t)=-k_{1}x(t)y(t) \end{matrix} \right . [/math] Por otra parte, definimos el problema de valor inicial asociado a este sistema, que es: [math] \left \{ \begin{matrix} y'(t)=k_{1}x(t)y(t) \\ x'(t)=-k_{1}x(t)y(t) \\ x(0)=1 \\ y(0)=0.01 \end{matrix} \right . [/math] siendo [math] k=1\frac{mol}{s} [/math]

2.3.1 Método de Euler

El método de Euler, se basa en la fórmula expuesta en el apartado 2.2.1. En este caso, al ser el sistema de ecuaciones no lineal, no podemos aplicar el método usando la técnica de la matriz explicada en las sesiones de numérico, siendo necesario por lo tanto aplicar el método en cuestión a cada ecuación por separado. El código del programa es el siguiente:

Gráfica por el método de Euler con paso h=0.1

%Datos:

%Concentración A: 1mol/L,
%Concentración B: 0.01mol/L,
%k=1mol/s, h=0.1, t=[0,10] (s)

%-----------------------------------------------------------------------

k=1;
t0=0; tN=10;
h=0.1; y0=0.01; x0=1;
%Definimos el vector del tiempo, con un paso de 0.1.
t=t0:h:tN;
%La constante C, del P.V.I., es la suma de las concentraciones (Principio
%de Conservación de la Masa).
% C=0.1+y0; ----> C=1.01;
%Preasignamos los vectores x e y, que serán nuestras soluciones de la 
%concentración A y B, respectivamente.
y=zeros(1,length(t));
x=y;
%Introducimos los valores inciales de cada concentración en la primera
%posición de los vectores.
y(1)=y0;
x(1)=x0;
%Definimos ambas funciones, para aplicarlas a continuación en el Método de
%Euler.
F=inline('x*y','t','y','x');
G=inline('-x*y','t','y','x');
%Recorremos, con un bucle, los vectores,y en ellos almacenamos los
%resultados obtenidos con Euler.
for i=1:length(t)-1   
    %Euler:
    y(i+1)=y(i)+h*F(t(i),y(i),x(i));
    x(i+1)=x(i)+h*G(t(i),y(i),x(i));
end
%Dibujamos ambas gráficas.
plot(t,y,'r');
hold on
plot(t,x,'g');
hold off
title('Gráfica Concentración-Tiempo');
xlabel('tiempo (segundos)');
ylabel('concentración (mol/L)');
legend('Concentración de B','Concentración de A','Location','Best');

Sin embargo, podemos usar un método alternativo. Podemos crear en primer lugar un archivo .m, en concreto una función, que nos servirá para escribir nuestras ecuaciones de manera más sencilla. El problema lo vamos a abordar, por tanto, vectorialmente. Este método es más cómodo cuando tratamos con sistemas de muchas ecuaciones.

% Empezamos creando el archivo sys_Euler_C2.m , donde definimos nuestro
% sistema:
function syst = sys_Euler_C2(t,y)
dy1=y(2)*y(1);
dy2=-y(2)*y(1);
syst=[dy1;dy2];
end

Gráfica por el método de Euler con paso h=0.1

%Y creamos nuestro programa, haciendo una llamada a la función anterior.
%Definiciones previas.
k=1;
t0=0; tN=10;
h=0.1; y0=0.01; x0=1;
%Definimos el vector del tiempo, con un paso de 0.1.
t=t0:h:tN;
%Creamos el vector donde tendremos nuestras condiciones iniciales.
yini=[y0;x0];
%Preasignamos la matriz Y, que será nuestra solución de la 
%concentración A y B,en la primera y segunda fila, respectivamente.
Y=zeros(2,length(t));
%Damos las condiciones iniciales a la primera columna.
Y(:,1)=yini;
%Recorremos, con un bucle, la matriz (por columnas),y almacenamos los
%resultados obtenidos con Euler.
for k=1:length(t)-1   
     Y(:,k+1)=Y(:,k)+h*sys_Euler_C2(t(k),Y(:,k));
end
%Dibujamos ambas gráficas.
plot(t,Y(1,:),'r');
hold on
plot(t,Y(2,:),'g');
hold off
title('Gráfica Concentración-Tiempo');
xlabel('tiempo (segundos)');
ylabel('concentración (mol/L)');
legend('Concentración de B','Concentración de A','Location','Best');

Por supuesto, la gráfica es exactamente la misma.

2.3.2 Método de Runge-Kutta

Ahora vamos a aplicar el método de Runge-Kutta de orden 4 a la ecuación expuesta anteriormente. Este método, al ser de orden superior al método de Euler, nos ofrecerá una mayor aproximación a la solución real mediante el cálculo numérico de esta. El código MatLab del programa es el que se muestra a continuación.

Gráfica por el método de Runge-Kutta con paso h=0.1

%Datos:

%Concentración A: 1mol/L,
%Concentración B: 0.01mol/L,
%k=1mol/s, h=0.1, t=[0,10] (s)

%-----------------------------------------------------------------------

k=1;
t0=0; tN=10;
h=0.1; y0=0.01; x0=1;
%Definimos el vector del tiempo, con un paso de 0.1.
t=t0:h:tN;
%La constante C, del P.V.I., es la suma de las concentraciones (Principio
%de Conservación de la Masa).
% C=0.1+y0; ----> C=1.01;
%Preasignamos los vectores x e y, que serán nuestras soluciones de la 
%concentración A y B, respectivamente.
y=zeros(1,length(t));
x=y;
%Introducimos los valores inciales de cada concentración en la primera
%posición de los vectores.
y(1)=y0;
x(1)=x0;
%Definimos ambas funciones, para aplicarlas a continuación en el Método de
%Runge-Kutta.
F=inline('x*y','t','y','x');
G=inline('-x*y','t','y','x');
%Recorremos, con un bucle, los vectores, y en ellos almacenamos los
%resultados obtenidos con RK4.
for k=1:length(t)-1
    %Runge-Kutta de orden 4 (RK4):
    K1_y=F(t(k),y(k),x(k));
    K1_x=G(t(k),y(k),x(k));
    K2_y=F(t(k)+(h/2),y(k)+h*K1_y/2,x(k)+h*K1_x/2);
    K2_x=G(t(k)+(h/2),y(k)+h*K1_y/2,x(k)+h*K1_x/2);
    K3_y=F(t(k)+(h/2),y(k)+h*K2_y/2,x(k)+h*K2_x/2);
    K3_x=G(t(k)+(h/2),y(k)+h*K2_y/2,x(k)+h*K2_x/2);
    K4_y=F(t(k)+h,y(k)+K3_y*h,x(k)+K3_x*h);
    K4_x=G(t(k)+h,y(k)+K3_y*h,x(k)+K3_x*h);
    
    y(k+1)=y(k)+(h/6)*(K1_y+2*K2_y+2*K3_y+K4_y);
    x(k+1)=x(k)+(h/6)*(K1_x+2*K2_x+2*K3_x+K4_x);
end
%Dibujamos ambas gráficas.
plot(t,y,'r');
hold on
plot(t,x,'g');
hold off
title('Gráfica Concentración-Tiempo');
xlabel('tiempo (segundos)');
ylabel('concentración (mol/L)');
legend('Concentración de B','Concentración de A','Location','Best');

Podemos observar como la representación gráfica vuelve a coincidir con los apartados anteriores, ya que se trata del mismo problema resuelto por caminos diferentes.

Un aspecto importante a la hora de aplicar este método a sistemas de ecuaciones diferenciales no lineales es que, al depender K2 de K1, K3 de K2, y K4 de K3, es necerario definir primero las K1 de ambas incógnitas, luego las K2 y así sucesivamente. Si no, a la hora de aplicar el bucle, si definimos primero las Ki de la incógnita x(t) y luego las Kj de la incógnita y(t), los valores de las Kj de y(t) serían los de la anterior iteración del bucle. realizándose por lo tanto un cálculo erróneo.

2.4 Conclusión

Una vez resuelto el problema de la reacción de autocátálisis [math] A + B \rightarrow _{k1} B [/math], procedemos a interpretar los resultados.

Tal y como se observó al despejar y(t) para pasar de un sistema de ecuaciones diferenciales a una sola ecuación, la concentración tanto de A como de B se rige por funciones logísticas del tipo: [math]\frac{dy}{dt}=ry\left(1 - \frac{y}{K}\right)[/math] En nuestro problema, las constantes r y K toman los siguientes valores: [math] r=k_{1} \cdot K= 1.01 [/math] : [math] K=c=1.01 [/math] Esta función logística se caracteriza, tanto como para x(t) como para y(t), por tener asíntotas horizontales en 0 y en 1.01. Además la función y(t) es siempre creciente, mientras que la función x(t) es decreciente. También podemos observar una simetría entre ambas funciones, siendo la recta que define el eje de simetría x=0.505, lo que también nos dice que los valores de ambas funciones siempre suman 1.01.

Toda esta información que nos proporciona el gráfico tiene una interpretación química clara. Las asíndotas nos muestran los valores máximos y mínimos que pueden presentar ambas concetraciones. Por otra parte, el crecimiento o decrecimiento de ambas funciones nos muestra como la concentracion de x(t) al principio de la reacción es muy alta, mientras que la de y(t) es mínima. Según transcurre el tiempo, las moléculas de A se van transformando en moléculas de B, por lo que la concentración de y(t) va aumentando a costa de la disminución de x(t). Por último, la simetría se las funciones nos muestra que y(t) se produce al mismo ritmo que disminuye x(t), sumando ambos valores 1.01 en cualquier instante, lo que verifica el principio de la conservación de la masa.

3 Segunda reacción

3.1 Interpretación

Consideramos la reacción consecutiva propuesta por Lotka en 1920.

[math] A + X \rightarrow _{k1} 2X [/math]

[math] X+ Y \rightarrow _{k2} 2Y [/math]

[math] Y \rightarrow _{k3} B [/math]

Tomando A, X, B, e Y en la resolución de las ecuaciones diferenciales como las concentraciones de las diferentes sustancias. Esta reacción consecutiva describe la transformación de A para producir B, estando controladas la velocidad y la mezcla de este proceso por las reacciones autocatalíticas en las que participan X e Y. Basándonos de nuevo en el principio de de conservación de la masa, partimos de la ecuación de que la suma de las concentraciones de todas las sustancias ha de ser constante: [math] A + x + y + B = cte [/math] Derivando esta ecuación deducimos que: [math] A' + x' + y' + B' = 0 [/math] Por otra parte, la ley de acción de masas nos indica que la velocidad de una reacción es proporcional al producto de las concentraciones de los reactivos. De esta forma, podemos obtener el resto de ecuaciones, teniendo en cuenta que tanto la sustancia x como la sustancia y no solo se producen, sino que también se consumen, por lo que habrá que restar a la parte consumida la parte producida Todo esto nos da lugar al siguiente sistema de ecuaciones: [math] \left \{ \begin{matrix} x' = k1Ax − k2xy \\ y' = k2xy − k3y \\ B' = k3y \\ A' + x' + y' + B' = 0 \end{matrix} \right . [/math] Por lo que podemos ver, a modo de mayor aclaración a lo explicado en el párrafo anterior, en la reacción autocatalítica [math] A + X \rightarrow _{k1} 2X [/math] se forma X, ya que es positivo, con una velocidad de reacción k1. De la misma manera observamos que en la reacción [math] X+ Y \rightarrow _{k2} 2Y [/math] X se consume con una velocidad de reacción k2.

La reacción [math] A' + x' + y' + B' = 0 [/math] está en función de las otras, por lo que sustituyendo llegamos a que [math]A'=-k1Ax [/math]

En conclusión vamos a resolver el problema de valor inicial tomando k1=k2=2k3=0.1 y las condiciones iniciales propuestas: [math] \left \{ \begin{matrix} x' = 0.1Ax − 0.1xy \\ y' = 0.1xy − 0.05y \\ B' = 0.05y \\ A'=-0.1Ax \\ A(0)=5 \\ X(0)=0.0005 \\ Y(0)=0.00001 \\ B(0)=0 \end{matrix} \right . [/math]

3.2 Método de Euler

Gráfica método de Euler h=0.01

Gráfica método de Euler h=0.001

%Datos del problema de valor inicial
t0=0;
tN=200;
A0=5;
X0=5*(10^-4);
Y0=10^(-5);
B0=0;
k1=0.1;
k2=0.1;
k3=0.05;
h=input('Introduzca un valor para h:');

%Vector de tiempo
t=t0:h:tN;

%Definir el número de intervalos
N=(tN-t0)/h;

%Vectores vacíos para nuestra solución
A=zeros(1,N+1);
X=zeros(1,N+1);
Y=zeros(1,N+1);
B=zeros(1,N+1);

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

%Resolvemos el sistema

for i=1:N
X(i+1)=X(i)+h.*(k1.*A(i).*X(i)-k2.*X(i).*Y(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.*X(i).*A(i));
end

%Gráfica de resultados
hold on
plot(t,X)
plot(t,Y,'g')
plot(t,B,'r')
plot(t,A,'y')
hold off
legend('Concentración X','Concentración de Y','Concentración de B','Concentración de A','Location','best');

3.3 Método de Heun

Gráfica método de Heun

%Datos del problema de valor inicial
t0=0;
tN=200;
A0=5;
X0=5*(10^-4);
Y0=10^(-5);
B0=0;
k1=0.1;
k2=0.1;
k3=0.05;
h=0.01;

%Vector de tiempo
t=t0:h:tN;

%Definir el número de intervalos
N=(tN-t0)/h;

%Vectores vacíos para nuestra solución
A=zeros(1,N+1);
X=zeros(1,N+1);
Y=zeros(1,N+1);
B=zeros(1,N+1);

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

%Definimos las constantes
for i=1:N
K1X=k1.*A(i).*X(i)-k2.*X(i).*Y(i);
K1Y=k2.*X(i).*Y(i)-k3.*Y(i);
K1B=k3.*Y(i);
K1A=-k1.*X(i).*A(i);
K2X=k1.*(A(i)+K1X.*h).*(X(i)+K1X.*h)-k2.*(X(i)+K1X.*h).*(Y(i)+K1X.*h);
K2Y=k2.*(X(i)+K1Y.*h).*(Y(i)+K1Y.*h)-k3.*(Y(i)+K1Y.*h);
K2B=k3.*(Y(i)+K1B.*h);
K2A=-k1.*(X(i)+K1A.*h).*(A(i)+K1A.*h);

X(i+1)=X(i)+0.5*h.*(K1X+K2X);
Y(i+1)=Y(i)+0.5*h.*(K1Y+K2Y);
B(i+1)=B(i)+0.5*h.*(K1B+K2B);
A(i+1)=A(i)+0.5*h.*(K1A+K2A);
end

%Gráficas
hold on
plot(t,X);
plot(t,Y,'g');
plot(t,B,'r');
plot(t,A,'y');
hold off
legend('Concentración X','Concentración de Y','Concentración de B','Concentración de A','Location','best');

3.4 Conclusión