Diferencia entre revisiones de «Reacciones complejas (Grupo 3-C)»

Revisión actual del 17:17 9 mar 2015

Trabajo realizado por:

• María Victoria Manget Sánchez Sacristán, 1021.
• Álvaro Ramírez Fernández de la Puente, 907.
• Ventura Rubí Zaldívar, 1218.
• Iván Díez Berjano 1111.
• Marta Ruiz Martínez, 1219.
• Alicia Vilalta Duce, 876.

1 Enunciado

Se considera una reacción química irreversible en una solución bien mezclada. Supondremos que la reacción ocurre para un volumen y temperatura constantes. Al inicio se encuentran dos reactivos A y B que van formando un producto C en lo que se conoce como una reacción bimolecular, es decir, una molécula de A y una de B producen una de C,

                                                  A + B → C

Supondremos también que se satisface la ley de acción de masas que establece que la velocidad de reacción es proporcional al producto de las concentraciones de los reactivos. Se pide:

1. Comprobar que la concentración del producto C a lo largo del tiempo puede obtenerse a partir de la ecuación diferencial

                                 [math]y'(t)=k_{1}[/math]([math]a_{0}[/math] −[math]y(t)[/math])([math]b_{0}[/math] −[math]y(t)[/math]), [math]t \gt 0[/math],

para algunas constantes [math]a_{0}[/math], [math]b_{0}[/math] y [math]k_{1}[/math]. Interpretar dichas constantes y enunciar un problema de valor inicial adecuado para calular la concentración de C a lo largo del tiempo. ¿Tiene solución única?

2. ¿Cómo cambiaría la ecuación diferencial si el proceso fuera reversible?

3. Suponiendo [math]a_{0}[/math] = 3 mol/l, [math]b_{0}[/math] = 1 mol/l y [math]k_{1}[/math] = 1 mol/s resolver el PVI por el método de Euler, eligiendo un paso h = 0.1, en los primeros 2 segundos. Dibujar en una misma gráfica las concentraciones de los dos reactivos y el producto. Interpretar las gráficas.

4. Resolver numéricamente el sistema para un tiempo suficientemente grande como para establecer una aproximación de los límites de las concentraciones cuando t → ∞. ¿Es lógico lo que se observa?

5. Resolver el PVI numéricamente usando el método del Trapecio y el método de Runge-Kutta de cuarto orden con el mismo paso de tiempo.

6. Supongamos ahora que se produce una reacción consecutiva de la forma

                                    A + B →[math]k_{1}[/math] C →[math]k_{2}[/math] D 

Plantear ahora un sistema de dos ecuaciones diferenciales para las concentraciones de C y D, basado en la ley de acción de masas. Obtener el PVI asociado.

7. Si [math]k_{2}[/math] = 5, lo que supone que la segunda reacción es mucho más rápida que la primera, resolver numéricamente el problema en los primeros 10 segundos con los métodos de Euler y Runge-Kutta. Dibujar las gráficas de las concentraciones en una misma figura e interpretar los resultados. ¿Se puede usar el paso h = 0.3?

8. Si [math]k_{2}[/math] = 1/5, lo que supone que la segunda reacción es mucho más lenta que la primera, resolver numéricamente el problema en los primeros 10 segundos con los métodos de Euler y Runge-Kutta. Dibujar las gráficas de las concentraciones en una misma figura e interpretar los resultados.

2 Planteamiento de nuestro problema y PVI asociado

En química, una reacción compleja (o reacción compuesta) es aquella que se produce, a nivel molecular, a través de varias etapas o reacciones elementales. Una reacción compleja se describe y explica a través de un mecanismo de reacción (la secuencia de etapas elementales por la que los reactivos pasan a productos). Por ejemplo, una reacción con al menos un intermedio y como mínimo dos etapas (o reacciones) elementales es una reacción compleja.

Dos sustancias A y B reaccionan formando una sustancia C. Interesa conocer la cantidad de sustancia C formada en el instante t. Para resolver el problema anterior nos apoyamos en la siguiente ley de la química: «La rapidez de una reacción química en la que se forma un compuesto C a partir de dos sustancias A y B es proporcional al producto de las cantidades de A y de B que no han reaccionado.» Denotamos por:

• [math]y(t)[/math]=‘Cantidad de C en el instante t’.
• [math]a_{0}[/math]=‘Cantidad inicial de A’.
• [math]αy(t)[/math]=‘Cantidad de A usada en el instante t’.
• [math]b_{0}[/math]=‘Cantidad inicial de B’.
• [math]βy(t)[/math]=‘Cantidad de B usada en el instante t’.
• [math]k_{1}[/math]=‘Constante de proporcionalidad’.

Siendo [math]α, β=1[/math] ya que la reacción es bimolecular. Lo que quiere decir que se necesitan una molécula de A y una molécula de B para formar una molécula de C. En el caso de que nos diesen otra relación entre las variables, para formar [math]y(t)[/math] gramos del compuesto C se necesitarían M partes de la sustancia A y N partes de la sustancia B, esto es [math]\frac{M}{M+N}y(t)[/math] de A y [math]\frac{N}{M+N}y(t)[/math] de B; donde [math]\frac{M}{M+N}[/math] lo hemos definido como [math]α[/math] y [math]\frac{N}{M+N}[/math] como [math]β[/math].

‘‘La Ley de Conservación de la Masa de Lavoisier garantiza que la cantidad de sustancia C creada en el instante [math]t[/math] coincide con la suma de las cantidades usadas de A y B’’. Formulando la ley anterior para [math]y(t)[/math] encontramos que esta función verifica la ecuación diferencial ordinaria de variables separadas [math]y'(t)=k_{1}[/math]([math]a_{0}[/math] −[math]y(t)[/math])([math]b_{0}[/math] −[math]y(t)[/math]), [math]t \gt 0[/math], con condición inicial [math]y(0)=0[/math] (pues en el instante [math]t=0[/math] no se ha formado aún sustancia C).

En nuestro caso el problema de Cauchy o valor inicial sería::

[math]\left. y'(t)=k_{1}(a_{0}−y(t))(b_{0}−y(t)), t \gt 0 \atop y(0)=0 \right\}[/math]

Vemos ahora si nuestro problema tiene o no solución única y si ésta existe a partir del Teorema de Existencia y Unicidad:

La función [math]f (t,y)= k_{1}[/math]([math]a_{0}[/math]−[math]y(t)[/math])([math]b_{0}[/math]−[math]y(t)[/math]) es continua en el intervalo [math] I=(0,∞)∩B((0,0),r\gt0)[/math], por lo que ya aseguramos la existencia de al menos una solución. Por otro lado, vemos como [math]\textstyle\frac{df}{dy}(t,y)[/math]=2[math]k_{1}[/math][math]y[/math]−[math]k_{1}[/math][math]a_{0}[/math] −[math]k_{1}[/math][math]b_{0}[/math], también es continua en todo [math]{ℝ}^2[/math] por ser polinómica y lo es entonces en ([math]t_{0}[/math],[math]y_{0}[/math]), es decir, en [math](0,0)[/math], y más concretamente en cualquier bola con [math]s\gt0[/math].

Por lo que diremos que el problema de valor inicial posee una única solución.

3 Reacción reversible

Se llama reacción reversible a la reacción química en la cual los productos de la reacción vuelven a combinarse para generar los reactivos. Este tipo de reacción se representa con una doble flecha, donde la flecha indica el sentido de la reacción. Esta ecuación representa una reacción directa (hacia la derecha) que ocurre simultáneamente con una reacción inversa (hacia la izquierda):

El sistema muestra que A y B desaparecen por la reacción directa, pero se crean en la reacción inversa. Por lo tanto

[math]\textstyle\frac{d[A]}{dt}+\textstyle\frac{d[B]}{dt}[/math]=(−[math]k_{1}[A]−k_{1}[B])+k_{-1}[/math][C]

[math]\textstyle\frac{d[A]}{dt}+\textstyle\frac{d[B]}{dt}[/math]=-[math]\textstyle\frac{dy}{dt}[/math]

[math]\textstyle\frac{dy}{dt}[/math]=([math]k_{1}[A]+k_{1}[B])-k_{-1}[/math][C]

[math][A][/math]=[math]a_{0}[/math]−[math]y[/math]; [math][B][/math]=[math]b_{0}[/math]−[math]y[/math]; [math][C][/math]=[math]c_{0}[/math]+[math]y[/math]; donde [math]y[/math] se define como el grado de avance en que se forma C.

[math]y'(t)=k_{1}[/math]([math]a_{0}[/math] −[math]y(t)[/math])([math]b_{0}[/math] −[math]y(t)[/math])−[math]k_{-1}[/math]([math]c_{0}[/math] +[math]y(t)[/math])

Como la condición inicial que hemos definido en el instante [math]y(0)=0[/math], ya que consideramos que no hay cantidad de C al principio de la reacción, [math]c_{0}=0[/math]. Por lo tanto nuestra ecuación diferencial quedaría

[math]y'(t)=k_{1}[/math]([math]a_{0}[/math] −[math]y(t)[/math])([math]b_{0}[/math] −[math]y(t)[/math])−[math]k_{-1}[/math]([math]y(t)[/math])

Planteamos entonces nuestro problema de valor inicial::

[math]\left. y'(t)=k_{1}(a_{0}−y(t))(b_{0}−y(t))−k_{-1}(y(t)), t \gt 0 \atop y(0)=0 \right\}[/math]

Siendo [math]k_{-1}[/math] la constante de proporcionalidad en el proceso inverso. Para que la reacción sea reversible, [math]k_{-1}[/math] tiene que ser mayor que 0 pero al mismo tiempo inferior que 1, ya que si no se estaría descomponiendo más compuesto C del que realmente se habría formado y eso no podría ser.

Procedemos a obtener una solución numérica mediante el método de Euler dándole un valor a [math]k_{-1}=0.1[/math]:

Método de Euler

Gráfica comparativa

Gráfica del error

%PVI por el método de Euler
clear all; clc; clf;

%Datos del problema 
t0=0;; 
tN=input('introducir el valor de t final:');%El problema se resuelve para tN=2 pero usaremos tN=6 
                                            %para observar cómo se estabilizan las concentraciones
y0=0; k1=1;
k2=0.1; %Sería lo que hemos definido teóricamente como k(-1)
a0=3; b0=1;

%Discretización temporal
h=0.1; %Paso del problema
t=t0:h:tN;
N=(tN-t0)/h; %N intervalos, N+1 valores

%Preparación vector solución
y=zeros(1,length(t)); %También válido y=zeros(1,length(t))
y(1)=y0; %Inicialización
for i=1:length(t)-1
    y(i+1)=y(i)+h*((3-y(i))*(1-y(i))-0.1*y(i));
end

%Mostrar tabla de resultados
[t',y']

%Definimos
ya=zeros(1,N+1); 
yb=zeros(1,N+1); 














%Solución Exacta
ye=(30.*(41+sqrt(481)).*(exp((sqrt(481)/10).*t)-1))./(((1081+(41*sqrt(481))).*(exp((sqrt(481)/10).*t)))-600)

%Concentración de los reactivos
ya=3-ye; %Resultado obtenido de hacer: ya=3-y(t); y(t)=ye
yb=1-ye; %Resultado obtenido de hacer: yb=1-y(t); y(t)=ye

%Solución de la gráfica
figure(1)
hold on
plot(t,y)
plot(t,ya,'r')
plot(t,yb,'g')
legend('Concentración de C','Concentración de A','Concentración de B','location','best');
title('Solución mediante el método de Euler');
hold off

%Error cometido
ye=zeros(1,length(t));
for i=1:N+1
    ye(i)=(30*(41+sqrt(481))*(exp((sqrt(481)/10)*t(i))-1))/(((1081+(41*sqrt(481)))*(exp((sqrt(481)/10)*t(i))))-600)
end

error_euler=abs(y-ye) %Devuelve el error cometido en cada punto
error=max(error_euler) %Devuelve el error máximo cometido que es: 0.0652

%Comparación de la solución real y el error cometido por Euler
figure(3)
hold on
plot(t,y)
plot(t,ye,'r')
legend('Euler','Real','location','best');
title('Gráfica comparativa');
hold off

%Gráfica del error
hold on
figure(4)
plot(t,error_euler)
legend('Error cometido','location','best');
title('Gráfica del error');
hold off

%Cálculo de los límites de las funciones de las concentraciones
limye=0.95 %Para el cálculo del límite de nuestra función y(t) nos apoyaremos en otro programa
ya=a0-limye; disp(ya) %Nos dará que la concentración A se estabiliza cuando alcanza el valor 2.05
yb=b0-limye; disp(yb) %Nos dará que la concentración B se estabiliza cuando alcanza el valor 0.05

Esta reacción alcanza eventualmente una situación de equilibrio en la que las velocidades de las reacciones directa e inversa se igualan y, por lo tanto, las concentraciones de ambos componentes permanecen estables. El cociente [math]k_{1}[/math]/[math]k_{-1}[/math] proporciona entonces el valor de la constante de equilibrio [math]K_{c}[/math]=[math][C]_{eq}[/math]/[math][A]_{eq}[B]_{eq}[/math].

4 Resolución numérica del PVI por el método de Euler

Lo resolveremos también para h=0.01 para ver cómo influye el valor del paso en la aproximación utilizando el método de Euler.

Método de Euler con h=0.1

Solución real con h=0.1

Gráfica comparativa con h=0.1

Gráfica del error con h=0.1

%PVI por el método de Euler
clear all; clc; clf;

%Datos del problema 
t0=0; tN=2; y0=0; k1=1; a0=3; b0=1;

%Discretización temporal
h=input('Introducimos valor del paso:'); %Paso del problema h=0.1, h=0.01
t=t0:h:tN;
N=(tN-t0)/h; %N intervalos, N+1 valores

%Preparación vector solución
y=zeros(1,N+1); %También valido y=zeros(1,length(t))
y(1)=y0; %Inicialización
for i=1:N       
  y(i+1) =y(i)+h*(k1*((a0-y(i))*(b0-y(i))));
end

%Mostrar tabla de resultados
[t',y']

%Definimos 
ya=zeros(1,N+1);
yb=zeros(1,N+1);

%Concentración de los reactivos
ya=(-6*exp(2*t))./(1-(3*exp(2*t))); %Resultado obtenido de hacer: ya=3-y(t); y(t)=ye
yb=(-2./(1-3*exp(2*t))); %Resultado obtenido de hacer: yb=1-y(t); y(t)=ye

%Solución de la gráfica
figure(1)
hold on
plot(t,y,'linewidth',2)
plot(t,ya,'r','linewidth',2)
plot(t,yb,'g','linewidth',2)
legend('Concentración de C','Concentración de A','Concentración de B','location','best');
title('Solución mediante el método de Euler');
hold off

%Comparación con la solución real
ye=(3-3*exp(2.*t))./(1-3*exp(2.*t));
figure(2)
hold on
plot(t,ye,'linewidth',2)
plot(t,ya,'r','linewidth',2)
plot(t,yb,'g','linewidth',2)
legend('Concentración de C','Concentración de A','Concentración de B','location','best');
title('Solución real');
hold off

%Error cometido
ye=zeros(1,length(t));
for i=1:N+1
  ye(i)=(3 - 3*exp(2.*t(i)))./(1 - 3*exp(2.*t(i)));
end

error_euler=abs(y-ye) %Devuelve el error cometido en cada punto
error=max(error_euler) %Devuelve el error maximo cometido que es: 0.065104

%Comparación de la solución real y el error cometido por Euler 
figure(3)
hold on
plot(t,y,'o-')
plot(t,ye,'ro-')
plot(t,0*t,'k*-')
legend('Euler','Real','Aprox.puntos','location','best');
title('Gráfica comparativa');
hold off

%Gráfica del error
figure(4)
hold on
plot(t,error_euler,'o-')
plot(t,0*t,'k*-')
legend('Error cometido','Aprox.puntos','location','best');
title('Gráfica del error');
hold off

Como podemos ver en la imagen la concentración de los reactivos y productos se vuelve asintótica tendiendo la de A a 2 mol/l, la de B a 0 mol/l y la de C a 1 mol/l. Lo que por la estequiometría de la reacción es lógico, ya que el producto C dejará de formarse cuando se agote uno de los dos reactivos que lo forman; en nuestro caso disminuye progresivamente la concentración de B, al que se denomina reactivo limitante, y a la vez se estabilizan las concentraciones.

Para un valor de paso h=0.01 quedaría:

Método de Euler con h=0.01

Solución real con h=0.01

Gráfica comparativa con h=0.01

Gráfica del error con h=0.01

Podemos comprobar que a medida que el paso, h, disminuye, el método de Euler se hace cada vez más preciso.

5 Resolución numérica para t → ∞

Vamos a calcular ahora cuánta cantidad de C se ha formado trancurrido un tiempo largo (t=∞):

Lo que quiere decir que transcurrido dicho tiempo, la concentración de B se agotará mientras que la de A será 2. Utilizaremos a continuación Matlab para realizar dichos cálculos mediante el uso de un programa y obtendremos dos gráficas en las que veremos los resultados con mayor claridad.

Cálculo del valor del tiempo

Aproximación del PVI para t->inf

%Primero calculamos el valor más aproximado del tiempo en el que las funciones de nuestro problema
%de valor inicial darán un valor lo más aproximado posible al que resulta de calcular el límite:

%NO SABEMOS EL TIEMPO FINAL. HAY QUE CALCULARLO
clear all;clf;clc;

%Datos del problema 
t0=0; y0=0; k1=1; h=0.1; a0=3; b0=1;

t(1)=t0;
y(1)=y0;
i=1;

%El límite de y(t) es 1, el objetivo es alcanzar un punto que diste de 1/1000000 (una parte por 
%millón) del límite, luego si estamos a una distancia de 1 menor que 1*(1/1000000) hemos logrado
%el objetivo.

lim=1;

while abs(lim-y(i))>1/1000000*lim
    y(i+1) =y(i)+h*(k1*((a0-y(i))*(b0-y(i)))); 
    t(i+1)=t(i)+h;
    i=i+1;
end

%Tabla de Resultados
[t',y']

disp('Tiempo final:')
disp(t(end)) %Obtenemos un t=6

hold on
figure(1)
plot(t,y)
title('Cálculo del valor del tiempo')

%Ahora sabemos que para cualquier t>6 (en realidad el mínimo valor que tomaría t es: 5.90030501871171741990900999periodo) nuestra 
%aproximación al valor que resulta de calcular el límite es muy precisa; entonces procedemos a resolver el PVI:

%PVI por el método de Euler
%Datos del problema 
t0=0; y0=0; k1=1; a0=3; b0=1;
tN=t(end); %Calculado anteriormente, pues ahora el intervalo es [0,6]

%Discretización temporal
h=0.01; %Paso del problema
t=t0:h:tN;
N=round((tN-t0)/h); %N intervalos, N+1 valores

%Preparación vector solución
y=zeros(1,length(t)); %También válido y=zeros(1,length(t))
y(1)=y0; %Inicialización
for i=1:(length(t)-1)      
  y(i+1) =y(i)+h*(k1*((a0-y(i))*(b0-y(i))));
end

%Mostrar tabla de resultados
[t',y']

%Definimos
ya=zeros(1,N+1);
yb=zeros(1,N+1);

%Concentración de los reactivos
ya=(-6*exp(2*t))./(1-(3*exp(2*t))); %Resultado obtenido de hacer: ya=3-y(t); y(t)=ye
yb=(-2./(1-3*exp(2*t))); %Resultado obtenido de hacer: yb=1-y(t); y(t)=ye

%Solución de la gráfica
figure(2)
hold on
plot(t,y)
plot(t,ya,'r')
plot(t,yb,'g')
legend('Concentración de C','Concentración de A','Concentración de B','location','best');
title('Aproximación del PVI para t->inf');
hold off

Como ya hemos mencionado en los apartados anteriores, es lógico lo que se observa.

6 Métodos del Trapecio y de Runge-Kutta

Al igual que el método de Euler, el método de Runge-Kutta se compone de un sistema explícito, por lo que no hay que hacer ninguna manipulación al esquema numérico. En el caso del método del Trapecio, el esquema numérico mediante el cual calculamos es implícito, por lo que habremos de realizar los siguientes cálculos:

Introduciendo el siguiente código Matlab implementamos el método:

Gráfica de los tres métodos de aproximación

Gráfica de los errores

%PVI por el método de Euler
clear all; clc; clf;

%Datos del problema 
t0=0; tN=2; y0=0; k1=1; a0=3; b0=1;

%Discretización temporal
h=0.1; %Paso del problema
t=t0:h:tN;
N=(tN-t0)/h; %N intervalos, N+1 valores

%Preparación vector solución para trapecio
y=zeros(1,N+1); %También válido y=zeros(1,length(t))
y(1)=y0; %Inicialización

%Preparación vector solución para RK4
z=zeros(1,N+1);
z(1)=y0;

%Preparación vector solución para Euler
r=zeros(1,N+1);
r(1)=y0;

for i=1:N
  %Trapecio
  y(i+1)=(-h)\((-2*h-1)+sqrt((((2*h)+1)^2)-2*h*(y(i)+(h/2)*(6-(4*y(i))+(y(i)^2)))));
  %RK4
  K1=(3-z(i))*(1-z(i));
  K2=(3-(z(i)+1/2*K1*h))*(1-(z(i)+1/2*K1*h));
  K3=(3-(z(i)+1/2*K2*h))*(1-(z(i)+1/2*K2*h));
  K4=(3-(z(i)+K3*h))*(1-(z(i)+K3*h));
  z(i+1)=z(i)+h/6*(K1+2*K2+2*K3+K4);
  %Extra Euler
  r(i+1) =r(i)+h*(k1*((a0-r(i))*(b0-r(i))));
end

%Mostrar tabla de resultados
[t',r',y',z']

%Solución Real
ye=(3-3*exp(2.*t))./(1-3*exp(2.*t));

%Gráfica del ejercicio
figure(1)
hold on
plot(t,y,'r')
plot(t,z,'g')
plot(t,r)
plot(t,ye,'k')
legend('trapecio','Runge Kutta','Euler','Solucion Real','Location','best')
title('Gráfica de los tres métodos de aproximación')
hold off

%Error cometido
ye=zeros(1,length(t));
for i=1:N+1
  ye(i)=(3 - 3*exp(2.*t(i)))./(1 - 3*exp(2.*t(i)));
end

error_euler= abs(r-ye);
error_trapecio= abs(y-ye);
error_RK4= abs(z-ye);

%Tabla de Errores cometidos en cada punto por cada método
[t',error_euler',error_trapecio',error_RK4']

error1=max(error_euler) %Devuelve el error máximo cometido por Euler que es: 0.065104
error2=max(error_trapecio)%Devuelve el error máximo cometido por Trapecio que es: 0.0042206
error3=max(error_RK4)%Devuelve el error máximo cometido por RK4 que es: 2.4027e-005

%Gráfica del error
figure(2)
hold on
plot(t,error_euler,'o-')
plot(t,error_trapecio,'ro-')
plot(t,error_RK4,'go-')
plot(t,0*t,'k*-')
legend('Error Euler','Error Trapecio','Error Runge Kutta-4','Aprox.puntos','location','best');
title('Gráfica de los errores');
hold off

Añadimos además la resolución de nuestra ecuación por el método de Euler para comparar los tres métodos en una sola gráfica. Podemos observar que mientras el error cometido por el método de Euler se dispara, el error cometido por el método del Trapecio disminuye considerablemente. En cambio, el método de Runge-Kutta aproxima con mucha exactitud, cometiendo un error muy pequeño (nulo en comparación con los otros dos).

7 Reacción consecutiva

Una reacción consecutiva o sucesiva es aquella donde el producto de una primera reacción es el reactivo de la siguiente. Nos regiremos de nuevo por la ley de acción de masas para determinar el sistema de ecuaciones diferenciales que modeliza el comportamiento de esta reacción química. La velocidad de reacción de este nuevo proceso químico consecutivo al anterior también será proporcional a la concentración de su reactivo, en este caso de C, y la velocidad de reacción del proceso anterior se verá reducida por la velocidad de formación del nuevo compuesto D. Así podríamos modelizar el comportamiento de la reacción anteriormente descrita con el siguiente sistema de ecuaciones diferenciales::

[math] \left . \begin{matrix} y'_1(t)={k_{1}}({a_{0}}-y_1(t))({b_{0}}-y_1(t))\\ y'_2(t)={k_{2}}(y_1(t)-y_2(t))\\ y_1(0)=0\\y_2(0)=0 \end{matrix} \right \} [/math]

Donde:

• [math]y'_{1}(t)[/math]=‘Velocidad con la que se forma C a través del tiempo’.
• [math]y'_{2}(t)[/math]=‘Velocidad con la que se forma D a través del tiempo’.
• [math]k_1[/math]=‘Constante de proporcionalidad de la primera reacción’.
• [math]k_2[/math]=‘Constante de proporcionalidad de la segunda reacción’.
• [math]y_{2}(t)[/math]=‘Concentración del producto D’.
• [math]y_{1}(t)[/math]-[math]y_{2}(t)[/math]=‘Concentración de C’.

Podemos observar que se trata de un sistema de ecuaciones diferenciales no lineales de primer orden. En los siguientes apartados lo resolveremos numéricamente de forma vectorial.

7.1 Resolución numérica por Trapecio y RK4 cuando k2=5

Solución mediante el método de Euler para h=0.1

Solución mediante el método de Runge-Kutta para h=0.1

Gráfico comparativo para compuesto C y D para h=0.1

%Sistema de ecuaciones no lineal por Trapecio y RK4
clear all; clc; clf;

%Datos del problema 
t0=0; tN=10; y0=[0;0]; k1=1; k2=5; a0=3; b0=1;

%Discretización temporal
h=input('Introducimos valor del paso:'); %Paso del problema h=0.1, h=0.3
t=t0:h:tN;
N=round((tN-t0)/h); %N intervalos, N+1 valores

%Método de Euler 
%Preparación vector solución para un tiempo t 
y=zeros(2,length(t));
y(:,1)=y0; %Inicialización
for i=1:length(t)-1
  y(:,i+1) =y(:,i)+h*[k1*((a0-(y(1,i)))*(b0-(y(1,i))));k2*(y(1,i)-y(2,i))];
end

%Cálculo de la concentración real de C tras formarse simultáneamente D
yC=y(1,:)-y(2,:);

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

%Método de Runge-Kutta 4
%Preparación vector solución para un tiempo t 
z=zeros(2,length(t));
z(:,1)=y0;%Inicializamos

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

