Ver el código de «Reacción Autocatálisis A2»


== ''' Introducción''' ==
A continuación se presenta un estudio sobre las reacciones de autocatálisis, basadas en el principio de conservación de la materia y en la Ley de Acción de masas.

•	'''Principio de conservación de la materia''': En una reacción química ordinaria, la masa permanece constante, es decir, la masa consumida de los reactivos es igual a la masa obtenida de los productos.

•	'''Ley de Acción de masas''': establece que para una reacción química reversible en equilibrio a una temperatura constante, una relación determinada de concentraciones de reactivos y productos, tienen un valor constante.
 
El enunciado del estudio a realizar será: 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.
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

'''A+B->k12B'''

 
 
== '''Deducción Problema de valor inicial''' ==
Lo primero es denominar x (t) e y (t) como las concentraciones de los reactivos que aparecen  en la ecuación de arriba. Para encontrar la ecuación en función de una única variable y (t), se parte de la ley de conservación de la masa, obteniendo como resultado 

'''X (t)+y (t)=cte.=c''' 

Como ya hemos explicado anteriormente podemos obtener a partir de la Ley de Acción de masas que: 

'''y’(t)=k1*x(t)*y(t)''' 

Sustituyendo la variable x (t) de la primera ecuación en la segunda obtenemos la ecuación diferencial deseada: 

'''y′ (t)=k1∗ (c−y(t))∗y(t)'''

Si se parte de que la concentración inicial en t=0 es y0, obtendremos el problema de valor inicial pedido.

 

== '''Resolución Numérica PVI''' ==
 

Suponiendo que las concentraciones iniciales de A y B son 1 mol/l y 0.01 mol/l respectivamente, y k = 1 mol/s, resolver el PVI por el método de Euler, eligiendo un paso h = 0.1, en los primeros 10 segundos. Dibujar en una misma gráfica las concentraciones de las dos sustancias. Interpretar las gráficas.

•	'''Método 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\[ PVI = \begin{cases} y_0 \\ y_{n+1} = y_n + h*f(t_n,y_n) \\ \end{cases}\quad \] 

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




  <small>%DATOS DEL PROBLEMA
  %Concentración inicial sustancia A = 1 mol/l 
  %Concentración inicial sustancia B = 0.01 mol/l
  %k = 1 mol/s
  %h = 0.1 en los primeros diez segundos t=[0,10] s
  k=1;
  h=0.1;
  t0=0;
  tN=10;
  y0=0.01;
  %Definimos el vector tiempo con h=0.1
  t=t0:h:tN;
  %La constante del P.V.I. es la suma de las concentraciones %como resultado del Principio de conservación de la masa.
  %C=0.1+y0=1.01;
  y=zeros(1, length(t));
  %La variable y será la solución de la concentración B.
  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');
  for i=1:length(t)-1
  y(i+1)=y(i)+h*f(t(i),y(i)); %Aplicamos Euler
  end
  %Calculamos la concentración de A a partir de las %concentraciones de B y de C, como la diferencia entre %ambas (siendo la concentración de B "y").
  x=1.01-y;
  %Dibujamos ambas gráficas
  plot(t,y,'k x');
  hold on
  plot (t,x,'g x'); 
  hold off
  title('Gráfica Concentración-Tiempo');
  xlabel('tiempo(s)');
  ylabel('Concentración (mol/l)');
  legend('Concentración de B', 'Concentración de A', 'Location','Best');</small>





•	'''Método 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\[ 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 \] 

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




  <small>
  dy=k1*(c-y)*y   x0=1   y0=0.01   k1=1   h=0.1   [0,10]   c=1.01
  % Datos del problema
  k=1;
  c=1.01;
  t0=0;
  tN=10; 
  y0=0.01;
  h=0.1; % tamaño de paso
  t=t0:h:tN; 
  N=(tN-t0)/h; % número de intervalos
  y=zeros(1,length(t)); % matriz de 1 por N+1
  y(1)=y0;
  % RK4
  for i=1:N
     K1=k*(c-y(i))*y(i);
     K2=k*(c-y(i)+0.5*K1*h)*(y(i)+0.5*K1*h);
     K3=k*(c-y(i)+0.5*K2*h)*(y(i)+0.5*K2*h);
     K4=k*(c-y(i)+K3*h)*(y(i)+K3*h);
     y(i+1)=y(i)+h/6*(K1+2*K2+2*K3+K4); 
   end
   x=1.01-y;
   % gráfica
   hold on
   plot(t,y,'k x')
   plot(t,x,'g x')
   hold off
   legend('Concentración de B','Concentración de A','location','best')

</small>



