Diferencia entre revisiones de «Modelo de mezclas-Grupo18»

Revisión actual del 19:11 1 may 2016

El problema que se nos plantea es el siguiente:

Consideramos tres pantanos A, B y C . Los pantanos B y C tienen 200 Hm3 de agua cada uno mientas que el A tiene 400Hm3 . Por un accidente, en un momento dado empieza a entrar agua contaminada al pantano A a razón de 5Hm3/día . Las mediciones indican que el porcentaje de contaminante que está entrando es de 3kg/Hm3. Se decide bombear agua de una pantano hacia otro o hacia un río quedando la situación como sigue

•  A recibe de B 2Hm3/día y B de A 3Hm3/día.
•  B recibe de C 2Hm3/día y C de B 1Hm3/día.
•  C recibe de A 4Hm3/día y A no recibe agua de C.
• De B sale a un río 2Hm3/día y de C 3Hm3/día. 

Si xA(t), xB(t) y xC(t) designa las cantidades en kilos de contaminante en cada pantano cuando han pasado t días desde que que se empezó a bombear agua de unos pantanos hacia otros,

1 Encuentre el problema de valor inicial que satisfacen xA(t), xB(t) y xC(t). 


Para este problema hemos usado el modelo de mezclas,que se basa en la ley de conservación de la masa, y es el siguiente:

En general, si disolvemos Q kg (o cualquier unidad de masa) de un soluto en V litros (o alguna otra unidad de volumen) de un solvente, obtenemos una solución con una concentración C = Q (kg/V) de soluto (leído C kilogramos por litro y entendido como C kilogramos de soluto por cada litro de solución). Ahora supongamos que inicialmente (en t = 0)) tenemos en un tanque una cantidad V0 de solución donde hay disuelta una cantidad Q0 de un soluto. Supongamos también que, a partir de t = 0, se deja entrar otra solución al tanque con una rapidez Re (flujo de entrada) y con una concentración Ce (concentración de entrada) del mismo soluto y que, al mismo tiempo, se deja salir del tanque la nueva solución (considerada uniforme por mezclado) con una rapidez Rs (flujo de salida) y una concentración Cs (concentración de salida) del mismo soluto. Observamos:

 • Si Re = Rs, entonces la cantidad V0 de solución se mantiene constante al paso del tiempo.
 • Si Re≠ Rs, entonces la cantidad V de solución será función del tiempo t, V = V(t). Aún más, Si Re > Rs, entonces V(t) > V0 (ademas es creciente);mientras que si Re < Rs, entonces V (t) > V0 (ademas es decreciente).
 • En general la cantidad Q de soluto es el tanque será función del tiempo, Q = Q(t).
 • Igualmente, la concentración C del soluto en el tanque será función del tiempo t; y variará  según que Re = Rs o que Re ̸= Rs, esto es,

 • Un problema que es de interés en esta clase de procesos consiste en determinar la cantidad Q(t) de soluto en el tanque en cualquier instante t ≥ 0. Para resolver este problema procedemos como sigue:

Consideramos primero la rapidez con que cambia la cantidad de soluto Q(t) en el tanque, la cual está dada por la rapidez con que entra el soluto menos la rapidez con la que sale

Q′(t) = (rapidez con que entra el soluto) − (rapidez con que sale el soluto)

Si tenemos en cuenta que la rapidez con que entra el soluto es Re · Ce y la rapidez con que sale el soluto es Rs · Cs, el modelo de este proceso queda

􏰌 Sabiendo esto, el problema a resolver queda como sigue:

 1. x’A (t) =  15 - 7/400* xA (t)  + 2/200 * xB (t)
 2. x’B (t) =  3/400*xA (t)  -5/200* xB (t) + 2/200* xC (t)
 3. x’C (t) =  4/400*xA (t)  + 1/200*xB (t) -5/200*  xC (t)
  
  xA(0), xB(0) , xC(0)= 0; t= [0,T]

2 Resuelva numéricamente el problema de valor inicial en los primeros 60 días desde e el bombeo de agua, con el método del trapecio y con h = 0.1 y h = 1. 


Usamos un programa genérico (a excepción de las matrices A y T) y en la ventana de comandos ejecutamos todas las condiciones

ne=input('numero de ecuaciones');%----2
t0=input('introduce instante inicial');%----0
tN=input('introduce instante final');%-----60
y0=input('inttroduce vector columna valor inicial');%-----[0;0;0]
funcion=input('introduce las funciones de cada ecuacion coom vector columna','s');%-----[15-7/400*y(1)+ 0.01*y(2);3/400*y(1)-5/200*y(2)+1/100*y(3); 1/100*y(1)+1/200*y(2)-5/200*y(3)]
f=inline(funcion,'t','y');
h=input('introduce tamaÒo de paso:')%;-----0.1 ó 1
N=round((tN-t0)/h);
t=linspace(t0,tN,N+1);
y=zeros(ne,N+1);%euler
y(:,1)=y0;
ytr=zeros(ne,N+1);%trapecio
ytr(:,1)=y0;
I=eye(ne);
 A=[-7/400 0.01 0;3/400 -5/200 1/100; 1/100 1/200 -5/200];
