Reacciones complejas - Grupo 1 A

Trabajo realizado por estudiantes
Título	Reacciones complejas. Grupo 1A
Asignatura	Ecuaciones Diferenciales
Curso	Curso 2014-15
Autores	María Ramírez Ignacio Posada Antonio López-Mateos Pablo Bueno
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 Introducción

El problema que vamos a tratar es un problema que aborda el tema de Reacciones Complejas. Para ello consideramos, según el enunciado dado, una reacción química irreversible en una solución bien mezclada. Tal y como se nos indica inicialmente en los datos del problema, suponemos también que dicha reacción dada ocurre a un volumen y temperatura constantes. La reacción química con la que trabajaremos responde a la estructura que a continuación se presenta:

                                                      A + B → C

Ésta nos indica que tenemos inicialmente dos reactivos A y B que van formando un producto C en lo que se conoce como una reacción bimolecular, es decir, una molécula de A y una de B producen una de C, tal y como expresamos en el esquema de la reacción.

Otro dato también a suponer es que dicha reacción química satisface la Ley de Acción de Masas, ley que establece que la velocidad de reacción es proporcional al producto de las concentraciones de los reactivos.

Una vez conocidos todos estos datos sobre nuestro problema, nos disponemos a ir resolviéndolo respondiendo a las preguntas que se nos han planteado y que a continuación en los siguientes apartados desarrollaremos.

4 Resolución Numérica por el Método de Euler

Ahora, nos disponemos a resolver numéricamente el problema de valor inicial descrito en el apartado anterior, suponiendo a0 = 3 mol/l, b0 = 1 mol/l y k1 = 1 mol/s. Utilizaremos, en este caso, el Método de Euler como método numérico para la resolución del problema, eligiendo un paso h = 0.1 en los primeros 2 segundos. El código Matlab resultante tras la aplicación numérica de este método a nuestro problema de valor inicial es el siguiente:

%PVI por Euler con una h = 0.1
clear all

%condiciones dadas
y0 = 0;
t0 = 0;
tN = 2;
h = 0.1;
k = 1;
%calculamos N
N = (tN - t0)/h;
%definimos la variable t
t = linspace(t0,tN,N+1);
%funcion y
y = zeros(1,N+1);

y(1)= y0;

%solucion obtenida
for i=1:N
    y(i+1)=y(i)+ h*(k*(3 - y(i))*(1 - y(i)));
end

%Concentracion de los reactivos
ya = zeros(1,N+1);
yb = zeros(1,N+1);

ya = 3.- y;
yb = 1.- y; 


%solucion grafica
hold on
plot(t,y,'linewidth',2)
plot(t,ya,'g','linewidth',2)
plot(t,yb,'r','linewidth',2)
legend('Concentracion de C','Concentracion de A','Concentracion de B','Location','best');
hold off

Una vez establecido dicho código Matlab, dibujamos en una misma gráfica tanto las concentraciones de los dos reactivos como la del producto, en función de su evolución en el tiempo, quedando las mismas de la siguiente manera:

Gráfica de las dos concentraciones y el producto. Utilizando el método de Euler y con una h = 0.1

5 Cuando t tiende a infinito

Gráfica de las dos concentraciones y el producto. Utilizando el método de Euler y con una h = 0.1

6 Resolucion por trapecio y runge kutta

6.1 Trapecio

%PVI por Trapecio con una h = 0.1
clear all

%condiciones dadas
y0 = 0;
t0 = 0;
tN = 2;
h = 0.1;
k = 1;
a = 3;
b = 1;
%calculamos N
N = (tN - t0)/h;
%definimos la variable t
t = linspace(t0,tN,N+1);
%funcion y
y = zeros(1,N+1);

y(1)= y0;

%solucion obtenida
for i=1:N
    %Definimos tres variables para reducir la ecuación
    Q = 1 + (h*k*(a+b))/2;
    W = -(h*k)/2;
    E = -y(i)- h*k*(1/2)*((a - y(i))*(b - y(i)) + a*b);
    y(i+1)=(-Q + sqrt(Q^2- (4*W*E)))/(2*W);
end

%Concentracion de los reactivos
ya = zeros(1,N+1);
yb = zeros(1,N+1);

ya = 3.- y;
yb = 1.- y; 


%solucion grafica
hold on
plot(t,y,'linewidth',2)
plot(t,ya,'g','linewidth',2)
plot(t,yb,'r','linewidth',2)
legend('Concentracion de C','Concentracion de A','Concentracion de B','Location','best');
hold off

Gráfica de las dos concentraciones y el producto. Utilizando el método del Trapecio y con una h = 0.1

6.2 Runge Kutta

%PVI por Runge Kutta orden 4 con una h = 0.1
clear all

%condiciones dadas
y0 = 0;
t0 = 0;
tN = 2;
h = 0.1;
k = 1;
a = 3;
b = 1;
%calculamos N
N = (tN - t0)/h;
%definimos la variable t
t = linspace(t0,tN,N+1);
%funcion y
y = zeros(1,N+1);

y(1)= y0;

%solucion obtenida
for i=1:N
    K1 = k*(a - y(i))*(b - y(i));
    K2 = k*(a-(y(i)+K1*h/2))*(b - (y(i)+K1*h/2));
    K3 = k*(a-(y(i)+K2*h/2))*(b - (y(i)+K2*h/2));
    K4 = k*(a-(y(i)+(K3*h)))*(b - (y(i)+K3*h));
    y(i+1)= y(i) + h/6*(K1+(2*K2)+(2*K3)+K4);
end

%Concentracion de los reactivos
ya = zeros(1,N+1);
yb = zeros(1,N+1);

ya = 3.- y;
yb = 1.- y; 


%solucion grafica
hold on
plot(t,y,'linewidth',2)
plot(t,ya,'g','linewidth',2)
plot(t,yb,'r','linewidth',2)
legend('Concentracion de C','Concentracion de A','Concentracion de B','Location','best');
hold off

Gráfica de las dos concentraciones y el producto. Utilizando el método de Runge Kutta de orden 4

7 Apartado 6

%Establecemos los valores iniciales
t0=0;
tN=10;
y0=0;
%Discretizamos
N=1000;
h=(tN-t0)/N;
%Vector de tiempos
t=t0:h:tN;
%Creamos el vector y, donde almacenaremos las soluciones
y1=zeros(1,N+1);
y2=zeros(1,N+1);
%Damos el valor para el primer elemento
y1(1)=y0;
y2(1)=y0;
%Creamos la variable yy para operar
yy1=y0;
yy2=y0;
%Creamos el bucle que nos da soluciones de y
for n=1:N
  yy1=yy1+h*(1*(3-yy1)*(1-yy1)-(1/5)*(yy1-yy2));
  yy2=yy2+h*((1/5)*(yy1-yy2));
  y1(n+1)=yy1;
  y2(n+1)=yy2;
end

%Dibujamos la solucione
a=3.-y1;
b=1.-y1;
c=y1-y2;
d=y2;
plot(t,c)
hold on
plot(t,a,'g')
hold on
plot(t,b,'y')
hold on
plot(t,d,'r')