Ver el código de «Tratamiento de residuos mediante análisis compartimental (Grupo 15B)»


=Introducción=

Tenemos tres pantanos A, B y C con 300 Hm3 de agua cada uno y con capacidad máxima de 500Hm3. Los pantanos A y C pueden recibir agua procedente de un rio y también pueden desalojarla por canales artificiales, pero tanto la entrada como salida de agua está regulada. Al pantano B puede fluir agua del rio, pero este pantano no tiene canal de desalojo de agua. Se produce un vertido tóxico al rio y al activar la entrada y salida a los pantanos del agua procedente del río, al pantano B empieza a entrar agua contaminada a razón de 2Hm3/día con 10kg/Hm de contaminante, mientras que:

hacia el pantano A entra agua contaminada a razón de 2Hm3/día con 5kg/Hm3 de contaminante y el agua contaminada sale por el canal a razón de 3Hm3/día;
hacia el pantano C entra agua contaminada a razón de 2Hm3/día con 0,05kg/Hm3 de contaminante y el agua contaminada sale por el canal a razón de 1Hm3=día;

Cuando el pantano B se ha llenado, las autoridades cortan el acceso de agua al pantano B y filtran la entrada de agua a los pantanos A y C, de tal manera que deja de entrar agua contaminada.

En este momento, si designamos por xA(0), xB(0), xC(0) y VA, VB y VC la cantidad de contaminante y el volumen de agua respectivamente que hay en los pantanos A, B y C, determine estas cantidades.

=Condiciones iniciales=

Buscamos definir el aumento de la masa del contaminante respecto al tiempo en cada uno de los pantanos, incomunicados entre si. Para ello, calcularemos el desarrollo de la concentración, definida como la masa del soluto entre el volumen de la disolución.

[[Archivo:Situacion0analitico.JPG|732px|thumb|center|Situación inicial para el cálculo analítico]]

La cual aplicada al desarrollo del volumen respecto al tiempo nos permite obtener las ecuaciones diferenciales correspondientes.

==Pantano A:==

Dados los datos iniciales del pantano:

:<math>X_{entrada}= 5 Kg/Hm^3</math>

:<math>V_{entrada}= 2 Hm^3/día</math>

:<math>V_{salida}= 3 Hm^3/día</math>

:<math>V_{entrada soluto} (t)= 2 Hm^3/día \times 5 Kg/Hm^3=10 Kg/día</math>

:<math>V_{inicial}=300 Hm^3</math>

Podemos desarrollar las ecuaciones que definen las variables respecto al tiempo:


:<math>V(t)=V(0)+(V_{entrada}-V_{salida})t=300+(2-3)t</math>

:<math>C(t)=X(t)/(300+(2-3)t)</math>

:<math>V_{salida soluto} (t)=3 Hm^3/día \times Q(t)/(300-t)</math>

Y así, el incremento diferencial de masa respecto al tiempo será:

:<math>X_A’=10-3X_A(t)/(300-t)</math>

Esta ecuación diferencial lineal se resuelve analíticamente, obteniendo una solución tal que, con la condición inicial de <math>Q(0)=0</math> podemos reducir el haz a una solución final:

:<math>X_A(t)= 5(300-t)-5,55 \times 10^{-5} \times(300-t)</math>

==Pantano B==

En este caso el problema que se plantea es más sencillo, ya que este pantano no tiene salida antes de alcanzar las condiciones iniciales. Así pues:

:<math>X_{entrada} = 2 Hm^3/día</math>

:<math>V_{max}= 500Hm^3</math>

:<math>V_{inicial}= 300 Hm^3</math>

:<math>V_{salida}= 0 Hm^3</math>

:<math>X_{entrada}=10 Kg/Hm^3</math>

Al no vaciarse, tendrá un incremento directo de volumen:

:<math> \Delta V= 500-300= 200 Hm^3</math>

Con lo que calculamos el tiempo que tarda en llenarse:

:<math>t= 200 Hm^3/ 2 Hm^3/día = 100 días</math>

==Pantano C==

Este caso es idéntico al Pantano A:

:<math>V_{entrada}= 2 Hm^3/día</math>

:<math>V_{max}= 500 Hm^3</math>

:<math>V_{inicial}= 300 Hm^3</math>

:<math>X_{entrada}= 0,05 Kg/Hm^3</math>

:<math>V_{salida}=1 Hm^3/día</math>

:<math>V_{entrada soluto}= 0,05 Kg/Hm^3 \times 2 Hm^3/día= 0,1 Kg/día</math>

:<math>V(t)=V(0)+(V_{entrada}-V_{salida})t=300+(2-3)t</math>

:<math>C(t)=X(t)/(300+(2-1)t)</math>

:<math>V_{salida}= 300 Hm^3/díax X_C/(300+t)</math>

La ecuación será:

:<math>X_C'= 0,1-3X_C/(300+t)</math>

Siendo su solución en  <math>Q(0)=0</math> la siguiente:

:<math>X_C(t)= (300+t)/40-300^4/40(300+t)^3</math>

==Obtención de los Valores Iniciales==

Obtenidas las ecuaciones que definen los 3 pantanos, basta con aplicar la condición de t=100 días para obtener las condiciones iniciales de volumen de disolución y masa del soluto en cada uno de ellos:

Pantano A:

:<math>X_A=556 Kg</math>

:<math>V_A=200 Hm^3</math>

Pantano B:

:<math>X_B=2000 Kg</math>

:<math>V_B=500 Hm^3</math>

Pantano C:

:<math>X_C=6, 8359 Kg</math>

:<math>V_C=400 Hm^3</math>

=Resolución numérica=
==Planteamiento==

Para intentar eliminar el contaminante de los pantanos, se bombea agua entre estos, quedando la situación como sigue:
al pantano A sigue entrando agua desde el rio a razón de 2Hm3/día, pero sin contaminar, y sale por el canal artificial 3Hm3/día de agua contaminada. Del pantano B al A entra agua a razón de 3Hm3/día y de A hacia B a razón de 2Hm3/día;
al pantano C sigue entrando agua desde el rio a razón de2Hm3/día, pero sin contaminar y sale por el canal artificial 1Hm3/día de agua contaminada. Del pantano B recibe 1Hm3/día y sale hacia B 2Hm3/día;

Resolvemos numéricamente este problema de valor inicial con N: 50, 100, 300 (medimos 
en días) utilizando el método de Euler, método trapezoidal y método de Runge-Kuta, con una longitud de paso h = 1 y h = 0:1. 

Con la nueva dinámica que producen los trasbases, pasamos a trabajar con los siguientes datos, además de los datos iniciales anteriormente obtenidos:

[[Archivo:Situacion0numerico.JPG|777px|thumb|center|Situación inicial para el cálculo numérico]]

A partir de ahora los pantanos no trabajan de forma independiente, sino que el desarrollo de cada uno es dependiente de los otros dos. Así pues, trabajaremos con el siguiente sistema de ecuaciones diferenciales:

[[Archivo:Sistema.JPG|252px|thumb|center|Sistema a resolver por los distintos métodos numéricos]]
==Métodos==
A continuación se resolverá el problema mediante los métodos de Euler,Trapezoidal y Runge-Kuta. Estos métodos pueden ser utilizados para cualquier intervalo de tiempo, por grande que sea y variando su h se conseguirá mayor o menor precisión en los resultados, a menos h, mayor precisión.
==Método de Euler==
El código que se presenta está resuelto para el caso particular de h=0.1 y pasados 50 días, para el resto de las gráficas el método es el mismo, variando h=1 y tf=100 o tf=300

{{matlab|codigo=

%Ejercicio de mezclas por el método de Euler
t0=0;
tf=50;
y0=[556 2000 6.8359]';
A=[-1/40 3/500 0; 1/100 -1/125 1/200; 0 1/500 -3/400]; 
%matriz de coeficientes
%Datos de la discretizacion
h=0.1;
N=(tf-t0)/h; 
%Construir de los vectores de tiempos y solucion aproximada
t=t0:h:tf;
y=zeros(3,N+1); 
%Hacemos una matriz con todas las soluciones conjuntas
%Inicialización
y(:,1)=y0; 
%la primera columna de la matriz de soluciones vale el valor de y0
yy=y0;
%Iteraciones y(tn)>>y(tn+1)
for n=1:N
  yy=yy+h*(A*yy);  
  y(:,n+1)=yy;
 end
clf
plot(t,y)
xlabel('Tiempo(días)')
ylabel('Cantidad soluto (kg)')
title('Metodo Euler h=0,1 N=50')
legend('Pantano A','Pantano B', 'Pantano C')
}}

[[Archivo:Grafico_Euler_h01_n50.jpg|400px|thumb|left|Gráfico para h=0.1 y n=50]]
[[Archivo:Grafico_Euler_h1_n50.jpg|400px|thumb|center|Gráfico para h=1 y n=50]]
[[Archivo:Grafico_Euler_h01_n100.jpg|400px|thumb|left|Gráfico para h=0.1 y n=100]]
[[Archivo:Grafico_Euler_h1_n100.jpg|400px|thumb|center|Gráfico para h=1 y n=100]]
[[Archivo:Grafico_Euler_h01_n300.jpg|400px|thumb|left|Gráfico para h=0.1 y n=300]]
[[Archivo:Grafico_Euler_h1_n300.jpg|400px|thumb|center|Gráfico para h=1 y n=300]]


_

==Método Trapezoidal== 
El código que se presenta está resuelto para el caso particular de h=0.1 y pasados 50 días, para el resto de las gráficas el método es el mismo, variando h=1 y tf=100 o tf=300.
{{matlab|codigo=
%Ejercicio de mezclas por el método del Trapecio
%Datos del problema
t0=0;
tf=50;
y0=[556 2000 6.8359]';
A=[-1/40 3/500 0; 1/100 -1/125 1/200; 0 1/500 -3/400];
%Matriz de coeficientes
%Datos de la discretizacion
h=0.1;
N=(tf-t0)/h; 
%Construir de los vectores de tiempos y solucion aproximada
t=t0:h:tf;
y=zeros(3,N+1); 
%Hacemos una matriz con todas las soluciones conjuntas
%Inicializacion
y(:,1)=y0; 
%La primera columna de la matriz de soluciones vale el valor de y0
yy=y0;
%Iteraciones y(tn)>>y(tn+1)
for n=1:N
  yy=(eye(3)-h/2*A)\((eye(3)+h/2*A)*yy);
  y(:,n+1)=yy;
 end
clf
plot(t,y)
xlabel('Tiempo(días)')
ylabel('Cantidad soluto (kg)')
title('Metodo Trapecio con h=0.1 y N=50')
legend('Pantano A','Pantano B','Pantano C')
legend('Pantano A','Pantano B','Pantano C')
}}
[[Archivo:Grafico_Trapecio_h01_n50.jpg|400px|thumb|left|Gráfico para h=0.1 y n=50]]
[[Archivo:Grafico_Trapecio_h1_n50.jpg|400px|thumb|center|Gráfico para h=1 y n=50]]
[[Archivo:Grafico_Trapecio_h01_n100.jpg|400px|thumb|left|Gráfico para h=0.1 y n=100]]
[[Archivo:Grafico_Trapecio_h1_n100.jpg|400px|thumb|center|Gráfico para h=1 y n=100]]
[[Archivo:Grafico_Trapecio_h01_n300.jpg|400px|thumb|left|Gráfico para h=0.1 y n=300]]
[[Archivo:Grafico_Trapecio_h1_n300.jpg|400px|thumb|center|Gráfico para h=1 y n=300]]

