Revisión actual del 22:19 6 mar 2015

1 Introducción y explicación del PVI

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

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

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

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

Se nos plantea la siguiente ecuación diferencial en el enunciado que, supuestamente, corresponde con la variación de la concentración de C a lo largo del tiempo:

                                                            [math]
   \left . 
 \begin{matrix} 
y'(t) = k_1 · (a_0 - y(t)) · (b_0 - y(t)),\\
  t\gt0
 \end{matrix}
   \right \}
[/math]

Tal y como podemos comprobar, si que podemos obtener la variación de la concentración del producto C a lo largo del tiempo a partir de dicha ecuación diferencial. Esto es así debido a que nuestra reacción química se trata de una reacción bimolecular y satisface a su vez la ley de acción de masas, conceptos que ya se han explicado.

Una vez conocidos todos estos datos, planteamos ahora nuestro Problema de Valor Inicial con el que trabajaremos en los siguientes apartados, que es el siguiente:

                                                            [math]
   \left . 
      \begin{matrix} 
                                                                y'(t) = k_1 · (a_0 - y(t)) · (b_0 - y(t)) \\
       A(t) = a_0 - y(t)\\B(t) = b_0 - y(t)
 \end{matrix}
   \right \}
[/math]

donde [math]a_0[/math]=concentración inicial del reactivo A, [math]b_0[/math]=concentración inicial del reactivo B y [math]k_1[/math]= constante de proporcionalidad de la velocidad de variación de la concentración del producto C.

Como hemos visto, la concentración de A [math](A(t) = a_0 - y(t))[/math] y la concentración de B [math](B(t) = b_0 - y(t))[/math] son inversamente proporcional a la concentración de C [math](C(t)=y(t))[/math]: las concentraciones de los dos reactivos A y B parten de su concentración inicial [math]a_0[/math] y [math] b_0[/math] van disminuyendo a medida que la concentración del producto C aumenta con el tiempo.

1.1 ¿Tiene solución única?

Podemos ver que la función es continua de 0 a ∞, por lo que ya aseguramos la existencia de al menos una solución. Observando la derivada parcial.

 ∂f∂y=k(2y-a-b)

No hay problemas de continuidad ya que podemos ver que es un polinomio, por tanto podemos afirmar que tiene una única solución.

1.2 ¿Qué sucedería si el proceso fuese reversible?

Hasta ahora hemos planteado el PVI suponiendo el proceso irreversible. En cambio, si lo suponemos reversible, únicamente tendríamos que añadir un componente negativo en la velocidad de reacción de la transformación de C a A+B , proporcional a la concentración de c [math]k_2*C[/math], destinada a la obtención de A y B.

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

Concentración reactivos k1=1 k2=2

Concentración reactivos k1=1 k2=3

En las gráficas vemos cómo las reactivos acaban estabilizándose. Siendo k1 la constante de proporcionalidad de la transformación de A+B a C y k2 la de transformación de C a A+B.

Si [math]k_2[/math]=∞ ó [math]k_1[/math]=0, no habría reacción. En el primer caso la transformación para obtener A y B es más rápida que la transformación en la que se pierde A+B. En el segundo nunca se llegaría a ceder nada de A y B.

Si [math]k_2[/math]=0 ó [math]k_1[/math]=∞, la reacción no sería reversible. En el primer caso no se obtendria nada de A y B por lo que sólo se obtendria C. En el segundo caso la obtención de C es mucho mayor que la de A y B.

2 Resolución Numérica del problema asociado

Ahora, nos disponemos a resolver numéricamente el problema de valor inicial descrito en el apartado anterior, suponiendo:

           [math]
          \left . 
      \begin{matrix} 
a_0 = 3 mol/l \\  
b_0 = 1 mol/l \\
k_1 = 1 mol/s 
     \end{matrix}
        \right \}
           [/math]

