Diferencia entre revisiones de «Reacciones complejas GRUPO 1A»

El código de Matlab para el método del trapecio es el siguiente:

6.1.1 Código en Matlab

%Trapecio
clear all
%Datos del problema
t0=0; %Instante inicial
tN=2; %Instante final
h=0.1; %Tamaño de paso
k1=1; %valor de la constante específica de la velocidad de reacción
y0=0;
a0=3; %Concentración inicial de A
b0=1; %Concentración inicial de B
%Discretización
N=round((tN-t0)/h);
%Creamos el vector t
t=t0:h:tN;
%Preparación del vector solución y defino primer valor
y=zeros(1,N+1);
y(1)=y0;
%Aplicación del esquema numérico
for i=1:N
y(i+1)=(-h)\((-2*h-1)+sqrt((((2*h)+1)^2)-2*h*(y(i)+(h/2)*(6-(4*y(i))+(y(i)^2)))));
end
%grafico
hold on
plot(t,y,'b','linewidth',2);
a=3.-y;
b=1.-y;
plot(t,(a0-y),'y','linewidth',2);
plot(t,(b0-y),'m','linewidth',2);
legend('C','A','B')
xlabel('tiempo(s)');
ylabel('concentracion[mol/l]');
hold off

6.1.2 Gráfica del Método del Trapecio

Vemos que esta gráfica es prácticamente igual a la del método de Euler, pero sabemos que esta gráfica nos da una aproximación mayor.

6.2 Método de Runge Kutta

El método de Runge Kutta, de cuarto orden, al contrario del método del Trapecio, se trata de un método explícito, igual que el de Euler, por lo que no es necesario hacer ninguna manipulación al esquema numérico, por lo que procederemos directamente. El código Matlab utilizado para el método de Runge-Kutta es el siguiente:

6.2.1 Código en Matlab

%RUNGEKUTTA
clear all
%Datos del problema
t0=0; %Instante inicial
tN=2; %Instante final
h=0.1; %Tamaño de paso
k1=1; %valor de la constante específica de la velocidad de reacción 
y0=0;
a0=3; %Concentración inicial de A
b0=1; %Concentración inicial de B
%Discretización
N=round((tN-t0)/h);
%Creamos el vector t
t=t0:h:tN;
%Preparación del vector solución y defino primer valor
y=zeros(1,N+1);
y(1)=y0;
%Aplicación del esquema numérico
for i=1:N
K1=k1*(a0-y(i))*(b0-y(i));
K2=k1*(a0-(y(i)+(h/2)*K1))*(b0-(y(i)+(h/2)*K1));
K3=k1*(a0-(y(i)+(h/2)*K2))*(b0-(y(i)+(h/2)*K2));
K4=k1*(a0-(y(i)+h*K3))*(b0-(y(i)+h*K3));
y(i+1)=y(i)+(h/6)*(K1+2*K2+2*K3+K4);
end
%grafico
hold on
plot(t,y,'b','linewidth',2);
a=3.-y;
b=1.-y;
plot(t,(a0-y),'y','linewidth',2);
plot(t,(b0-y),'m','linewidth',2);
legend('C','A','B')
xlabel('tiempo(s)');
ylabel('concentracion[mol/l]');
hold off

6.2.2 Gráfica del Método de RUNGE KUTTA

Una vez que hemos obtenido las tres gráficas de los métodos de Euler, Trapecio y Runge-Kutta, podemos observar que todos tienen una gran similitud, pero el método de Euler, se separará más de los demás debido a que es un método de orden 1 y su precisión es más reducida que el resto. Los otros dos métodos serán más similares entre ellos puesto que tienen un mayor orden que el de Euler, el método del trapecio es de segundo orden y el de Runge-Kutta es de cuarto orden, por lo que el método que tendrá una mayor precisión será el de Runge-Kutta, gracias a su orden.

7 Reacción consecutiva

Ahora, suponemos que hay una reacción consecutiva de la siguiente forma:

                  A + B → C → D

Esta reacción, al igual que la reacción anterior, se puede expresar en forma de problema de valor inicial. Para ello, según explica el enunciado, usaremos de nuevo la ley de acción de masas, pero esta vez lo haremos dos veces consecutivas, una para de pasar de A y B a C, y otra para pasar de C a D. Al imponerse la ley de acción de masas, suponemos que la velocidad con la que se forma compuesto C es directamente proporcional a las concentraciones de los compuestos A y B, y de igual forma la de formación de D lo es a la concentración de C que se ha formado.

