Diferencia entre revisiones de «Reacciones con Autocatálisis (A5)»

De MateWiki
Saltar a: navegación, buscar
(Introducción)
 
(No se muestran 39 ediciones intermedias del mismo usuario)
Línea 1: Línea 1:
 
{{ TrabajoED | Reacciones con Autocatálisis. Grupo 5A | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED14/15|Curso 2014-15]] |Javier Blanco Villarroel<br />Marta Cavero Guillén<br />Alba Bringas Gil<br />Irene Bendala Sugrañes<br />Paula Botella Andreu}}
 
{{ TrabajoED | Reacciones con Autocatálisis. Grupo 5A | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED14/15|Curso 2014-15]] |Javier Blanco Villarroel<br />Marta Cavero Guillén<br />Alba Bringas Gil<br />Irene Bendala Sugrañes<br />Paula Botella Andreu}}
 
+
[[Categoría:Ecuaciones Diferenciales]]
 +
[[Categoría:ED14/15]]
 +
[[Categoría:Trabajos 2014-15]]
 
<br />
 
<br />
  
 
=== Introducción ===
 
=== Introducción ===
 
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.
 
En este ejercicio analizaremos el caso particular en el que A se transforma en B pero suponiendo que la presencia de B hace de efecto catalítico en la reacción. Escribiremos este proceso como una reacción bimolecular, A + B →(k1) 2B, en la que B hace al mismo tiempo papel de reactivo y producto. Se pide:
 
 
=== Apartados 1 2 y 3 ===
 
  
 
Consideraremos una reacción química irreversible en una solución bien mezclada. <br />
 
Consideraremos una reacción química irreversible en una solución bien mezclada. <br />
Línea 31: Línea 25:
 
en la que B hace al mismo tiempo papel de reactivo y producto.<br />
 
en la que B hace al mismo tiempo papel de reactivo y producto.<br />
  
 +
 +
=== Apartados 1 2 y 3 ===
  
 
Llamaremos:<br />
 
Llamaremos:<br />
Línea 216: Línea 212:
 
ylabel('Concentración (mol/L)','FontSize', 11);
 
ylabel('Concentración (mol/L)','FontSize', 11);
 
legend('Concentración de B','Concentración de A','location','best')}}<br />
 
legend('Concentración de B','Concentración de A','location','best')}}<br />
 
(interpretación de las gráficas)
 
  
  
 
=== Apartado 4 ===
 
=== Apartado 4 ===
  
En este apartado resolveremos el problema anterior tratándolo como un sistema de dos variables.
+
En este apartado resolveremos el problema anterior tratándolo como un sistema de dos variables. Para ello emplearemos los métodos de Euler y Runge-Kutta.
 
+
[[Archivo:Grafica24bis.jpg|500px|thumb|best|Gráfica que relaciona relaciona la variación de la concentración de los dos compuestos de la mezcla con el tiempo.]]
+
 
+
  
 +
[[Archivo:Grafica24bis.jpg|500px|thumb|right|Gráfica que relaciona relaciona la variación de la concentración de los dos compuestos de la mezcla con el tiempo.]]
 
{{matlab | codigo=
 
{{matlab | codigo=
  
Línea 233: Línea 225:
 
h=0.1; N=(tf-t0)/h;
 
h=0.1; N=(tf-t0)/h;
 
t=t0:h:tf;
 
t=t0:h:tf;
w=zeros(2,N+1); %las filas es el número de  incógnitas en el sistema
+
w=zeros(2,N+1);  
 +
%las filas de la matriz w coinciden número de  incógnitas en el sistema
 
w(:,1)=w0;
 
w(:,1)=w0;
 
ww=w0;
 
ww=w0;
Línea 248: Línea 241:
  
 
subplot(1,2,1)
 
subplot(1,2,1)
plot(t,w,'*'); %MATLAB dibuja una grágica para cada columna, NO poner color para que nos ponga cada una de un colo
+
plot(t,w,'*');  
 +
%MATLAB dibuja una gráfica para cada columna
 +
%NO poner color para que nos ponga cada una de un color
 
legend('Concentración de A','Concentración de B', 'location','best')
 
legend('Concentración de A','Concentración de B', 'location','best')
 
xlabel('Tiempo (s)','FontSize', 11);
 
xlabel('Tiempo (s)','FontSize', 11);
Línea 278: Línea 273:
 
     zz=zz+(h/6)*(k1+2*k2+2*k3+k4);
 
     zz=zz+(h/6)*(k1+2*k2+2*k3+k4);
 
     z(:,n+1)=zz;
 
     z(:,n+1)=zz;
 +
  
 
end
 
end
Línea 290: Línea 286:
  
 
}}
 
}}
 +
  
 
   
 
   
  
