Modelo de distribución de un contaminante en una red de pantanos

De MateWiki
Saltar a: navegación, buscar
Trabajo realizado por estudiantes
Título Modelo de distribución de un contaminante en una red de pantanos. Grupo 14
Asignatura Ecuaciones Diferenciales
Curso Curso 2013-14
Autores Manuel Campayo Fernández, Alejandro Arturo Serrano Leo, David Suárez Cuesta, Arturo Rodríguez Gárate, Jorge Martín Sebastián, Joaquín Sánchez Molina
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura


1 Situación inicial del problema

Disponemos de tres pantanos, a los que llamaremos a partir de ahora A, B y C. Estos pantanos poseen unas capacidades de 400, 200 y 200 [math]Hm^3[/math] respectivamente.

Un fallo en una tubería hace que esta comience a verter agua contaminada sobre el pantano A. Con motivo de examinar este fallo con mayor detenimiento, se llama a un grupo de ingenieros, que determina que la tubería está vertiendo agua contaminada en el pantano A a razón de 5 [math]Hm^3/día[/math]. Un análisis químico posterior muestra que la concentración de contaminante que está entrando es de 3 [math]Kg/Hm^3[/math].

Entonces, se toma la decisión de bombear agua de un pantano hacia otro o hacia un río, tal y como se decribe en la siguiente figura: Planteamientotrabajo1edos.jpg

2 Problema de valor inicial

En primer lugar, se definen $X_{A}(t)$, $X_{B}(t)$ y $X_{C}(t)$ como las cantidades (en kilogramos) de contaminante que se encuentran en cada pantano pasado un tiempo determinado t desde que comenzó el bombeo.

Para plantear el problema de valor inicial que modelice la variación de contaminante que existe en cada uno de los pantanos al producirse el movimiento de este entre los distintos pantanos y el río, se ha acordado obtener un sistema basado en la siguiente ecuación:: [math] \frac{dx_i}{dt} = velocidad(entrada) – velocidad(salida)[/math] Siendo i = Pantanos A, B y C