El problema de valor inicial es un sistema de dos ecuaciones en el que la primera muestra la producción de C y la segunda muestra la producción de D, que está relacionada con la cantidad de producto C en función del tiempo. Los valores en t=0 son las concentraciones iniciales de los productos, es decir, nulos. El sistema sería el siguiente:

  y'1(t) =k1(a0-y(t))(b0-y(t))
  y'2(t) =k2(y1(t)-y2(t))

En el que :

  k1,como ya dijimos anteriormente, es la constante de proporcionalidad de la primera reacción
  k2 es la constante de proporcionalidad de la segunda reacción
  y'1(t) es la velocidad de la reacción química de C
  y'2(t) es la velocidad de la reacción química de D

Dado que los valores de las concentraciones en t=0 son nulos, es decir:

  y1 (t=0) = 0
  y2 (t=0) = 0   

8 Reacción consecutiva: obtención de la ecuación diferencial y PVI

Las reacciones irreversibles se pueden definir como aquellas que empiezan con un reactivo inicial y producen productos o intermediarios generalmente en una sola dirección. Las reacciones consecutivas son reacciones secuenciales irreversibles

Plantearemos  un sistema con  dos ecuaciones para obtener el PVI asociado. 

[math] \left . \begin{matrix} y'_1(t)={k_{1}}({a_{0}}-y_1(t))({b_{0}}-y_1(t))\\ y'_2(t)={k_{2}}(y_1(t)-y_2(t))\\ y_1(0)=0\\y_2(0)=0 \end{matrix} \right \} [/math] Observamos que la primera es la misma que obtuvimos anteriormente, ya que los compuestos A y B se siguen descomponiendo independientemente de lo que haya al otro lado de la igualdad, solo depende de la concentración de ambos compuestos. Y la segunda ecuación del sistema es de carácter parecido solo que, la concentración de C mas D es la antigua concentración de C cuando no era consecutiva la reacción. . La velocidad de reacción de este nuevo proceso químico consecutivo al anterior también será proporcional a la concentración de su reactivo, en este caso de C, y la velocidad de reacción del proceso anterior se verá reducida por la velocidad de formación del nuevo compuesto D.

k2 hace referencia a una nueva constante específica de velocidad de reacción, independiente de la primera ya que trata otra reacción química totalmente distinta, la de disgregación de C, que no tiene nada que ver con A y B.
===Resolución e interpretación del sistema por Euler y Runge-Kutta con k2=5===

8.1 Resolución numérica por Trapecio cuando k2=5

%Sistema de ecuaciones no lineal por Trapecio y RK4
%Datos del problema 
T0=0;
 TN=10; 
Y0=[0;0];
 k1=1;
 k2=5; 
a0=3;
 b0=1;

%Discretizamos para h=0.1
h=0.1
T=T0:h:TN;
N=round((TN-T0)/h);

%EULER
y=zeros(2,length(T));
y(:,1)=Y0;
for i=1:length(T)-1
  y(:,i+1) =y(:,i)+h*[k1*((a0-(y(1,i)))*(b0-(y(1,i))));k2*(y(1,i)-y(2,i))];
end

%Cálculamos la concentración real de C 
yC=y(1,:)-y(2,:);
%Representación Euler
figure(1)
hold on
plot(T,yC,'linewidth',2)
plot(T,y(2,:),'y','linewidth',2)
legend('Concentración de C','Concentración de D','location','best');
title('Solución mediante el método de Euler');
hold off

=== Resolución numérica por RK4 cuando k2=5 ===
%RK 4
z=zeros(2,length(T));
z(:,1)=Y0;%Inicializamos