Se denomina reacción autocatalítica aquella en la que uno de los productos actúa como catalizador, la velocidad de reacción presenta un máximo cuando las concentraciones de A y B son iguales, mientras que es igual a cero cuando la concentración de A o de B es nula. Esto se debe a que en una reacción autocatalítica si comenzamos con una cantidad pequeña de A la velocidad de reacción aumentaría a medida que se vaya formando más R. En el otro extremo, cuando haya desaparecido todo el componente A la velocidad ha de tender a cero. En la gráfica podemos observa como al principio la variación de la concentración de ambas sustancia es muy pequeña, pero a medida que se empieza a generarse B la velocidad va aumentando, hasta hacerse máxima en el punto que las concentraciones de A y B se igualan, y vuelve a descender hasta hacerse nula cuando desaparece todo el A.
+
En esta reacción biomolecular, tenemos inicialmente A en concentración <math>1\frac{mol}{l}</math> y B en concentración <math>0.1\frac{mol}{l}</math>, en ella B actúa como catalizador, lo que implica que a medida que su concentración va aumentando, la velocidad de la reacción aumenta de manera exponencial; esto se aprecia en la pendiente de las curvas de ambos reactivos, las cuales, al ser la representación de la derivada de las concentraciones, nos indican cómo de rápido varían dichas concentraciones (la velocidad con la que varían), esta pendientes son inicialmente prácticamente horizontales, y van aumentando a medida que se forma B.
 +
 
 +
También cabe destacar que las pendiente de las curvas de las concentraciones de A y B son iguales pero de sentido contrario para cada instante, lo cual nos indica que la velocidad con la que se pierde A es la misma con la que se genera B. Por otro parte, esto le da el carácter simétrico a la gráfica, ya que, teniendo en cuenta la ley de concentración de masas, todo lo que desaparece de A se "transforma" en B.
 +
 
 +
La máxima velocidad de la reacción se alcanzará a los cinco segundos del inicio de la reacción, instante en el que las concentraciones de A y B son iguales; a partir de ese momento la velocidad de la reacción empieza a disminuir, es decir, la pendiente de las curvas de las concentraciones se va acercando a la horizontal, haciéndose nula en el instante en que la concentración de A ha desaparecido por completo (a los diez segundo de iniciarse la reacción), momento en el que obtenemos la concentración final de B (<math>1\frac{mol}{l}</math>), siendo este el final de la reacción autocatalítica.
 +
 
 +
== Apartado 5 ==
 +
En este apartado se dispone de una reacción consecutiva propuesta por Lotka en 1920:<br />
 +
<math> A + X \rightarrow 2X </math><br />
 +
<math> X + Y \rightarrow 2Y </math><br />
 +
<math> Y \rightarrow B </math>
 +
 
 +
con constantes de reacción k1, k2 y k3 respectivamente.
 +
Donde A,B,X e Y son sustancias diferentes.<br />
 +
Partiendo de la ley de acción de masas y de la ley de conservación de la masa llegamos a 4 ecuaciones diferenciales que describen las concentraciones de los elementos en función del tiempo.
 +
 
 +
 
 +
La primera ecuacón diferencial se deduce fácilmente de la ley de conservación de la masa:
 +
<math> A+B+X+Y=cte  \rightarrow  A'+B'+X'+Y'=0 </math>
 +
 
 +
 
 +
La segunda ED se deduce también de forma inmediata de la ley de acción de masas, ya que la velocidad de formación de formación de B (velocidad de reacción) es proporcional a la concentración de Y, siendo k3 la constante de proporcionalidad:
 +
<math> B'= k3*Y </math>
 +
 
 +
 
 +
La tercera ecuación tiene una deducción ligeramente más compleja pero sabemos que:<br />
 +
 
 +
<math> Y formada = k2*X*Y </math><br />
 +
<math> Y gastada = B formada = k3*Y </math><br />
 +
 
 +
 
 +
Por lo tanto podemos deducir que la variación de Y con el tiempo es la diferencia entre ambas:
 +
 
 +
<math> Y' = k2*X*Y - K3*Y </math>
 +
 
 +
 
 +
 
 +
La cuarta y última ecuación ecuación se obtiene de una forma muy parecida a la anterior pues sabemos que:<br />
 +
 
 +
<math> X formada = k1*A*X</math><br />
 +
<math> X gastada = k2*X*Y</math><br />
 +
 
 +
 
 +
Por lo que de nuevo la variación de X se expresa como la difeencia de ambas:
 +
 
 +
<math> X'= k1*A*X - K2*X*Y</math>
 +
 
 +
 
 +
En el apartado 6 procedemos a resolver este sistema de 4 ecuaciones.
  
 
== Apartado 6 ==
 
== Apartado 6 ==
[[Archivo:Apartado61.png|500px|thumb|best]]
+
En este apartado plantearemos el problema de valor inicial asociado al sistema anterior y lo resolveremos con el método de Euler con tamaño de paso h=0.01 y h=0.001. Para ello hemos tomado el tiempo t en un intervalo (0,200), los siguientes valores k1=k2=2k3=0.1 y concentraciones iniciales [A]=5<math>\frac{mol}{l}</math>, [X]=5*10<sup>-4</sup><math>\frac{mol}{l}</math>, [Y]=10<sup>-5</sup><math>\frac{mol}{l}</math>, [B]=0<math>\frac{mol}{l}</math>. 
 +
[[Archivo:Apartado61.png|500px|thumb|Tamaño de paso h=0.01]]
 