%Cálculo de la concentración real de C tras formarse simultáneamente D
zC=z(1,:)-z(2,:);

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

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

Apreciamos cómo los dos compuestos C y D se forman simultáneamente, siendo la velocidad de formación de D más rápida. En un determinado punto vemos que el compuesto C deja de formarse ya que B (reactivo limitante) se ha agotado, y a partir de este momento, el reactivo C comienza a consumirse hasta agotarse, mientras que D sigue formandose hasta alcanzar un equilibrio. Cabe observar también cómo en el primer instante, tanto mediante el método de Euler como para el de Runge-Kutta, cuando empieza a formarse C el compuesto D no se forma hasta pasado un determinado tiempo. Esto es erróneo puesto que C y D se forman simultáneamente, pero no se trata de un error cometido por el paso sino por un error que cometen los métodos: el bucle, para el primer valor de i, calcula que la concentración de D da cero y no es hasta que se toma el segundo valor de i que nos proporciona un valor determinado.

Para un valor de paso h=0.3 quedaría:

Solución mediante el método de Euler para h=0.3

Solución mediante el método de Runge-Kutta para h=0.3

Gráfico comparativo para compuesto C y D para h=0.3

Como ya vimos anteriormente a medida que el paso, h, aumenta, el método de Euler va cometiendo más error y se vuelve más impreciso. Usa como ordenada, el punto perteneciente a la recta tangente del anterior punto, por lo que la diferencia entre la ordenada real y la aproximada es notable. Se observa que al aplicar Euler explícito con tamaño del paso decreciente, la solución mejora aproximándose a la solución exacta. Esta es la idea de convergencia del método. Por otra parte se ve que para valores del paso, superiores a una cierta cantidad, la solución tiene un carácter oscilatorio. En el caso de Runge-Kutta, vemos que la aproximación que se realiza es igualmente errónea, salvo que para este método se observa un carácter oscilatorio menos pronunciado, al ser Runge-Kutta un método de aproximación más exacto (Gráfico comparativo). Con lo que concluimos descartando la opción de utilizar un valor de paso h=0.3 por conllevar un gran error.

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

Solución mediante el método de Euler para h=0.1

Solución mediante el método de Runge-Kutta para h=0.1

Gráfico comparativo para compuesto C y D para h=0.1

%Sistema de ecuaciones no lineal por Trapecio y RK4
clear all; clc; clf;

%Datos del problema 
t0=0; tN=10; y0=[0;0]; k1=1; k2=1/5; a0=3; b0=1;

%Discretización temporal
h=0.1; %Paso del problema 
t=t0:h:tN;
N=round((tN-t0)/h); %N intervalos, N+1 valores

%Método de Euler 
%Preparación vector solución para un tiempo t 
y=zeros(2,length(t));
y(:,1)=y0; %Inicialización
for i=1:length(t)-1
  y(:,i+1) =y(:,i)+h*[k1*((a0-(y(1,i)))*(b0-(y(1,i))));k2*(y(1,i)-y(2,i))];
end

%Cálculo de la concentración real de C tras formarse simultáneamente D
yC=y(1,:)-y(2,:);

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

%Método de Runge-Kutta 4
%Preparación vector solución para un tiempo t 
z=zeros(2,length(t));
z(:,1)=y0;%Inicializamos

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

%Cálculo de la concentración real de C tras formarse simultáneamente D
zC=z(1,:)-z(2,:);

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

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

En este caso, podemos apreciar como la curva de la concentración de D tarda mucho más en formarse y alcanzar a C. Esto es debido a que la segunda reacción es mucho más lenta que la primera, por lo tanto, el reactivo C tardará más tiempo en agotarse y por ello, para el intervalo de tiempo [0,10] vemos como aún existe C a diferencia del apartado anterior. Mediante el método de Euler, como ya hemos visto acumulamos error. En cambio, Runge-Kutta nos da una aproximación más exacta, con un error prácticamente nulo.