for i=1:length(T)-1
    K1=[k1*(a0-z(1,i))*(b0-z(1,i));k2*(z(1,i)-z(2,i))];
    K2=[k1*(a0-(z(1,i)+(h/2)*K1(1)))*(b0-(z(1,i)+(h/2)*K1(1)));k2*((z(1,i)+(h/2)*K1(2))-(z(2,i)+(h/2)*K1(2)))];
    K3=[k1*(a0-(z(1,i)+(h/2)*K2(1)))*(b0-(z(1,i)+(h/2)*K2(1)));k2*((z(1,i)+(h/2)*K2(2))-(z(2,i)+(h/2)*K2(2)))];
    K4=[k1*(a0-(z(1,i)+h*K3(1)))*(b0-(z(1,i)+h*K3(1)));k2*((z(1,i)+h*K3(2))-(z(2,i)+h*K3(2)))];
    z(:,i+1)=z(:,i)+h/6*(K1+2*K2+2*K3+K4);
end

%Cálculamos  la concentración real de C  
zC=z(1,:)-z(2,:);

%Representación  RK4
figure(2)
hold on
plot(T,zC,'linewidth',2)
plot(T,z(2,:),'y','linewidth',2)
legend('Concentración de C','Concentración de D','location','best');
title('Solución mediante el método de Runge-Kutta');
hold off
=== Grafica conjunta de Euler y RK4 cuando k2=5 ===
%Gráfico conjunta
figure(3)
hold on
subplot(1,2,1)
hold on
plot(T,yC,'b','linewidth',2)
plot(T,zC,'y--','linewidth',1.5)
legend('Concentración de C(Euler)','Concentración de C(RK4)','location','best');
title('Gráfico comparativo para compuesto C')
hold off
subplot(1,2,2)
hold on
plot(T,y(2,:),'b','linewidth',1.5)
plot(T,z(2,:),'y--','linewidth',1.5)
legend('Concentración de D(Euler)','Concentración de D(RK4)','location','best');
title('Gráfico comparativo para compuesto D')
hold off
hold off

Comentarios de las graficas

Se observa como el compuesto C comienza con gran velocidad. En el momento en que empieza a haber C se  generara D, que como tiene una velocidad de creación aún mayor, en cuanto pasa un poco de tiempo el compuesto D termina por hacer desaparecer al C y el compuesto D sigue formandose hasta alcanzar el equilibrio. Tambien se observa (ambos métodos) que en los primeros instantes el compuesto C se empieza a formar, mientras que el D no hasta que que pasa un tiempo debido a un error de los métodos. En conclusión vemos como se van formando ambos compuestos teniendo una mayor velocidad D

. %Datos del problema T0=0;

TN=10; 

Y0=[0;0];

k1=1;
k2=5; 

a0=3;

b0=1;

%Discretizamos para h=0.3 h=0.3 T=T0:h:TN; N=round((TN-T0)/h);

%EULER y=zeros(2,length(T)); y(:,1)=Y0; for i=1:length(T)-1

 y(:,i+1) =y(:,i)+h*[k1*((a0-(y(1,i)))*(b0-(y(1,i))));k2*(y(1,i)-y(2,i))];

end

%Cálculamos la concentración real de C yC=y(1,:)-y(2,:); %Representación Euler figure(1) hold on plot(T,yC,'linewidth',2) plot(T,y(2,:),'y','linewidth',2) legend('Concentración de C','Concentración de D','location','best'); title('Solución mediante el método de Euler'); hold off

8.2 Resolución numérica por RK4 cuando k2=5

%RK 4 z=zeros(2,length(T)); z(:,1)=Y0;%Inicializamos

for i=1:length(T)-1

   K1=[k1*(a0-z(1,i))*(b0-z(1,i));k2*(z(1,i)-z(2,i))];
   K2=[k1*(a0-(z(1,i)+(h/2)*K1(1)))*(b0-(z(1,i)+(h/2)*K1(1)));k2*((z(1,i)+(h/2)*K1(2))-(z(2,i)+(h/2)*K1(2)))];
   K3=[k1*(a0-(z(1,i)+(h/2)*K2(1)))*(b0-(z(1,i)+(h/2)*K2(1)));k2*((z(1,i)+(h/2)*K2(2))-(z(2,i)+(h/2)*K2(2)))];
   K4=[k1*(a0-(z(1,i)+h*K3(1)))*(b0-(z(1,i)+h*K3(1)));k2*((z(1,i)+h*K3(2))-(z(2,i)+h*K3(2)))];
   z(:,i+1)=z(:,i)+h/6*(K1+2*K2+2*K3+K4);

end

%Cálculamos la concentración real de C zC=z(1,:)-z(2,:);

