Diferencia entre revisiones de «Modelos de mezclas(Grupo 15)»

Revisión actual del 00:28 3 jun 2013

Dos pantanos A y B con 100 [math]Hm^3[/math] de agua cada uno están unidos por una presa que deja pasar agua de A a B. El pantano A recibe 3 [math]Hm^3/dia[/math] de agua limpia provniente de rios y el B 1.5 [math]Hm^3/dia[/math]. Para manatener estable la presa, de A a B se deja pasar una media de 3 [math]Hm^3/dia[/math] mientras que la presa al final de B desaloja 4.5 [math]Hm^3/dia[/math]. Se produce un vertido en el pantano A que deja 20 toneladas de un cierto contaminante. Suponemos las siguientes hipotesis:

• El contaminante está disuelto de forma homogénea en el agua de los pantanos.
• Al entrar o salir agua en un pantano, está se mezcla con el agua del pantano de forma inmediata creando un mezcla homogénea.
• La variación de contaminante en un lago es la diferencia ente el contaminate que entra y el que sale del lago.

1 PRIMER SISTEMA DE ECUACIONES

Llamamos [math]X_A(t)[/math] a la cantidad de contaminante del lago A y [math]X_B(t)[/math] a la cantidad de contaminante del lago B. Definimos la variación de contaminante en un lago como:

[math]X'_A(t)=velocidad(entrada)-velocidad(salida)[/math]

Conocemos la velocidad de entrada y salida de volumen de agua de cada pantano, está multiplicada por la concentración de contaminante en el instante t nos da la velocidad de entrada y salida de contaminante. Definimos como concentración de contaminante, siendo vol(t) el volumen del lago en un instante t:

[math]C(t)={X(t)\over vol(t)}[/math]

Aplicando las hipotesis y sabiendo que el contaminante se traspasa de lago a lago y que el agua de los rios es agua limpia, nos quedan las siguiente ecuaciones: \begin{array}{c} X'_A(t)=-{3\over 100}X_A \\ X'_B(t)={3\over 100}X_A-{4.5\over 100}X_B \end{array}

En forma matricial y considerando que [math]X_A(0)=20[/math] nos queda el siguiente problema de Cauchy:

[math]\begin{pmatrix} X'_A\\ X'_B\end{pmatrix}={1\over 100}\begin{pmatrix} -3 & 0 \\ 3 & -4.5 \end{pmatrix}\begin{pmatrix} X_A \\ X_B \end{pmatrix}[/math]

[math]\begin{pmatrix} X_A(0) \\ X_B(0) \end{pmatrix}=\begin{pmatrix} 20 \\ 0 \end{pmatrix}[/math]

• La solución general de este sistema es inmediata, por variables separadas se obtiene [math]X_A[/math] y sustituyendo en la otra ecuación se obtiene [math]X_B[/math] :

[math] X_A(t)=20e^{-3t/100}[/math]

[math] X_B(t)=40(e^{^-3t/100}-e^{-4.5t/100})[/math]

• Ahora,suponemos que hay un tercer pantano C unido a B por una segunda presa que suelta 6 [math]Hm^3/dia[/math], reciviendo 1.5 [math]Hm^3/dia[/math] de agua limpia de rios.

El sistema pasaría a depender también de la velocidad en C, ahora como no conocemos el volumen de C consideramos VolA=VolB=VolC= 100 [math] Hm^3 [/math]=cte, tampoco sabemos la velocidad de C, pero si se tiene que mantener el volumen constante tendrá que recivir 4.5 [math] Hm^3 [/math] de B (que es lo que suelta) y él mismo soltará 6 [math] Hm^3 [/math], quedando el sistema de la siguiente forma:

[math] \left\{\begin{matrix} X'_A=\frac{-3X_A}{100}\\ X'_B=\frac{3X_A}{100}-\frac{4.5X_B}{100}\\ X'_C=\frac{4.5X_B}{100}-\frac{6X_C}{100}\\X_A(t_{0})=20, \quad X_B(t_{0})=0\end{matrix}\right. [/math]

2 SEGUNDO SISTEMA DE ECUACIONES

Supongamos que se activa un plan de limpieza que consiste en bombear 1 [math]Hm^3/dia[/math] de agua del pantano B al A ajustando las cantidades de agua que dejan pasar las presas para mantener estables los niveles de los pantanos. El sistema de ecuaciones va a cambiar, ahora el pantano A recivirá agua contaminada de B y para que se mantenga estable, dejara pasar más agua a B reciviendo este más agua contaminada. Los cambios se pueden ver en las siguientes ecuaciones:

[math] X'_A(t)=-{4\over 100}X_A+{1\over 100}X_B [/math]

[math] X'_B(t)={4\over 100}X_A-{5.5\over 100}X_B [/math]

En forma matricial tenemos el siguiente problema de Cauchy:

[math] \begin{pmatrix} X'_A \\ X'_B \end{pmatrix}={1\over 100}\begin{pmatrix} -4 & 1 \\ 4 & -5.5 \end{pmatrix}\begin{pmatrix} X_A \\ X_B \end{pmatrix} [/math]

[math] \begin{pmatrix} X_A(0) \\ X_B(0) \end{pmatrix}=\begin{pmatrix} 20 \\ 0 \end{pmatrix} [/math]

Para solucionarlo es útil pasarlo a una ecuación diferencial de segundo orden, siendo esta de coeficientes constantes y homogenea, resolviendola sacando soluciones con el polinomio caracteristico, considerando que ambas forman un espacio vectorial de dimension 2 y sustituyendo las condiciones iniciales obtenemos la solución general, una vez obtenida sustituimos en el sistema inicial y sacamos la concentración del otro pantano.

[math] \begin{matrix}X''_A+{19\over 200}X'_A+{9\over 5000}X_A=0 \\ X_A(0)=20 \\ X'_A(0)=-{4\over 100}X_A(0)+{1\over 100}X_B(0)=-{4\over 5}\end{matrix}[/math]

• Las soluciones son:

[math] X_A(t)=(10/\sqrt{73})e^{-19t/400}\left [ (\sqrt{73}+35)e^{\sqrt{73}t/400}+(\sqrt{73}-35)e^{-\sqrt{73}t/400}\right ] [/math]

[math] X_B(t)=(320/\sqrt{73})e^{-19t/400}sinh(\sqrt{73}/400)[/math]

3 RESOLUCIÓN MEDIANTE EL MÉTODO DE EULER

El método de Euler es un método númerico explicito de primer orden.Resolvemos los dos sistemas con MATLAB mediante el método de Euler y damos una aproximación a las soluciones generales obtenidas en los aparatados anteriores.

PRIMER SISTEMA

clear all
format long
% PRIMER SISTEMA
t0=0; tN=6*100/4.5; % Extremos del intervalo temporal a estudiar
N=100; h=(tN-t0)/N; % Número y amplitud de subintervalos
t=t0:h:tN; % Valores de t
x0=20;y0=0; % Asignacion de valores iniciales x->XA y->XB
x(1)=x0;y(1)=y0; % Inicialización de los vectores
%Algoritmo
for n=1:N
    x(n+1)=x(n)+h*(-3/100)*x(n); %Aplicación del método de Euler 
    y(n+1)=y(n)+h*(3/100)*x(n)-h*(4.5/100)*y(n); 
end
% Dibujo de gráficas
plot(t,x,'b',t,y,'r')

SEGUNDO SISTEMA

clear all
format long
% SEGUNDO SISTEMA
t0=0; tN=6*100/4.5; % Extremos del intervalo temporal a estudiar
N=100; h=(tN-t0)/N; % Número y amplitud de subintervalos
t=t0:h:tN; % Valores de t
x0=20;y0=0; % Asignacion de valores iniciales x->XA y->XB
x(1)=x0;y(1)=y0; % Inicialización de los vectores
%Algoritmo
for n=1:N
    x(n+1)=x(n)+h*(-4/100)*x(n)+h*(1/100)*y(n); %Aplicación del método de Euler 
    y(n+1)=y(n)+h*(4/100)*x(n)-h*(5.5/100)*y(n); 
end
% Dibujo de gráficas
plot(t,x,'b',t,y,'r')

Gráficamente, las soluciones son:

centro

centro

Se observa que ambas soluciones son del tipo exponencial y se ve que para un tiempo grande la contaminación acaba desapareiendo. Tambien se observa que ambos pantanos llegan a tener la misma contaminación y que para el sistema 2, en el que se activa un plan de limpieza, la grafica del pantano B se desplaza hacia la izquierda lo que indica que tardaría menos en desaparecer toda la contaminación.

• La diferencia de tiempo que tarda en desaparecer la mitad del contaminante inicial (20 toneladas) en el pantano A cuando se activa el sistema de limpieza será la diferencia de valores medidos en ambas graficas, es decir, t2-t1=23-18=5 dias y para la tercera parte es t4-t3=13-10=3 dias.

4 RESOLUCIÓN MEDIANTE EL MÉTODO DE RUNGE-KUTTA

El método de Runge-Kutta es un método implicito de cuarto orden con lo que nos dára una mejor aproximación que el método de Euler al tener mayor orden, pero al ser implicito la programación será más difícl. Resolvemos los dos sistemas en MATLAB mediante el método de Runge-Kutta:

PRIMER SISTEMA

clear all
format long
t0=0; tN=6*100/4.5; 
N=100; h=(tN-t0)/N; 
t=t0:h:tN; 
x0=20;y0=0;
x(1)=x0;y(1)=y0; 
for n=1:N
% Aplicación del método de Runge-Kutta
    p1=(-3/100)*x(n);
    q1=(3/100)*x(n)-(4.5/100)*y(n);
    p2=(-3/100)*(x(n)+h*p1/2);
    q2=(3/100)*(x(n)+h*p1/2)-(4.5/100)*(y(n)+h*q1/2);
    p3=(-3/100)*(x(n)+h*p2/2);
    q3=(3/100)*(x(n)+h*p2/2)-(4.5/100)*(y(n)+h*q2/2);
    p4=(-3/100)*(x(n)+h*p3);
    q4=(3/100)*(x(n)+h*p3)-(4.5/100)*(y(n)+h*q3);
    x(n+1)=x(n)+(h/6)*(p1+2*p2+2*p3+p4);
    y(n+1)=y(n)+(h/6)*(q1+2*q2+2*q3+q4);
end
plot(t,x,'b',t,y,'r')

SEGUNDO SISTEMA

clear all
format long
t0=0; tN=6*100/4.5; 
N=100; h=(tN-t0)/N; 
t=t0:h:tN; 
x0=20;y0=0;
x(1)=x0;y(1)=y0; 
for n=1:N
% Aplicación del método de Runge-Kutta
    p1=(-4/100)*x(n)+(1/100)*y(n);
    q1=(4/100)*x(n)-(5.5/100)*y(n);
    p2=(-4/100)*(x(n)+h*p1/2)+(1/100)*(y(n)+h*q1/2);
    q2=(4/100)*(x(n)+h*p1/2)-(5.5/100)*(y(n)+h*q1/2);
    p3=(-4/100)*(x(n)+h*p2/2)+(1/100)*(y(n)+h*q2/2);
    q3=(4/100)*(x(n)+h*p2/2)-(5.5/100)*(y(n)+h*q2/2);
    p4=(-4/100)*(x(n)+h*p3)+(1/100)*(y(n)+h*q3);
    q4=(4/100)*(x(n)+h*p3)-(5.5/100)*(y(n)+h*q3);
    x(n+1)=x(n)+(h/6)*(p1+2*p2+2*p3+p4);
    y(n+1)=y(n)+(h/6)*(q1+2*q2+2*q3+q4);
end
plot(t,x,'b',t,y,'r')

Gráficamente, las soluciones son:

centro

centro

• Para comparar el método de Runge-Kutta con el de Euler, se pueden comparar los gráficos de ambos métodos o bien como en nuestro caso, hemos cojido el 70 valor de la solucion x (por cojer uno) del Runge-Kutta y resulta que nos da el mismo valor que haciendolo con la solucion exacta(x=1,2658) y haciendo lo mismo con el de Euler nos da un valor diferente(x=1,1960) con lo que se verifica que el método de Runge-Kutta da una mejor aproximación de la solución.

5 CONTAMINANTE INICIAL DESCONOCIDO

Si no conocemos la cantidad de contaminante inicial pero sabemos tras unos dias se redujo el contaminante a sólo una tonelada en A y dos en B,¿Cuánto contaminante se estima que se vertió inicialmente? En este caso desconocemos el contaminante inicial,pero para el tiempo tN=0 y otro tiempo inicial por ejmplo t0=100, conocemos los valores del contaminante [math](X_A,X_B)=(1,2)[/math], con esto aplicando el método de Euler podemos estimar el contaminante inicial.

CÓDIGO EN MATLAB

t0=100; tN=0; 
x0=[1 2]'; %Valores iniciales para t0
N=100; h=(tN-t0)/N; %Número y amplitud de subintervalos
A=[-3/100,0;3/100,-4.5/100]; %Matriz de coeficientes del primer sistema
x=x0; 
xA(1)=x(1); %Iniciamos el vector cantidad de contaminante en A
xB(1)=x(2); %Iniciamos el vector cantidad de contaminante en B
for n=1:N
    x=x+h*(A*x); %Aplicamos el método de Euler
    xA(n+1)=x(1);
    xB(n+1)=x(2); 
end
x=t0:h:tN;
%Dibujo
plot(x,xA,'b',x,xB,'r')

Gráficamente la solución es:

centro

• De la gráfica se puede estimar la cantidad de contaminate inicial y sera XA(0) y XB(0).