•	'''Método Runge Kutta'''

Resolver el PVI numéricamente usando el método del trapecio y el método de Runge Kutta de cuarto orden con el mismo paso de tiempo.
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 \(k_i\), ahorrándonos tener que escribir innumerables veces la correspondiente expresión matemática.




<small>
  % dy=k1*(c-y)*y   x0=1   y0=0.01   k1=1   h=0.1   [0,10]   c=1.01
  % Datos del problema
  k=1;
  c=1.01;
  t0=0;
  tN=10; 
  y0=0.01;
  h=0.1; % tamaño de paso
  t=t0:h:tN; 
  N=(tN-t0)/h; % número de intervalos
  y=zeros(1,length(t)); % matriz de 1 por N+1
  y(1)=y0;
  % Trapecio
  for i=1:N
  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
  x=1.01-y;
  % gráfico
  hold on
  plot(t,y,'k x')
  plot(t,x,'g x')
  hold off
  legend('Concentración de B','Concentración de A','location','best')


</small>

== '''SISTEMA DE ECUACIONES''' == 

Otra forma de plantear la ecuación anteriormente estudiada es mediante un sistema de ecuaciones:

Este sistema es\[ \left \{ \begin{matrix} y'(t)=k_{1}x(t)y(t) \\ x'(t)=-k_{1}x(t)y(t) \end{matrix} \right . \] Por otra parte, definimos el problema de valor inicial asociado a este sistema, que es\[ \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 . \] siendo \( k=1\frac{mol}{s} \)

•	Aplicamos el método de '''Euler''', previamente:


 <small>%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
</small>

•	Aplicamos el método de '''Runge-Kutta''', que también hemos usado con anterioridad:

  <small>%Método de Runge-Kutta</small>
  <small>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

</small>
 
Podemos observar como la representación gráfica coincide con la obtenida en el apartado anterior, ya que se trata del mismo problema resuelto por métodos diferentes.
Un aspecto importante a la hora de aplicar Runge-Kutta a sistemas de ecuaciones diferenciales no lineales es que, al depender K2 de K1, K3 de K2, y K4 de K3, es necesario definir primero las Ki de ambas incógnitas antes que las Kj sindo i<j, de lo contrario, a la hora de aplicar el bucle, los valores de las Ki utilizados serían los de la anterior iteración del bucle, obteniendo, por lo tanto, resultados erróneos.

 
== '''CONCLUSIONES''' == 

En esta reacción biomolecular, tenemos inicialmente A en concentración \(1\frac{mol}{l}\) y B en concentración \(0.1\frac{mol}{l}\), en ella B actúa como catalizador, lo que implica que a medida que su concentración va aumentando, la velocidad de la reacción aumenta de manera exponencial, hasta llegar a un punto en que la cantidad de sustancia que queda por reaccionar es demasiado pequeña y la velocidad de reacción empieza a disminuir. 
La máxima velocidad de la reacción se alcanzará a los 4,5610 segundos del inicio de la reacción, instante en el que las concentraciones de A y B son iguales (0.5050 mol/L) a partir de ese momento, como ya hemos mencionado, la velocidad de la reacción empieza a disminuir, es decir, la pendiente de las curvas de las concentraciones se va acercando a la horizontal, haciéndose nula en el instante en que la concentración de A ha desaparecido prácticamente por completo (a los diez segundo de iniciarse la reacción), siendo éste el final de la reacción autocatalítica.
Podemos observar que nuestras funciones tienen asíntotas horizontales en 0 y en 1.01. Estas asíntotas nos muestran los valores máximos y mínimos que pueden presentar ambas concentraciones. 
Por otra parte, la función y(t) es siempre creciente, mientras que la función x(t) es siempre decreciente. El crecimiento o decrecimiento de las funciones nos muestra como la concentración de A al principio es muy alta, mientras que la de B es mínima y, 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 esta última aumenta a costa de la disminución de la otra (cumple la ley de concentración de masa). 
También podemos observar una simetría entre ambas funciones. La recta x=4.5610 es una recta de simetría, al igual que la y=0.5050. La simetría de las funciones nos muestra que B se produce al mismo ritmo que disminuye A, sumando ambos valores 1.01 en cualquier instante, lo que verifica el principio de la conservación de masa.

== '''Deducción e interpretación de las ecuaciones diferenciales a partir de la reacción consecutiva de Lotka.''' ==


Partimos de la reacción consecutiva propuesta por Lotka:

''A + X → k₁*2X''

''X + Y → k₂*2Y''

''Y → k₃*B'''

Con contantes: k₁, k₂ y k₃.

Donde A, X, B, Y son sustancias distintas. Observamos que las dos primeras ecuaciones son autocatalíticas, la reacción transcurre consumiendo A y originándose B, como se ve en la reacción global 

''A + B → k₁*2B.''

Para deducir las ecuaciones diferenciales nos basamos en los dos principios físicos ya nombrados, ley de acción de masas y principio de conservación de la masa.

Las concentración de las cuatro sustancias depende del tiempo, B(t),A(t),X(t),Y(t) pero las simplificamos a A,X,B,Y para un mejor manejo.

Para deducir las tres primeras ecuaciones nos apoyamos en la ley de acción de masas, que dice: la velocidad de reacción es proporcional al producto de la concentración de los reactivos.

Vemos que X e Y intervienen en las reacciones como producto y reactivo por tanto nos da lugar a X’=X₁’+X₂’ al igual que Y Y’=Y₁’+Y₂’ de ahí obtenemos:

''X’ = k₁*A*X - k₂*X*Y''

Los signos se establecen debido a que en la primera ecuación se produce X en cambio, en la segunda se consume para producir Y.

''Y’ = k₂*X*Y - k₃*Y''

Los signos siguen el mismo criterio que en el caso anterior donde en la segunda reacción se produce Y y en la tercera se va consumiendo.

La tercera ecuación a deducir es más sencilla puesto que B solo interviene en la última reacción como producto:

''B’ = k₃*Y''

Para la última ecuación diferencial a deducir partimos del 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):

''A+B+X+Y= cte.''

Lo que derivando nos da la siguiente ecuación:
''
X’+Y’+A’+B’= 0''

Quedando así deducidas las cuatro ecuaciones diferenciales pedidas, que se emplearan en los siguientes apartados.


== '''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,'c')
 plot(t, X, 'r')
 plot(t,Y, 'g')
 plot(t, B, 'k')
 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

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:
Como no se separan ambos métodos podemos afirmar que esestable.


 
== '''Resolución por Heun''' ==

 %Intorducimos datos del tiempo y calculamos su vector
 t0=0; tf=200;
 h=0.01; N=(tf-t0)/h;
 t=t0:h:tf;
 % introducimos datos aportados por el enunciado
 A0=5; X0=5*(10^-4); Y0=10^(-5); B0=0;
 %creamos las matrices A B X Y que tendrán guardados los puntos solución
 %para cada instante de tiempo
 A=zeros(1,N+1);
 X=zeros(1,N+1);
 Y=zeros(1,N+1);
 B=zeros(1,N+1);
 % introducimos los valores iniciales dentro de las matrices
 A(1)=A0; X(1)=X0; Y(1)=Y0; B(1)=B0;
 aa=A0; xx=X0; yy=Y0; bb=B0;
 % introducimos los valores de las constantes dadas en el enunciado
 K1=0.1; K2=K1; K3=(K1/2);
 
 % creamos un bucle for que nos va a ir solucionando el PVI por el método de
 % Heun para cada valor de t y metiendo este valor en cada uno de los
 % lugares correspondientes a cada matriz
 for n=1:N
    K1A=-K1*aa*xx;
    K2A=K1A+K1A*h;
    K1X=K1*aa*xx-K2*xx*yy;
    K2X=K1X+K1X*h;
    K1Y=K2*xx*yy-K3*yy;
    K2Y=K1Y+K1Y*h;
    K1B=K3*yy;
    K2B=K1B+K1B*h;
    aa=aa+0.5*h*(K1A+K2A);
    xx=xx+0.5*h*(K1X+K2X);
    yy=yy+0.5*h*(K1Y+K2Y);
    bb=bb+0.5*h*(K1B+K2B);
    A(n+1)=aa;
    X(n+1)=xx;
    Y(n+1)=yy;
    B(n+1)=bb;
    end
 % una vez tenemos todas las matrices completas dibujamos, para ello usa bos
 % la función plot. Para poder dibujar varias gráficas dentro de una misma
 % pestaña usamos el comando hold on // hold off
 plot(t,A,'*r');
 hold on
 plot(t,X,'*g');
 plot(t,Y,'*y');
 plot(t,B,'*');
 hold off
 xlabel('Tiempo(s)','FontSize',11);
 ylabel('concentración (mol/L)','Fontsize',11);
 title('CONCENTRACIÓN - TIEMPO (Heun)','FontName','Berlin sans FB','FontSize',11);
 legend('concentración de A','concentración de X','contenctración de Y','concentración de B','location','best');
 %xlabel, ylabel, title y legend son comandos que nos permiten escribir en
 %la gráfica qué es cada eje, un titulo y la legenda para la mejor
 %comprensión de la gráfica.