Utilizaremos, en este caso, el Método de Euler como método numérico para la resolución del problema, eligiendo un paso [math]h[/math]=0.1 en los primeros 2 segundos. El código Matlab resultante tras la aplicación numérica de este método a nuestro problema de valor inicial es el siguiente:

%PVI por Euler con una h = 0.1
clear all

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

y(1)= y0;

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

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

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


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

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

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

Observamos ahora en la gráfica la variación de las tres concentraciones según el paso del tiempo. Como era de esperar, las concentraciones de A y B disminuirán exponencialmente, mientras que la concentración de C, por el contrario, irá aumentando de manera logarítmica. Esto es debido a que los reactivos necesitan irse consumiendo para poder dar lugar a C, que va creciendo y se va formando más cantidad del mismo. El punto inicial de las tres concentraciones para un tiempo [math]t[/math]=0 representado en la gráfica nos indica la concentración de cada elemento que tenemos en el instante inicial, antes de que comience a darse lugar la reacción. Este valor coincide con los datos que se nos planteaban en el problema de valor inicial:

      [math]
          \left . 
      \begin{matrix} 
a_0 = 3 mol/l \\  
b_0 = 1 mol/l \\
k_1 = 1 mol/s 
     \end{matrix}
        \right \}
           [/math]

Las gráficas de las funciones de las concentraciones tienden a valores diferentes a medida que se avanza en el tiempo. Observamos que la gráfica que representa la variación de la concentración de A con el tiempo (color verde) tiende a 2, la de B (color rojo) tiende a 0 y la de C (color azul) tiende a 1. Esto es así debido a tratarse de una reacción bimolecular, que, como ya se ha explicado anteriormente, significa que una molécula de A y una de B producen una de C. Por tanto, esa es la razón por la que la concentración de A pasa de [math]a_0[/math]=3 mol/l a [math]A(t)[/math]=2 mol/l en función del tiempo t transcurrido, y así ocurre también con B y C, disminuyendo B de [math]b_0[/math]=1 mol/l a [math]B(t)[/math]=0 mol/l y, por el contrario, aumentando C de [math]c_0[/math]=0 mol/l a [math]C(t)[/math]=1 mol/l.

Cuando t tiende a infinito:

Una vez analizada cómo es la variación de las concentraciones en función del transcurso del tiempo, resolveremos numéricamente el mismo sistema para un tiempo lo suficientemente grande como para ser capaces de establecer una aproximación de los límites de las concentraciones cuando t → ∞. Para ello, lo único a modificar en el programa anterior es el dato del tiempo, que lo aumentaremos razonablemente para lograr ver hacia dónde tienden las gráficas con el paso del tiempo, mientras que el resto de datos se mantienen.

Ahora, la gráfica del Método de Euler para un tiempo t → ∞ queda de la siguiente manera:

Gráfica con un tiempo mayor para observar los límites de las concentraciones. Utilizando el método de Euler

Como ya se podía intuir en el anterior gráfico a este, las gráficas de las funciones de las concentraciones tienen sus límites para un tiempo t → ∞ determinados según la reacción bimolecular en la que están participando. Por tanto, tal y como ya se ha comentado, la concentración de A tiende a 2 mol/l, la de B tiende a 0 mol/l y la de C tiende a 1 mol/l, cuando la reacción tiende a un tiempo t → ∞.

2.1 Método del Trapecio

Resolvemos ahora nuestro mismo problema de valor inicial por otros dos métodos numéricos de aproximación, que son el Método del Trapecio y el Método de Runge Kutta. De esta manera, podemos comparar la solución al problema obtenida por los tres métodos numéricos empleados, debido a que, al ser métodos aproximados, las soluciones serán más precisas en función del método empleado. Mantenemos todos los datos iguales, utilizando también el mismo paso de tiempo.

Tendremos, por un lado, el problema de valor inicial resuelto por el Método del Trapecio.

Al ser un método de aproximación implícito, tenemos que despejar primero la variable [math]y_{n+1}[/math] que queremos calcular, que es la que necesitamos para operar el método numéricamente. Aplicando el método:

[math]y_{n+1} = y_n + \frac{h}{2}·(f(t_n,y_n)+f(t_{n+1},y_{n+1}))[/math]

Esta variable, una vez despejada, quedaría de la siguiente forma:

[math]y_{n+1} = \frac{-Q+\sqrt{Q^2-4·W·E}}{2·W}[/math]

siendo:

           [math]
          \left . 
      \begin{matrix} 
    Q = 1 + (h·k·(a+b))/2 \\
    W = -(h·k)/2 \\
    E = - y_n - (1/2)·h·k·((a_0 - y_n)·(b_0 - y_n) + a·b)
     \end{matrix}
        \right \}
           [/math]

Utilizando esta ecuación, el código Matlab del Método del Trapecio será:

%PVI por Trapecio con una h = 0.1
clear all

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

y(1)= y0;

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

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

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


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

Acompañado de su gráfica:

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

2.2 Método de Runge-Kutta de Orden 4

Y, por otro lado, el problema de valor inicial resuelto empleando el Método de Runge Kutta, quedando así:

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

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

y(1)= y0;

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

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

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


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

Acompañado también de su gráfica:

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

3 Segunda reacción. A partir de C obtenemos otro compuesto D

Suponemos ahora, según el enunciado dado, que está ocurriendo una reacción consecutiva que tiene la siguiente estructura:

                                                    A + B → C → D

La forma en la que esta reacción química funciona es muy similar a la reacción anterior con la que hemos trabajado: es igualmente una reacción bimolecular, luego sigue ocurriendo que una molécula de A y una de B van a producir una de C, pero ahora hay que considerar también el hecho de que este compuesto C se irá descomponiendo en la misma proporción en la que lo estaba haciendo anteriormente, dando lugar, por tanto, al compuesto D.

Para poder calcular la velocidad con la que se forma el compuesto D, tenemos también en cuenta la conservación de la Ley de Acción de masas en esta reacción consecutiva que estamos analizando.

La [math]y(t)[/math] calculada anteriormente, es decir, la concentración del compuesto C en función del tiempo, ahora pasaría a ser la suma de las concentraciones y(t)= yreal(t)+z(t) es decir: y(t)= [C] + [D]

                                  yreal(t)=y(t)- z(t)
 

Por tanto, después de analizar nuestra nueva reacción, definimos el nuevo problema de valor inicial asociado a esta reacción, que es:

[math] \left . \begin{matrix} yreal'(t)=y'(t)- z'(t)\\ z'(t)=k·yreal(t)\\ yreal'(t)=k_1·(a_0-y(t))·(b_0-y(t))-k_2·(y(t)- z(t))\\ a_0=3\\b_0=1\\y(0)=0\\z(0)=0 \end{matrix} \right \} [/math]

Será con este PVI con el que trabajaremos a partir de ahora.

3.1 Resolución numérica

Para la resolución numérica de este sistema, vamos a proponer dos casos diferentes con los que trabajar. Suponemos los mismos datos iniciales que hemos utilizado en la primera reacción, que son, a saber [math]a_0[/math]=3 , [math]b_0[/math]=1 y [math]k_1[/math]=1, para los primeros 10 segundos de duración de la reacción.

1) Utilizaremos, en el primer caso un valor de [math]k_2[/math]=5, lo que supone que [math]k_2[/math] es mayor que [math]k_1[/math], luego la segunda reacción es mucho más rápida que la primera.

2) En el segundo caso, utilizaremos un valor de [math]k_2[/math]=1/5, por lo que [math]k_2[/math] es menor que [math]k_1[/math], luego la segunda reacción es mucho más lenta que la primera.

Emplearemos dos Métodos Numéricos diferentes para resolver ambos casos, siendo estos el Método de Euler y el Método de Runge-Kutta. Esto nos va a permitir comparar los resultados obtenidos en los distintos procedimientos y así establecer una solución lo más precisa posible.

3.1.1 Euler