%Representación RK4 figure(2) hold on plot(T,zC,'linewidth',2) plot(T,z(2,:),'y','linewidth',2) legend('Concentración de C','Concentración de D','location','best'); title('Solución mediante el método de Runge-Kutta'); hold off

8.3 Grafica conjunta de Euler y RK4 cuando k2=5

%Gráfico conjunta figure(3) hold on subplot(1,2,1) hold on plot(T,yC,'b','linewidth',1.5) plot(T,zC,'y--','linewidth',1.5) legend('Concentración de C(Euler)','Concentración de C(RK4)','location','best'); title('Gráfico comparativo para compuesto C') hold off subplot(1,2,2) hold on plot(T,y(2,:),'b','linewidth',1.5) plot(T,z(2,:),'y--','linewidth',1.5) legend('Concentración de D(Euler)','Concentración de D(RK4)','location','best'); title('Gráfico comparativo para compuesto D') hold off hold off}} Para h=0.3 observamos que la gráfica se desequilibra y deducimos de ello que es mejor utilizar intervalos pequeños en este caso 0.1 ya que como hemos visto con 0.3 hace que la gráfica se empiece a separar de la solución real

8.4 Resolución numérica por Trapecio y RK4 cuando k2=1/5

%Sistema de ecuaciones no lineal por Trapecio y RK4
%Datos del problema 
t0=0; 
tN=10;
 y0=[0;0];
 k1=1;
 k2=1/5; 
a0=3;
 b0=1;

%Discretización 
h=0.1; 
t=t0:h:tN;
N=round((tN-t0)/h); 

%EULER 
y=zeros(2,length(t));
y(:,1)=y0; %Inicialización
for i=1:length(t)-1
  y(:,i+1) =y(:,i)+h*[k1*((a0-(y(1,i)))*(b0-(y(1,i))));k2*(y(1,i)-y(2,i))];
end

%Cálculamos  la concentración real de C 
yC=y(1,:)-y(2,:);

%Representación (Euler)
figure(1)
hold on
plot(t,yC,'linewidth',2)
plot(t,y(2,:),'y','linewidth',2)
legend('Concentración de C','Concentración de D','location','best');
title('Solución mediante el método de Euler');
hold off

%RK4  
z=zeros(2,length(t));
z(:,1)=y0;%Inicializamos

for i=1:length(t)-1
    K1=[k1*(a0-z(1,i))*(b0-z(1,i));k2*(z(1,i)-z(2,i))];
    K2=[k1*(a0-(z(1,i)+(h/2)*K1(1)))*(b0-(z(1,i)+(h/2)*K1(1)));k2*((z(1,i)+(h/2)*K1(2))-(z(2,i)+(h/2)*K1(2)))];
    K3=[k1*(a0-(z(1,i)+(h/2)*K2(1)))*(b0-(z(1,i)+(h/2)*K2(1)));k2*((z(1,i)+(h/2)*K2(2))-(z(2,i)+(h/2)*K2(2)))];
    K4=[k1*(a0-(z(1,i)+h*K3(1)))*(b0-(z(1,i)+h*K3(1)));k2*((z(1,i)+h*K3(2))-(z(2,i)+h*K3(2)))];
    z(:,i+1)=z(:,i)+h/6*(K1+2*K2+2*K3+K4);
end

%Cálculamos  la concentración real de C 
zC=z(1,:)-z(2,:);

%Representación De  (RK4)
figure(2)
hold on
plot(t,zC,'linewidth',2)
plot(t,z(2,:),'y','linewidth',2)
legend('Concentración de C','Concentración de D','location','best');
title('Solución mediante el método de Runge-Kutta');
hold off

%Gráfico conjunto
figure(3)
hold on
subplot(1,2,1)
hold on
plot(t,yC,'b','linewidth',1.5)
plot(t,zC,'y--','linewidth',1.5)
legend('Concentración de C(Euler)','Concentración de C(RK4)','location','best');
title('Gráfico comparativo para compuesto C')
hold off
subplot(1,2,2)
hold on
plot(t,y(2,:),'b','linewidth',1.5)
plot(t,z(2,:),'y--','linewidth',1.5)
legend('Concentración de D(Euler)','Concentración de D(RK4)','location','best');
title('Gráfico comparativo para compuesto D')
hold off
hold off