{{matlab|codigo=
 
{{matlab|codigo=
 
t0=0;  
 
t0=0;  
Línea 341: Línea 387:
 
}}
 
}}
 
----------------------------
 
----------------------------
 +
[[Archivo:grafica11.png|500px|thumb|Tamaño de paso h=0.001]]
 
{{matlab|codigo=
 
{{matlab|codigo=
 
t0=0;  
 
t0=0;  
Línea 383: Línea 430:
 
title('Concentración - Tiempo (Euler)','FontName','Berlin sans FB','FontSize', 20);
 
title('Concentración - Tiempo (Euler)','FontName','Berlin sans FB','FontSize', 20);
 
legend('Concentración de A','Concentración de X', 'Concentración de Y', 'Concentración de B', 'location','east')
 
legend('Concentración de A','Concentración de X', 'Concentración de Y', 'Concentración de B', 'location','east')
 +
}}
 +
Ahora vamos a estudiar la estabilidad del PVI, al cambiar el paso de h=0.01 a h=0.001 a simple vista no se aprecia diferencia, para comparar mejor las gráficas resultantes hemos dibujado juntas las gráficas obtenidas para la concentración de Y con h=0.01 y h=0.001. De esta manera vemos que las dos quedan superpuestas por lo que todo indica que el PVI es estable.
 +
[[Archivo:grafica22.png|500px|thumb|Comparación de la concentración de Y para h=0.01 y h=0.001]]
 +
{{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=K1; K3=K1/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,Y,'-y');
 +
 +
t01=0;
 +
tf1=200;
 +
h1=0.01;
 +
N1=(tf1-t01)/h1;
 +
t1=t01:h1:tf1;
 +
A01=5;
 +
X01=5*(10^-4);
 +
Y01=10^(-5);
 +
B01=0;
 +
 +
A1=zeros(1,N1+1);
 +
X1=zeros(1,N1+1);
 +
Y1=zeros(1,N1+1);
 +
B1=zeros(1,N1+1);
 +
 +
A1(1)=A01;
 +
X1(1)=X01;
 +
Y1(1)=Y01;
 +
B1(1)=B01;
 +
 +
 +
K11=0.1; K12=K1; K13=K11/2;
 +
 +
for i=1:N1
 +
X1(i+1)=X1(i)+h1*(K11*A1(i)*X1(i)-K12*X1(i)*Y1(i));
 +
Y1(i+1)=Y1(i)+h1*(K12*X1(i)*Y1(i)-K13*Y1(i));
 +
B1(i+1)=B1(i)+h1*(K13*Y1(i));
 +
A1(i+1)= A1(i)+h1*(-K11*A1(i)*X1(i));
 +
end
 +
 +
plot(t1,Y1,'-r');
 +
 +
hold off
 +
 +
xlabel('Tiempo (s)','FontSize', 11);
 +
ylabel('Concentración (mol/L)','FontSize', 11);
 +
title('Concentración - Tiempo (Euler)','FontName','Berlin sans FB','FontSize', 20);
 +
legend('Concentración de Y con paso h=0.01','Concentración de Y con paso h=0.001', 'location','east')
 +
}}
 +
Se dice que un PVI es estable si pequeñas perturbaciones de la función o de las condiciones iniciales afectan poco a la solución. Para asegurarnos de que el problema es estable hemos modificado en pequeña medida las concentraciones iniciales, hemos tomado las siguientes concentraciones: [A]=5.00001, [B]=10<sup>-5</sup><math>\frac{mol}{l}</math> [X]=5.1*10<sup>-4</sup><math>\frac{mol}{l}</math> [Y]=20<sup>-5</sup><math>\frac{mol}{l}</math>.                                                                                                                                                                        En la gráfica obtenida para las concentraciones iniciales utilizadas inicialmente junto con la gráfica de las nuevas concentraciones observamos que no se producen cambios bruscos, podríamos decir que estos cambios afectan poco a la solución, luego afirmamos que el PVI es estable.
 +
[[Archivo:grafica321.png|500px|thumb|Concentración de X,Y,A,B para distintas concentraciones iniciales]]
 +
{{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=K1; K3=K1/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,'-r');
 +
plot(t,X,'-g');
 +
plot(t,Y,'-y');
 +
plot(t,B,'-');
 +
 +
t0=0;
 +
tf=200;
 +
h=0.001;
 +
N=(tf-t0)/h;
 +
t=t0:h:tf;
 +
A10=5.00001;
 +
X10=5.1*(10^-4);
 +
Y10=20^(-5);
 +
B10=0.00001;
 +
 +
A1=zeros(1,N+1);
 +
X1=zeros(1,N+1);
 +
Y1=zeros(1,N+1);
 +
B1=zeros(1,N+1);
 +
 +
A1(1)=A10;
 +
X1(1)=X10;
 +
Y1(1)=Y10;
 +
B1(1)=B10;
 +
 +
 +
K1=0.1; K2=K1; K3=K1/2;
 +
 +
for i=1:N
 +
X1(i+1)=X1(i)+h*(K1*A1(i)*X1(i)-K2*X1(i)*Y1(i));
 +
Y1(i+1)=Y1(i)+h*(K2*X1(i)*Y1(i)-K3*Y1(i));
 +
B1(i+1)=B1(i)+h*(K3*Y1(i));
 +
A1(i+1)= A1(i)+h*(-K1*A1(i)*X1(i));
 +
end
 +
 +
plot(t,A1,'-k');
 +
plot(t,X1,'-k');
 +
plot(t,Y1,'-k');
 +
plot(t,B1,'-k');
 +
hold off
 +
 +
xlabel('Tiempo (s)','FontSize', 11);
 +
ylabel('Concentración (mol/L)','FontSize', 11);
 +
title('Concentración - Tiempo (Euler)','FontName','Berlin sans FB','FontSize', 20);
 
}}
 
}}
  
 
== Apartado 7 ==
 
