MateWiki - Contribuciones del usuario [es]

Reacciones con autocatálisis, Grupo C26

2015-03-06T20:23:00Z

Alejandropcr:

''' 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]]

Reacciones con autocatálisis, Grupo C26

2015-03-06T20:18:32Z

Alejandropcr:

''' 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]]

'''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]]

{{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]]

Archivo:Cod7.png

2015-03-06T20:16:50Z

Alejandropcr:

Reacciones con autocatálisis, Grupo C26

2015-03-06T20:16:31Z

Alejandropcr:

''' 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]]

'''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]]

[[Archivo:cod7.png]]
[[Categoría:Ecuaciones Diferenciales]]
[[Categoría:ED14/15]]
[[Categoría:Trabajos 2014-15]]

Reacciones con autocatálisis, Grupo C26

2015-03-06T20:09:17Z

Alejandropcr:

''' 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]]

'''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]]

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

Archivo:Cod2.png

2015-03-06T20:06:36Z

Alejandropcr:

Reacciones con autocatálisis, Grupo C26

2015-03-06T20:06:15Z

Alejandropcr:

''' 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]]

'''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]]

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

Reacciones con autocatálisis, Grupo C26

2015-03-06T20:05:16Z

Alejandropcr:

''' 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');}}

'''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]]

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

Archivo:Cooood1.png

2015-03-06T20:00:57Z

Alejandropcr:

Reacciones con autocatálisis, Grupo C26

2015-03-06T20:00:08Z

Alejandropcr:

''' 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]]

'''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]]

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

Reacciones con autocatálisis, Grupo C26

2015-03-06T19:58:24Z

Alejandropcr:

''' 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:cod1.png]]

'''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]]

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

Reacciones con autocatálisis, Grupo C26

2015-03-06T19:56:13Z

Alejandropcr:

''' 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:Cod1.png]]

'''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]]

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

Reacciones con autocatálisis, Grupo C26

2015-03-06T19:54:08Z

Alejandropcr:

''' 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:cod1.png]]

'''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]]

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

Archivo:Cod1 (1).png

2015-03-06T19:52:29Z

Alejandropcr:

Reacciones con autocatálisis, Grupo C26

2015-03-06T19:47:16Z

Alejandropcr:

''' 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}}

'''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]]

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

Reacciones con autocatálisis, Grupo C26

2015-03-06T19:45:39Z

Alejandropcr:

''' 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=1. %Método de Euler
2. clear all
3.
4. %Definimos el vector t, N y h
5. t0=0;
6. tN=10;
7. h=0.1;
8. t=t0:h:tN
9. N=(tN-t0)/h;
10.
11. %Creamos los vectores de las variables dependientes
12. x=zeros(1,N+1);
13. y=zeros(1,N+1);
14. x(1)=1;
15. y(1)=0.01;
16. %aplicamos el método de Euler
17. for i=1:N
18. x(i+1)=x(i)+h*(-x(i).*y(i));
19. y(i+1)=y(i)+h*(y(i).*x(i));
20. end
21.
22. %Representamos las soluciones en gráficas
23. hold on
24. plot(t,x,);
25. plot(t,y,'r');
26. grid on
27. legend('Sustancia A: x(t)','Sustancia B: y(t)')
28. xlabel('Tiempo (s)')
29. ylabel('Concentración (mol/l)')
30. hold off}}

'''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]]

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

Reacciones con autocatálisis, Grupo C26

2015-03-06T19:08:32Z

Alejandropcr:

''' 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]]

'''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]]

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

Archivo:Cdg5.png

2015-03-06T19:08:15Z

Alejandropcr:

Reacciones con autocatálisis, Grupo C26

2015-03-06T19:06:25Z

Alejandropcr:

''' 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]]

'''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:cdg4.png]]

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

Archivo:Cdg4.png

2015-03-06T19:06:03Z

Alejandropcr:

Reacciones con autocatálisis, Grupo C26

2015-03-06T19:05:52Z

Alejandropcr:

''' 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]]

'''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.

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

Archivo:Cdg3.png

2015-03-06T19:04:18Z

Alejandropcr:

Reacciones con autocatálisis, Grupo C26

2015-03-06T18:58:26Z

Alejandropcr:

''' 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]]

'''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}}

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

Reacciones con autocatálisis, Grupo C26

2015-03-06T18:57:18Z

Alejandropcr:

Archivo:Cdg2.png

2015-03-06T18:52:47Z

Alejandropcr:

Reacciones con autocatálisis, Grupo C26

2015-03-06T18:52:23Z

Alejandropcr:

Reacciones con autocatálisis, Grupo C26

2015-03-06T18:41:09Z

Alejandropcr:

Reacciones con autocatálisis, Grupo C26

2015-03-06T18:40:22Z

Alejandropcr:

Archivo:Cdg1.png

2015-03-06T18:36:28Z

Alejandropcr: Código usando método de Euler con tamaño de paso h=0.01

Código usando método de Euler con tamaño de paso h=0.01

Reacciones con autocatálisis, Grupo C26

2015-03-06T18:30:02Z

Alejandropcr:

Reacciones con autocatálisis, Grupo C26

2015-03-06T18:29:20Z

Alejandropcr:

Reacciones con autocatálisis, Grupo C26

2015-03-06T18:24:09Z

Alejandropcr:

Reacciones con autocatálisis, Grupo C26

2015-03-06T17:40:31Z

Alejandropcr:

Archivo:RUNGEKUTTA.jpg

2015-03-06T17:38:18Z

Alejandropcr: metodo numérico de Runge kutta

metodo numérico de Runge kutta

Reacciones con autocatálisis, Grupo C26

2015-03-06T17:36:18Z

Alejandropcr:

Archivo:TRAPECIO.jpg

2015-03-06T17:27:42Z

Alejandropcr:

Reacciones con autocatálisis, Grupo C26

2015-03-06T17:26:31Z

Alejandropcr:

Reacciones con autocatálisis, Grupo C26

2015-03-06T17:18:09Z

Alejandropcr:

Reacciones con autocatálisis, Grupo C26

2015-03-06T17:07:31Z

Alejandropcr:

Archivo:EULER.png

2015-03-06T17:06:13Z

Alejandropcr: metodo de Euler

metodo de Euler

Reacciones con autocatálisis, Grupo C26

2015-03-06T16:55:34Z

Alejandropcr:

Reacciones con autocatálisis, Grupo C26

2015-03-06T16:43:51Z

Alejandropcr:

Reacciones con autocatálisis, Grupo C26

2015-03-06T16:42:32Z

Alejandropcr:

Reacciones con autocatálisis, Grupo C26

2015-03-06T16:41:30Z

Alejandropcr:

Reacciones con autocatálisis, Grupo C26

2015-03-06T16:34:48Z

Alejandropcr:

Reacciones con autocatálisis, Grupo C26

2015-03-06T15:28:54Z

Alejandropcr:

Reacciones con autocatálisis, Grupo C26

2015-03-06T15:07:45Z

Alejandropcr:

Reacciones con autocatálisis, Grupo C26

2015-03-06T15:03:11Z

Alejandropcr:

Reacciones con autocatálisis, Grupo C26

2015-03-06T14:49:49Z

Alejandropcr:

Reacciones con autocatálisis, Grupo C26

2015-03-06T14:48:24Z

Alejandropcr: