Ver el código de «Reacciones con autocatálisis, Grupo C26»

''' 1 INTRODUCCION'''

Este trabajo estudia las concentraciones  de reactivos y productos a lo largo del tiempo en una reaccíon química  bimolecular irreversible, con un volumen y una temperatura constantes  y  con la particularidad de que  uno de los componentes (B) hace de efecto catalítico y actuará tanto de reactivo como de producto,es decir, partiendo de un contenido de A y B iniciales, reaccionaran  para producir B que a su vez participara de nuevo como reactivo,  por lo que es de la siguiente forma:

A +B →k2B

Supondremos que en nuestra reacción se satisface la '''Ley de Accion de masas''', que establece que la velocidad de reaccion es proporcional al producto de la concentracion de los reactivos, y la '''ley de de conservacion de masa''' ,por la  que definimos que la masa de reactivos mas productos es constante en el tiempo
Para la resolución del problema se usan los metodos de cálculo  de Euler, del Trapecio , de Runge-Kutta 

''' 2 INTERPRETACION DEL PROBLEMA Y EXPRESION NUMÉRICA'''

comenzamos nombrando las variables y expresando el problema numéricamente:

contenido de A= x(t) [mol/l]

contenido de B= y(t) [mol/l]
 
'''Ley de Acción de Masas''': por definición es: 

'''y'(t)=K*x(t)*y(t)'''  '''(I)'''


y'(t)= velocidad de creación del producto [mol/(l*s)]
k = constante de proporcionalidad mol/s

'''Ley de conservación de masas''' 

masa reactivos + masa productos = cte

'''x(t)+y(t) =(C) (II)'''

derivamos esta expresión y se obtiene:

'''x'(t)+y'(t)=0   x'(t)=-y'(t)  (III)'''

La velocidad de desaparición de los reactivos es igual a la de aparicion de los productos(que es un concepto muy intuitivo)
Por lo tanto la formula (III) se explica a partir derivar respecto al tiempo de La ley de conservación de masas.
ahora sustitumos el contenido de A , x(t) de (II) en (I)  
 
y'(t)=k*y(t)*(C-y(t)) → '''y'(t)=k*C*y(t)-K*(y(t)) (IV)'''
 
además conocemos las concentraciones iniciales(en t=0) y el valor de la k.
A(0)=1mol/l
B(0)=0,01mol/l
k=mol/s