== Apartado 7 ==
 +
 +
En este último apartado resolveremos el sistema de reacciones encadenadas mediante el método de Heun.
  
 
[[Archivo:Grafica27bisbis.jpg|500px|thumb|best|Gráfica que relaciona relaciona la variación de la concentración de los cuatro compuestos de la mezcla con el tiempo.]]
 
[[Archivo:Grafica27bisbis.jpg|500px|thumb|best|Gráfica que relaciona relaciona la variación de la concentración de los cuatro compuestos de la mezcla con el tiempo.]]
Línea 447: Línea 650:
 
title('Concentración - Tiempo (Heun)','FontName','Berlin sans FB','FontSize', 20);
 
title('Concentración - Tiempo (Heun)','FontName','Berlin sans FB','FontSize', 20);
 
legend('Concentración de A','Concentración de X', 'Concentración de Y', 'Concentración de B', 'location','best')
 
legend('Concentración de A','Concentración de X', 'Concentración de Y', 'Concentración de B', 'location','best')
 
  
 
}}
 
}}
  
En la gráfica podemos ver como originalmente tenemos una alta concentración de A, que va transformándose, al principio lentamente y luego más rápido en X (alcanzando la máxima velocidad en torno a los 20 segundos cuando las concentraciones de A y X se igualan). A continuación, al compuesto X se transforma en Y, en una reacción de apenas 20 segundos, no llega a formarse una alta concentración de Y, en seguida empieza a decrecer su presencia transformándose en B, compuesto final, lo que se aprecia en la horizontalidad de la gráfica a la que llegamos una vez el contenido de Y ha desaparecido por completo.
+
En esta ocasión tenemos una seria de reacciones bimoleculares encadenadas, y mediante esta gráfica intentamos ilustrar las conexiones entre ellas.
 +
 
 +
En la primera reacción interviene A con una concentración inicial de <math>5\frac{mol}{l}</math> y X con una concentración inicial de 5*10<sup>-4</sup> <math>\frac{mol}{l}</math>, siendo X el catalizador de la reacción en este caso. La reacción alcanza su máxima velocidad a lo 20 segundos y finaliza, al consumirse todo a los 30 segundos aproximadamente, dando lugar a una concentración final de X de <math>5\frac{mol}{l}</math>.
 +
 
 +
A continuación, los <math>5\frac{mol}{l}</math> de X reaccionan con una concentración inicial de Y de 10<sup>-4</sup> que actúa como catalizador de esta segunda reacción. En este caso, la duración de la misma es mucho menor, de apenas 20 segundos, la concentración de X desaparece rápidamente por lo que sólo da tiempo a que se formen <math>3,5\frac{mol}{l}</math> de Y. La velocidad máxima de la reacción se alcanza a los 10 segundos de iniciarse.
 +
 
 +
Finalmente, en la última reacción de la serie, la concentración resultante de Y se transforma en B, cuya concentración inicial es cero, y que actúa como catalizador. Esta reacción, más lenta que la anterior, alcanza su velocidad máxima a los 20 segundos de iniciarse (aproximadamente) y finaliza a los 150 segundos con una concentración de B de <math>5\frac{mol}{l}</math>, que coincide con la concentración inicial de A, muestra de que el principio de conservación de la masa se cumple. Como en los casos anteriores la reacción, y con ella la serie, finaliza al igualarse la concentración de Y a cero.

Revisión actual del 10:11 6 mar 2015

Trabajo realizado por estudiantes
Título Reacciones con Autocatálisis. Grupo 5A
Asignatura Ecuaciones Diferenciales
Curso Curso 2014-15
Autores Javier Blanco Villarroel
Marta Cavero Guillén
Alba Bringas Gil
Irene Bendala Sugrañes
Paula Botella Andreu
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura


1 Introducción

Consideraremos 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,


[math] A + B \rightarrow C [/math]

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. En nuestro caso analizaremos el caso particular en el que A se transforma en B pero suponiendo que la presencia de B hace de efecto catalítico en la reacción. Escribiremos este proceso como una reacción bimolecular


[math] A + B \rightarrow k . 2B [/math]

en la que B hace al mismo tiempo papel de reactivo y producto.


2 Apartados 1 2 y 3

Llamaremos:

[math] x= A+B[/math]
[math] y= 2B[/math]


Despejando el sistema obtenemos que:

[math] B=\frac{y}{2}[/math]
[math] A=x-\frac{y}{2}[/math]


A partir del principio de conservación de masas sabemos que [math] x+y= cte [/math] por lo que [math] x'+ y' = 0[/math] , quedando demostrada la primera ecuación.

A partir de la ley de acción de masas obtenemos que la velocidad de reacción [math] v= k.A.B[/math] y además sabemos que [math]v= y' = -x'[/math]
Del sistema anterior obtenemos que [math] AB=\frac{y}{2}. (x-\frac{y}{2})[/math], sustituyendo:


[math]y' = k.\frac{y}{2}.(x-\frac{y}{2}) [/math] y como [math] x=cte-y[/math] podemos decir que [math]y' = k.\frac{y}{2}.(cte-y-\frac{y}{2}) [/math]


[math]\frac{k}{2} = k[/math] puesto que es una constante.


Seguimos despejando [math] y'=k.y(\frac{2}{3}cte-y) = k.y(cte-y) \rightarrow y'= k.x.y [/math]


Integrando la primera ecuación y sustituyendo el valor de x(t) en la segunda, obtenemos que [math] y'=-k.y²[/math] por lo que el PVI asociado será:

  • [math] y'=-k.y²[/math]
  • [math] y(x₀)= y₀[/math]


Para que tenga una solución única f(x,y) tiene que ser continua en (x0,y0) [math] \rightarrow [/math]Punto de condición inicial.
Además [math] \frac{\partial f}{\partial f} (x,y)[/math] tiene que ser continua en (x0,y0).


Como [math] \frac{\partial f}{\partial f}(x,y) = -2kyy' \rightarrow [/math] Es continua para todo valor (x0,y0).
Y [math]f(x,y)= -k.y²[/math] también lo es para todo valor (x0,y0), podemos afirmar que tiene una solución única.


Suponiendo que las concentraciones iniciales de A y B son [math]1\frac{mol}{l}[/math] y [math]0,01\frac{mol}{l}[/math] respectivamente, y [math]k=1\frac{mol}{l}[/math], resolvemos el PVI por el método de Euler, eligiendo un paso h = 0,1, en los primeros 10 segundos.


Gráficaque relaciona la concentración y el tiempo hecha mediante el método de Euler.
%Ejercicio 2.2 Euler
t0=0
tn=10
n=100
h=(tn-t0)/n;
t=t0:h:tn;

y0=1
y=zeros(1,n+1);
y(1)=y0;
yy=y0;

k=1

x0=0.01
x=zeros(1,n+1);
x(1)=x0;
xx=x0;

for i=1:n
   yy=yy+h*(-k*xx.*yy);
   xx=xx+h*k*(xx.*yy);
   y(i+1)=yy;
   x(i+1)=xx;
   end

plot(t,y,'g','linewidth',2) 
title('Concentración - Tiempo','FontName','Berlin sans FB','FontSize', 20);
xlabel('Tiempo (s)','FontSize', 11);
ylabel('Concentración (mol/L)','FontSize', 11);
legend('Concentración de B','Concentración de A','location','east')



Ahora lo resolveremos por el método del trapecio:


Gráfica que relaciona la concentración y el tiempo hecha mediante el método del trapecio.
%Ejercicio 2.3 trapecio
t0=0
tn=10
n=100
h=(tn-t0)/n;
t=t0:h:tn;

y0=1
y=zeros(1,n+1);
y(1)=y0;
yy=y0;

k=1

x0=0.01;
x=zeros(1,n+1);
x(1)=x0;
xx=x0;

for i=1:n
   yy=(((2*yy)/(-k*h))+xx.*yy)/((-2/(k*h))+(-k*xx));
   xx=(((2*xx)/(k*h))+xx.*yy)/((2/(k*h))-(k*yy));
   y(i+1)=yy;
   x(i+1)=xx;
   end

plot(t,x,'k','linewidth',2)
hold on
plot(t,y,'c','linewidth',2) 