%A=[0 1;-1 4*t(k+1)];%si A depende de t pongo ese termino como t(k+1) y lo meto dentro del bucle
for k=1:N
   %y(:,k+1)=y(:,k)+h*f(t(k),y(:,k));%euler
   T=[15;0;0];
    ytr(:,k+1)=(I-h/2*A)\(ytr(:,k)+h/2*f(t(k),ytr(:,k))+h/2*T);
end
%PREGUNTA 3
tt=input('Introduce el tiempo en el que queremos conocer la cantidad de contaminante: ');
J=round((tt-t0)/h);
ytr(:,J+1)

%PREGUNTA 4
a=find(ytr(2,:)>400);
diapantano2=(a(1)+1)*h

b=find(ytr(3,:)>400);
diapantano3=(b(1)+1)*h

%PREGUNTA 5
c=find(ytr(1,:)>1000);
diapantano1=(c(1)+1)*h

%PREGUNTA 6
if ytr(2,601)>ytr(3,601)
    disp('El pantano B tiene más contaminante a los 600 días');
end

if ytr(3,601)>ytr(2,601)
    disp('El pantano C tiene más contaminante a los 600 días');
end
hold on
plot(t,ytr(1,:),'b')
plot(t,ytr(2,:),'r')
plot(t,ytr(3,:),'g')
legend('PANTANO A','PANTANO B','PANTANO C');
hold off

El programa es el mismo para h=0,1 y h=1

Cantidad de contaminante a los 60 días con tamaño de paso h=0.1

Tamaño de paso h=1

La gráfica refleja la situación planteada en el sistema de ecuaciones. Al pantano A es al que le entra una mayor cantidad de contaminante, mientras que en los pantanos B y C entra una cantidad menor. Pero a su vez, el aumento de la cantidad en A hace que también aumente las cantidades en B y c (los pantanos están conectados). Esto hace que cambie la pendiente de las curvas en un cierto punto. Como B recibe menos de A que C, que es el pantano con más cantidad de contaminante, conforme pasa el tiempo aumenta la cantidad de C respecto de B.

3 Obtenga una aproximación de xA(t), xB(t) y xC(t) cuando han pasado justo 60 días del bombeo. 


Si vemos en el programa anterior el vector YTR, el último termino correspondiente a cada fila será el siguiente:

     xA (60)= 577.6 Kg  
     xB (60)= 133.3 Kg
     xC (60)= 114.9 Kg

4 A partir de cuanto tiempo podemos asegurar que los pantanos B y C tienen más de 400 kg. de contaminannte. 


Metemos en el programa un tiempo largo para asegurarnos que llegamos a los 400 Kg. Una vez ejecutado, buscamos en la variable ytr (una fila por pantano) el valor 400 y esa columna nos dirá en que ‘paso’ se ha ejecutado. Si vamos al vector t (tiempo) en dicho paso, nos da el tiempo.

También podemos hacerlo como postpreoceso en el propio programa.

Pantano B: 184 días

Pantano C: 174 días

5 A partir de cuanto tiempo podemos asegurar que el pantano A tiene más de 1 tonelada de contaminante. 


Aplicamos el mismo criterio que antes pero buscando el valor 1000 en la primera fila Pantano A: 192 días

6 A lo largo de los primeros 600 días desde el bombeo, ¿qué pantano tiene más contaminación, el B o el C?. 


En un gráfico correspondiente a 600 días se observa que tiene más contaminación el C, aunque por poco. La gráfica se amplía para observar la diferencia

Grafica ampliada de la diferencia entre B y C

Si lo hacemos con el postproceso, lo podemos hacer con el comando disp e if, que nos verifica lo que nos ha mostrado la gráfica

7 ¿A qué tienden las concentraciones de contaminación en los pantanos A,B y C para tiempos largos?

En este caso se mete un tiempo muy grande, lo suficiente para que se estabilicen las cantidades. Hemos tomado un tiempo de 3000 días.

Pantano A: 1200 Kg

Pantano B: 600 Kg

Pantano C: 600 Kg

Cantidades estables de contaminante en los pantanos

8 Segunda situación

 Suponemos ahora una situación distinta a la tratada. A partir de día 60, se ha logrado evitar que entre la contaminación al pantano A, pero con el fin de activar un plan de limpieza en los pantanos, sigue la entrada de agua al pantano A sin contaminar a la misma velocidad, y el mismo bombeo de agua que en los primeros 60 días. Si seguimos llamando , por comodidad, xA(t), xB(t) y xC(t) las cantidades en kilos de contaminante en cada pantano cuando han pasado t días desde que empezó a entrar agua no contaminada,

8.1 Calcule xA(0), xB(0) y xC(0). 


Las condiciones iniciales serán las del primer caso cuando han pasado 60 días

 xA(0)= 577.6 Kg  
 xB(0)= 133.3 Kg
 xC(0)= 114.9 Kg

8.2 Problema de valor inicial que satisface xA(t), xB(t) y xC(t). 


Es el mismo problema que antes pero quitando el contaminante que entra en A:

         x’A (t) =  - 7/400* xA (t)  + 2/200 * xB (t)
         x’B (t) =  3/400*xA (t)  -5/200* xB (t) + 2/200* xC (t)
         x’C (t) =  4/400*xA (t)  + 1/200*xB (t) -5/200*  xC (t)
      
     xA(0)= 577.6 Kg  ; xB(0)= 133.3 Kg;  xC(0)= 114.9 Kg ; t= [0,T]

8.3 Resuelva numéricamente el problema de valor inicial utilizando el método de Runge-Kuta y la longitud de paso h = 1 en


(a) primer año desde que dejó de entrar contaminante a los pantanos,


(b) dos primeros años desde que dejó de entrar contaminante a los pantanos,

(c) tres primeros años desde que dejó de entrar contaminante a los pantanos. 


Lo único que cambia es el tiempo que metemos, pero como lo ejecutamos en la ventana de comandos, el programa es igual en los 3 casos

clear, close all
ne=input('numero de ecuaciones');%-------3
t0=input('Introduzca instante inicial del problema');%--------0
tN=input('Introduzca instante final del problema');%-------365 ó 730 ó 1095 
y0=input('Introduzca valor inicial');---------[577.6; 133.3 ; 114.9 ]
fun=input('Introduce forma analítica de la función','s');%-------[15-7/400*y(1)+ 0.01*y(2);3/400*y(1)-5/200*y(2)+1/100*y(3); 1/100*y(1)+1/200*y(2)-5/200*y(3)]
f=inline(fun,'t','y');
%N=input('Introduce número de subintervalos');
%h=(tN-t0)/N;
%o tambien se puede
h=input('Introduce tamaño de paso');--------1
N=round((tN-t0)/h);
t=linspace(t0,tN,N+1);
y=zeros(ne,N+1);%euler
y(:,1)=y0;
for i=1:N
 k1=f(t(i),y(:,i));
 k2=f(t(i)+(1/2)*h,y(:,i)+(1/2)*k1*h);
 k3=f(t(i)+(1/2)*h,y(:,i)+(1/2)*k2*h);
 k4=f(t(i)+h,y(:,i)+(1/2)*k3*h);
 y(:,i+1)=y(:,i)+(h/6)*(k1+2*k2+2*k3+k4);
end


hold on
plot(t,y(1,:),'b')
plot(t,y(2,:),'r')
plot(t,y(3,:),'g')
legend('PANTANO A','PANTANO B','PANTANO C');
hold off

%PREGUNTA 5
%Apartado1
MaxPantano2=max(y(2,:))
MaxPantano3=max(y(3,:))
%Apartado3
a=find(y(1,:)<4);
b=find(y(2,:)<2);
c=find(y(3,:)<2);
lim=[a(1),b(1),c(1)];
minimo=disp((max(lim)+1)*h)


hold on
plot(t,y(1,:),'b')
plot(t,y(2,:),'r')
plot(t,y(3,:),'g')
legend('PANTANO A','PANTANO B','PANTANO C');
hold off

Cantidad de contaminante después de un año

Cantidad de contaminante después de dos años

Cantidad de contaminante después de tres años ]

8.4 ¿Puede resolver numéricamente en [0,T∞] cuando T∞ es tan grande como uno quiera?. ¿Hay alguna relación con el tamaño de h > 0?.

Se puede resolver aunque el programa MATLAB tardará más en ejecutar el programa cuanto mayor sea T. Asimismo, cuanto mayor sea h, menos puntos tendrá que analizar el programa, con lo que el tardará menos en ejecutarse. En cuanto al error, al no haber prácticamente pendiente (el final de las gráficas se puede asociar a rectas horizontales), si aumentamos mucho T, el error no va a depender de h.

8.5 Parece razonable de que xB(t) y xC(t) empezarán siendo creciente, para luego pasar a ser decrecientes, mientras que xA(t) siempre va a ser decreciente.

• Dé un valor aproximado del valor máximo de xB(t) y xC(t). ¿Existe alguna manera razonable de medir, o tener una idea, del error que estamos cometiendo?. 


Obtenemos los máximos en el postproceso con MATLAB y los resultados son:

Cantidad máxima de contaminante en el pantano B: 164.2075 Kg

Cantidad máxima de contaminante en el pantano C: 174.9483 Kg

• Si t es muy grande, ¿a qué valores se aproximan xA(t), xB(t) y xC(t)?. 
 Tiende a 0, que es el objetivo

• Las autoridades considerarán que ha desaparecido la contaminación de los pantanos, cuando 
las cantidades xA (t), xB (t) y xC (t) sean menores que 0.01 Kg/Hm3.

¿Puede dar una estimación del número de días a partir de los cuales, podemos pensar que ha desaparecido la contaminación?

Para este caso también metemos un tiempo grande, para asegurarnos que pasamos dicho valor. Aplicamos el mismo criterio que para el apartado anterior, mediante postproceso con el comando 'find'

El resultado es que transcurren 642 días