por lo que sustituyendo estos valores  en (II) para deducir la (C)  y en (IV) planteamos el problema de valor inicial
 
 '''{y'(t)= 1,01*y(t)-y²(t)
 '''{ y(0)=0,01'''                      '''(t,y) ∈ [0,10]x R''''''
'''UNICIDAD DE SOLUCIÓN'''

Aplicamos a continuación el Teorema de Existencia y Unicidad:

Sea f(t,y)=y′(t)⇒f(t,y)=1.01y(t)−y(t), y δf(t,y)/δx=1.01−2y(t).

y  B[(t0,y0),r], r>0,en(t0,y0)=(0,0.01)

f(t,y) es continua en B∩D⇒∃ y(t)/y′(t)=1,01y(t)−y²(t)

y su derivada  respecto de ''y'' →f'(t,y) es continua en la intersección con el dominio 

Por lo que  hemos demostrado que sí tiene solución y ésta es única.

'''3 RESOLUCION NUMÉRICA DEL PROBLEMA DE VALOR INICIAL'''

ahora procedemos a resolver el P.V.I dado anteriormente por loe métodos numéricos de Euler, del Trapecio y de Runge-kutta. Estudiamos el problema a lo largo de los 10 primeros segundos.

'''3.1 METODO DE EULER'''
{{matlab|codigo=
%Método de Euler

%Datos iniciales
t0=0;
y0=0.01;
h=0.1;
t0=0;
tN=10;
 
%Cálculo de subintervalos equidistantes
N=(tN-t0)/h;
t=t0:h:tN;
 
%vector solución:
y=zeros(1,N+1);
y(1)=y0;
 
%aplicamos el metodo de 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}}

[[Archivo:EULER.png]]

Con esto vemos que si bien al principio la veocidad de reaccion  (o de creacion de producto)es lenta porque hay poco contenido en B y mucho en A, va progresando hasta su punto máximo cuando A Y B son iguales  y luego vuelve a disminuir porque se va agotando A , y esto ocurre por la Ley de acción de Masas y'(x)=kx(t)y(t) 

'''3.2 METODO DEL TRAPECIO'''
ahora calculamos el P.V.I con el metodo del trapecio con el mismo paso del tiempo que anteriormente.
%Método del Trapecio

{{Matlab|codigo=

%Datos iniciales
t0=0;
y0=0.01;
h=0.1;
t0=0;
tN=10;
 
%calculo de subintervalos equidistantes:
N=(tN-t0)/h;
t=t0:h:tN;
 
%vector solución:
y=zeros(1,N+1);
y(1)=y0;
%metodo del trapecio:
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'''

[[Archivo:TRAPECIO.jpg]]
'''3.3 METODO DE RUNGE KUTTA DE ORDEN CUATRO'''

El método de Runge Kutta es mas complicado por lo que hay que usar funciones ya en Matlab para que sea mas sencillo de expresarlo

{{Matlab|codigo=%Método de Runge-Kutta
 
%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;
 
t=t0:h:tN;
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
 
x=1.01-y;
 
hold on
plot(t,x,'b')
plot(t,y,'r')
legend('A: x(t)','B: y(t)')
xlabel('Time(s)')
ylabel('Concentración (mol/l)')
grid on}}
[[Archivo:RUNGEKUTTA.jpg]]

'''4. RESOLUCION EN FORMA DE SISTEMA DE ECUACIONES'''

'''4.1 METODO DE EULER'''

Aplicamos el método de Euler

{{Matlab|codigo=%Método de Euler
clear all

%Definimos el vector t, N y h
t0=0;
tN=10;
h=0.1;
t=t0:h:tN
N=(tN-t0)/h;

%Creamos los vectores de las variables dependientes
x=zeros(1,N+1);
y=zeros(1,N+1);
x(1)=1;
y(1)=0.01;
%aplicamos el método de Euler
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 las soluciones en gráficas
hold on
plot(t,x,);
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}}

[[Archivo:cooood1.png]]

'''4.2 METODO DE RUNGE-KUTTA'''

El método Runge-Kutta de orden 4 es de orden superior al de Euler, con lo cual nos ofrecerá una mayor aproximación a la solución real. El código del programa es el que se muestra a continuación.

{{Matlab|codigo=%metemos los datos de k, h y creamos el vector t
k=1;
t0=0;
tN=10;
h=0.1;
y0=0.01;
x0=1;
t=t0:h:tN;
%creamos los vectores x e y.
y=zeros(1,length(t));
x=y;
%Introducimos los valores iniciales de las concentraciones
y(1)=y0;
x(1)=x0;
%Definimos las funciones
F=inline('x*y','t','y','x');
G=inline('-x*y','t','y','x');
%Aplicamos el método Runge Kutta de orden 4
for k=1:length(t)-1
    K1y=F(t(k),y(k),x(k));
    K1x=G(t(k),y(k),x(k));
    K2y=F(t(k)+(h/2),y(k)+h*K1y/2,x(k)+h*K1x/2);
    K2x=G(t(k)+(h/2),y(k)+h*K1y/2,x(k)+h*K1x/2);
    K3y=F(t(k)+(h/2),y(k)+h*K2y/2,x(k)+h*K2x/2);
    K3x=G(t(k)+(h/2),y(k)+h*K2y/2,x(k)+h*K2x/2);
    K4y=F(t(k)+h,y(k)+K3y*h,x(k)+K3x*h);
    K4x=G(t(k)+h,y(k)+K3y*h,x(k)+K3x*h);
    y(k+1)=y(k)+(h/6)*(K1y+2*K2y+2*K3y+K4y);
    x(k+1)=x(k)+(h/6)*(K1x+2*K2x+2*K3x+K4x);
end
%Dibujamos ambas gráficas.
hold on
plot(t,y);
plot(t,x,'r');
hold off
legend('Concentración de B','Concentración de A','Location','Best');}}


[[Archivo:cod2.png]]

'''4.3 CONCLUSION'''

Como podemos observar se obtienen los mismos resultados resolviéndolo mediante una EDO o por un sistema. La velocidad de aparición de B es igual a la velocidad en la que desaparece A. A acaba desapareciendo casi por completo y B se acaba acercando a 1mol.
Alrededor de los 4.5 segundos las concentraciones de ambos son iguales.


'''5. REACCIÓN DE LOTKA. INTERPRETACIÓN'''

Consideramos a partir de ahora la reacción consecutiva propuesta por Lotka en 1920:
 '''A + X --> {k1} 2X
 '''X + Y --> {k2} 2Y'''                      
 '''Y --> {k3} B'''''
Donde A, X, Y, D son sustancias distintas e interpretaremos como la concentración de dichas sustancias.  La reacción consume A para producir B y la velocidad y mezcla de las reacciones están controladas por X e Y.
Según la ley de la conservación de la masa: "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 por los productos". En nuestro caso, al no tratarse de una reacción nuclear, podemos aplicar esta ley de la siguiente manera: la suma de las concentraciones de las sustancias implicadas debe de permanencer constante.
            
 '''A + X + Y + B = cte'''''

Sin embargo, debemos recurrir también a la ley de acción de masas: "En 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". En este punto, es crítico ver que las sustancias X e Y no solo se consumen, sino que las primeras dos reaccionas también las producen. Por ello mismo, restaremos la porción producida a la porción inicial o consumida.

A partir de esto podemos calcular el sistema de ecuaciones que corresponde a nuestro problema:

 '''X' = k1*A*X - k2*X*Y
 '''Y' = k2*X*Y - k3*Y'''                      
 '''B' = k3*Y'''''

Si derivamos la expresión obtenida a partir de la ley de conservación de la masa, obtenemos:
               
 '''A' + X' + Y' + B' = 0'''''

Expresión en función de X', Y' y B', las cuales ya tenemos, y podemos por tanto sustituir en nuestro problema:
 '''A' + k1*A*X - k2*X*Y + k2*X*Y - k3*Y + k3*Y = 0
 '''A' + k1*A*X = 0'''                      
 '''A' = -k1*A*X'''''

Donde k1 es la velocidad de reacción.

Por ello, en los próximos apartados deberemos resolver el siguiente problema de valor inicial con los datos que se proporcionan, que son:

 ''k1 = k2 = 2*k3 = 0.1
 '''Y' = 0.1*X*Y - 0.1/2*Y'''
 '''B' = 0.1/2*Y'''                     
 '''A' = -0.1*A*X'''''

con valores iniciales: 

A(0) = 5; X(0)=0.0005; Y(0)=0.00001; B(0)=0

'''6. PVI ASOCIADO AL PROBLEMA ANTERIOR'''


Para los valores de k1, k2, k3, y de las concentraciones dados, planteo el PVI usando el método de Euler para los tamaños de paso h=0.01 y h=0.001

{{Matlab|codigo=
t0=0;                                                                                                              
tf=200;
h=0.01; 
N=(tf-t0)/h;
t=t0:h:tf;
a0=5; 
x0=5*(10^-4); 
y0=10^(-5); 
b0=0; 
A=zeros(1,N+1); 
X=zeros(1,N+1);
Y=zeros(1,N+1);
B=zeros(1,N+1); 
A(1)=a0; 
X(1)=x0; 
Y(1)=y0; 
B(1)=b0; 
k1=0.1; 
k2=0.1; 
k3=0.1/2;
 
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*A(i)*X(i));
end
 
hold on
plot(t,A,'-b'); 
plot(t,X,'-r');
plot(t,Y,'-y');
plot(t,B,'-g');
hold off
xlabel('Tiempo (s)','FontSize', 11);
ylabel('Concentraciones (mol/L)','FontSize', 11);
title('Tamaño de paso 0.01','FontName','Berlin sans FB','FontSize', 20);
legend('Concentracion A','Concentracion X', 'Concentracion Y', 'Concentracion B', 'location','east')}}

[[Archivo:cdg1.png]] 

{{Matlab|codigo=
t0=0; 
tf=200;
h=0.001; 
N=(tf-t0)/h;
t=t0:h:tf;
a0=5; 
x0=5*(10^-4); 
y0=10^(-5); 
b0=0; 
A=zeros(1,N+1); 
X=zeros(1,N+1);
Y=zeros(1,N+1);
B=zeros(1,N+1); 
A(1)=a0; 
X(1)=x0; 
Y(1)=y0; 
B(1)=b0; 
k1=0.1; 
k2=0.1; 
k3=0.1/2;
 
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*A(i)*X(i));
end
 
hold on
plot(t,A,'-b'); 
plot(t,X,'-r');
plot(t,Y,'-y');
plot(t,B,'-g');
hold off
xlabel('Tiempo (s)','FontSize', 11);
ylabel('Concentraciones (mol/L)','FontSize', 11);
title('Tamaño de paso 0.001','FontName','Berlin sans FB','FontSize', 20);
legend('Concentracion A','Concentracion X', 'Concentracion Y', 'Concentracion B', 'location','east')}}

[[Archivo:cdg2.png]] 

'''¿ES ESTABLE?'''

Para comprobar si es o no estable, cambiaremos ligeramente las concentraciones iniciales, realizaremos una nueva gráfica (todo ello con un mismo tamaño de paso) y comprobaremos si la diferencia entre ambas es casi inapreciable. Si es así, es estable.

{{Matlab|codigo=
t0=0; 
tf=200;
h=0.01; 
N=(tf-t0)/h;
t=t0:h:tf;
a0=5.001; 
x0=5.001*(10^-4); 
y0=10.001^(-5); 
b0=0.001; 
A=zeros(1,N+1); 
X=zeros(1,N+1);
Y=zeros(1,N+1);
B=zeros(1,N+1); 
A(1)=a0; 
X(1)=x0; 
Y(1)=y0; 
B(1)=b0; 
k1=0.1; 
k2=0.1; 
k3=0.1/2;
 
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*A(i)*X(i));
end

a02=5; 
x02=5*(10^-4); 
y02=10^(-5); 
b02=0; 
A2=zeros(1,N+1); 
X2=zeros(1,N+1);
Y2=zeros(1,N+1);
B2=zeros(1,N+1); 
A2(1)=a0; 
X2(1)=x0; 
Y2(1)=y0; 
B2(1)=b0; 
 
for i=1:N
X2(i+1)=X2(i)+h*(k1*A2(i)*X2(i)-k2*X2(i)*Y2(i));
Y2(i+1)=Y2(i)+h*(k2*X2(i)*Y2(i)-k3*Y2(i));
B2(i+1)=B2(i)+h*(k3*Y2(i));
A2(i+1)=A2(i)+h*(-k1*A2(i)*X2(i));
end
hold on
plot(t,A,'-b'); 
plot(t,X,'-b');
plot(t,Y,'-b');
plot(t,B,'-b');
plot(t,A2,'-r'); 
plot(t,X2,'-r');
plot(t,Y2,'-r');
plot(t,B2,'-r');
hold off}}