title('Concentración - Tiempo (Trapecio)','FontName','Berlin sans FB','FontSize', 20);
xlabel('Tiempo (s)','FontSize', 11);
ylabel('Concentración (mol/L)','FontSize', 11);
legend('Concentración de B','Concentración de A','location','best'


Y por último por el método de Runge-Kutta:


Gráfica que relaciona la concentración y el tiempo hecha mediante el método de Runge-Kutta.
%Ejercicio 2.3 Runge Kutta
t0=0
tn=10
n=100
h=(tn-t0)/n;
t=t0:h:tn;

y0=1
y=zeros(1,n+1);
y(1)=y0;
yy=y0;

k=1

x0=0.01
x=zeros(1,n+1);
x(1)=x0;
xx=x0;

for i =1:n

 k1A=(-k*xx.*yy);
 k1B=(k*xx.*yy);
 y1=yy+h*k1A/2;
 x1=xx+h*k1B/2
 
 k2A=(-k*x1.*y1);
 k2B=(k*x1.*y1);
 y2=yy+h*k2A/2;
 x2=xx+h*k2B/2;
 
 k3A=(-k*x2.*y2);
 k3B=(k*x2.*y2);
 y3=yy+h*k3A/2;
 x3=xx+h*k3B/2;
 
 k4A=(-k*x3.*y3);
 k4B=(k*x3.*y3);
 
 yy=yy+((h/6)*(k1A+2*k2A+2*k3A+k4A));
 xx=xx+((h/6)*(k1B+2*k2B+2*k3B+k4B));
 y(i+1)=yy;
 x(i+1)=xx;
 end

figure(1)   
plot(t,x,'y','linewidth',2)
hold on
plot(t,y,'m','linewidth',2) 

title('Concentración - Tiempo (Runge-Kutta)','FontName','Berlin sans FB','FontSize', 20);
xlabel('Tiempo (s)','FontSize', 11);
ylabel('Concentración (mol/L)','FontSize', 11);
legend('Concentración de B','Concentración de A','location','best')



3 Apartado 4

En este apartado resolveremos el problema anterior tratándolo como un sistema de dos variables. Para ello emplearemos los métodos de Euler y Runge-Kutta.

Gráfica que relaciona relaciona la variación de la concentración de los dos compuestos de la mezcla con el tiempo.
t0=0; tf=10;
w0=[1;0.01];
h=0.1; N=(tf-t0)/h;
t=t0:h:tf;
w=zeros(2,N+1); 
%las filas de la matriz w coinciden número de  incógnitas en el sistema
w(:,1)=w0;
ww=w0;

%Euler

for n=1:N
    
    F=[-1*ww(1)*ww(2);1*ww(1)*ww(2)];
    
    ww=ww+h*F;
    w(:,n+1)=ww;
end

subplot(1,2,1)
plot(t,w,'*'); 
%MATLAB dibuja una gráfica para cada columna
%NO poner color para que nos ponga cada una de un color
legend('Concentración de A','Concentración de B', 'location','best')
xlabel('Tiempo (s)','FontSize', 11);
ylabel('Concentración (mol/L)','FontSize', 11);
title('Concentración - Tiempo (Euler)','FontName','Berlin sans FB','FontSize', 10);

%Runge Kutta

z0=[1;0.01]; %para volver a empezar
z=zeros(2,N+1); %las filas es el número de  incógnitas en el sistema
z(:,1)=z0;
zz=z0;

for n=1:N
    
    F=[-1*zz(1)*zz(2);1*zz(1)*zz(2)];
     
    k1=F;
    
    z1=F+(1/2)*k1*h;
    k2=z1;
    
    z2=F+(1/2)*k2*h;
    k3=z2;
    
    z3=F+k3*h;
    k4=z3;
    
    zz=zz+(h/6)*(k1+2*k2+2*k3+k4);
    z(:,n+1)=zz;


end

subplot(1,2,2)
plot(t,z,'o')
legend('Concentración de A','Concentración de B', 'location','best')
xlabel('Tiempo (s)','FontSize', 11);
ylabel('Concentración (mol/L)','FontSize', 11);
title('Concentración - Tiempo (Runge-Kutta)','FontName','Berlin sans FB','FontSize', 10);



En esta reacción biomolecular, tenemos inicialmente A en concentración [math]1\frac{mol}{l}[/math] y B en concentración [math]0.1\frac{mol}{l}[/math], en ella B actúa como catalizador, lo que implica que a medida que su concentración va aumentando, la velocidad de la reacción aumenta de manera exponencial; esto se aprecia en la pendiente de las curvas de ambos reactivos, las cuales, al ser la representación de la derivada de las concentraciones, nos indican cómo de rápido varían dichas concentraciones (la velocidad con la que varían), esta pendientes son inicialmente prácticamente horizontales, y van aumentando a medida que se forma B.

También cabe destacar que las pendiente de las curvas de las concentraciones de A y B son iguales pero de sentido contrario para cada instante, lo cual nos indica que la velocidad con la que se pierde A es la misma con la que se genera B. Por otro parte, esto le da el carácter simétrico a la gráfica, ya que, teniendo en cuenta la ley de concentración de masas, todo lo que desaparece de A se "transforma" en B.

La máxima velocidad de la reacción se alcanzará a los cinco segundos del inicio de la reacción, instante en el que las concentraciones de A y B son iguales; a partir de ese momento la velocidad de la reacción empieza a disminuir, es decir, la pendiente de las curvas de las concentraciones se va acercando a la horizontal, haciéndose nula en el instante en que la concentración de A ha desaparecido por completo (a los diez segundo de iniciarse la reacción), momento en el que obtenemos la concentración final de B ([math]1\frac{mol}{l}[/math]), siendo este el final de la reacción autocatalítica.

4 Apartado 5

En este apartado se dispone de una reacción consecutiva propuesta por Lotka en 1920:
[math] A + X \rightarrow 2X [/math]
[math] X + Y \rightarrow 2Y [/math]
[math] Y \rightarrow B [/math]

con constantes de reacción k1, k2 y k3 respectivamente. Donde A,B,X e Y son sustancias diferentes.
Partiendo de la ley de acción de masas y de la ley de conservación de la masa llegamos a 4 ecuaciones diferenciales que describen las concentraciones de los elementos en función del tiempo.


La primera ecuacón diferencial se deduce fácilmente de la ley de conservación de la masa: [math] A+B+X+Y=cte \rightarrow A'+B'+X'+Y'=0 [/math]


La segunda ED se deduce también de forma inmediata de la ley de acción de masas, ya que la velocidad de formación de formación de B (velocidad de reacción) es proporcional a la concentración de Y, siendo k3 la constante de proporcionalidad: [math] B'= k3*Y [/math]


La tercera ecuación tiene una deducción ligeramente más compleja pero sabemos que:

[math] Y formada = k2*X*Y [/math]
[math] Y gastada = B formada = k3*Y [/math]


Por lo tanto podemos deducir que la variación de Y con el tiempo es la diferencia entre ambas:

[math] Y' = k2*X*Y - K3*Y [/math]


La cuarta y última ecuación ecuación se obtiene de una forma muy parecida a la anterior pues sabemos que:

[math] X formada = k1*A*X[/math]
[math] X gastada = k2*X*Y[/math]


Por lo que de nuevo la variación de X se expresa como la difeencia de ambas:

[math] X'= k1*A*X - K2*X*Y[/math]


En el apartado 6 procedemos a resolver este sistema de 4 ecuaciones.

5 Apartado 6

En este apartado plantearemos el problema de valor inicial asociado al sistema anterior y lo resolveremos con el método de Euler con tamaño de paso h=0.01 y h=0.001. Para ello hemos tomado el tiempo t en un intervalo (0,200), los siguientes valores k1=k2=2k3=0.1 y concentraciones iniciales [A]=5[math]\frac{mol}{l}[/math], [X]=5*10-4[math]\frac{mol}{l}[/math], [Y]=10-5[math]\frac{mol}{l}[/math], [B]=0[math]\frac{mol}{l}[/math].

Tamaño de paso h=0.01
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=K1; K3=K1/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,'-r'); 
plot(t,X,'-g');
plot(t,Y,'-y');
plot(t,B,'-');
hold off