A la hora de crear el código Matlab con el método de Euler, ponemos una entrada de datos (input) tanto en el valor de la constante [math]k_2[/math], para poder realizar con el mismo código los dos casos descritos anteriormente, como en el tamaño de longitud de paso [math]h[/math], ya que vamos a utilizar varios valores también de la misma para comparar las distintas soluciones obtenidas.

El código numérico que resuelve nuestro sistema con el método de Euler es el siguiente:

%Establecemos los valores iniciales
t0=0;
tN=10;
y0=0;
k1=1;
a=3;
b=1;
k2=input('introduce k2: ');
%Discretizamos
h=input('introduce longitud de paso: ');
N=round((tN-t0)/h);
t=t0:h:tN;
%Creamos el vector y, donde almacenaremos las soluciones
y1=zeros(1,N+1);
yd=zeros(1,N+1);

y1(1)=y0;
yd(1)=y0;

yy1=y0;
yyd=y0;

%Utilizamos el 'for' para crear los vectores solución
for n=1:N
  yy1=yy1+h*(k1*(a-yy1)*(b-yy1)-k2*(yy1-yyd));
  yyd=yyd+h*(k2*(yy1-yyd));
  y1(n+1)=yy1;
  yd(n+1)=yyd;
end

%Definimos las variables de las concentraciones para dibujarlas
ya=3.-y1;
yb=1.-y1;
yc=y1-yd;

hold on

plot(t,yc,'b','linewidth',2)
plot(t,ya,'g','linewidth',2)
plot(t,yb,'r','linewidth',2)
plot(t,yd,'k','linewidth',2)
legend('Concentracion de C','Concentracion de A','Concentracion de B','Concentracion de D','Location','best');
hold off

Una vez hemos descrito el código, nos disponemos a pintar en gráficas las distintas evoluciones de las concentraciones de los elementos a lo largo del tiempo según nuestra reacción química consecutiva.

La primera gráfica que representamos es la correspondiente al primer caso que hemos descrito anteriormente, es decir, a un valor de [math]k_2[/math]=5 y a una longitud de paso [math]h[/math]=0.1, quedando de la siguiente manera:

Gráfica de las dos concentraciones A y B, el producto y la concentración D originada a partir de C. Utilizando el método de Euler.h=0.1

Como se puede observar, C y D empiezan a crecer mientras que A y B descienden a medida que avanza el tiempo hasta un límite. D siempre crecerá tendiendo a 1 debido a la estructura bimolecular de la reacción. Todo el compuesto de C se transforma para dar lugar a D. El compuesto C en este caso, al transformarse luego en el D, comienza creciendo hasta un máximo de 0.3 mol/l(en la gráfica no se aprecia muy bien con h=0.1)para luego decrecer hasta 0, ya que todo el compuesto se transforma en D. Este hecho es debido a que la velocidad de reacción de D es mayor que la velocidad de reacción de C, que era lo que suponía el utilizar una [math]k_2[/math] superior a la [math]k_1[/math] de la primera reacción, tal y como se ha explicado previamente.

A continuación, vamos a estudiar qué ocurriría si, utilizando la misma [math]k_2[/math] anterior, decidimos utilizar una longitud de paso [math]h[/math]=0.3. Introduciendo nuestros nuevos datos en el input y dibujando la gráfica solución del sistema, observamos lo que ocurre con la misma:

Gráfica de las dos concentraciones A y B, el producto y la concentración D originada a partir de C. Utilizando el método de Euler.h=0.3

Tal y como se puede observar en la gráfica, es un tamaño de paso demasiado grande para resolver este problema, ya que sitúa a D por encima de 1,a A por debajo de 2 y a B y C con valores negativos. Cabe destacar que, a pesar de que transcurridos los 10 segundos de duración de la reacción se puede observar que los límites de las gráficas de las concentraciones tienden a los mismos valores con los que lo hacía para el tamaño de longitud de paso anterior, no es un tamaño de h que permita ver claramente cuál es la solución del problema.

Tras este pequeño análisis, vamos a resolver ahora el segundo caso inicialmente descrito, es decir, emplearemos ahora una [math]k_2[/math]=1/5 y mantendremos la misma longitud de paso que hemos utilizado en la primera gráfica, [math]h[/math]=0.1.

La gráfica nos quedaría así:

Gráfica de las dos concentraciones A y B, el producto y la concentración D originada a partir de C. Utilizando el método de Euler.h=0.1

Al ser la velocidad de la reacción C 5 veces mas rápida que la de D, observamos que C tiene un máximo bastante mayor que D, en este caso de valor aproximado de 0.9 ( siendo 1 el límite al que tiende la concentración de C), y también se observa, en comparación con el primer caso, que la gráfica de la concentración de C decrece con menos pendiente hacia 0.

3.1.2 Runge-Kutta de Orden 4

Ahora nos disponemos a desarrollar el problema con el segundo método numérico anteriormente propuesto, el Método de Runge-Kutta. Igual que ocurría para el método de Euler, en este caso también ponemos una entrada de datos (input) tanto en el valor de la constante [math]k_2[/math] como en el tamaño de longitud de paso [math]h[/math], ya que también vamos a utilizar varios valores diferentes de las mismas.

El código numérico que resuelve nuestro PVI con el método de Runge-Kutta es el siguiente:

%C como reactivo y D como producto
t0=0;
tn=10;
h=input('introduce longitud de paso: ');
N=round((tn-t0)/h);
t=t0:h:tn;
y1=zeros(1,N+1);
yd=zeros(1,N+1);
k1=1;
k2=input('introduce k2: ');
a=3;
b=1;
for n=1:N
    
    %runge kutta 
    kc1=k1*(a-y1(n))*(b-y1(n));
    kc2=k1*(a-y1(n)-h/2*kc1)*(b-y1(n)-h/2*kc1);
    kc3=k1*(a-y1(n)-h/2*kc2)*(b-y1(n)-h/2*kc2);
    kc4=k1*(a-y1(n)-h*kc3)*(b-y1(n)-h*kc3);
    y1(n+1)=y1(n)+(h/6)*(kc1+2*kc2+2*kc3+kc4);
    kd1=k2*(y1(n)-yd(n));
    kd2=k2*(y1(n)-yd(n)-h/2*kd1);
    kd3=k2*(y1(n)-yd(n)-h/2*kd2);
    kd4=k2*(y1(n)-yd(n)-h*kd3);
    yd(n+1)=yd(n)+(h/6)*(kd1+2*kd2+2*kd3+kd4);
 
end
%concentraciones
ya=3.-y1;
yb=1.-y1;
yc=y1-yd;

hold on

plot(t,yc,'b','linewidth',2)
plot(t,ya,'g','linewidth',2)
plot(t,yb,'r','linewidth',2)
plot(t,yd,'k','linewidth',2)
legend('Concentracion de C','Concentracion de A','Concentracion de B','Concentracion de D','Location','best');
hold off

El Método de Runge-Kutta es una forma numérica de resolver problemas de valor inicial que trabaja de la misma manera que lo hace el Método de Euler, dando soluciones aproximadas al problema de valor inicial propuesto. Por tanto, las gráficas que nos van a salir con este método van a ser muy parecidas a las obtenidas con el Método de Euler, pero sabiendo de antemano que el Método de Runge-Kutta consigue hacer una mejor aproximación de las soluciones de la que puede hacer Euler.

Una vez conocido esto, nos disponemos a pintar las mismas gráficas representadas con el método anterior pero utilizando ahora Runge-Kutta:

Primer caso: valor de [math]k_2[/math]=5 y longitud de paso [math]h[/math]=0.1, quedando de la siguiente manera:

Gráfica de las dos concentraciones A y B, el producto y la concentración D originada a partir de C. Utilizando el método de Runge Kutta h=0.1 k2=5

Las gráficas evolucionan de la misma manera que lo hacían antes, dando sin embargo una solución más aproximada que la primera obtenida.

Estudiamos también qué ocurriría si decidimos utilizar una longitud de paso [math]h[/math]=0.3, con las mismas condiciones con las que lo hemos hecho con el anterior método pero usando ahora Runge-Kutta. Representamos la gráfica obtenida:

Gráfica de las dos concentraciones A y B, el producto y la concentración D originada a partir de C. Utilizando el método de Runge Kutta h=0.3 k2=5

El error para este método es notable, ya que C llega a valer mas de 0.5 cuando no debería sobrepasar 0.3. Por lo tanto, este tamaño de longitud de paso [math]h[/math]=0.3 no nos sirve tampoco para el Método de Runge-Kutta. Aunque no nos vale para ninguno de los dos métodos, ya que obtenemos errores en ambos dos, si que es preciso comentar que error obtenido con el Método de Runge-Kutta para [math]h[/math]=0.3 es un error menor al obtenido con el Método de Euler para [math]h[/math]=0.3.

Ahora el segundo caso: una [math]k_2[/math]=1/5 y mantendremos la misma longitud de paso que hemos utilizado en la primera gráfica, [math]h[/math]=0.1. La gráfica sería:

Gráfica de las dos concentraciones A y B, el producto y la concentración D originada a partir de C. Utilizando el método de Runge Kutta h=0.1 k2=1/5

Las gráficas evolucionan igual, siendo los valores mucho mejor aproximados.

3.2 Discusión y análisis

Una vez resueltos los dos casos anteriormente planteados ([math]k_2=5[/math] y [math]k_2=1/5[/math]) por los métodos de Euler y de Runge-Kutta, hacemos una comparación de los resultados obtenidos:

Comparación de los dos métodos con k2 = 5

Comparación de los dos métodos con k2 = 1/5

1) Primer caso [math]k_2=5[/math] por Euler y Runge-Kutta:

En los primeros instantes en los que ocurre la reacción, las concentraciones de C y D aumentan. Al imponer que la reacción de D es más rápida que la de C, se observa como la concentración de C comienza a decrecer (hasta alcanzar el valor 0 molar) mucho antes de lo que lo hace la concentración de D (que se estabiliza después que C, alcanzando el valor de 1 molar). Esto ocurre cuando la reacción alcanza un tiempo t de 1 segundo.

1) Segundo caso [math]k_2=1/5[/math] por Euler y Runge-Kutta:

Imponemos ahora que la reacción C es más rápida que D. Ocurre algo parecido al primer caso, pero ahora las reacciones evolucionan de manera más lenta en el tiempo, aunque alcanzan los mismo límites( C se va acercando a 0 mientras que D se estabiliza para 1 molar).

Se observa por tanto que, aunque en los dos casos lo único que cambia es la velocidad de las dos reacciones, la solución obtenida nos da al final las mismas concentraciones finales, teniendo en cuenta que para el segundo caso, tardan más tiempo en obtenerse, debido a que la reacción tarda más en estabilizarse.

En cuanto a los métodos numéricos, los dos nos dan la misma solución, diferenciándose el uno del otro en que el Método de Runge-Kutta nos ofrece una mayor aproximación que el de Euler. De acuerdo en la variación de las dos concentraciones podemos ver que es considerablemente similar, teniendo una variación algo más significativa la concentración que alcanza C en los dos métodos. También cabe destacar que a mayor tamaño de la longitud de paso, peor es la aproximación obtenida.

Por último, decir que el error en los métodos con un tamaño de paso h=0.3 es grande y lógico, Esto se debe al representar un curva con pequeños segmentos de rectas. Cada segmento de recta tiene una pendiente igual a la de la curva en cada punto.Cuantos más segmentos más próxima será la solución , y un segmento de recta cada 0.3 son 3 veces menos que los que habría cada 0.1. Por esto se explica los valores tan disparatados que presentan las gráficas. Un caso hipotético para una aproximación perfecta seria por ejemplo, un tamaño de paso de h=0.1^∞ .