Archivo:H03euler.jpg

2016-05-03T14:56:42Z

Medimed:

Reacciones complejas GRUPO 1A

2016-05-03T14:56:17Z

Medimed: /* Gráfica comparativa con h=0.1 */

{{ TrabajoED | Reacciones complejas Grupo 1A | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED15/16|Curso 2015-16]] | Pablo Medina Higueras
Jesús Caballero Pozo

Jaime Delage Ramírez }}

==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,
<big>'''A + B → C'''</big>
Supondremos 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.

==Interpretación de las constantes y PVI para calcular la concentración de C==

En primer lugar, interpretaremos las constantes y variables de la siguiente ecuación diferencial,esta ecuación diferencial será la de la ley de masas:
<big>'''y'(t) = k1(a0 − y(t))(b0 − y(t)), siendo t>0'''</big>
donde:
'''<big>y'(t)</big>''' = Representa la velocidad de la reacción química, es decir, la velocidad con que se produce C.
'''<big>k1</big>''' = Constante de proporcionalidad.
'''<big>a0</big>''' = Concentración inicial de A.
'''<big>y(t)</big>''' = Evolución de la concentración de C a lo largo del tiempo.
'''<big>b0</big>''' = Concentración inicial de B.

Es conveniente mencionar que la concentración de A respecto al tiempo será de a0 - y(t), y para B será de b0 - y(t), tal y como nos indica el enunciado, tomaremos el tiempo como mayor que 0: t>0.

Además, sabemos que para t = 0, la concentración de C será nula, por lo que: y(t = 0) = 0.

Con todo lo deducido anteriormente, llegamos a la conclusión que estamos ante un Problema de Valor Inicial (PVI) o de Cauchy:
<big>'''y'(t) = k1(a0 − y(t))(b0 − y(t)), siendo t>0'''</big>
<big>'''y(0) = 0'''</big>

Ahora procederemos a ver si nuestro problema tiene solución única o no, a través del Teorema de existencia y Unicidad. Observamos que la función f(t,y) = y'(t) es continua en el intervalo (I = (0,∞) ∩ B(0,0),r>0), por lo que admite al menos una solución. Por otro lado, observamos la derivada parcial:
<big>'''∂f/∂y=k(2y-a-b)'''</big>
y vemos que no hay problemas de continuidad, porque la derivada parcial nos resulta un polinomio, así que podremos afirmar que tiene una única solución.

== Proceso reversible==

Hasta ahora hemos supuesto que la reacción era irreversible, tal y como nos indicaba el enunciado, esto significaba que A y B se juntaban para formar C. Ahora nos pide ver cómo cambiaría la ecuación diferencial si el proceso fuera reversible; que sea reversible, quiere decirnos que además de que A y B se junten para formar C, también tendremos C descomponiéndose para formar A y B, por lo que la ecuación diferencial se vería afectada por este cambio. Necesitaremos una nueva constante de velocidad de reacción (k2) dado que suponemos que se sigue satisfaciendo la ley de masas. Este nuevo término disminuirá la cantidad de producto, por lo tanto el término aparecerá restando en nuestra ecuación.

Por lo que la ecuación diferencial quedará de la siguiente manera:
<big>'''y'(t) = k1(a0 − y(t))(b0 − y(t)) - k2 * y(t)'''</big>

==Resolución por el método de Euler==

Ahora, resolveremos el problema de valor inicial con el método de Euler de primer orden (método explícito). Del enunciado, sabemos que <big>'''f(t,y)=k1(a0-y(t))(b0-y(t))'''</big> , por lo tanto aplicamos la fórmula del método de Euler:

<big>'''yn+1 = yn + h * f(tn, yn)'''</big>

Cuando despejamos, obtenemos : <big>'''yn+1 = yn + h * k1(a0-y(t))(b0-y(t))'''</big>

Según los datos del enunciado, tenemos que:

'''<big>a0</big>''' = 3 [mol/L]
'''<big>b0</big>''' = 1 [mol/L]
'''<big>k1</big>''' = 1 [mol/s]
'''<big>h</big>''' = 0.1 (salto)
'''<big>t</big>''' = [0,2] (segundos)

===Código en Matlab===

{{matlab|codigo=
%Euler
clear all
%Datos enunciado
a0=3;
b0=1;
t0=0;
tN=2;
y0=0;
h=0.1;
k1=1;
%Número subintervalos
N=round((tN-t0)/h);
%Variable independiente t
t=t0:h:tN;
%Vector y
y=zeros(1,N+1);
y(1)=y0;
%Bucle
for i=1:N
y(i+1)=y(i)+h*k1*(a0-y(i))*(b0-y(i));
end
%Gráfica
hold on
plot(t,y,'b','linewidth', 3);
plot(t,(a0-y),'y','linewidth',3);
plot(t,(b0-y),'m','linewidth',3);
legend('C','A','B');
xlabel('tiempo(segundos)');
ylabel('concentracion[mol/l]');
hold off

}}

===Gráfica método de Euler===

[[Archivo:trabajoeceuler1.jpg|marco|centro]]
Como se puede observar, la concentración de C aumenta en la proporción que disminuyen los reactivos A y B, como es lógico. También se puede observar que la velocidad disminuye en el tiempo hasta ser 0 cuando se agota el reactivo limitante, B.

==Resolución del PVI por el método de Euler cuando '''t→∞'''==

Una vez hecho el método de Euler de primer orden para un intervalo de 2 segundos, es fácil observar que las concentraciones convergen y presuponemos que se mantendrán constantes, pero hay que comprobarlo para el límite '''t→∞'''. Es suficiente con realizar el mismo programa que el anterior pero modificando el tiempo hasta un valor suficientemente alto, comparable a '''∞''' para ver que se mantienen constantes las concentraciones, tomaremos 15 segundos debido a que la reacción es muy rápida. Esto se debe a que la velocidad de la reacción se anula cuando se acaba el contenido del reactivo limitante, que en este caso es B.

===Código en Matlab===

{{matlab|codigo=
%Euler
clear all
%Datos enunciado
a0=3;
b0=1;
t0=0;
tN=15;
y0=0;
h=0.1;
k1=1;
%Número subintervalos
N=round((tN-t0)/h);
%Variable independiente t
t=t0:h:tN;
%Vector y
y=zeros(1,N+1);
y(1)=y0;
%Bucle
for i=1:N
y(i+1)=y(i)+h*k1*(a0-y(i))*(b0-y(i));
end
%Gráfica
hold on
plot(t,y,'b','linewidth',3);
plot(t,(a0-y),'y','linewidth',3);
plot(t,(b0-y),'m','linewidth',3);
legend('C','A','B');
xlabel('tiempo(segundos)');
ylabel('concentracion[mol/l]');
hold off

}}

===Gráfica método de Euler (t→∞)===

[[Archivo:trabajoeceulerinf.jpg|marco|centro]]

Tal y como habíamos pensado, las concentraciones tienden a A=2 mol/l; B=0 mol/l; C=1 mol/l. Esto se debe a que el compuesto B es el limitante de la reacción, por lo cuál cuando se acaba, el resto de concentraciones se mantienen constantes.

==Resolución del PVI por el Método del Trapecio y por Runge Kutta de cuarto orden==

Una vez que hemos obtenido el PVI por el Método de Euler, queremos obtener el PVI por el método del Trapecio y por Runge Kutta, hay que tener en cuenta que vamos a utilizar los mismos datos que para el método de Euler, por lo que:

'''<big>a0</big>''' = 3 [mol/L]
'''<big>b0</big>''' = 1 [mol/L]
'''<big>k1</big>''' = 1 [mol/s]
'''<big>h</big>''' = 0.1 (salto)
'''<big>t</big>''' = [0,2] (segundos)

===Método del Trapecio===

Hay que tener en cuenta que el método del trapecio es un método implícito, a diferencia del método de Euler que era un método explícito. Como siempre partimos de: <big>'''f(t,y)=k1(a0-y(t))(b0-y(t))'''</big> , y tendremos que aplicar la fórmula del trapecio:
'''<big>yn+1 = yn + h/2*(f(tn,yn) + f(tn+1,yn+1)

El código de Matlab para el método del trapecio es el siguiente:

==== Código en Matlab====

{{matlab|codigo=
%Trapecio
clear all
%Datos del problema
t0=0; %Instante inicial
tN=2; %Instante final
h=0.1; %Tamaño de paso
k1=1; %valor de la constante específica de la velocidad de reacción
y0=0;
a0=3; %Concentración inicial de A
b0=1; %Concentración inicial de B
%Discretización
N=round((tN-t0)/h);
%Creamos el vector t
t=t0:h:tN;
%Preparación del vector solución y defino primer valor
y=zeros(1,N+1);
y(1)=y0;
%Aplicación del esquema numérico
for i=1:N
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)))));
end
%grafico
hold on
plot(t,y,'b','linewidth',2);
a=3.-y;
b=1.-y;
plot(t,(a0-y),'y','linewidth',2);
plot(t,(b0-y),'m','linewidth',2);
legend('C','A','B')
xlabel('tiempo(s)');
ylabel('concentracion[mol/l]');
hold off
}}

====Gráfica del Método del Trapecio====

[[Archivo:Graficatrapecioy.jpg|centro|600px]]

Vemos que esta gráfica es prácticamente igual a la del método de Euler, pero sabemos que esta gráfica nos da una aproximación mayor.

===Método de Runge Kutta===

El método de Runge Kutta, de cuarto orden, al contrario del método del Trapecio, se trata de un método explícito, igual que el de Euler, por lo que no es necesario hacer ninguna manipulación al esquema numérico, por lo que procederemos directamente. El código Matlab utilizado para el método de Runge-Kutta es el siguiente:

==== Código en Matlab ====

{{matlab|codigo=
%RUNGEKUTTA
clear all
%Datos del problema
t0=0; %Instante inicial
tN=2; %Instante final
h=0.1; %Tamaño de paso
k1=1; %valor de la constante específica de la velocidad de reacción
y0=0;
a0=3; %Concentración inicial de A
b0=1; %Concentración inicial de B
%Discretización
N=round((tN-t0)/h);
%Creamos el vector t
t=t0:h:tN;
%Preparación del vector solución y defino primer valor
y=zeros(1,N+1);
y(1)=y0;
%Aplicación del esquema numérico
for i=1:N
K1=k1*(a0-y(i))*(b0-y(i));
K2=k1*(a0-(y(i)+(h/2)*K1))*(b0-(y(i)+(h/2)*K1));
K3=k1*(a0-(y(i)+(h/2)*K2))*(b0-(y(i)+(h/2)*K2));
K4=k1*(a0-(y(i)+h*K3))*(b0-(y(i)+h*K3));
y(i+1)=y(i)+(h/6)*(K1+2*K2+2*K3+K4);
end
%grafico
hold on
plot(t,y,'b','linewidth',2);
a=3.-y;
b=1.-y;
plot(t,(a0-y),'y','linewidth',2);
plot(t,(b0-y),'m','linewidth',2);
legend('C','A','B')
xlabel('tiempo(s)');
ylabel('concentracion[mol/l]');
hold off
}}

==== Gráfica del Método de RUNGE KUTTA ====

[[Archivo:PARTE6TRABAJO.jpg|centro|600px]]

Una vez que hemos obtenido las tres gráficas de los métodos de Euler, Trapecio y Runge-Kutta, podemos observar que todos tienen una gran similitud, pero el método de Euler, se separará más de los demás debido a que es un método de orden 1 y su precisión es más reducida que el resto. Los otros dos métodos serán más similares entre ellos puesto que tienen un mayor orden que el de Euler, el método del trapecio es de segundo orden y el de Runge-Kutta es de cuarto orden, por lo que el método que tendrá una mayor precisión será el de Runge-Kutta, gracias a su orden.

==Reacción consecutiva==

Ahora, suponemos que hay una reacción consecutiva de la siguiente forma:

<big>'''A + B → C → D'''</big>

Esta reacción, al igual que la reacción anterior, se puede expresar en forma de problema de valor inicial. Para ello, según explica el enunciado, usaremos de nuevo la ley de acción de masas, pero esta vez lo haremos dos veces consecutivas, una para de pasar de A y B a C, y otra para pasar de C a D. Al imponerse la ley de acción de masas, suponemos que la velocidad con la que se forma compuesto C es directamente proporcional a las concentraciones de los compuestos A y B, y de igual forma la de formación de D lo es a la concentración de C que se ha formado.

El problema de valor inicial es un sistema de dos ecuaciones en el que la primera muestra la producción de C y la segunda muestra la producción de D, que está relacionada con la cantidad de producto C en función del tiempo.
Los valores en t=0 son las concentraciones iniciales de los productos, es decir, nulos. El sistema sería el siguiente:

<big>'''y'1(t) =k1(a0-y(t))(b0-y(t))'''</big>
<big>'''y'2(t) =k2(y1(t)-y2(t))'''</big>

En el que :

<big>'''k1'''</big>,como ya dijimos anteriormente, es la constante de proporcionalidad de la primera reacción
<big>'''k2'''</big> es la constante de proporcionalidad de la segunda reacción
<big>'''y'1(t)'''</big> es la velocidad de la reacción química de C
<big>'''y'2(t)'''</big> es la velocidad de la reacción química de D

Dado que los valores de las concentraciones en t=0 son nulos, es decir:
<big>'''y1 (t=0) = 0'''</big>
<big>'''y2 (t=0) = 0'''</big>

== Reacción consecutiva: obtención de la ecuación diferencial y PVI ==
[[Archivo:Reaccion consecutiva.jpg]]
Las reacciones irreversibles se pueden definir como aquellas que empiezan con un reactivo inicial y producen productos o intermediarios generalmente en una sola dirección. Las reacciones consecutivas son reacciones secuenciales irreversibles