xlabel('Tiempo (s)','FontSize', 11);
ylabel('Concentración (mol/L)','FontSize', 11);
title('Concentración - Tiempo (Euler)','FontName','Berlin sans FB','FontSize', 20);
legend('Concentración de A','Concentración de X', 'Concentración de Y', 'Concentración de B', 'location','east')

Tamaño de paso h=0.001
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=K1; K3=K1/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,'-r'); 
plot(t,X,'-g');
plot(t,Y,'-y');
plot(t,B,'-');
hold off

xlabel('Tiempo (s)','FontSize', 11);
ylabel('Concentración (mol/L)','FontSize', 11);
title('Concentración - Tiempo (Euler)','FontName','Berlin sans FB','FontSize', 20);
legend('Concentración de A','Concentración de X', 'Concentración de Y', 'Concentración de B', 'location','east')

Ahora vamos a estudiar la estabilidad del PVI, al cambiar el paso de h=0.01 a h=0.001 a simple vista no se aprecia diferencia, para comparar mejor las gráficas resultantes hemos dibujado juntas las gráficas obtenidas para la concentración de Y con h=0.01 y h=0.001. De esta manera vemos que las dos quedan superpuestas por lo que todo indica que el PVI es estable.

Comparación de la concentración de Y para h=0.01 y h=0.001
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=K1; K3=K1/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,Y,'-y');

t01=0; 
tf1=200;
h1=0.01; 
N1=(tf1-t01)/h1;
t1=t01:h1:tf1;
A01=5; 
X01=5*(10^-4); 
Y01=10^(-5); 
B01=0;

A1=zeros(1,N1+1); 
X1=zeros(1,N1+1);
Y1=zeros(1,N1+1);
B1=zeros(1,N1+1);

A1(1)=A01; 
X1(1)=X01; 
Y1(1)=Y01; 
B1(1)=B01;


K11=0.1; K12=K1; K13=K11/2;

for i=1:N1
X1(i+1)=X1(i)+h1*(K11*A1(i)*X1(i)-K12*X1(i)*Y1(i));
Y1(i+1)=Y1(i)+h1*(K12*X1(i)*Y1(i)-K13*Y1(i));
B1(i+1)=B1(i)+h1*(K13*Y1(i));
A1(i+1)= A1(i)+h1*(-K11*A1(i)*X1(i));
end

plot(t1,Y1,'-r');

hold off

xlabel('Tiempo (s)','FontSize', 11);
ylabel('Concentración (mol/L)','FontSize', 11);
title('Concentración - Tiempo (Euler)','FontName','Berlin sans FB','FontSize', 20);
legend('Concentración de Y con paso h=0.01','Concentración de Y con paso h=0.001', 'location','east')

Se dice que un PVI es estable si pequeñas perturbaciones de la función o de las condiciones iniciales afectan poco a la solución. Para asegurarnos de que el problema es estable hemos modificado en pequeña medida las concentraciones iniciales, hemos tomado las siguientes concentraciones: [A]=5.00001, [B]=10-5[math]\frac{mol}{l}[/math] [X]=5.1*10-4[math]\frac{mol}{l}[/math] [Y]=20-5[math]\frac{mol}{l}[/math]. En la gráfica obtenida para las concentraciones iniciales utilizadas inicialmente junto con la gráfica de las nuevas concentraciones observamos que no se producen cambios bruscos, podríamos decir que estos cambios afectan poco a la solución, luego afirmamos que el PVI es estable.