[[Archivo:cdg3.png]]

Se puede ver que la diferencia es apenas apreciable, hay que realizar un zoom importante para ver con claridad la separación de las gráficas.

[[Archivo:cdg5.png]]

'''7. METODO DE HEUN'''

El método de Heun es un método explícito. Primero definimos unas constantes y luego aplicamos el método a nuestro sistema

{{Matlab|codigo=%Introducimos los datos
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;
%Creamos los vectores tiempo y los vectores vacios
t=t0:h:tN;
N=(tN-t0)/h;
A=zeros(1,N+1);
X=zeros(1,N+1);
Y=zeros(1,N+1);
B=zeros(1,N+1);
 
%Introducimos los valores iniciales
A(1)=A0;
X(1)=X0;
Y(1)=Y0;
B(1)=B0;
 
for i=1:N
%Definimos las constantes
K1A=-k1.*X(i).*A(i);
K1X=k1.*A(i).*X(i)-k2.*X(i).*Y(i);
K1Y=k2.*X(i).*Y(i)-k3.*Y(i);
K1B=k3.*Y(i);
K2A=-k1.*(X(i)+K1A.*h).*(A(i)+K1A.*h);
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);
%Método de heun
A(i+1)=A(i)+0.5*h.*(K1A+K2A);
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);
end

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

[[Archivo:cod7.png]]

[[Categoría:Ecuaciones Diferenciales]]
[[Categoría:ED14/15]]
[[Categoría:Trabajos 2014-15]]