===Resolución e interpretación del sistema por Euler y Runge-Kutta con k2=5===
=== Resolución numérica por Trapecio cuando k2=5 y h=0.1 ===
{{matlab|codigo=
%Sistema de ecuaciones no lineal por Trapecio y RK4
%Datos del problema
T0=0;
TN=10;
Y0=[0;0];
k1=1;
k2=5;
a0=3;
b0=1;

%Discretizamos para h=0.1
h=0.1
T=T0:h:TN;
N=round((TN-T0)/h);

%EULER
y=zeros(2,length(T));
y(:,1)=Y0;
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álculamos la concentración real de C
yC=y(1,:)-y(2,:);
%Representación Euler
figure(1)
hold on
plot(T,yC,'linewidth',2)
plot(T,y(2,:),'y','linewidth',2)
legend('Concentración de C','Concentración de D','location','best');
title('Solución mediante el método de Euler');
hold off

%RK 4
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álculamos la concentración real de C
zC=z(1,:)-z(2,:);

%Representación RK4
figure(2)
hold on
plot(T,zC,'linewidth',2)
plot(T,z(2,:),'y','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 conjunta
figure(3)
hold on
subplot(1,2,1)
hold on
plot(T,yC,'b','linewidth',2)
plot(T,zC,'y--','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,:),'y--','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
}}

===Gráfica por el metodo de Euler con h=0.1===
[[Archivo:Ecuaciones4b.jpg|centro|600px]]
===Gráfica RK4 con h=0.1===
[[Archivo:Ecuaciones4a.jpg|centro|600px]]
===Gráfica comparativa con h=0.1===
[[Archivo:Figura1234.jpg|centro|600px]]

Comentarios de las graficas
Se observa como el compuesto C comienza con gran velocidad. En el momento en que empieza a haber C se generara D, que como tiene una velocidad de creación aún mayor, en cuanto pasa un poco de tiempo el compuesto D termina por hacer desaparecer al C y el compuesto D sigue formandose hasta alcanzar el equilibrio. Tambien se observa (ambos métodos) que en los primeros instantes el compuesto C se empieza a formar, mientras que el D no hasta que que pasa un tiempo debido a un error de los métodos. En conclusión vemos como se van formando ambos compuestos teniendo una mayor velocidad D

=== Resolución numérica por Trapecio cuando k2=5 y h=0.3 ===
{{matlab|codigo=
%Sistema de ecuaciones no lineal por Trapecio y RK4
%Datos del problema
T0=0;
TN=10;
Y0=[0;0];
k1=1;
k2=5;
a0=3;
b0=1;

%Discretizamos para h=0.3
h=0.3
T=T0:h:TN;
N=round((TN-T0)/h);

%EULER
y=zeros(2,length(T));
y(:,1)=Y0;
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álculamos la concentración real de C
yC=y(1,:)-y(2,:);
%Representación Euler
figure(1)
hold on
plot(T,yC,'linewidth',2)
plot(T,y(2,:),'y','linewidth',2)
legend('Concentración de C','Concentración de D','location','best');
title('Solución mediante el método de Euler');
hold off

%RK 4
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álculamos la concentración real de C
zC=z(1,:)-z(2,:);

%Representación RK4
figure(2)
hold on
plot(T,zC,'linewidth',2)
plot(T,z(2,:),'y','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 conjunta
figure(3)
hold on
subplot(1,2,1)
hold on
plot(T,yC,'b','linewidth',2)
plot(T,zC,'y--','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,:),'y--','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
}}

===Gráfica por el metodode Euler con h=0.3====
[[Archivo:Ecuaciones5b.jpg|centro|600px]]
===Gráfica de RK4 con h=0.3====
[[Archivo:Ecuaciones5a.jpg|centro|600px]]
===Gráfica comparativa con h=0.3====
[[Archivo:Ecuaciones50.31.jpg|centro|600px]]

Para h=0.3 observamos que la gráfica se desequilibra y deducimos de ello que es mejor utilizar intervalos pequeños en este caso 0.1 ya que como hemos visto con 0.3 hace que la gráfica se empiece a separar de la solución real

=== Resolución numérica por Trapecio y RK4 cuando k2=1/5 ===
{{matlab|codigo=
%Sistema de ecuaciones no lineal por Trapecio y RK4
%Datos del problema
t0=0;
tN=10;
y0=[0;0];
k1=1;
k2=1/5;
a0=3;
b0=1;

%Discretización
h=0.1;
t=t0:h:tN;
N=round((tN-t0)/h);

%EULER
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álculamos la concentración real de C
yC=y(1,:)-y(2,:);

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

%RK4
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álculamos la concentración real de C
zC=z(1,:)-z(2,:);

%Representación De (RK4)
figure(2)
hold on
plot(t,zC,'linewidth',2)
plot(t,z(2,:),'y','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 conjunto
figure(3)
hold on
subplot(1,2,1)
hold on
plot(t,yC,'b','linewidth',1.5)
plot(t,zC,'y--','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,:),'y--','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
}}
===Gráfica por el método de Euler===
[[Archivo:Ecuacones3.jpg|centro|600px]]
===Gráfica RK-4===
[[Archivo:Ecuaciones2.jpg|centro|600px]]
===Gráfico comparativo===
[[Archivo:Ecuaciones1.jpg|centro|600px]]

En estos graficos podemos apreciar que la segunda reacción es más lenta que la primera y como la curva de la concentración de D tarda más en formarse y por tanto en alcanzar a C, es decir, el reactivo C tarda más en acabarse a diferencia del apartado anterior. Con R-K obtenemos una aproximación mayor y un error bastante menor que el obtenido mediante el método de Euler.

[[Categoría:Ecuaciones Diferentiales]]
[[Categoría:ED15/16]]
[[Categoría:Trabajos 2015-16]]

Archivo:Figura1234.jpg

2016-05-03T14:55:34Z

Medimed:

Reacciones complejas GRUPO 1A

2016-05-03T14:49:35Z

Medimed: /* Resolución numérica por Trapecio cuando k2=5 y h=0.3 */

{{ TrabajoED | Reacciones complejas Grupo 1A | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED15/16|Curso 2015-16]] | Pablo Medina Higueras
Jesús Caballero Pozo

Jaime Delage Ramírez }}

==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,
<big>'''A + B → C'''</big>
Supondremos 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.

==Interpretación de las constantes y PVI para calcular la concentración de C==

En primer lugar, interpretaremos las constantes y variables de la siguiente ecuación diferencial,esta ecuación diferencial será la de la ley de masas:
<big>'''y'(t) = k1(a0 − y(t))(b0 − y(t)), siendo t>0'''</big>
donde:
'''<big>y'(t)</big>''' = Representa la velocidad de la reacción química, es decir, la velocidad con que se produce C.
'''<big>k1</big>''' = Constante de proporcionalidad.
'''<big>a0</big>''' = Concentración inicial de A.
'''<big>y(t)</big>''' = Evolución de la concentración de C a lo largo del tiempo.
'''<big>b0</big>''' = Concentración inicial de B.

Es conveniente mencionar que la concentración de A respecto al tiempo será de a0 - y(t), y para B será de b0 - y(t), tal y como nos indica el enunciado, tomaremos el tiempo como mayor que 0: t>0.

Además, sabemos que para t = 0, la concentración de C será nula, por lo que: y(t = 0) = 0.

Con todo lo deducido anteriormente, llegamos a la conclusión que estamos ante un Problema de Valor Inicial (PVI) o de Cauchy:
<big>'''y'(t) = k1(a0 − y(t))(b0 − y(t)), siendo t>0'''</big>
<big>'''y(0) = 0'''</big>

Ahora procederemos a ver si nuestro problema tiene solución única o no, a través del Teorema de existencia y Unicidad. Observamos que la función f(t,y) = y'(t) es continua en el intervalo (I = (0,∞) ∩ B(0,0),r>0), por lo que admite al menos una solución. Por otro lado, observamos la derivada parcial:
<big>'''∂f/∂y=k(2y-a-b)'''</big>
y vemos que no hay problemas de continuidad, porque la derivada parcial nos resulta un polinomio, así que podremos afirmar que tiene una única solución.

== Proceso reversible==

Hasta ahora hemos supuesto que la reacción era irreversible, tal y como nos indicaba el enunciado, esto significaba que A y B se juntaban para formar C. Ahora nos pide ver cómo cambiaría la ecuación diferencial si el proceso fuera reversible; que sea reversible, quiere decirnos que además de que A y B se junten para formar C, también tendremos C descomponiéndose para formar A y B, por lo que la ecuación diferencial se vería afectada por este cambio. Necesitaremos una nueva constante de velocidad de reacción (k2) dado que suponemos que se sigue satisfaciendo la ley de masas. Este nuevo término disminuirá la cantidad de producto, por lo tanto el término aparecerá restando en nuestra ecuación.

Por lo que la ecuación diferencial quedará de la siguiente manera:
<big>'''y'(t) = k1(a0 − y(t))(b0 − y(t)) - k2 * y(t)'''</big>

==Resolución por el método de Euler==

Ahora, resolveremos el problema de valor inicial con el método de Euler de primer orden (método explícito). Del enunciado, sabemos que <big>'''f(t,y)=k1(a0-y(t))(b0-y(t))'''</big> , por lo tanto aplicamos la fórmula del método de Euler:

<big>'''yn+1 = yn + h * f(tn, yn)'''</big>

Cuando despejamos, obtenemos : <big>'''yn+1 = yn + h * k1(a0-y(t))(b0-y(t))'''</big>

Según los datos del enunciado, tenemos que:

'''<big>a0</big>''' = 3 [mol/L]
'''<big>b0</big>''' = 1 [mol/L]
'''<big>k1</big>''' = 1 [mol/s]
'''<big>h</big>''' = 0.1 (salto)
'''<big>t</big>''' = [0,2] (segundos)

===Código en Matlab===

{{matlab|codigo=
%Euler
clear all
%Datos enunciado
a0=3;
b0=1;
t0=0;
tN=2;
y0=0;
h=0.1;
k1=1;
%Número subintervalos
N=round((tN-t0)/h);
%Variable independiente t
t=t0:h:tN;
%Vector y
y=zeros(1,N+1);
y(1)=y0;
%Bucle
for i=1:N
y(i+1)=y(i)+h*k1*(a0-y(i))*(b0-y(i));
end
%Gráfica
hold on
plot(t,y,'b','linewidth', 3);
plot(t,(a0-y),'y','linewidth',3);
plot(t,(b0-y),'m','linewidth',3);
legend('C','A','B');
xlabel('tiempo(segundos)');
ylabel('concentracion[mol/l]');
hold off

}}

===Gráfica método de Euler===

[[Archivo:trabajoeceuler1.jpg|marco|centro]]
Como se puede observar, la concentración de C aumenta en la proporción que disminuyen los reactivos A y B, como es lógico. También se puede observar que la velocidad disminuye en el tiempo hasta ser 0 cuando se agota el reactivo limitante, B.

==Resolución del PVI por el método de Euler cuando '''t→∞'''==

Una vez hecho el método de Euler de primer orden para un intervalo de 2 segundos, es fácil observar que las concentraciones convergen y presuponemos que se mantendrán constantes, pero hay que comprobarlo para el límite '''t→∞'''. Es suficiente con realizar el mismo programa que el anterior pero modificando el tiempo hasta un valor suficientemente alto, comparable a '''∞''' para ver que se mantienen constantes las concentraciones, tomaremos 15 segundos debido a que la reacción es muy rápida. Esto se debe a que la velocidad de la reacción se anula cuando se acaba el contenido del reactivo limitante, que en este caso es B.

===Código en Matlab===

{{matlab|codigo=
%Euler
clear all
%Datos enunciado
a0=3;
b0=1;
t0=0;
tN=15;
y0=0;
h=0.1;
k1=1;
%Número subintervalos
N=round((tN-t0)/h);
%Variable independiente t
t=t0:h:tN;
%Vector y
y=zeros(1,N+1);
y(1)=y0;
%Bucle
for i=1:N
y(i+1)=y(i)+h*k1*(a0-y(i))*(b0-y(i));
end
%Gráfica
hold on
plot(t,y,'b','linewidth',3);
plot(t,(a0-y),'y','linewidth',3);
plot(t,(b0-y),'m','linewidth',3);
legend('C','A','B');
xlabel('tiempo(segundos)');
ylabel('concentracion[mol/l]');
hold off

}}

===Gráfica método de Euler (t→∞)===

[[Archivo:trabajoeceulerinf.jpg|marco|centro]]

Tal y como habíamos pensado, las concentraciones tienden a A=2 mol/l; B=0 mol/l; C=1 mol/l. Esto se debe a que el compuesto B es el limitante de la reacción, por lo cuál cuando se acaba, el resto de concentraciones se mantienen constantes.

==Resolución del PVI por el Método del Trapecio y por Runge Kutta de cuarto orden==

Una vez que hemos obtenido el PVI por el Método de Euler, queremos obtener el PVI por el método del Trapecio y por Runge Kutta, hay que tener en cuenta que vamos a utilizar los mismos datos que para el método de Euler, por lo que:

'''<big>a0</big>''' = 3 [mol/L]
'''<big>b0</big>''' = 1 [mol/L]
'''<big>k1</big>''' = 1 [mol/s]
'''<big>h</big>''' = 0.1 (salto)
'''<big>t</big>''' = [0,2] (segundos)

===Método del Trapecio===

Hay que tener en cuenta que el método del trapecio es un método implícito, a diferencia del método de Euler que era un método explícito. Como siempre partimos de: <big>'''f(t,y)=k1(a0-y(t))(b0-y(t))'''</big> , y tendremos que aplicar la fórmula del trapecio:
'''<big>yn+1 = yn + h/2*(f(tn,yn) + f(tn+1,yn+1)

El código de Matlab para el método del trapecio es el siguiente:

==== Código en Matlab====

{{matlab|codigo=
%Trapecio
clear all
%Datos del problema
t0=0; %Instante inicial
tN=2; %Instante final
h=0.1; %Tamaño de paso
k1=1; %valor de la constante específica de la velocidad de reacción
y0=0;
a0=3; %Concentración inicial de A
b0=1; %Concentración inicial de B
%Discretización
N=round((tN-t0)/h);
%Creamos el vector t
t=t0:h:tN;
%Preparación del vector solución y defino primer valor
y=zeros(1,N+1);
y(1)=y0;
%Aplicación del esquema numérico
for i=1:N
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)))));
end
%grafico
hold on
plot(t,y,'b','linewidth',2);
a=3.-y;
b=1.-y;
plot(t,(a0-y),'y','linewidth',2);
plot(t,(b0-y),'m','linewidth',2);
legend('C','A','B')
xlabel('tiempo(s)');
ylabel('concentracion[mol/l]');
hold off
}}

====Gráfica del Método del Trapecio====

[[Archivo:Graficatrapecioy.jpg|centro|600px]]

Vemos que esta gráfica es prácticamente igual a la del método de Euler, pero sabemos que esta gráfica nos da una aproximación mayor.

===Método de Runge Kutta===

El método de Runge Kutta, de cuarto orden, al contrario del método del Trapecio, se trata de un método explícito, igual que el de Euler, por lo que no es necesario hacer ninguna manipulación al esquema numérico, por lo que procederemos directamente. El código Matlab utilizado para el método de Runge-Kutta es el siguiente:

==== Código en Matlab ====

{{matlab|codigo=
%RUNGEKUTTA
clear all
%Datos del problema
t0=0; %Instante inicial
tN=2; %Instante final
h=0.1; %Tamaño de paso
k1=1; %valor de la constante específica de la velocidad de reacción
y0=0;
a0=3; %Concentración inicial de A
b0=1; %Concentración inicial de B
%Discretización
N=round((tN-t0)/h);
%Creamos el vector t
t=t0:h:tN;
%Preparación del vector solución y defino primer valor
y=zeros(1,N+1);
y(1)=y0;
%Aplicación del esquema numérico
for i=1:N
K1=k1*(a0-y(i))*(b0-y(i));
K2=k1*(a0-(y(i)+(h/2)*K1))*(b0-(y(i)+(h/2)*K1));
K3=k1*(a0-(y(i)+(h/2)*K2))*(b0-(y(i)+(h/2)*K2));
K4=k1*(a0-(y(i)+h*K3))*(b0-(y(i)+h*K3));
y(i+1)=y(i)+(h/6)*(K1+2*K2+2*K3+K4);
end
%grafico
hold on
plot(t,y,'b','linewidth',2);
a=3.-y;
b=1.-y;
plot(t,(a0-y),'y','linewidth',2);
plot(t,(b0-y),'m','linewidth',2);
legend('C','A','B')
xlabel('tiempo(s)');
ylabel('concentracion[mol/l]');
hold off
}}

==== Gráfica del Método de RUNGE KUTTA ====

[[Archivo:PARTE6TRABAJO.jpg|centro|600px]]

Una vez que hemos obtenido las tres gráficas de los métodos de Euler, Trapecio y Runge-Kutta, podemos observar que todos tienen una gran similitud, pero el método de Euler, se separará más de los demás debido a que es un método de orden 1 y su precisión es más reducida que el resto. Los otros dos métodos serán más similares entre ellos puesto que tienen un mayor orden que el de Euler, el método del trapecio es de segundo orden y el de Runge-Kutta es de cuarto orden, por lo que el método que tendrá una mayor precisión será el de Runge-Kutta, gracias a su orden.

==Reacción consecutiva==

Ahora, suponemos que hay una reacción consecutiva de la siguiente forma:

<big>'''A + B → C → D'''</big>

Esta reacción, al igual que la reacción anterior, se puede expresar en forma de problema de valor inicial. Para ello, según explica el enunciado, usaremos de nuevo la ley de acción de masas, pero esta vez lo haremos dos veces consecutivas, una para de pasar de A y B a C, y otra para pasar de C a D. Al imponerse la ley de acción de masas, suponemos que la velocidad con la que se forma compuesto C es directamente proporcional a las concentraciones de los compuestos A y B, y de igual forma la de formación de D lo es a la concentración de C que se ha formado.

El problema de valor inicial es un sistema de dos ecuaciones en el que la primera muestra la producción de C y la segunda muestra la producción de D, que está relacionada con la cantidad de producto C en función del tiempo.
Los valores en t=0 son las concentraciones iniciales de los productos, es decir, nulos. El sistema sería el siguiente:

<big>'''y'1(t) =k1(a0-y(t))(b0-y(t))'''</big>
<big>'''y'2(t) =k2(y1(t)-y2(t))'''</big>

En el que :

<big>'''k1'''</big>,como ya dijimos anteriormente, es la constante de proporcionalidad de la primera reacción
<big>'''k2'''</big> es la constante de proporcionalidad de la segunda reacción
<big>'''y'1(t)'''</big> es la velocidad de la reacción química de C
<big>'''y'2(t)'''</big> es la velocidad de la reacción química de D

Dado que los valores de las concentraciones en t=0 son nulos, es decir:
<big>'''y1 (t=0) = 0'''</big>
<big>'''y2 (t=0) = 0'''</big>

== Reacción consecutiva: obtención de la ecuación diferencial y PVI ==
[[Archivo:Reaccion consecutiva.jpg]]
Las reacciones irreversibles se pueden definir como aquellas que empiezan con un reactivo inicial y producen productos o intermediarios generalmente en una sola dirección. Las reacciones consecutivas son reacciones secuenciales irreversibles

===Resolución e interpretación del sistema por Euler y Runge-Kutta con k2=5===
=== Resolución numérica por Trapecio cuando k2=5 y h=0.1 ===
{{matlab|codigo=
%Sistema de ecuaciones no lineal por Trapecio y RK4
%Datos del problema
T0=0;
TN=10;
Y0=[0;0];
k1=1;
k2=5;
a0=3;
b0=1;

%Discretizamos para h=0.1
h=0.1
T=T0:h:TN;
N=round((TN-T0)/h);

%EULER
y=zeros(2,length(T));
y(:,1)=Y0;
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álculamos la concentración real de C
yC=y(1,:)-y(2,:);
%Representación Euler
figure(1)
hold on
plot(T,yC,'linewidth',2)
plot(T,y(2,:),'y','linewidth',2)
legend('Concentración de C','Concentración de D','location','best');
title('Solución mediante el método de Euler');
hold off

%RK 4
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álculamos la concentración real de C
zC=z(1,:)-z(2,:);

%Representación RK4
figure(2)
hold on
plot(T,zC,'linewidth',2)
plot(T,z(2,:),'y','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 conjunta
figure(3)
hold on
subplot(1,2,1)
hold on
plot(T,yC,'b','linewidth',2)
plot(T,zC,'y--','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,:),'y--','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
}}

===Gráfica por el metodo de Euler con h=0.1===
[[Archivo:Ecuaciones4b.jpg|centro|600px]]
===Gráfica RK4 con h=0.1===
[[Archivo:Ecuaciones4a.jpg|centro|600px]]
===Gráfica comparativa con h=0.1===
[[Archivo:Ecuaciones40.11.jpg|centro|600px]]

Comentarios de las graficas
Se observa como el compuesto C comienza con gran velocidad. En el momento en que empieza a haber C se generara D, que como tiene una velocidad de creación aún mayor, en cuanto pasa un poco de tiempo el compuesto D termina por hacer desaparecer al C y el compuesto D sigue formandose hasta alcanzar el equilibrio. Tambien se observa (ambos métodos) que en los primeros instantes el compuesto C se empieza a formar, mientras que el D no hasta que que pasa un tiempo debido a un error de los métodos. En conclusión vemos como se van formando ambos compuestos teniendo una mayor velocidad D

=== Resolución numérica por Trapecio cuando k2=5 y h=0.3 ===
{{matlab|codigo=
%Sistema de ecuaciones no lineal por Trapecio y RK4
%Datos del problema
T0=0;
TN=10;
Y0=[0;0];
k1=1;
k2=5;
a0=3;
b0=1;

%Discretizamos para h=0.3
h=0.3
T=T0:h:TN;
N=round((TN-T0)/h);

%EULER
y=zeros(2,length(T));
y(:,1)=Y0;
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álculamos la concentración real de C
yC=y(1,:)-y(2,:);
%Representación Euler
figure(1)
hold on
plot(T,yC,'linewidth',2)
plot(T,y(2,:),'y','linewidth',2)
legend('Concentración de C','Concentración de D','location','best');
title('Solución mediante el método de Euler');
hold off

%RK 4
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álculamos la concentración real de C
zC=z(1,:)-z(2,:);

%Representación RK4
figure(2)
hold on
plot(T,zC,'linewidth',2)
plot(T,z(2,:),'y','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 conjunta
figure(3)
hold on
subplot(1,2,1)
hold on
plot(T,yC,'b','linewidth',2)
plot(T,zC,'y--','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,:),'y--','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
}}

===Gráfica por el metodode Euler con h=0.3====
[[Archivo:Ecuaciones5b.jpg|centro|600px]]
===Gráfica de RK4 con h=0.3====
[[Archivo:Ecuaciones5a.jpg|centro|600px]]
===Gráfica comparativa con h=0.3====
[[Archivo:Ecuaciones50.31.jpg|centro|600px]]

Para h=0.3 observamos que la gráfica se desequilibra y deducimos de ello que es mejor utilizar intervalos pequeños en este caso 0.1 ya que como hemos visto con 0.3 hace que la gráfica se empiece a separar de la solución real

=== Resolución numérica por Trapecio y RK4 cuando k2=1/5 ===
{{matlab|codigo=
%Sistema de ecuaciones no lineal por Trapecio y RK4
%Datos del problema
t0=0;
tN=10;
y0=[0;0];
k1=1;
k2=1/5;
a0=3;
b0=1;

%Discretización
h=0.1;
t=t0:h:tN;
N=round((tN-t0)/h);

%EULER
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álculamos la concentración real de C
yC=y(1,:)-y(2,:);

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

%RK4
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álculamos la concentración real de C
zC=z(1,:)-z(2,:);

%Representación De (RK4)
figure(2)
hold on
plot(t,zC,'linewidth',2)
plot(t,z(2,:),'y','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 conjunto
figure(3)
hold on
subplot(1,2,1)
hold on
plot(t,yC,'b','linewidth',1.5)
plot(t,zC,'y--','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,:),'y--','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
}}
===Gráfica por el método de Euler===
[[Archivo:Ecuacones3.jpg|centro|600px]]
===Gráfica RK-4===
[[Archivo:Ecuaciones2.jpg|centro|600px]]
===Gráfico comparativo===
[[Archivo:Ecuaciones1.jpg|centro|600px]]

En estos graficos podemos apreciar que la segunda reacción es más lenta que la primera y como la curva de la concentración de D tarda más en formarse y por tanto en alcanzar a C, es decir, el reactivo C tarda más en acabarse a diferencia del apartado anterior. Con R-K obtenemos una aproximación mayor y un error bastante menor que el obtenido mediante el método de Euler.

[[Categoría:Ecuaciones Diferentiales]]
[[Categoría:ED15/16]]
[[Categoría:Trabajos 2015-16]]

Reacciones complejas GRUPO 1A

2016-05-03T14:48:50Z

Medimed: /* Resolución numérica por Trapecio cuando k2=5 y h=0.1 */

{{ TrabajoED | Reacciones complejas Grupo 1A | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED15/16|Curso 2015-16]] | Pablo Medina Higueras
Jesús Caballero Pozo

Jaime Delage Ramírez }}

==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,
<big>'''A + B → C'''</big>
Supondremos 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.

==Interpretación de las constantes y PVI para calcular la concentración de C==

En primer lugar, interpretaremos las constantes y variables de la siguiente ecuación diferencial,esta ecuación diferencial será la de la ley de masas:
<big>'''y'(t) = k1(a0 − y(t))(b0 − y(t)), siendo t>0'''</big>
donde:
'''<big>y'(t)</big>''' = Representa la velocidad de la reacción química, es decir, la velocidad con que se produce C.
'''<big>k1</big>''' = Constante de proporcionalidad.
'''<big>a0</big>''' = Concentración inicial de A.
'''<big>y(t)</big>''' = Evolución de la concentración de C a lo largo del tiempo.
'''<big>b0</big>''' = Concentración inicial de B.

Es conveniente mencionar que la concentración de A respecto al tiempo será de a0 - y(t), y para B será de b0 - y(t), tal y como nos indica el enunciado, tomaremos el tiempo como mayor que 0: t>0.

Además, sabemos que para t = 0, la concentración de C será nula, por lo que: y(t = 0) = 0.

Con todo lo deducido anteriormente, llegamos a la conclusión que estamos ante un Problema de Valor Inicial (PVI) o de Cauchy:
<big>'''y'(t) = k1(a0 − y(t))(b0 − y(t)), siendo t>0'''</big>
<big>'''y(0) = 0'''</big>

Ahora procederemos a ver si nuestro problema tiene solución única o no, a través del Teorema de existencia y Unicidad. Observamos que la función f(t,y) = y'(t) es continua en el intervalo (I = (0,∞) ∩ B(0,0),r>0), por lo que admite al menos una solución. Por otro lado, observamos la derivada parcial:
<big>'''∂f/∂y=k(2y-a-b)'''</big>
y vemos que no hay problemas de continuidad, porque la derivada parcial nos resulta un polinomio, así que podremos afirmar que tiene una única solución.

== Proceso reversible==

Hasta ahora hemos supuesto que la reacción era irreversible, tal y como nos indicaba el enunciado, esto significaba que A y B se juntaban para formar C. Ahora nos pide ver cómo cambiaría la ecuación diferencial si el proceso fuera reversible; que sea reversible, quiere decirnos que además de que A y B se junten para formar C, también tendremos C descomponiéndose para formar A y B, por lo que la ecuación diferencial se vería afectada por este cambio. Necesitaremos una nueva constante de velocidad de reacción (k2) dado que suponemos que se sigue satisfaciendo la ley de masas. Este nuevo término disminuirá la cantidad de producto, por lo tanto el término aparecerá restando en nuestra ecuación.

Por lo que la ecuación diferencial quedará de la siguiente manera:
<big>'''y'(t) = k1(a0 − y(t))(b0 − y(t)) - k2 * y(t)'''</big>

==Resolución por el método de Euler==

Ahora, resolveremos el problema de valor inicial con el método de Euler de primer orden (método explícito). Del enunciado, sabemos que <big>'''f(t,y)=k1(a0-y(t))(b0-y(t))'''</big> , por lo tanto aplicamos la fórmula del método de Euler:

<big>'''yn+1 = yn + h * f(tn, yn)'''</big>

Cuando despejamos, obtenemos : <big>'''yn+1 = yn + h * k1(a0-y(t))(b0-y(t))'''</big>

Según los datos del enunciado, tenemos que:

'''<big>a0</big>''' = 3 [mol/L]
'''<big>b0</big>''' = 1 [mol/L]
'''<big>k1</big>''' = 1 [mol/s]
'''<big>h</big>''' = 0.1 (salto)
'''<big>t</big>''' = [0,2] (segundos)

===Código en Matlab===

{{matlab|codigo=
%Euler
clear all
%Datos enunciado
a0=3;
b0=1;
t0=0;
tN=2;
y0=0;
h=0.1;
k1=1;
%Número subintervalos
N=round((tN-t0)/h);
%Variable independiente t
t=t0:h:tN;
%Vector y
y=zeros(1,N+1);
y(1)=y0;
%Bucle
for i=1:N
y(i+1)=y(i)+h*k1*(a0-y(i))*(b0-y(i));
end
%Gráfica
hold on
plot(t,y,'b','linewidth', 3);
plot(t,(a0-y),'y','linewidth',3);
plot(t,(b0-y),'m','linewidth',3);
legend('C','A','B');
xlabel('tiempo(segundos)');
ylabel('concentracion[mol/l]');
hold off

}}

===Gráfica método de Euler===

[[Archivo:trabajoeceuler1.jpg|marco|centro]]
Como se puede observar, la concentración de C aumenta en la proporción que disminuyen los reactivos A y B, como es lógico. También se puede observar que la velocidad disminuye en el tiempo hasta ser 0 cuando se agota el reactivo limitante, B.

==Resolución del PVI por el método de Euler cuando '''t→∞'''==

Una vez hecho el método de Euler de primer orden para un intervalo de 2 segundos, es fácil observar que las concentraciones convergen y presuponemos que se mantendrán constantes, pero hay que comprobarlo para el límite '''t→∞'''. Es suficiente con realizar el mismo programa que el anterior pero modificando el tiempo hasta un valor suficientemente alto, comparable a '''∞''' para ver que se mantienen constantes las concentraciones, tomaremos 15 segundos debido a que la reacción es muy rápida. Esto se debe a que la velocidad de la reacción se anula cuando se acaba el contenido del reactivo limitante, que en este caso es B.

===Código en Matlab===

{{matlab|codigo=
%Euler
clear all
%Datos enunciado
a0=3;
b0=1;
t0=0;
tN=15;
y0=0;
h=0.1;
k1=1;
%Número subintervalos
N=round((tN-t0)/h);
%Variable independiente t
t=t0:h:tN;
%Vector y
y=zeros(1,N+1);
y(1)=y0;
%Bucle
for i=1:N
y(i+1)=y(i)+h*k1*(a0-y(i))*(b0-y(i));
end
%Gráfica
hold on
plot(t,y,'b','linewidth',3);
plot(t,(a0-y),'y','linewidth',3);
plot(t,(b0-y),'m','linewidth',3);
legend('C','A','B');
xlabel('tiempo(segundos)');
ylabel('concentracion[mol/l]');
hold off

}}

===Gráfica método de Euler (t→∞)===

[[Archivo:trabajoeceulerinf.jpg|marco|centro]]

Tal y como habíamos pensado, las concentraciones tienden a A=2 mol/l; B=0 mol/l; C=1 mol/l. Esto se debe a que el compuesto B es el limitante de la reacción, por lo cuál cuando se acaba, el resto de concentraciones se mantienen constantes.

==Resolución del PVI por el Método del Trapecio y por Runge Kutta de cuarto orden==

Una vez que hemos obtenido el PVI por el Método de Euler, queremos obtener el PVI por el método del Trapecio y por Runge Kutta, hay que tener en cuenta que vamos a utilizar los mismos datos que para el método de Euler, por lo que:

'''<big>a0</big>''' = 3 [mol/L]
'''<big>b0</big>''' = 1 [mol/L]
'''<big>k1</big>''' = 1 [mol/s]
'''<big>h</big>''' = 0.1 (salto)
'''<big>t</big>''' = [0,2] (segundos)

===Método del Trapecio===

Hay que tener en cuenta que el método del trapecio es un método implícito, a diferencia del método de Euler que era un método explícito. Como siempre partimos de: <big>'''f(t,y)=k1(a0-y(t))(b0-y(t))'''</big> , y tendremos que aplicar la fórmula del trapecio:
'''<big>yn+1 = yn + h/2*(f(tn,yn) + f(tn+1,yn+1)

El código de Matlab para el método del trapecio es el siguiente:

==== Código en Matlab====

{{matlab|codigo=
%Trapecio
clear all
%Datos del problema
t0=0; %Instante inicial
tN=2; %Instante final
h=0.1; %Tamaño de paso
k1=1; %valor de la constante específica de la velocidad de reacción
y0=0;
a0=3; %Concentración inicial de A
b0=1; %Concentración inicial de B
%Discretización
N=round((tN-t0)/h);
%Creamos el vector t
t=t0:h:tN;
%Preparación del vector solución y defino primer valor
y=zeros(1,N+1);
y(1)=y0;
%Aplicación del esquema numérico
for i=1:N
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)))));
end
%grafico
hold on
plot(t,y,'b','linewidth',2);
a=3.-y;
b=1.-y;
plot(t,(a0-y),'y','linewidth',2);
plot(t,(b0-y),'m','linewidth',2);
legend('C','A','B')
xlabel('tiempo(s)');
ylabel('concentracion[mol/l]');
hold off
}}

====Gráfica del Método del Trapecio====

[[Archivo:Graficatrapecioy.jpg|centro|600px]]

Vemos que esta gráfica es prácticamente igual a la del método de Euler, pero sabemos que esta gráfica nos da una aproximación mayor.

===Método de Runge Kutta===

El método de Runge Kutta, de cuarto orden, al contrario del método del Trapecio, se trata de un método explícito, igual que el de Euler, por lo que no es necesario hacer ninguna manipulación al esquema numérico, por lo que procederemos directamente. El código Matlab utilizado para el método de Runge-Kutta es el siguiente:

==== Código en Matlab ====

{{matlab|codigo=
%RUNGEKUTTA
clear all
%Datos del problema
t0=0; %Instante inicial
tN=2; %Instante final
h=0.1; %Tamaño de paso
k1=1; %valor de la constante específica de la velocidad de reacción
y0=0;
a0=3; %Concentración inicial de A
b0=1; %Concentración inicial de B
%Discretización
N=round((tN-t0)/h);
%Creamos el vector t
t=t0:h:tN;
%Preparación del vector solución y defino primer valor
y=zeros(1,N+1);
y(1)=y0;
%Aplicación del esquema numérico
for i=1:N
K1=k1*(a0-y(i))*(b0-y(i));
K2=k1*(a0-(y(i)+(h/2)*K1))*(b0-(y(i)+(h/2)*K1));
K3=k1*(a0-(y(i)+(h/2)*K2))*(b0-(y(i)+(h/2)*K2));
K4=k1*(a0-(y(i)+h*K3))*(b0-(y(i)+h*K3));
y(i+1)=y(i)+(h/6)*(K1+2*K2+2*K3+K4);
end
%grafico
hold on
plot(t,y,'b','linewidth',2);
a=3.-y;
b=1.-y;
plot(t,(a0-y),'y','linewidth',2);
plot(t,(b0-y),'m','linewidth',2);
legend('C','A','B')
xlabel('tiempo(s)');
ylabel('concentracion[mol/l]');
hold off
}}

==== Gráfica del Método de RUNGE KUTTA ====

[[Archivo:PARTE6TRABAJO.jpg|centro|600px]]

Una vez que hemos obtenido las tres gráficas de los métodos de Euler, Trapecio y Runge-Kutta, podemos observar que todos tienen una gran similitud, pero el método de Euler, se separará más de los demás debido a que es un método de orden 1 y su precisión es más reducida que el resto. Los otros dos métodos serán más similares entre ellos puesto que tienen un mayor orden que el de Euler, el método del trapecio es de segundo orden y el de Runge-Kutta es de cuarto orden, por lo que el método que tendrá una mayor precisión será el de Runge-Kutta, gracias a su orden.

==Reacción consecutiva==

Ahora, suponemos que hay una reacción consecutiva de la siguiente forma:

<big>'''A + B → C → D'''</big>

Esta reacción, al igual que la reacción anterior, se puede expresar en forma de problema de valor inicial. Para ello, según explica el enunciado, usaremos de nuevo la ley de acción de masas, pero esta vez lo haremos dos veces consecutivas, una para de pasar de A y B a C, y otra para pasar de C a D. Al imponerse la ley de acción de masas, suponemos que la velocidad con la que se forma compuesto C es directamente proporcional a las concentraciones de los compuestos A y B, y de igual forma la de formación de D lo es a la concentración de C que se ha formado.

El problema de valor inicial es un sistema de dos ecuaciones en el que la primera muestra la producción de C y la segunda muestra la producción de D, que está relacionada con la cantidad de producto C en función del tiempo.
Los valores en t=0 son las concentraciones iniciales de los productos, es decir, nulos. El sistema sería el siguiente:

<big>'''y'1(t) =k1(a0-y(t))(b0-y(t))'''</big>
<big>'''y'2(t) =k2(y1(t)-y2(t))'''</big>

En el que :

<big>'''k1'''</big>,como ya dijimos anteriormente, es la constante de proporcionalidad de la primera reacción
<big>'''k2'''</big> es la constante de proporcionalidad de la segunda reacción
<big>'''y'1(t)'''</big> es la velocidad de la reacción química de C
<big>'''y'2(t)'''</big> es la velocidad de la reacción química de D

Dado que los valores de las concentraciones en t=0 son nulos, es decir:
<big>'''y1 (t=0) = 0'''</big>
<big>'''y2 (t=0) = 0'''</big>

== Reacción consecutiva: obtención de la ecuación diferencial y PVI ==
[[Archivo:Reaccion consecutiva.jpg]]
Las reacciones irreversibles se pueden definir como aquellas que empiezan con un reactivo inicial y producen productos o intermediarios generalmente en una sola dirección. Las reacciones consecutivas son reacciones secuenciales irreversibles

===Resolución e interpretación del sistema por Euler y Runge-Kutta con k2=5===
=== Resolución numérica por Trapecio cuando k2=5 y h=0.1 ===
{{matlab|codigo=
%Sistema de ecuaciones no lineal por Trapecio y RK4
%Datos del problema
T0=0;
TN=10;
Y0=[0;0];
k1=1;
k2=5;
a0=3;
b0=1;

%Discretizamos para h=0.1
h=0.1
T=T0:h:TN;
N=round((TN-T0)/h);

%EULER
y=zeros(2,length(T));
y(:,1)=Y0;
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álculamos la concentración real de C
yC=y(1,:)-y(2,:);
%Representación Euler
figure(1)
hold on
plot(T,yC,'linewidth',2)
plot(T,y(2,:),'y','linewidth',2)
legend('Concentración de C','Concentración de D','location','best');
title('Solución mediante el método de Euler');
hold off

%RK 4
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álculamos la concentración real de C
zC=z(1,:)-z(2,:);

%Representación RK4
figure(2)
hold on
plot(T,zC,'linewidth',2)
plot(T,z(2,:),'y','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 conjunta
figure(3)
hold on
subplot(1,2,1)
hold on
plot(T,yC,'b','linewidth',2)
plot(T,zC,'y--','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,:),'y--','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
}}

===Gráfica por el metodo de Euler con h=0.1===
[[Archivo:Ecuaciones4b.jpg|centro|600px]]
===Gráfica RK4 con h=0.1===
[[Archivo:Ecuaciones4a.jpg|centro|600px]]
===Gráfica comparativa con h=0.1===
[[Archivo:Ecuaciones40.11.jpg|centro|600px]]

Comentarios de las graficas
Se observa como el compuesto C comienza con gran velocidad. En el momento en que empieza a haber C se generara D, que como tiene una velocidad de creación aún mayor, en cuanto pasa un poco de tiempo el compuesto D termina por hacer desaparecer al C y el compuesto D sigue formandose hasta alcanzar el equilibrio. Tambien se observa (ambos métodos) que en los primeros instantes el compuesto C se empieza a formar, mientras que el D no hasta que que pasa un tiempo debido a un error de los métodos. En conclusión vemos como se van formando ambos compuestos teniendo una mayor velocidad D

=== Resolución numérica por Trapecio cuando k2=5 y h=0.3 ===
{{matlab|codigo=
%Sistema de ecuaciones no lineal por Trapecio y RK4
%Datos del problema
T0=0;
TN=10;
Y0=[0;0];
k1=1;
k2=5;
a0=3;
b0=1;

%Discretizamos para h=0.3
h=0.3
T=T0:h:TN;
N=round((TN-T0)/h);

%EULER
y=zeros(2,length(T));
y(:,1)=Y0;
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álculamos la concentración real de C
yC=y(1,:)-y(2,:);
%Representación Euler
figure(1)
hold on
plot(T,yC,'linewidth',2)
plot(T,y(2,:),'y','linewidth',2)
legend('Concentración de C','Concentración de D','location','best');
title('Solución mediante el método de Euler');
hold off

%RK 4
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álculamos la concentración real de C
zC=z(1,:)-z(2,:);

%Representación RK4
figure(2)
hold on
plot(T,zC,'linewidth',2)
plot(T,z(2,:),'y','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
=== Grafica conjunta de Euler y RK4 cuando k2=5 ===
%Gráfico conjunta
figure(3)
hold on
subplot(1,2,1)
hold on
plot(T,yC,'b','linewidth',2)
plot(T,zC,'y--','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,:),'y--','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
}}
===Gráfica por el metodode Euler con h=0.3====
[[Archivo:Ecuaciones5b.jpg|centro|600px]]
===Gráfica de RK4 con h=0.3====
[[Archivo:Ecuaciones5a.jpg|centro|600px]]
===Gráfica comparativa con h=0.3====
[[Archivo:Ecuaciones50.31.jpg|centro|600px]]

Para h=0.3 observamos que la gráfica se desequilibra y deducimos de ello que es mejor utilizar intervalos pequeños en este caso 0.1 ya que como hemos visto con 0.3 hace que la gráfica se empiece a separar de la solución real

=== Resolución numérica por Trapecio y RK4 cuando k2=1/5 ===
{{matlab|codigo=
%Sistema de ecuaciones no lineal por Trapecio y RK4
%Datos del problema
t0=0;
tN=10;
y0=[0;0];
k1=1;
k2=1/5;
a0=3;
b0=1;

%Discretización
h=0.1;
t=t0:h:tN;
N=round((tN-t0)/h);

%EULER
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álculamos la concentración real de C
yC=y(1,:)-y(2,:);

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

%RK4
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álculamos la concentración real de C
zC=z(1,:)-z(2,:);

%Representación De (RK4)
figure(2)
hold on
plot(t,zC,'linewidth',2)
plot(t,z(2,:),'y','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 conjunto
figure(3)
hold on
subplot(1,2,1)
hold on
plot(t,yC,'b','linewidth',1.5)
plot(t,zC,'y--','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,:),'y--','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
}}
===Gráfica por el método de Euler===
[[Archivo:Ecuacones3.jpg|centro|600px]]
===Gráfica RK-4===
[[Archivo:Ecuaciones2.jpg|centro|600px]]
===Gráfico comparativo===
[[Archivo:Ecuaciones1.jpg|centro|600px]]

En estos graficos podemos apreciar que la segunda reacción es más lenta que la primera y como la curva de la concentración de D tarda más en formarse y por tanto en alcanzar a C, es decir, el reactivo C tarda más en acabarse a diferencia del apartado anterior. Con R-K obtenemos una aproximación mayor y un error bastante menor que el obtenido mediante el método de Euler.

[[Categoría:Ecuaciones Diferentiales]]
[[Categoría:ED15/16]]
[[Categoría:Trabajos 2015-16]]

Reacciones complejas GRUPO 1A

2016-05-03T14:27:23Z