Diferencia entre revisiones de «Modelos de mezclas (grupo 5)»

Revisión actual del 19:23 2 feb 2014

1 Introducción

Tenemos dos pantanos A y B, con 100 [math]Hm^3[/math] de agua cada uno, que están unidos por una presa que deja pasar agua de A a B. El pantano A recibe 3 [math]Hm^3/día[/math] de agua limpia proveniente de ríos y el B 1.5 [math]Hm^3/día[/math]. Para mantener el nivel de los pantanos estable la presa de A a B deja pasar una media de 3 [math]Hm^3/día[/math] mientras que la presa al final de B desaloja 4.5 [math]Hm^3/día[/math]. Se produce un vertido tóxico en el pantano A que deja 20 toneladas de un cierto contaminante. Supongamos que se dan las siguientes hipótesis::

• El contaminante está disuelto de forma homogénea en el agua de los pantanos.
• Al entrar o salir agua en un pantano, esta se mezcla con el agua del pantano de forma inmediata creando una mezcla homogénea.
• La variación de contaminante en un lago es la diferencia entre el contaminante que entra y sale en el lago, es decir si denotamos por [math] x_A (t) [/math] la cantidad de contaminante en el lago A, se tiene:: [math] \frac{dx_A}{dt} = velocidad de entrada – velocidad de salida[/math]

2 Planteamiento

2.1 Sin plan de limpieza

Sean [math] x_A [/math] y [math] x_B [/math] las cantidades de contaminante en cada uno de los lagos, la cantidad de contaminante en un instante de tiempo determinado se obtiene como diferencia entre la velocidad de entrada y salida de agua contaminada. En el lago A sólo entra agua limpia mientras que la velocidad de salida se obtiene como el producto de la concentración de contaminante en el lago A [math] \frac{X_A}{100} [/math] por el caudal de agua que sale 3 [math]^3/día [/math]. Así mismo, en el lago B la velocidad de entrada de agua contaminada es la misma que la velocidad de salida de agua contaminada del lago A, ya que el lago A vierte el agua contaminada en B, y la velocidad de salida del agua de B es el producto de su concentración [math] \frac{X_B}{100} [/math] por el caudal de agua que sale 4.5 [math] Hm^3/día [/math].

Para modelizar el problema, se plantea el siguiente sistema dos ecuaciones diferenciales de primer orden:

[math] \left\{\begin{matrix}\frac{\mathrm{d}X_A}{\mathrm{d} t}=-\frac{3X_A}{100}\\\frac{\mathrm{d} X_B}{\mathrm{d} t}=\frac{3X_A}{100}-\frac{4.5X_B}{100}\\X_A(t_{0})=20, \quad X_B(t_{0})=0\end{matrix}\right. [/math]

2.2 Con plan de limpieza

A continuación se propone activar un plan de limpieza basado en bombear 1[math] Hm^3/día [/math] de agua del pantano B al A, por lo que para mantener el volumen de agua entre los lagos constante se tiene que modificar el flujo entre ambos.

De este modo, si al lago A le llegan 3[math]Hm^3/día[/math] provenientes de los ríos y 1[math]Hm^3/día[/math] proveniente de B, el lago A deberá verter a B 4[math]Hm^3/día [/math] en lugar de los '3 [math] Hm^3/día [/math] que estaba vertiendo antes de activar el plan de limpieza. El lago B recibe 1.5 [math] Hm^3/día [/math] provenientes de los ríos y, como se ha visto, 4 [math] Hm^3/día [/math] provenientes del lago A. Tras activar el plan de limpieza bombea 1 [math] Hm^3/día [/math] a A, por lo que el vertido que desaloja B aguas abajo se mantendrá constante igual a 4.5 [math] Hm^3/día [/math] para mantener el volumen de agua. Así se plantea un sistema de ecuaciones diferenciales similar al anterior:: [math] \left\{\begin{matrix}\frac{\mathrm{d}X_A}{\mathrm{d}t}=-\frac{4X_A}{100}+\frac{1X_B}{100}\\\frac{\mathrm{d}X_B}{\mathrm{d} t}=\frac{4X_A}{100}-\frac{5.5X_B}{100}\\X_A(t_{0})=20, \quad X_B(t_{0})=0\end{matrix}\right. [/math]

2.3 Unión de un tercer pantano

Si hubiese un tercer pantano, al que se le llamará C, unido al anterior que desalojase 6 [math]Hm^3/día [/math] y recibiese 1.5 [math]Hm^3/día [/math], la velocidad de entrada de agua sería igual a la de salida de B, es decir, 4.5 [math] Hm^3/día [/math] ,mientras que la de salida serían estos 4.5 [math]Hm^3/día [/math] más 1.5 [math] Hm^3/día [/math] provenientes de los ríos, dando lugar a los 6 [math]Hm^3/día [/math] para que se mantenga constante el volumen. Así, el sistema de ecuaciones diferenciales quedará de la siguiente forma::

[math] \left\{\begin{matrix}\frac{\mathrm{d}X_A}{\mathrm{d}t}=\frac{-3X_A}{100}\\\frac{\mathrm{d}X_B}{\mathrm{d} t}=\frac{3X_A}{100}-\frac{4.5X_B}{100}\\\frac{\mathrm{d}X_C}{\mathrm{d} t}=\frac{4.5X_B}{100}-\frac{6X_C}{100}\\X_A(t_{0})=20, \quad X_B(t_{0})=0\end{matrix}\right. [/math]

3 Resolución numérica

A partir de este momento nos olvidamos del pantano C. Para resolver estos sistemas se utilizará un programade análisis numérico (en este caso Matlab®).

3.1 Método de Euler

Se utilizan los siguientes códigos::

Sin plan de limpieza:

clear all

N=100000; %número de iteraciones
t0=0;tN=325; %tiempo inicial y final
h=(tN-t0)/N; %paso: espacio entre dos iteraciones
x0=[20 0]'; %condiciones iniciales
A=[-0.03 0;0.03 -0.045];
x=x0;
xA(1)=x(1);
xB(1)=x(2);

for n=1:N    
  x=x+h*A*x;    
  xA(n+1)=x(1);    
  xB(n+1)=x(2);
end
y=t0:h:tN;

%dibujamos la gráfica
hold on
plot(y,xA,'b') %contaminante en pantano A: azul
plot(y,xB,'r') %contaminante en pantano B: rojo
title('Sin plan de limpieza')
xlabel('Tiempo (días)')
ylabel('Contaminante (T)')
hold off

Con plan de limpieza:

clear all

N=100000; %número de iteraciones
t0=0;tN=325; %tiempo inicial y final
h=(tN-t0)/N; %paso: espacio entre dos iteraciones
x0=[20 0]'; %condiciones iniciales
B=[-0.04 0.01;0.04 -0.055];
X=x0;
Xa(1)=X(1);
Xb(1)=X(2);

for n=1:N    
  X=X+h*B*X;    
  Xa(n+1)=X(1);    
  Xb(n+1)=X(2);
end
y=t0:h:tN;

%dibujamos la gráfica
hold on
plot(y,Xa,'g') %contaminante en pantano A: verde
plot(y,Xb,'m') %contaminante en pantano B: magenta
title('Con plan de limpieza')
xlabel('Tiempo (días)')
ylabel('Contaminante (T)')
hold off

Si queremos ver la diferencia de tiempo que tarda en desaparecer la mitad de contaminante en el pantano A activando el plan de limpieza o sin activarlo no tenemos más que observar las gráficas.

Trazando una horizontal desde 10 toneladas cortará a la gráfica sin limpieza en el tiempo [math] t_{mitad1}=23 días [/math] y a la gráfica con limpieza en [math] t_{mitad2}=18 días [/math], por lo que la diferencia de tiempos será [math] t_{mitad1}- t_{mitad2}=23-18=5 días[/math].

Si hacemos lo mismo para ver cuánto tarda en desparecer la tercera parte [math] t_{tercera1}-t_{tercera2}=13-10=3 días[/math].

3.2 Método de Runge-Kutta

Aplicando el método de Runge-Kutta de cuarto orden se utilizan el siguiente código:

clear all

N=10000; %número de interacciones
t0=0; tN=325;  %tiempo inicial y final
x0=[20 0]'; %condiciones iniciales
h=(tN-t0)/N; %paso: espacio entre dos iteraciones


%sin plan de limpieza
A=[-0.03 0;0.03 -0.045];
x=x0;
xA(1)=x(1);
xB(1)=x(2);
for n=1:N     
   k1=A*x;    
   k2=A*(x+(1/2*k1*h));    
   k3=A*(x+(1/2*k2*h));    
   k4=A*(x+(k3*h));    
   x=x+h/6*(k1+2*k2+2*k3+k4);    
   xA(n+1)=x(1);    
   xB(n+1)=x(2);
end


%Con plan de limpieza
B=[-0.04 0.01;0.04 -0.055];
X=x0;
Xa(1)=X(1);
Xb(1)=X(2);
for n=1:N     
   k1=B*X;    
   k2=B*(X+(1/2*k1*h));    
   k3=B*(X+(1/2*k2*h));    
   k4=B*(X+(k3*h));    
   X=X+h/6*(k1+2*k2+2*k3+k4);    
   Xa(n+1)=X(1);    
   Xb(n+1)=X(2);
end


y=t0:h:tN;


%dibujamos las dos gráficas
hold on
figure(1) %sin plan de limpieza
plot(y,xA,'g'); %contaminante pantano A: verde
plot(y,xB,'m'); %contaminante pantano B: magenta
title('Sin plan de limpieza')
xlabel('Tiempo (dias)')
ylabel('Contaminante (T)')
figure(2) %con plan de limpieza
plot(y,Xa,'b') %contaminante pantano A: azul
plot(y,Xb,'r') %contaminante pantano B: rojo
title('Con plan de limpieza')
xlabel('Tiempo (dias)')
ylabel('Contaminante (T)')
hold off

Ambas soluciones quedan reflejadas en las siguientes gráficas:

3.3 Comparación método Euler con método Runge-Kutta

Para comparar ambos métodos basta con comparar sus gráficas. Realizamos el procedimiento con ambos sistemas. La gráfica del primer sistema es la siguiente:

Como puede verse la diferencia entre ambos métodos en tan pequeña que las gráficas aparecen superpuestas. Para poder apreciar la diferencia entre ellas ampliamos la grafica cogiendo dos intervalos pequeños de tiempo. Runge-Kutta aparece representado en azul y Euler en verde.

Realizamos el mismo proceso con el segundo sistema y como es de esperar sucede lo mismo:

Una vez comparados ambos métodos se puede afirmar que la diferencia entre ellos es casi inapreciable.

4 Contaminante inicial desconocido

No conocemos la cantidad de contaminante inicial pero sabemos que tras unos días se redujo el contaminante la sólo una tonelada en A y dos en B. Si queremos estimar el contaminante que se vertió inicialmente lo podemos hacer a través del método de Euler con el siguiente código:

clear all

t0 = 0;
tN = -600;
N = 10000;
h = (tN-t0)/N;
x0 = [1;2];
A = [-0.03 0;0.03 -0.045];
x = x0;
xA(10001) = x(1);
xB(10001) = x(2);
for i=N+1:-1:2
    x = x+h*A*x;
    if x(2)<=0
        xA(1:i-1) = [];
        xB(1:i-1) = [];
        r = length(xA);
        break
    else
        r = N+1;
    end
    xA(i-1) = x(1);
    xB(i-1) = x(2);
end
t = tN:-h:t0;
tr = length(t);
if tr~=r
    t(1:tr-r) = [];
end
hold on
xlabel('Tiempo (días)')
ylabel('Contaminante (T)')
plot(t,xA,'b');
plot(t,xB,'r');
hold off

Obtendremos el vertido inicial en el mismo instante en el que el contaminante en el pantano B se anule. Comopuede verse en la gráfica el contaminante en B nunca se hace cero ya que crece exponencialmente hacia la izquierda hasta un tiempo [math] t=-\infty [/math]. Este resultado no es lógico, por lo que modificamos mínimamente los datos que nos dan sobre la cantidad de contaminante en ambos pantanos tras unos días para ver si podemos conseguir un resultado aproximado que sea lógico.

Cambiamos la cantidad de contaminante en B de 2 toneladas a 1.8:

Cantidad de contaminante en B=1.999 toneladas:

Como puede comprobarse según aumentamos la cantidad de contaminante en B la cantidad de contaminante inicialse va incrementando hasta que, utilizando la condición dada de 2 toneladas de contaminante en B, tiende a [math] t=+\infty [/math]. Por lo tanto no se puede concretar ninguna cantidad inicial.

5 Autores

• Raquel Gómez Vázquez, 187
• Eva Périz Rodríguez, 334
• Daniel Prieto Alonso, 341
• Adrián Salas Calvo, 385