Las velocidades de entrada y de salida de contaminante se obtienen multiplicando el caudal que se mueve de un pantano a otro ([math]Hm^3/día[/math]) por la concentración de contaminante que posee dicho volumen ([math]X_i(t)/Volumen del pantano i[/math]). De esta manera, el problema de valor inicial buscado es el que se muestra a continuación:: [math] \left\{\begin{matrix} X'_A(t)=\frac{-7X_A(t)}{400}+\frac{X_B(t)}{100}+15 \\ X'_B(t)=\frac{3X_A(t)}{400}-\frac{X_B(t)}{40}+\frac{X_C(t)}{100}\\ X'_C(t)=\frac{X_A(t)}{100}+\frac{X_B(t)}{200}-\frac{X_C(t)}{40}\\X_A(t_{0})=0, \quad X_B(t_{0})=0, \quad X_C(t_{0})=0 \end{matrix}\right. [/math]

Que expresado de forma matricial queda:

\begin{equation*}

\begin{bmatrix}

X’_A(t)\\

X’_B(t)\\

X’_C(t)

\end{bmatrix}

=

\begin{bmatrix}

-7/400 & 1/100 & 0 \\

3/400 & -1/40 & 1/100 \\

1/100 & 1/200 & -1/40

\end{bmatrix}

\begin{bmatrix}

X_A \\

X_B \\

X_C

\end{bmatrix}

+

\begin{bmatrix}

15\\

0\\

0

\end{bmatrix}

\end{equation*}

3 Resolución númerica del problema

Para obtener la solución del problema de valor inicial que se acaba de plantear, se utiliza un método numérico, concretamente el del trapecio.

Se quiere obtener dicha solución en los primeros 60 días desde el bombeo, considerando un tamaño del intervalo (h) de 0,1 y 1.

3.1 Resolución numérica para h=0,1

El código utilizado para resolverlo es: 

t0=0; tn=60;
y0=[0,0,0]'; 
A=[-7/400,1/100,0;3/400,-1/40,1/100;1/100,1/200,-1/40];%Matriz de coeficientes
%Datos para la discretización
h=0.1;
N=(tn-t0)/h;
%Vectores de tiempos y de solución aproximada
t=t0:h:tn 
y=zeros(3,N+1);
%Inicializamos el bucle
y(:,1)=y0;
yy=y0;
for n=1:N
  g=[15, 0, 0]' 
  gn1=[15, 0, 0]'
  yy=(eye(3)-h/2*A)\((eye(3)+h/2*A)*yy + h/2*(g+gn1)); 
  y(:,n+1)=yy;
end
figure(1)
plot(t,y,'-')
legend('xa(t)','xb(t)','xc(t)')
title('Cantidad de contaminante en pantano')
grid on


Por tanto, la gráfica que se obtiene es: Aproximación numérica para h=0,1

3.2 Resolución numérica para h=1

El código utilizado para hallar la solución es: 

t0=0; tn=60;
y0=[0,0,0]'; 
A=[-7/400,1/100,0;3/400,-1/40,1/100;1/100,1/200,-1/40];%Matriz de coeficientes
%Datos para la discretización
h=1;
N=(tn-t0)/h;
%Vectores de tiempos y de solución aproximada
t=t0:h:tn 
y=zeros(3,N+1);
%Inicializamos el bucle
y(:,1)=y0;
yy=y0;
for n=1:N
  g=[15, 0, 0]' 
  gn1=[15, 0, 0]'
  yy=(eye(3)-h/2*A)\((eye(3)+h/2*A)*yy + h/2*(g+gn1)); 
  y(:,n+1)=yy;
end
figure(1)
plot(t,y,'-')
legend('xa(t)','xb(t)','xc(t)')
title('Cantidad de contaminante en pantano')
grid on


Y la gráfica obtenida: centro

4 Conclusiones extraídas de la solución

Ahora, se plantean una serie de interrogantes en torno a la solución obtenida en el problema, con el fin de entenderla mejor.

4.1 Aproximación de la solución a los 60 días

Las aproximaciónes de $X_{A}(t)$, $X_{B}(t)$ y $X_{C}(t)$ cuando han pasado justo 60 días del bombeo se han sacado del caso en el que h era 0.1, ya que de esta forma el resultado es mucho más preciso.

Las aproximaciones obtenidas son : $X_{A}(60)=577.588$; $X_{B}(60)=114.934$ y $X_{C}(60)=133.263$

4.2 Tiempo en el que XB(t) y XC(t) sobrepasarán los 400 Kg

Para ver cuándo los pantanos B y C tienen más de 400 Kg de contaminante en su interior se observa el vector que almacena las soluciones del problema. Cuando éste sea mayor que 400 Kg se mirará en qué subintervalo ocurre, obteniendo de esta manera el tiempo.

Sin embargo, se tendrá que ampliar el intervalo temporal de estudio, puesto que analizando en los primeros 60 días no conseguimos alcanzar aún los 400 kg de contaminante en los pantanos B y C. Por consiguiente, se analizará lo que ocurre en los 1000 días después de iniciar el vertido.

De esta manera, se procede a obtener la gráfica con el intervalo ampliado y con un paso de h=0,1, quedando como se muestra a continuación: Grafica1000dias.jpg De esta forma, se obtiene que la cantidad de 400 kg de contaminante se alcanza para los pantanos B y C en:

Pantano B: en la posición del vector de tiempos 1839. Entonces, como se ha tomado h=0,1, el tiempo en el que ocurre es a los 183,8 días (la cantidad de contaminante es de 400,143).

Pantano C: se saca que ocurre en el día 173,1 (actuando de forma análoga al pantano B).

4.3 Tiempo en el que XA(t) sobrepasará 1 tonelada

Para calcular esta información se trabaja análogamente a como se ha trabajado en el epígrafe anterior, de manera que el pantano A alcanza la tonelada de contaminante cuando pasan 191,8 días desde el primer vertido.

4.4 Pantano con mayor contaminación a lo largo de los 600 primeros días

Si se observa el vector de resultados que proporciona la gráfica con intervalo de 1000 días y h=0,1, concluimos que el pantano B tiene 592,753 kg de contaminante, mientras que el C posee 593,389 kg cuando han pasado 600 días desde que se empezó con el bombeo. Entonces, el pantano con más contaminación es el C.

4.5 Tendencias de las concentraciones de contaminación para tiempos largos

Para obtener la concentración de contaminante en cada uno de los pantanos basta con dividir entre el volumen. Modificamos el código tal y como indicamos a continuación, de manera que nos genere el gráfico de la concentración de contaminante a lo largo del tiempo. 

t0=0; tn=1000;
y0=[0,0,0]'; 
A=[-7/400,1/100,0;3/400,-1/40,1/100;1/100,1/200,-1/40];%Matriz de coeficientes
%Datos para la discretización
h=0.1;
N=(tn-t0)/h;
%Vectores de tiempos y de solución aproximada
t=t0:h:tn 
y=zeros(3,N+1);
%Inicializamos el bucle
y(:,1)=y0;
yy=y0;
for n=1:N
  g=[15, 0, 0]' 
  gn1=[15, 0, 0]'
  yy=(eye(3)-h/2*A)\((eye(3)+h/2*A)*yy + h/2*(g+gn1)); 
  y(:,n+1)=yy;
end
figure(1)
plot(t,y,'-')
legend('xa(t)','xb(t)','xc(t)')
title('Cantidad de contaminante en pantano')
grid on

%Para las concentraciones (el volumen es constante)
figure(2)
y(1,:)=y(1,:)/400
y(2:3,:)=y(2:3,:)/200
plot(t,y,'-')
grid on
title ('Concentraciones en kg/hm3')


La gráfica queda: Graficaconcentraciones.jpg

Como bien ilustra el gráfico, para tiempos largos las concentraciones tienden a 3 [math]Kg/Hm^3[/math].

5 Nueva situación de la red de pantanos

Tras dedicar grandes esfuerzos para eliminar el vertido de contaminante, el día 60 desde que se inició el mismo se coloca una depuradora que limpia el flujo de agua que entra al pantano A.

Teniendo en cuenta este hecho, volvemos a calcular la concentración de contaminante en cada pantano. Seguiremos designándolas por comodidad como $X_{A}(t)$, $X_{B}(t)$ y $X_{C}(t)$. Podríamos plantear el problema de valor inicial de la siguiente forma:


[math] \left\{\begin{matrix} X'_A(t)=\frac{-7X_A(t)}{400}+\frac{X_B(t)}{100} \\ X'_B(t)=\frac{3X_A(t)}{400}-\frac{X_B(t)}{40}+\frac{X_C(t)}{100}\\ X'_C(t)=\frac{X_A(t)}{100}+\frac{X_B(t)}{200}-\frac{X_C(t)}{40}\\X_A(t_{0})=577,588, \quad X_B(t_{0})=114,934, \quad X_C(t_{0})=133,263 \end{matrix}\right. [/math]

Y en forma matricial tendríamos:

\begin{equation*}

\begin{bmatrix}

X’_{A} (t)\\

X’_{B}(t)\\

X’_{C}(t) \end{bmatrix}

=

\begin{bmatrix}

-7/400 & 1/100 & 0 \\

3/400 & -1/40 & 1/100 \\

1/100 & 1/200 & -1/40

\end{bmatrix}

\begin{bmatrix}

X_{A} (t)\\

X_{B}(t)\\

X_{C}(t) \end{bmatrix}

\end{equation*}

5.1 Cantidad de contaminante tras 1 año de vertido

Para realizar esta aproximación debemos partir de los datos de cantidad de contaminante que tenemos a los 60 días de vertido, que son los siguientes:

$X_{A}(60)=577,588$ $X_{B}(60)=114,934$ $X_{C}(60)=133,263$

Una vez tomados los datos, escribimos el código correspondiente: 

clf;
%Aplicaremos el método de Runge-Kutta de orden 4.
%Primero definimos los datos del problema.
t0=0; tn=365;
x0=[577.588,114.934,133.263]'; 
%Aquí no hay matriz A, porque es un sistema no lineal
%Datos para la discretización
h=1;
N=(tn-t0)/h
%Vectores de tiempos y de solución aproximada
t=t0:h:tn; %Es muy sencillo de resolver
x=zeros(3,N+1); %Es una matriz de dos filas y N+1 columnas
%Inicializamos el bucle
x(:,1)=x0;
xx=x0; 
for n=1:N
  k1=[-7/400*xx(1)+xx(2)/100,3/400*xx(1)-xx(2)/40+xx(3)/100,xx(1)/100+xx(2)/200-xx(3)/40]';
  tp=t(n)+h/2; xp=xx+1/2*h*k1;
  k2=[-7/400*xp(1)+xp(2)/100,3/400*xp(1)-xp(2)/40+xp(3)/100,xp(1)/100+xp(2)/200-xp(3)/40]';
  xp=xx+1/2*h*k2;
  k3=[-7/400*xp(1)+xp(2)/100,3/400*xp(1)-xp(2)/40+xp(3)/100,xp(1)/100+xp(2)/200-xp(3)/40]';
  xp=xx+1/2*h*k3;
  k4=[-7/400*xp(1)+xp(2)/100,3/400*xp(1)-xp(2)/40+xp(3)/100,xp(1)/100+xp(2)/200-xp(3)/40]';
  xx=xx+h/6*(k1+2*k2+2*k3+k4);
  x(:,n+1)=xx;
 end
plot(t,x,'-')
legend('xa(t)','xb(t)','xc(t)')
grid on


Siendo la gráfica correspondiente:

Contaminante en cada pantano

5.2 Cantidad de contaminante tras 2 años de vertido

Para conocer una aproximación visual a la cantidad de contaminante utilizamos el siguiente código: 

clf;
%Aplicaremos el método de Runge-Kutta de orden 4.
%Primero definimos los datos del problema.
t0=0; tn=730;
x0=[577.588,114.934,133.263]'; 
%Aquí no hay matriz A, porque es un sistema no lineal
%Datos para la discretización
h=1;
N=(tn-t0)/h
%Vectores de tiempos y de solución aproximada
t=t0:h:tn; %Es muy sencillo de resolver
x=zeros(3,N+1); %Es una matriz de dos filas y N+1 columnas
%Inicializamos el bucle
x(:,1)=x0;
xx=x0; 
for n=1:N
  k1=[-7/400*xx(1)+xx(2)/100,3/400*xx(1)-xx(2)/40+xx(3)/100,xx(1)/100+xx(2)/200-xx(3)/40]';
  tp=t(n)+h/2; xp=xx+1/2*h*k1;
  k2=[-7/400*xp(1)+xp(2)/100,3/400*xp(1)-xp(2)/40+xp(3)/100,xp(1)/100+xp(2)/200-xp(3)/40]';
  xp=xx+1/2*h*k2;
  k3=[-7/400*xp(1)+xp(2)/100,3/400*xp(1)-xp(2)/40+xp(3)/100,xp(1)/100+xp(2)/200-xp(3)/40]';
  xp=xx+1/2*h*k3;
  k4=[-7/400*xp(1)+xp(2)/100,3/400*xp(1)-xp(2)/40+xp(3)/100,xp(1)/100+xp(2)/200-xp(3)/40]';
  xx=xx+h/6*(k1+2*k2+2*k3+k4);
  x(:,n+1)=xx;
 end
plot(t,x,'-')
legend('xA(t)','xB(t)','xC(t)')
grid on


Contaminante en cada pantano

5.3 Cantidad de contaminante tras 3 años de vertido

En este caso, utilizamos el siguiente código para hacernos una idea de la cantidad de contaminante que tenemos: 

clf;
%Aplicaremos el método de Runge-Kutta de orden 4.
%Primero definimos los datos del problema.
t0=0; tn=1095;
x0=[577.588,114.934,133.263]'; 
%Aquí no hay matriz A, porque es un sistema no lineal
%Datos para la discretización
h=1;
N=(tn-t0)/h
%Vectores de tiempos y de solución aproximada
t=t0:h:tn; %Es muy sencillo de resolver
x=zeros(3,N+1); %Es una matriz de dos filas y N+1 columnas
%Inicializamos el bucle
x(:,1)=x0;
xx=x0; 
for n=1:N
  k1=[-7/400*xx(1)+xx(2)/100,3/400*xx(1)-xx(2)/40+xx(3)/100,xx(1)/100+xx(2)/200-xx(3)/40]';
  tp=t(n)+h/2; xp=xx+1/2*h*k1;
  k2=[-7/400*xp(1)+xp(2)/100,3/400*xp(1)-xp(2)/40+xp(3)/100,xp(1)/100+xp(2)/200-xp(3)/40]';
  xp=xx+1/2*h*k2;
  k3=[-7/400*xp(1)+xp(2)/100,3/400*xp(1)-xp(2)/40+xp(3)/100,xp(1)/100+xp(2)/200-xp(3)/40]';
  xp=xx+1/2*h*k3;
  k4=[-7/400*xp(1)+xp(2)/100,3/400*xp(1)-xp(2)/40+xp(3)/100,xp(1)/100+xp(2)/200-xp(3)/40]';
  xx=xx+h/6*(k1+2*k2+2*k3+k4);
  x(:,n+1)=xx;
 end
plot(t,x,'-')
legend('xa(t)','xb(t)','xc(t)')
grid on


Contaminante en cada pantano

5.4 Cantidad de contaminante en otros casos

En caso de que nos interesara calcular la cantidad de contaminante en cada pantano para un tiempo que tienda a infinito, debemos tener en cuenta que a partir de alrededor de 700 días, las tres variables $X_{A}(t)$, $X_{B}(t)$ y $X_{C}(t)$ tienden a cero, por lo cual es básicamente innecesario programar el método numérico para un tiempo mayor a éste, sin importar el paso que tomemos para las particiones.


5.5 Máxima cantidad de contaminante en cada pantano

Como podemos ver fácilmente en las gráficas, $X_{A}(t)$ es siempre decreciente, siendo su máximo el que toma en el valor t=0, es decir, 577,588. Sin embargo, tanto la $X_{B}(t)$ como la $X_{C}(t)$ aumentan su valor desde que acaba el vertido, llegando, aproximadamente, a 165 para $X_{B}(t)$ y 175 para $X_{C}(t)$. Los valores exactos son los siguientes:

En el caso del pantano B, lo calculamos mediante el siguiente código: 

t0=0; tn=600;
y0=[577.588,114.934,133.263]'; 
A=[-7/400,1/100,0;3/400,-1/40,1/100;1/100,1/200,-1/40];%Matriz de coeficientes
%Datos para la discretización
h=0.1;
N=(tn-t0)/h;
%Vectores de tiempos y de solución aproximada
t=t0:h:tn 
y=zeros(3,N+1);
%Inicializamos el bucle
y(:,1)=y0;
yy=y0;
for n=1:N
  yy=(eye(3)-h/2*A)\((eye(3)+h/2*A)*yy); 
  y(:,n+1)=yy;
  if 3*y(1,n+1)/400-y(2,n+1)/40+y(3,n+1)/100<=0.0000001
  break
  end
end
n*h
y(2,n+1)
figure(1)
plot(t,y,'-')
legend('xa(t)','xb(t)','xc(t)')
grid on


Obtenemos que el máximo de $X_{B}(t)$ se ha producido a los 43 días y su máximo ha sido 164,05 kg.

En el caso del pantano C, siguiendo el mismo procedimiento: 

t0=0; tn=600;
y0=[577.588,114.934,133.263]'; 
A=[-7/400,1/100,0;3/400,-1/40,1/100;1/100,1/200,-1/40];%Matriz de coeficientes
%Datos para la discretización
h=0.01;
N=(tn-t0)/h;
%Vectores de tiempos y de solución aproximada
t=t0:h:tn 
y=zeros(3,N+1);
%Inicializamos el bucle
y(:,1)=y0;
yy=y0;
for n=1:N
  yy=(eye(3)-h/2*A)\((eye(3)+h/2*A)*yy); 
  y(:,n+1)=yy;
  if y(1,n+1)/100+y(2,n+1)/200-y(3,n+1)/40<=0.0000001
  break
  end
end
n*h
y(3,n+1)
figure(1)
plot(t,y,'-')
legend('xa(t)','xb(t)','xc(t)')
grid on


Resultando que el máximo de $X_{C}(t)$ se alcanza a los 34 días, tomando un valor de 174,67 kg.

A la hora de calcular el error en ambos casos tenemos que tener en cuenta que en el método de Runge-Kutta de orden 4, el error local de truncamiento es del orden de [math]h^5[/math], mientras que el error total acumulado es del orden de [math]h^4[/math]. Habiendo tomado un h=0,01, el error pedido es de [math]1*10^-8[/math], lo cual es despreciable frente a las cantidades de las que estamos hablando. Aunque no conozcamos el error exacto en estas medidas, sí podemos realizar una comparativa entre los distintos métodos de cálculo. En el caso del método del trapecio y el método de Runge-Kutta de orden 4:

Diferencia entre el método del trapecio y Runge-Kutta de orden 4

Podemos observar que en el caso más desfavorable, la diferencia llega a ser de 0,25 kg, lo cual es pequeño para las cantidades que manejamos, de varios cientos de kg.

5.6 Potabilización del agua

Las autoridades consideran que las aguas del pantano dejan de estar contaminadas cuando la concentración de contaminante en el agua es menor a 0,01 [math]Kg/Hm^3[/math]. Para acabar con los cortes de agua y volver a usar el agua de los pantanos para uso humano, debemos calcular cuánto tiempo deberemos esperar para poder hacerlo. Usamos el siguiente programa: 


clf;
%Aplicaremos el método de Runge-Kutta de orden 4.
%Primero definimos los datos del problema.
t0=0; tn=1500;
x0=[577.588,114.934,133.263]'; 
%Aquí no hay matriz A, porque es un sistema no lineal
%Datos para la discretización
h=1;
N=(tn-t0)/h
%Vectores de tiempos y de solución aproximada
t=t0:h:tn; %Es muy sencillo de resolver
x=zeros(3,N+1); %Es una matriz de dos filas y N+1 columnas
%Inicializamos el bucle
x(:,1)=x0;
xx=x0; 
for n=1:N
  k1=[-7/400*xx(1)+xx(2)/100,3/400*xx(1)-xx(2)/40+xx(3)/100,xx(1)/100+xx(2)/200-xx(3)/40]';
  tp=t(n)+h/2; xp=xx+1/2*h*k1;
  k2=[-7/400*xp(1)+xp(2)/100,3/400*xp(1)-xp(2)/40+xp(3)/100,xp(1)/100+xp(2)/200-xp(3)/40]';
  xp=xx+1/2*h*k2;
  k3=[-7/400*xp(1)+xp(2)/100,3/400*xp(1)-xp(2)/40+xp(3)/100,xp(1)/100+xp(2)/200-xp(3)/40]';
  xp=xx+1/2*h*k3;
  k4=[-7/400*xp(1)+xp(2)/100,3/400*xp(1)-xp(2)/40+xp(3)/100,xp(1)/100+xp(2)/200-xp(3)/40]';
  xx=xx+h/6*(k1+2*k2+2*k3+k4);
  x(:,n+1)=xx;
end
x(1,:)=x(1,:)/400
x(2,:)=x(2,:)/200
x(2:3,:)=x(2:3,:)/200
plot(t,x,'-')
legend('xa(t)','xb(t)','xc(t)')
grid on


Contaminante en cada pantano

Siendo en la gráfica las concentraciones de contaminante en cada pantano. Si además queremos conocer a partir de qué día se puede beber el agua de cada pantano, podemos plantearlo de la siguiente manera:

Para el pantano A, utilizaremos el siguiente código: 

clf;
%Aplicaremos el método de Runge-Kutta de orden 4.
%Primero definimos los datos del problema.
t0=0; tn=1500;
x0=[577.588,114.934,133.263]'; 
%Aquí no hay matriz A, porque es un sistema no lineal
%Datos para la discretización
h=1;
N=(tn-t0)/h
%Vectores de tiempos y de solución aproximada
t=t0:h:tn; %Es muy sencillo de resolver
x=zeros(3,N+1); %Es una matriz de dos filas y N+1 columnas
%Inicializamos el bucle
x(:,1)=x0;
xx=x0; 
for n=1:N
  k1=[-7/400*xx(1)+xx(2)/100,3/400*xx(1)-xx(2)/40+xx(3)/100,xx(1)/100+xx(2)/200-xx(3)/40]';
  tp=t(n)+h/2; xp=xx+1/2*h*k1;
  k2=[-7/400*xp(1)+xp(2)/100,3/400*xp(1)-xp(2)/40+xp(3)/100,xp(1)/100+xp(2)/200-xp(3)/40]';
  xp=xx+1/2*h*k2;
  k3=[-7/400*xp(1)+xp(2)/100,3/400*xp(1)-xp(2)/40+xp(3)/100,xp(1)/100+xp(2)/200-xp(3)/40]';
  xp=xx+1/2*h*k3;
  k4=[-7/400*xp(1)+xp(2)/100,3/400*xp(1)-xp(2)/40+xp(3)/100,xp(1)/100+xp(2)/200-xp(3)/40]';
  xx=xx+h/6*(k1+2*k2+2*k3+k4);
  x(:,n+1)=xx;
  if xx(1)/400<=0.01
 break
  end
end
n+1
plot(t,x,'-')
legend('xA(t)','xB(t)','xC(t)')
grid on


Obtenemos el valor de 561 días desde el inicio del vertido.

En el caso del pantano B: 

clf;
%Aplicaremos el método de Runge-Kutta de orden 4.
%Primero definimos los datos del problema.
t0=0; tn=1500;
x0=[577.588,114.934,133.263]'; 
%Aquí no hay matriz A, porque es un sistema no lineal
%Datos para la discretización
h=1;
N=(tn-t0)/h
%Vectores de tiempos y de solución aproximada
t=t0:h:tn; %Es muy sencillo de resolver
x=zeros(3,N+1); %Es una matriz de dos filas y N+1 columnas
%Inicializamos el bucle
x(:,1)=x0;
xx=x0; 
for n=1:N
  k1=[-7/400*xx(1)+xx(2)/100,3/400*xx(1)-xx(2)/40+xx(3)/100,xx(1)/100+xx(2)/200-xx(3)/40]';
  tp=t(n)+h/2; xp=xx+1/2*h*k1;
  k2=[-7/400*xp(1)+xp(2)/100,3/400*xp(1)-xp(2)/40+xp(3)/100,xp(1)/100+xp(2)/200-xp(3)/40]';
  xp=xx+1/2*h*k2;
  k3=[-7/400*xp(1)+xp(2)/100,3/400*xp(1)-xp(2)/40+xp(3)/100,xp(1)/100+xp(2)/200-xp(3)/40]';
  xp=xx+1/2*h*k3;
  k4=[-7/400*xp(1)+xp(2)/100,3/400*xp(1)-xp(2)/40+xp(3)/100,xp(1)/100+xp(2)/200-xp(3)/40]';
  xx=xx+h/6*(k1+2*k2+2*k3+k4);
  x(:,n+1)=xx;
  if xx(2)/200<=0.01
 break
  end
end
n+1
plot(t,x,'-')
legend('xA(t)','xB(t)','xC(t)')
grid on


En cuyo caso obtenemos la cifra de 641 días.

Y para el pantano C: 

clf;
%Aplicaremos el método de Runge-Kutta de orden 4.
%Primero definimos los datos del problema.
t0=0; tn=1500;
x0=[577.588,114.934,133.263]'; 
%Aquí no hay matriz A, porque es un sistema no lineal
%Datos para la discretización
h=1;
N=(tn-t0)/h
%Vectores de tiempos y de solución aproximada
t=t0:h:tn; %Es muy sencillo de resolver
x=zeros(3,N+1); %Es una matriz de dos filas y N+1 columnas
%Inicializamos el bucle
x(:,1)=x0;
xx=x0; 
for n=1:N
  k1=[-7/400*xx(1)+xx(2)/100,3/400*xx(1)-xx(2)/40+xx(3)/100,xx(1)/100+xx(2)/200-xx(3)/40]';
  tp=t(n)+h/2; xp=xx+1/2*h*k1;
  k2=[-7/400*xp(1)+xp(2)/100,3/400*xp(1)-xp(2)/40+xp(3)/100,xp(1)/100+xp(2)/200-xp(3)/40]';
  xp=xx+1/2*h*k2;
  k3=[-7/400*xp(1)+xp(2)/100,3/400*xp(1)-xp(2)/40+xp(3)/100,xp(1)/100+xp(2)/200-xp(3)/40]';
  xp=xx+1/2*h*k3;
  k4=[-7/400*xp(1)+xp(2)/100,3/400*xp(1)-xp(2)/40+xp(3)/100,xp(1)/100+xp(2)/200-xp(3)/40]';
  xx=xx+h/6*(k1+2*k2+2*k3+k4);
  x(:,n+1)=xx;
  if xx(3)/200<=0.01
 break
  end
end
n+1
plot(t,x,'-')
legend('xA(t)','xB(t)','xC(t)')
grid on


Para este pantano, deberemos esperar 630 días para poder beber de sus aguas.

En conjunto, y tomando la condición más restrictiva, tendríamos que esperar 641 días para que el agua de todos los pantanos fuera potable.

Obtenido de «https://mat.caminos.upm.es/w/index.php?title=Modelo_de_distribución_de_un_contaminante_en_una_red_de_pantanos&oldid=10372»