Concentración de X,Y,A,B para distintas concentraciones iniciales
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=K1; K3=K1/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,'-r'); 
plot(t,X,'-g');
plot(t,Y,'-y');
plot(t,B,'-');

t0=0; 
tf=200;
h=0.001; 
N=(tf-t0)/h;
t=t0:h:tf;
A10=5.00001; 
X10=5.1*(10^-4); 
Y10=20^(-5); 
B10=0.00001;

A1=zeros(1,N+1); 
X1=zeros(1,N+1);
Y1=zeros(1,N+1);
B1=zeros(1,N+1);

A1(1)=A10; 
X1(1)=X10; 
Y1(1)=Y10; 
B1(1)=B10;


K1=0.1; K2=K1; K3=K1/2;

for i=1:N
X1(i+1)=X1(i)+h*(K1*A1(i)*X1(i)-K2*X1(i)*Y1(i));
Y1(i+1)=Y1(i)+h*(K2*X1(i)*Y1(i)-K3*Y1(i));
B1(i+1)=B1(i)+h*(K3*Y1(i));
A1(i+1)= A1(i)+h*(-K1*A1(i)*X1(i));
end

plot(t,A1,'-k'); 
plot(t,X1,'-k');
plot(t,Y1,'-k');
plot(t,B1,'-k');
hold off

xlabel('Tiempo (s)','FontSize', 11);
ylabel('Concentración (mol/L)','FontSize', 11);
title('Concentración - Tiempo (Euler)','FontName','Berlin sans FB','FontSize', 20);


6 Apartado 7

En este último apartado resolveremos el sistema de reacciones encadenadas mediante el método de Heun.

Gráfica que relaciona relaciona la variación de la concentración de los cuatro compuestos de la mezcla con el tiempo.


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;
aa=A0; xx=X0; yy=Y0; bb=B0;

K1=0.1; K2=K1; K3=K1/2;

for n=1:N
    
   K1A=-K1*aa*xx;
   K2A=K1A+K1A*h;
   
   K1X=K1*aa*xx-K2*xx*yy;
   K2X=K1X+K1X*h;
   
   K1Y=K2*xx*yy-K3*yy;
   K2Y=K1Y+K1Y*h;
   
   K1B=K3*yy;
   K2B=K1B+K1B*h;
    
    
   aa=aa+0.5*h*(K1A+K2A);
   xx=xx+0.5*h*(K1X+K2X);
   yy=yy+0.5*h*(K1Y+K2Y);
   bb=bb+0.5*h*(K1B+K2B);
   
   A(n+1)=aa;
   X(n+1)=xx;
   Y(n+1)=yy;
   B(n+1)=bb;

end


plot(t,A,'*r'); 
hold on
plot(t,X,'*g');
plot(t,Y,'*y');
plot(t,B,'*');
hold off

xlabel('Tiempo (s)','FontSize', 11);
ylabel('Concentración (mol/L)','FontSize', 11);
title('Concentración - Tiempo (Heun)','FontName','Berlin sans FB','FontSize', 20);
legend('Concentración de A','Concentración de X', 'Concentración de Y', 'Concentración de B', 'location','best')


En esta ocasión tenemos una seria de reacciones bimoleculares encadenadas, y mediante esta gráfica intentamos ilustrar las conexiones entre ellas.

En la primera reacción interviene A con una concentración inicial de [math]5\frac{mol}{l}[/math] y X con una concentración inicial de 5*10-4 [math]\frac{mol}{l}[/math], siendo X el catalizador de la reacción en este caso. La reacción alcanza su máxima velocidad a lo 20 segundos y finaliza, al consumirse todo a los 30 segundos aproximadamente, dando lugar a una concentración final de X de [math]5\frac{mol}{l}[/math].

A continuación, los [math]5\frac{mol}{l}[/math] de X reaccionan con una concentración inicial de Y de 10-4 que actúa como catalizador de esta segunda reacción. En este caso, la duración de la misma es mucho menor, de apenas 20 segundos, la concentración de X desaparece rápidamente por lo que sólo da tiempo a que se formen [math]3,5\frac{mol}{l}[/math] de Y. La velocidad máxima de la reacción se alcanza a los 10 segundos de iniciarse.

Finalmente, en la última reacción de la serie, la concentración resultante de Y se transforma en B, cuya concentración inicial es cero, y que actúa como catalizador. Esta reacción, más lenta que la anterior, alcanza su velocidad máxima a los 20 segundos de iniciarse (aproximadamente) y finaliza a los 150 segundos con una concentración de B de [math]5\frac{mol}{l}[/math], que coincide con la concentración inicial de A, muestra de que el principio de conservación de la masa se cumple. Como en los casos anteriores la reacción, y con ella la serie, finaliza al igualarse la concentración de Y a cero.