_
==Método de Runge-Kuta==
El código que se presenta está resuelto para el caso particular de h=0.1 y pasados 50 días, para el resto de las gráficas el método es el mismo, variando h=1 y tf=100 o tf=300.
{{matlab|codigo=
%Ejercicio de mezclas por el método Runge-Kuta
%Datos del problema
t0=0;
tf=50;
y0=[556 2000 6.8359]';
A=[-1/40 3/500 0; 1/100 -1/125 1/200; 0 1/500 -3/400];
%Datos de la discretizacion
h=0.1;
N=tf;
%Vectores de tiempos y aproximaciones
t=t0:h:tf;
y=zeros(3,N+1);
%metemos una matriz de ceros que iremos rellenando con el bucle 
%Inicializacion
y(:,1)=y0; 
%Para guardar la aproxmacion en cada tiempo vamos a usar yy es un valor numerico q va cambiando
yy=y0; 
%Ahora ya hacemos el bucle que vaya calculando el valor siguiente en cada instante las aproximaciones nos permite pasar del y(tn) al y(tn+1).
%Tendremos que hacer tantas iteraciones como intervalos, para eso tenemos N
for n=1:N
  %ahora depende del metodo numerico, antes de hacer la media ponderada tenemos q hacer la constante
  %Calculo k1=f(tn,yn)
    k1=[-1/40*yy(1)+3/500*yy(2),0.01*yy(1)-1/125*yy(2)+1/200*yy(3),1/500*yy(2)-3/400*yy(3)]';
  %Hemos metido el f(tn,yn) como una matriz en la q cada columna es una %de las dos ecuciacones del sistema y1 e y2 para k2 calculamos el argumento en el 
  %tiempo y en el valor de y
  yp=yy+0.5*h*k1;
  k2=[-1/40*yp(1)+3/500*yp(2),0.01*yp(1)-1/125*yp(2)+1/200*yp(3),1/500*yp(2)-3/400*yp(3)]';  
  %para k3 el argumento en el tiempo es el mismo q en k2 pero cambia el argumento en y
  yp=yy+1/2*h*k2;
  k3=[-1/40*yp(1)+3/500*yp(2),0.01*yp(1)-1/125*yp(2)+1/200*yp(3),1/500*yp(2)-3/400*yp(3)]';  %calculo de k3
   %para k4 cambiamos ambos argumentos
  yp=yy+h*k3;
  k4=[-1/40*yp(1)+3/500*yp(2),0.01*yp(1)-1/125*yp(2)+1/200*yp(3),1/500*yp(2)-3/400*yp(3)]';
    yy=yy+h/6*(k1+2*k2+2*k3+k4);
  y(:,n+1)=yy;
end

 
%Dibujo la solucion aproximada
figure(1)
clf
plot(t,y)
xlabel('Tiempo(días)')
ylabel('Cantidad soluto (kg)')
title('Metodo Runge Kutta 4 con h=0.1 y N=50')
legend('Pantano A','Pantano B','Pantano C')
}}

[[Archivo:Grafico_RK4_h01_n50.jpg|400px|thumb|left|Gráfico para h=0.1 y n=50]]
[[Archivo:Grafico_RK4_h1_n50.jpg|400px|thumb|center|Gráfico para h=1 y n=50]]
[[Archivo:Grafico_RK4_h01_n100.jpg|400px|thumb|left|Gráfico para h=0.1 y n=100]]
[[Archivo:Grafico_RK4_h1_n100.jpg|400px|thumb|center|Gráfico para h=1 y n=100]]
[[Archivo:Grafico_RK4_h01_n300.jpg|400px|thumb|left|Gráfico para h=0.1 y n=300]]
[[Archivo:Grafico_Trapecio_h1_n300.jpg|400px|thumb|center|Gráfico para h=1 y n=300]]


_
=Conclusiones=

Es razonable pensar que el pantano C, e incluso el pantano A vean aumentadas su cantidad de soluto durante los primeros días. Esto se debe a la considerable diferencia de contaminante entre el pantano B y los demás. Con ayuda de las gráficas se ve claramente como la cantidad de soluto del pantano C aumenta durante los primeros 200 días por culpa del agua contaminada que recibe del pantano B hasta llegar a un máximo aproximado de 250 Kg y como transcurrido un tiempo comienza a depurarse y recibir agua en mejor estado disminuyendo el soluto, hasta eliminar por completo su contaminante. 
Las gráficas tambien nos aclaran la duda que puede surgir con el pantano A, el cual no aumenta su cantidad de soluto y disminuye desde el primer día, por tanto su máximo se da en t=0, 556 Kg.

Aunque no se aprecie bien en las gráficas por su corto rango de tiempo, todas ellas tienden a 0. Por tanto cuanto mas tiempo pase, menos cantidad de soluto habrá

Asimismo, utilizando el comando 'find' matlab es capaz de calcular el numero de columnas y, por tanto, el numero de días para los que, con h=1, el contaminante desciende de 0,01Kg, hasta niveles despreciables.

- Para el primer pantano:

find(x(1,:)>0.01)

Dando un valor de t=3479 días.

- Para el segundo:

find(x(2,:)>0.01)

Dando un valor de t=3907 días.

- Y para el tercero:

find(x(3,:)>0.01)

Dando un valor de t=3641 días.

[[Categoría:Ecuaciones Diferentiales]]
 [[Categoría:ED13/14]]
 [[Categoría:Trabajos 2013-14]]