Diferencia entre revisiones de «Reacciones con Autocatálisis Grupo A17»
(→Método del Trapecio) |
|||
| Línea 103: | Línea 103: | ||
legend('Concentracion y','Concentracion x','location','best') | legend('Concentracion y','Concentracion x','location','best') | ||
}} | }} | ||
| + | [[Archivo:Graficapartado2b.jpg|marco|derecha|Método del Trapecio]] | ||
====Método del Trapecio==== | ====Método del Trapecio==== | ||
{{matlab|codigo= | {{matlab|codigo= | ||
| Línea 132: | Línea 133: | ||
legend('Concentracion y','Concentracion x','location','best') | legend('Concentracion y','Concentracion x','location','best') | ||
}} | }} | ||
| + | |||
====Método de Runge-Kutta==== | ====Método de Runge-Kutta==== | ||
{{matlab|codigo= | {{matlab|codigo= | ||
Revisión del 18:24 6 mar 2015
| Trabajo realizado por estudiantes | |
|---|---|
| Título | Reacciones con Autocatálisis. Grupo A17 |
| Asignatura | Ecuaciones Diferenciales |
| Curso | Curso 2014-15 |
| Autores | Daniel Diez Sanz, Jorge Fernández Mendoza, Guillermo Mella Martínez |
| Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura | |
Contenido
1 Introducción
En el presente informe se plantea un estudio de las concentraciones de los productos y
reactivos de una reacción química, basándose en la Ley de Acción de masas y el Principio de
conservación de la materia, basados en que la materia ni se crea ni se destruye, y por tanto se
deberá mantener constante entre los productos y los reactivos, así como que la velocidad de la
reacción es proporcional al producto de las concentraciones de los reactivos.
Para ello se usarán los métodos numéricos de Euler, Trapecio, Runge-Kutta y Heun en el
programa informático Matlab y se mostrarán las soluciones a través de gráficos en una misma
imagen que permitan comparar las evoluciones de las concentraciones de los elementos en el
tiempo.
A continuación se muestra el enunciado del ejercicio a realizar.
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.
A + B → C.
2 REACCIÓN 1
2.1 Interpretación del problema y deducción de las ecuaciones diferenciales con valor inicial
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 →k1 2B
En primer lugar se establece x(t) e y(t) como las concentraciones de los reactivos que aparecen
en la ecuación de arriba. Para deducir la ecuación en función de una única variable y(t), se
parte de la ley de conservación de la masa, que establece que la suma de concentraciones es
siempre constante.
x(t)+y(t)=cte=c
Asimismo se conoce la ley de conservación de masas que nos proporciona la velocidad de
reacción k1 para las concentraciones de los reactivos.
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 a la misma, se le añade que y(0)=y0, se dispondría de un PVI en función de las
concentraciones iniciales y la velocidad de reacción que se irán introduciendo en los
sucesivos apartados.
2.2 Resolución numérica del PVI
2.2.1 Método de Euler
% dy=k1*(c-y)*y x0=1 y0=0.01 k1=1 h=0.1 [0,10] c=1.01
% Datos del problema
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;
% Euler
for i=1:N
y(i+1)=y(i)+h*(y(i)*(1.01-y(i)));
end
x=1.01-y;
% grafico
hold on
plot(t,y,'r')
plot(t,x,'g')
hold off
legend('Concentracion y','Concentracion x','location','best')2.2.2 Método del Trapecio
% 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;
% grafico
hold on
plot(t,y,'r')
plot(t,x,'g')
hold off
legend('Concentracion y','Concentracion x','location','best')
2.2.3 Método de Runge-Kutta
% 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;
% grafico
hold on
plot(t,y,'r')
plot(t,x,'g')
hold off
legend('Concentracion y','Concentracion x','location','best')2.3 Resolver en forma de sistema por medio de Euler y RK4
% dy=k*x*y dx=-k*x*y
% Datos del problema
k=1;
t0=0;
tN=10;
y0=0.01;
x0=1;
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;
x=zeros(1,N+1);
x(1)=x0;
% Euler
for i=1:N
y(i+1)=y(i)+h*(k*y(i)*x(i));
x(i+1)=x(i)+(-h)*(k*y(i)*x(i));
end
% grafico
hold on
plot(t,y,'r')
plot(t,x,'g')
hold off
legend('Concentracion y','Concentracion x','location','best')% dy=k*x*y dx=-k*x*y
% Datos del problema
k=1;
t0=0;
tN=10;
y0=0.01;
x0=1;
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;
x=zeros(1,N+1);
x(1)=x0;
% RK4
for i=1:N
K1_y=k*x(i)*y(i);
K1_x=-k*x(i)*y(i);
K2_y=k*(x(i)+h/2*K1_x)*(y(i)+h/2*K1_y);
K2_x=-k*(x(i)+h/2*K1_x)*(y(i)+h/2*K1_y);
K3_y=k*(x(i)+h/2*K2_x)*(y(i)+h/2*K2_y);
K3_x=-k*(x(i)+h/2*K2_x)*(y(i)+h/2*K2_y);
K4_y=k*(x(i)+h*K3_x)*(y(i)+h*K3_y);
K4_x=-k*(x(i)+h*K3_x)*(y(i)+h*K3_y);
y(i+1)=y(i)+h/6*(K1_y+2*K2_y+2*K3_y+K4_y);
x(i+1)=x(i)+h/6*(K1_x+2*K2_x+2*K3_x+K4_x);
end
% resultados
[y',x',t']
% grafico
hold on
plot(t,y,'r')
plot(t,x,'g')
hold off
legend('Concentracion y','Concentracion x','location','best')3 REACCIÓN Nº2
3.1 Interpretación del problema y deducción de las ecuaciones diferenciales con valor inicial
En esta segunda parte del trabajo, se estudia la ecuación consecutiva de Lotka (1920)
A + X →k1 2X
X + Y →k2 2Y
Y →k3 B
Donde A, X, Y, B son sustancias distintas. Observamos que las dos primeras reacciones son
autocatalíticas. La reacción consume A para producir B mientras que X e Y dominan la
velocidad y mezcla en los estados intermedios. Siguiendo la estrategia de los apartados 1 y 2
deducir las siguientes ecuaciones diferenciales para las concentraciones interpretando
adecuadamente los términos:
X’(t) = k1Ax − k2xy,
y’(t)= k2xy − k3y
B’(t) = k3y
A’ + x’ + y’ + B’ = 0
En los sucesivos apartados se proporcionará un valor inicial x(0), y(0), A(0), B(0) que permitirá
resolver de manera completa las ecuaciones propuestas arriba.
3.2 Resolución numérica del sistema de ecuaciones diferenciales con valor inicial
3.2.1 Método de Euler
% dx=k1*A*x-k2*x*y dy=k2*x*y-k3*y dB=k3*y dA+dx+dy+dB=0 ==> dA=-k1*A*x
% Datos del problema
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; % tamaño de paso
t=t0:h:tN;
N=(tN-t0)/h; % número de intervalos
y=zeros(1,length(t)); % matrices de 1 por N+1
x=zeros(1,length(t));
A=zeros(1,length(t));
B=zeros(1,length(t));
y(1)=y0;
A(1)=A0;
B(1)=B0;
x(1)=x0;
% Euler
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
% grafico
hold on
plot(t,x)
plot(t,y,'r')
plot(t,B,'g')
plot(t,A,'c')
hold off
legend('Concentración X','Concentración de Y','Concentración de B','Concentración de A','Location','best');% dx=k1*A*x-k2*x*y dy=k2*x*y-k3*y dB=k3*y dA+dx+dy+dB=0 ==> dA=-k1*A*x
% Datos del problema
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.001; % tamaño de paso
t=t0:h:tN;
N=(tN-t0)/h; % número de intervalos
y=zeros(1,length(t)); % matrices de 1 por N+1
x=zeros(1,length(t));
A=zeros(1,length(t));
B=zeros(1,length(t));
y(1)=y0;
A(1)=A0;
B(1)=B0;
x(1)=x0;
% Euler
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
% grafico
hold on
plot(t,x)
plot(t,y,'r')
plot(t,B,'g')
plot(t,A,'c')
hold off
legend('Concentración X','Concentración de Y','Concentración de B','Concentración de A','Location','best');3.2.2 Método de Heun
% dx=k1*A*x-k2*x*y dy=k2*x*y-k3*y dB=k3*y dA+dx+dy+dB=0 ==> dA=-k1*A*x
% Datos del problema
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; % tamaño de paso
t=t0:h:tN;
N=(tN-t0)/h; % número de intervalos
y=zeros(1,length(t)); % matrices de 1 por N+1
x=zeros(1,length(t));
A=zeros(1,length(t));
B=zeros(1,length(t));
y(1)=y0;
A(1)=A0;
B(1)=B0;
x(1)=x0;
% Heun
for i=1:N
K1_x=k1.*A(i).*x(i)-k2.*x(i).*y(i);
K1_y=k2.*x(i).*y(i)-k3.*y(i);
K1_B=k3.*y(i);
K1_A=-k1.*x(i).*A(i);
K2_x=k1.*(A(i)+K1_A.*h).*(x(i)+K1_x.*h)-k2.*(x(i)+K1_x.*h).*(y(i)+K1_y.*h);
K2_y=k2.*(x(i)+K1_x.*h).*(y(i)+K1_y.*h)-k3.*(y(i)+K1_y.*h);
K2_B=k3.*(y(i)+K1_y.*h);
K2_A=-k1.*(x(i)+K1_x.*h).*(A(i)+K1_A.*h);
x(i+1)=x(i)+0.5*h.*(K1_x+K2_x);
y(i+1)=y(i)+0.5*h.*(K1_y+K2_y);
B(i+1)=B(i)+0.5*h.*(K1_B+K2_B);
A(i+1)=A(i)+0.5*h.*(K1_A+K2_A);
end
% grafico
hold on
plot(t,x)
plot(t,y,'r')
plot(t,B,'g')
plot(t,A,'c')
hold off
legend('Concentración X','Concentración de Y','Concentración de B','Concentración de A','Location','best');
