Trabajo realizado por estudiantes
Título	Difusión de una sustancia contaminante. Grupo 24-C
Asignatura	Ecuaciones Diferenciales
Curso	Curso 2014-15
Autores	Jose Antonio Martinez Montalvo 1494 Isaac Rebollo Palos 1522 Daniel Pascual Cobos 1690 Rodrigo Bellot Rodriguez 1270
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

Contenido

1 Introducción
2 Ecuación de difusión
- 2.1 Deducción
- 2.2 Conservación de la masa contaminante
3 Problema Propuesto

1 Introducción

El trabajo tiene como objetivo el estudio de la mezcla de dos sustancias, una de ellas contaminante, en un tubo largo.

Su planteamiento se realizará de forma similar al de la transmisión de calor a lo largo de una varilla, visto en clase.

Analizaremos el modelo que define nuestro problema (ecuación de difusión), varios puntos teóricos del mismo (la conservación de la masa contaminante, la solución estacionaria) y terminar por plantear un problema concreto el cual resolveremos mediante métodos numéricos, modificando según que condiciones para poder interpretar las distintas situaciones que pueden modelizar nuestro problema.

2 Ecuación de difusión

2.1 Deducción

Sea un tubo largo de longitud L, orientado en la dirección x y cuya sección es constante desde x=0 hasta x=L. Supondremos el tubo aislado tanto superficialmente como lateralmente.

Configuración inicial de las concentraciones de contaminante en el tubo.

Sean dos sustancias contenidas en el tubo, una de ellas contaminante. Denotaremos por u(x,t) la concentración de contaminante, la cual dependerá únicamente de la posición del tubo y del tiempo, manteniéndose constante en cada sección transversal del tubo. Medido en mol/m2s.

Sea F(x,t) el flujo de contaminante (análogo al flujo de calor), definiéndose como la cantidad de contaminante que atraviesa una sección transversal por unidad de tiempo y area, dependiendo del número de moles, siendo éste último proporcional a la masa de la sustancia. Debido a los aislantes laterales y superficiales, el flujo sera a lo largo del eje x.

Debido al principio de conservación de la masa la variación de la concentración de contaminante en cada posición del tubo en función del tiempo, es igual a la suma del flujo de contaminante a través de los extremos del tubo por unidad de tiempo, más la concentración de contaminante generada en el interior por unidad de tiempo. Suponemos igual a 0 las perdidas y/o ganancias en el interior del tubo (problema homogéneo). La ley de Flick (similar a Furier) determina que el flujo de difusión del contaminante es proporcional a la variación de concentración: [math]F(x,t)=-D\frac{\partial u}{\partial x}[/math]

siendo D el coeficiente de difusión (medido en m2/s), que dependerá de las propiedades químicas de los compuestos.

De esta forma y continuando con la semejanza a la transmisión de calor, la concentración de sustancia contaminante en un instante de tiempo t y en una sección transversal que dista x unidades del extremo izquierdo del tubo satisfará: [math]u_t(x,t)=-\frac{\partial }{\partial x}(-D u_x(x,t)) = D u_{xx}(x,t)[/math]

Pasando todo al termino de la izquierda obtendremos la ecuación de difusión: [math]u_t(x,t)- D u_{xx}(x,t)= 0[/math]

dónde D es el coeficiente de difusión anteriormente citado.

Para ello hemos definido [math] A [/math] como la variable de superficie, tomando una sección pequeña del tubo, designada por [math]\Delta x[/math]: La variación de la concentración de contaminante en función del tiempo es:

[math] A u(x,t) \Delta x [/math]

Derivamos respecto del tiempo:

[math] A u_t(x,t) \Delta x [/math]

A continuación, se considera nula la concentración de contaminante generada en el interior por unidad de tiempo debido a la ausencia de sumideros.

Suponemos que Δx>0 y la concentración de contaminante en un tiempo t es menor en x+ Δx que en x, entonces u(x+ Δx) – u(x,t) < 0, y como Δx es pequeño, se tiene que ux(x,t)<0 y el flujo de difusión del contaminante es positivo y va hacia la derecha en la dirección del eje x.

El ﬂujo de calor en un intervalo [math][x, x + ∆x ] [/math] viene dado dado por:

[math]F(x,t)A-F(x+\Delta x, t)A[/math]

Igualando los términos anteriores y dividiendo por [math]A \Delta x[/math] se obtiene:

[math]A u_t(x,t) \Delta x = F(x,t)A-F(x+\Delta x, t)A[/math]:

[math]u_t(x,t)= \frac {F(x,t)-F(x+\Delta x, t)}{\Delta x}[/math]

Haciendo que [math]\Delta x[/math] tienda a 0:

[math]u_t(x,t)= -F_x(x,t)[/math]

Y aplicando la ley de Fick:

[math]u_t(x,t)=-\frac{\partial }{\partial x}(-D u_x(x,t)) = D u_xx(x,t) [/math]

Así se obtiene la ecuación de difusión de la sustancia contaminante a lo largo del tubo anteriormente citada:

[math]u_t(x,t)- D u_xx(x,t)= 0 [/math]

2.2 Conservación de la masa contaminante

A partir de la ecuación u(x,t)se deduce que la masa total del contaminante se conserva a lo largo del tiempo.Para ello se integra:

[math]\frac{d}{dt}\displaystyle\int_{0}^{L}u(x,t)\,dx=\displaystyle\int_{0}^{L}u_{xx}(x,t)\,dx=u_x(x,t)|_0^L=0[/math]

clear all
clf

% Datos
L=7;
T=10;
q=1;

% Discretizacion espacial
h=0.1;
N=L/h;
x=0:h:L; 

% Discretizacion temporal
dt=h/4; % Tamaño de paso
t=0:dt:T; % Vector de tiempos

% Matriz K
K=2*diag(ones(1,N+1))-diag(ones(1,N),1)-diag(ones(1,N),-1);
K(1,2)=-2; % Condicion Neumann
K(N+1,N)=-2; % Condicion Neumann
K=q*K/h^2;

% Definimos F
F=zeros(N+1,1);

% Calculamos u0
u0=[zeros(1,1+N*5/7),3*ones(1,N*2/7)]';
% Guardamos la solucion para t=0 en la matriz solucion
sol(1,:)=u0'; 
U=u0; % Inicializacion

% Calculamos U_n+1 a partir de U_n con un bucle for
for j=1:length(t)-1 % Hasta length(t)-1 porque ya tenemos el valor para t=0
U=(eye(N+1)+dt*K)\(U+dt*F); % Euler Implicito
sol(j+1,:)= U'; % Guardamos la solucion
end

% Calculamos el area de la funcion para cada t (el area será la cantidad de contaminante)
contaminante=zeros(1,length(t));
for i=1:length(t)
f=sol(i,:);
area=trapz(x,f);
contaminante(i)=area;
end

% Obtenemos en la figura 1 la cantidad de contaminante con el tiempo,y en la
% figura 2 la concentracion en el punto medio de la barra a lo largo del tiempo
figure(1)
plot(t,contaminante)
xlabel('Tiempo')
ylabel('Cantidad de contaminante')
punto_medio=(length(x)-1)/2;
figure(2)
plot(t,sol(:,punto_medio))

En las gráficas obtenidas se observa cómo la masa total de contaminante se conserva a lo largo del tiempo, y como varía la concentracion de contaminante a lo largo del tiempo en el punto medio.

Conservación masa total de contaminante.

Concentracion en el punto medio.

3 Problema Propuesto

3.1 Sistema de ecuaciones

Sea el problema:

[math] (P)\left\{\begin{matrix}\\u_t-u_{xx}=0 ,x\in\ (0,7), t\gt0\\u_x(0,t)=0, t\gt0\\u_x(7,t)=0, t\gt0\\u(x,0)=u_0, x\in\ (0,7)\end{matrix}\right. [/math]

Sean las condiciones iniciales: :[math] u(x,0)=\left\{\begin{matrix}\\0, x≤5\\3, x\gt5\end{matrix}\right. [/math]

Pasamos a su interpretación física.

3.2 Modelización

El problema modeliza el proceso de difusión de una sustancia contaminante a lo largo de un tubo de longitud x=7, delgado, difusor, con superficie lateral aislada y con coeficiente de difusión D igual a 1.

En los puntos interiores del tubo no existentes ni fuentes ni sumideros de sustancia contaminante.

Tanto el extremo izquierdo como el extremo derecho del tubo se encuentran aislados de forma que el flujo de sustancia en ellos es igual a 0.

En el instante inicial, las secciones transversales que distan 5 o menos unidades de longitud del extremo izquierdo tendrán una concentración inicial de sustancia contaminante igual a 0. Las secciones que distan mas de 5 unidades de ese extremo tendrán una concentración de contaminante igual a 3.

3.3 Resolución

3.3.1 Método diferencias finitas

Se trata de un método numérico que nos proporciona una solución aproximada de distintos problemas de ecuaciones en derivadas parciales, aplicable a nuestro caso. Primero realizamos una discretización del espacio, tomando como longitud de paso h y desde x=0 hasta x=7, obteniendo un vector N con x+1 elementos.

Seguidamente aplicamos la ecuacion diferencial a cada nodo interior xn, aproximamos u_xx mediante: [math]u_{xx}(x,t)\simeq\frac{u(x_{n-1},t)-2u(x_n,t)+u(x_{n+1},t)}{h^2}=0[/math]

Aplicamos la notación [math] u_n(t) = u(x_n,t)[/math] resulta el sistema de N-1 ecuaciones: \begin{array}{c}u'_n(t)+\frac{-u_{n-1}(t)+2u_n(t)-u_{n+1}(t)}{h^2}=0\end{array} para n=1,2...N-1

Aplicamos las condiciones de contorno:

\begin{array}{c}u'_0(t)=0\\u'_n(t)=0\end{array} Resultando el sistema: \begin{array}{c}U'+KU=0=F\\U(0)=U^0\end{array}

Siendo [math]U^0[/math] las condiciones iniciales en el instante inicial se verifica :[math] u(x,0)=\left\{\begin{matrix}\\0, x≤5\\3, x\gt5\end{matrix}\right. [/math]

3.3.1.1 Trapecio

Se resuelve el sistema en primer lugar utilizando el método del trapecio tomando [math]\triangle t= \triangle x/4[/math] en [math]t\in\ [0,7][/math], con [math]\triangle x=0.1[/math]

El código en MATLAB del método sería:

clear all
clf
% u_t-qu_xx=0, x en (0,L)
% u_x(0,t)=0 Condicion Neumann
% u_x(L,t)=0 Condicion Neumann
% u(x,0)=u0(x)

% Datos
L=7;
T=7;
q=1;

% Discretizacion espacial
h=0.1;
N=L/h;
x=0:h:L; 

% Discretizacion temporal
dt=h/4; % Tamaño de paso
t=0:dt:T; % Vector de tiempos

% Matriz K
K=2*diag(ones(1,N+1))-diag(ones(1,N),1)-diag(ones(1,N),-1);
K(1,2)=-2; % Condicion Neumann
K(N+1,N)=-2; % Condicion Neumann
K=q*K/h^2;

% Definimos F
F=zeros(N+1,1);

% Calculamos u0
u0=[zeros(1,1+N*5/7),3*ones(1,N*2/7)]';
% Guardamos la solucion para t=0 en la matriz solucion
sol(1,:)=u0'; 
U=u0; % Inicializacion

% Calculamos U_n+1 a partir de U_n con un bucle for
for j=1:length(t)-1 % Voy hasta la longitud de t menos 1 porque ya conozco un valor
U=(eye(N+1)+(dt*K)/2)\((eye(N+1)-(K*dt)/2)*U); % Trapecio
sol(j+1,:)= U'; % Guardamos solucion
end
 % Dibujamos la solucion
[xx,tt]=meshgrid(x,t);
surf(xx,tt,sol);

Se obtiene la siguiente gráfica concentración / espacio-tiempo

Resolución del sistema de difusión de una sustancia contaminante por el método del trapecio.

3.3.1.2 Euler explícito

clear all
clf
% u_t-qu_xx=0, x en (0,L)
% u_x(0,t)=0 Condicion Neumann
% u_x(L,t)=0 Condicion Neumann
% u(x,0)=u0(x)

% Datos
L=7;
T=7;
q=1;

% Discretizacion espacial
h=0.1;
N=L/h;
x=0:h:L; 

% Discretizacion temporal
dt=h/4; % Tamaño de paso
t=0:dt:T; % Vector de tiempos

% Matriz K
K=2*diag(ones(1,N+1))-diag(ones(1,N),1)-diag(ones(1,N),-1);
K(1,2)=-2; % Condicion Neumann
K(N+1,N)=-2; % Condicion Neumann
K=q*K/h^2;

% Definimos F
F=zeros(N+1,1);

% Calculamos u0
u0=[zeros(1,1+N*5/7),3*ones(1,N*2/7)]';
% Guardamos la solucion para t=0 en la matriz solucion
sol(1,:)=u0'; 
U=u0; % Inicializacion

% Calculamos U_n+1 a partir de U_n con un bucle for
for j=1:length(t)-1 % Voy hasta la longitud de t menos 1 porque ya conozco un valor
U=U-dt*K*U; % metodo explicito
sol(j+1,:)= U'; % Guardamos solucion
end
% Dibujamos la solucion
[xx,tt]=meshgrid(x,t);
surf(xx,tt,sol);

La gráfica concentración / espacio-tiempo obtenida es:

Resolución del sistema de difusión de una sustancia contaminante por el método de euler explicito.

Se puede observar que el metodo es inestable para este paso de tiempo.

Para que el método de Euler explícito de diferencias finitas sea estable hace falta que se cumpla que dt/(h^2) sea menor que 0.5

Al usar un paso de tiempo dt=h/40 si que funciona, y nos sale esta gráfica:

Resolución del sistema de difusión de una sustancia contaminante por el método de euler explicito, con dt=h/40.

3.3.1.3 Euler implícito

clear all
clf
% u_t-qu_xx=0, x en (0,L)
% u_x(0,t)=0 Condicion Neumann
% u_x(L,t)=0 Condicion Neumann
% u(x,0)=u0(x)

% Datos
L=7;
T=7;
q=1;

% Discretizacion espacial
h=0.1;
N=L/h;
x=0:h:L; 

% Discretizacion temporal
dt=h/4; % Tamaño de paso
t=0:dt:T; % Vector de tiempos

% Matriz K
K=2*diag(ones(1,N+1))-diag(ones(1,N),1)-diag(ones(1,N),-1);
K(1,2)=-2; % Condicion Neumann
K(N+1,N)=-2; % Condicion Neumann
K=q*K/h^2;

% Definimos F
F=zeros(N+1,1);

% Calculamos u0
u0=[zeros(1,1+N*5/7),3*ones(1,N*2/7)]';
% Guardamos la solucion para t=0 en la matriz solucion
sol(1,:)=u0'; 
U=u0; % Inicializacion

% Calculamos U_n+1 a partir de U_n con un bucle for
for j=1:length(t)-1 % Hasta length(t)-1 porque ya tenemos el valor para t=0
U=(eye(N+1)+dt*K)\(U+dt*F); % Euler Implicito
sol(j+1,:)= U'; % Guardamos la solucion
end

% Dibujamos la solucion
[xx,tt]=meshgrid(x,t);
surf(xx,tt,sol);

Se obtiene la siguiente gráfica concentración / espacio-tiempo:

Resolución del sistema de difusión de una sustancia contaminante por el método de euler implícito.

3.3.1.4 Heun

clear all
clf
% u_t-qu_xx=0, x en (0,L)
% u_x(0,t)=0 Condicion Neumann
% u_x(L,t)=0 Condicion Neumann
% u(x,0)=u0(x)

% Datos
L=7;
T=7;
q=1;

% Discretizacion espacial
h=0.1;
N=L/h;
x=0:h:L; 

% Discretizacion temporal
dt=h/4; % Tamaño de paso
t=0:dt:T; % Vector de tiempos

% Matriz K
K=2*diag(ones(1,N+1))-diag(ones(1,N),1)-diag(ones(1,N),-1);
K(1,2)=-2; % Condicion Neumann
K(N+1,N)=-2; % Condicion Neumann
K=q*K/h^2;

% Definimos F
F=zeros(N+1,1);

% Calculamos u0
u0=[zeros(1,1+N*5/7),3*ones(1,N*2/7)]';
% Guardamos la solucion para t=0 en la matriz solucion
sol(1,:)=u0'; 
U=u0; % Inicializacion

% Calculamos U_n+1 a partir de U_n con un bucle for
for j=1:length(t)-1 % Voy hasta la longitud de t menos 1 porque ya conozco un valor
K1=K*U; % metodo heun
K2=K*(U-dt*K*U);
U=U-dt*0.5*(K1+K2);
sol(j+1,:)= U'; % Guardamos solucion
end
% Dibujamos la solucion
[xx,tt]=meshgrid(x,t);
surf(xx,tt,sol);

La gráfica concentración / espacio-tiempo obtenida es:

Resolución del sistema de difusión de una sustancia contaminante por el método de Heun.

Se observa que el método de Heun tampoco es estable para este paso de tiempo.

El paso de tiempo dt=h/4 nos sirve para resolver el problema con los métodos implícitos del Trapecio y de Euler implícito, mientras que Heun y Euler explícito que son métodos explícitos son inestables.

Para que funcione hace falta reducir tambien el tamaño de paso.

Cuando usamos dt=h/20 nos queda esta gráfica:

Resolución del sistema de difusión de una sustancia contaminante por el método de Heun con dt=h/20.

Y con dt=h/40 queda una gráfica mucho mas precisa:

Resolución del sistema de difusión de una sustancia contaminante por el método de Heun con dt=h/40.

3.3.2 Método de Furier

El método de Fourier nos proporciona, al igual que el método de diferencias finitas, una aproximación de la solución del problema:

[math] (P)\left\{\begin{matrix}\\u_t-u_{xx}=0 ,x\in\ (0,7), t\gt0\\u_x(0,t)=0, t\gt0\\u_x(7,t)=0, t\gt0\\u(x,0)=u_0, x\in\ (0,7)\end{matrix}\right. [/math]

Cuyas condiciones iniciales son: [math] u(x,0)=\left\{\begin{matrix}\\0, x≤5\\3, x\gt5\end{matrix}\right. [/math]

Dado que se trata de un sistema con condiciones de frontera homogéneas, el problema de autovalores asociado es:

[math] (PA)\left\{\begin{matrix}\\x''+μx=0 ,x\in\ (0,7)\\x'(0)=0\\x'(7)=0\end{matrix}\right. [/math]

Y los autovalores y autofunciones salvo una constante son:

[math]μ_k=\frac{\pi^2k^2}{49}[/math]: [math]φ_k=cos(kπ/7)x[/math] k=0,1,2,3...

Se aproxima el dato inicial y el segundo miembro, se elige Q y se toma:

[math]\displaystyle\sum_{k=0}^Q c_kφ_k(x)=0[/math]:

[math]\displaystyle\sum_{k=0}^Q a_kφ_k(x)=\left\{\begin{matrix}\\0, x≤5\\3, x\gt5\end{matrix}\right.[/math]

donde se deduce que c_k=0 y [math]a_k=\frac{\displaystyle\int_{0}^{7} u_0cos(kπ/7)x\, dx}{\displaystyle\int_{0}^{7} cos^2(kπ/7)x\, dx}[/math]

El problema aproximado queda de la forma:

[math] (P)\left\{\begin{matrix}\\u_t-u_{xx}=0 ,x\in\ (0,7), t\gt0\\u_x(0,t)=0, t\gt0\\u_x(7,t)=0, t\gt0\\u(x,0)=u_0=\displaystyle\sum_{k=0}^Q a_kφ_k(x), x\in\ (0,7)\end{matrix}\right. [/math]

La solución aproximada sería:

[math]u(x,t)=\displaystyle\sum_{k=0}^Q T_k(t)φ_k(x)[/math]

Sustituyendo en la ecuación de difusión de una sustancia contaminante e igualando los coeficicientes de Fourier se obtiene el sistema:

[math]\left\{\begin{matrix}\\T'_k(t)+μ_kT_k(t)=c_k=0\\T_k(0)=a_k\end{matrix}\right.[/math]

Resolviendo el problema de valor inicial, se deduce que [math] T_k(t)=Ce^{-μ_k t}[/math]

La solución aproximada queda de la forma [math]u(x,t)=\displaystyle\sum_{k=0}^Q Ce^{-μ_k t} cos(kπ/7)x[/math]

El valor de la constante C depende de la condición inicial,que es la que condiciona el valor de a_k.

El código en MATLAB del método de Fourier con 1,3,5,10 y 20 términos de la serie sería:

clear all

L=7;
T=0.5;
% Datos d
N=70; % Intervalos de x
h=L/N;
x=0:h:L;
dt=h/4;
t=0:dt:T;
% Aproximacion
sol=zeros(length(t),N+1);
u0=[zeros(1,1+N*5/7),3*ones(1,N*2/7)];
for k=0:1
phi=cos((k*pi/7)*x);
a=(trapz(x,u0.*phi))/(trapz(x,phi.^2));
T=a*exp(-(k^2)*(pi^2)*1/49*t);
sol=sol+T'*phi;
end
[xx,tt]=meshgrid(x,t);
subplot(3,2,1);
surf(xx,tt,sol)
 
for k=0:3
phi=cos((k*pi/7)*x);
a=(trapz(x,u0.*phi))/(trapz(x,phi.^2));
T=a*exp(-(k^2)*(pi^2)*1/49*t);
sol=sol+T'*phi;
end
[xx,tt]=meshgrid(x,t);
subplot(3,2,2);
surf(xx,tt,sol)
 
 
for k=0:5
phi=cos((k*pi/7)*x);
a=(trapz(x,u0.*phi))/(trapz(x,phi.^2));
T=a*exp(-(k^2)*(pi^2)*1/49*t);
sol=sol+T'*phi;
end
[xx,tt]=meshgrid(x,t);
subplot(3,2,3);
surf(xx,tt,sol)
 
for k=0:10
phi=cos((k*pi/7)*x);
a=(trapz(x,u0.*phi))/(trapz(x,phi.^2));
T=a*exp(-(k^2)*(pi^2)*1/49*t);
sol=sol+T'*phi;
end
[xx,tt]=meshgrid(x,t);
subplot(3,2,4);
surf(xx,tt,sol)
for k=0:20
phi=cos((k*pi/7)*x);
a=(trapz(x,u0.*phi))/(trapz(x,phi.^2));
T=a*exp(-(k^2)*(pi^2)*1/49*t);
sol=sol+T'*phi;
end
[xx,tt]=meshgrid(x,t);
subplot(3,2,5);
surf(xx,tt,sol)

Resolución del sistema de difusión de una sustancia contaminante por el método de Fourier.

Si cambiamos el valor de T=0.5 por T=7, podemos comprobar que la gráfica obtenida por Fourier es muy similar a las que obteníamos por el método de diferencias finitas.

clear all
clf
L=7;
T=7;
% Datos d
N=70; % Intervalos de x
h=L/N;
x=0:h:L;
dt=h/4;
t=0:dt:T;
% Aproximacion
sol=zeros(length(t),N+1);
u0=[zeros(1,1+N*5/7),3*ones(1,N*2/7)];
for k=0:20
phi=cos((k*pi/7)*x);
a=(trapz(x,u0.*phi))/(trapz(x,phi.^2));
T=a*exp(-(k^2)*(pi^2)*1/49*t);
sol=sol+T'*phi;
end
[xx,tt]=meshgrid(x,t);
surf(xx,tt,sol)

Resolución del sistema de difusión de una sustancia contaminante por el método de Fourier para 20 terminos con T=7.

3.4 Solución estacionaria

La solucion estacionaria debera ser nuestra cantidad de contaminante 5,85 dividido por 7, que es 0.83571

Se calcula la distancia entre la solución estacionaria y la solución u(x,t) para t=0,1,2,10:

clear all
clf
% u_t-qu_xx=0, x en (0,L)
% u_x(0,t)=0 Condicion Neumann
% u_x(L,t)=0 Condicion Neumann
% u(x,0)=u0(x)

% Datos
L=7;
T=10;
q=1;

% Discretizacion espacial
h=0.1;
N=L/h;
x=0:h:L; 

% Discretizacion temporal
dt=h/4; % Tamaño de paso
t=0:dt:T; % Vector de tiempos

% Matriz K
K=2*diag(ones(1,N+1))-diag(ones(1,N),1)-diag(ones(1,N),-1);
K(1,2)=-2; % Condicion Neumann
K(N+1,N)=-2; % Condicion Neumann
K=q*K/h^2;

% Definimos F
F=zeros(N+1,1);

% Calculamos u0
u0=[zeros(1,1+N*5/7),3*ones(1,N*2/7)]';
% Guardamos la solucion para t=0 en la matriz solucion
sol(1,:)=u0'; 
U=u0; % Inicializacion

% Calculamos U_n+1 a partir de U_n con un bucle for
for j=1:length(t)-1 % Voy hasta la longitud de t menos 1 porque ya conozco un valor
U=(eye(N+1)+(dt*K)/2)\((eye(N+1)-(K*dt)/2)*U); % Trapecio
sol(j+1,:)= U'; % Guardamos solucion
end
 % Dibujamos la solucion
[xx,tt]=meshgrid(x,t);
surf(xx,tt,sol);

dist0=(sqrt(sum((sol(1,:)-(5.85/7)).^2)/(N+1)))
dist1=(sqrt(sum((sol(41,:)-(5.85/7)).^2)/(N+1)))
dist2=(sqrt(sum((sol(81,:)-(5.85/7)).^2)/(N+1)))
dist10=(sqrt(sum((sol(401,:)-(5.85/7)).^2)/(N+1)))

De donde se obtiene:

dist0=1.3495

dist1=0.90614

dist2=0.71079

dist10=0.13935

Para tiempos grandes la cantidad de contaminante se ira acercando al valor estacionario de 5.85/7. Aqui tenemos una grafica donde los valores mas "planos" son de T=10 y T=20

Concentracion para distintos tiempos a lo largo de la barra.

vamos a calcular la media cuadrática de las distancias desde cada punto de la función al estado estacionario para averiguar el tiempo que tarda la concentración en alcanzar el estado estacionario con un error del 5%:

La media cuadrática para N valores [math][x_1,x_2,x_3,...,x_N][/math] de una variable discreta x viene dada por: [math]x=\sqrt{\mathstrut \frac{1}{N}\sum_{k=1}^N}x_k^2=(\frac{x_1^2+x_2^2+...+x_N^2}{N})^{1/2} [/math]

Vamos a obtener las soluciones de nuestra ecuacion de difusion hasta un tiempo mayor (T=100 por ejemplo), y una vez tengamos nuestras soluciones usaremos un bucle que se rompa cuando la media cuadratica de las distancias a la solución estacionaria sea menor que un 5% de esta solucion estacionaria.

clear all
clf
% u_t-qu_xx=0, x en (0,L)
% u_x(0,t)=0 Condicion Neumann
% u_x(L,t)=0 Condicion Neumann
% u(x,0)=u0(x)

% Datos
L=7;
T=20;
q=1;

% Discretizacion espacial
h=0.1;
N=L/h;
x=0:h:L; 

% Discretizacion temporal
dt=h/4; % Tamaño de paso
t=0:dt:T; % Vector de tiempos

% Matriz K
K=2*diag(ones(1,N+1))-diag(ones(1,N),1)-diag(ones(1,N),-1);
K(1,2)=-2; % Condicion Neumann
K(N+1,N)=-2; % Condicion Neumann
K=q*K/h^2;

% Definimos F
F=zeros(N+1,1);

% Calculamos u0
u0=[zeros(1,1+N*5/7),3*ones(1,N*2/7)]';
% Guardamos la solucion para t=0 en la matriz solucion
sol(1,:)=u0'; 
U=u0; % Inicializacion

% Calculamos U_n+1 a partir de U_n con un bucle for
for j=1:length(t)-1 % Voy hasta la longitud de t menos 1 porque ya conozco un valor
U=(eye(N+1)+(dt*K)/2)\((eye(N+1)-(K*dt)/2)*U); % Trapecio
sol(j+1,:)= U'; % Guardamos solucion
end
% Calculamos cual es el error del 5%
limite=0.05*(5.85/7);
% Buscamos para que tiempo se alcanza ese valor
for i=1:length(t)
dist=(sqrt(sum((sol(i,:)-(5.85/7)).^2)/(N+1)));
if dist<limite
tiempo=(i-1)*dt
break
end
end

El tiempo que obtenemos es T=16.

3.5 Variación de las condiciones de frontera

3.5.1 Condición limpiador

Pasado un tiempo T = 1, colocamos en el extremo izquierdo un limpiador que elimina el contaminante en ese extremo izquierdo, lo que modifica las condiciones de frontera del problema provocando que una condición de frontera tipo Neumann se convierta en una condición tipo Dirichlet, obteniendo el siguiente sistema:

[math] (P)\left\{\begin{matrix}\\u_t-u_{xx}=0 ,x\in\ (0,7), t\gt0\\u_x(0,t)=0, t\lt1\\u(0,t)=0, t\gt1\\u_x(7,t)=0, t\gt0\\u(x,0)=u_0, x\in\ (0,7)\end{matrix}\right. [/math]

Suponiendo que en el instante inicial se verifica :[math] u(x,0)=\left\{\begin{matrix}\\0, x≤5\\3, x\gt5\end{matrix}\right. [/math]

Se comprueba que está bien propuesto analizando la existencia de solución, la unicidad de ésta y su estabilidad,

Obtenemos una aproximación numérica de la solución utilizando el método de diferencias finitas.

El código en MATLAB del método sería:

clear all
clf
% u_t-qu_xx=0, x en (0,L)
% u_x(0,t)=0 Condicion Neumann para t<=1
% u(0,t)=0 Condicion Dirichlet para t>1
% u_x(L,t)=0 Condicion Neumann
% u(x,0)=u0(x)

% Datos
L=7;
T=7;
q=1;

% Discretizacion espacial
h=0.1;
N=L/h;
x=0:h:L; 
xint=h:h:L; % Vector de nodos interiores

% Discretizacion temporal
dt=h/4; % Tamaño de paso
t=0:dt:T; % Vector de tiempos

% Matriz K para t<=1
K=2*diag(ones(1,N+1))-diag(ones(1,N),1)-diag(ones(1,N),-1);
K(1,2)=-2; % Condicion Neumann
K(N+1,N)=-2; % Condicion Neumann
K=q*K/h^2;
% Matriz K para t>1
K2=2*diag(ones(1,N))-diag(ones(1,N-1),1)-diag(ones(1,N-1),-1);
K2(N,N-1)=-2;
K2=K2/h^2;
% Definimos F
F=zeros(N+1,1);
F2=zeros(N,1);
% Calculamos u0
u0=[zeros(1,1+N*5/7),3*ones(1,N*2/7)]';
% Guardamos la solucion para t=0 en la matriz solucion
sol(1,:)=u0'; 
U=u0; % Inicializacion

% Calculamos U_n+1 a partir de U_n con un bucle for
for j=1:((length(t)-1)/7) 
U=(eye(N+1)+dt*K)\(U+dt*F); % Euler Implicito
sol(j+1,:)= U'; % Guardamos la solucion
end
U(1)=[];
for j=((length(t)-1)/7):length(t)-1 
U=(eye(N)+dt*K2)\(U+dt*F2); % Euler Implicito
sol(j+1,:)= [0,U']; % Guardamos la solucion
end


% Dibujamos la solucion
[xx,tt]=meshgrid(x,t);
surf(xx,tt,sol);

La gráfica obtenida sería:

Resolución del sistema de difusión de una sustancia contaminante por el método del trapecio.

El valor estacionario de la concentración en el tubo será de cero, debido a que se ha colocado un limpiador en uno de sus extremos que irá eliminando el contaminante.

Para averiguar el tiempo que tarda la concentración en alcanzar este otro estado estacionario con un error del 5% vamos a calcular la media cuadrática de las concentraciones y ver cual coincide con el 5% de 5.85/7.

clear all
clf
% u_t-qu_xx=0, x en (0,L)
% u_x(0,t)=0 Condicion Neumann para t<=1
% u(0,t)=0 Condicion Dirichlet para t>1
% u_x(L,t)=0 Condicion Neumann
% u(x,0)=u0(x)
 
% Datos
L=7;
T=70;
q=1;
 
% Discretizacion espacial
h=0.1;
N=L/h;
x=0:h:L; 
xint=h:h:L; % Vector de nodos interiores
 
% Discretizacion temporal
dt=h/4; % Tamaño de paso
t=0:dt:T; % Vector de tiempos
 
% Matriz K para t<=1
K=2*diag(ones(1,N+1))-diag(ones(1,N),1)-diag(ones(1,N),-1);
K(1,2)=-2; % Condicion Neumann
K(N+1,N)=-2; % Condicion Neumann
K=q*K/h^2;
% Matriz K para t>1
K2=2*diag(ones(1,N))-diag(ones(1,N-1),1)-diag(ones(1,N-1),-1);
K2(N,N-1)=-2;
K2=K2/h^2;
% Definimos F
F=zeros(N+1,1);
F2=zeros(N,1);
% Calculamos u0
u0=[zeros(1,1+N*5/7),3*ones(1,N*2/7)]';
% Guardamos la solucion para t=0 en la matriz solucion
sol(1,:)=u0'; 
U=u0; % Inicializacion
 
% Calculamos U_n+1 a partir de U_n con un bucle for
for j=1:(1/dt) 
U=(eye(N+1)+dt*K)\(U+dt*F); % Euler Implicito
sol(j+1,:)= U'; % Guardamos la solucion
end
U(1)=[];
for j=(1/dt):length(t)-1 
U=(eye(N)+dt*K2)\(U+dt*F2); % Euler Implicito
sol(j+1,:)= [0,U']; % Guardamos la solucion
end
 
 % Calculamos cual es el error del 5%
limite=0.05*(5.85/7);
% Buscamos para que tiempo se alcanza ese valor
for i=1:length(t)
dist=(sqrt(sum(sol(i,:).^2)/(N+1)));
if dist<limite
tiempo=(i-1)*dt
break
end
end

De esta forma el tiempo que obtenemos es 65.775

3.5.2 Condiciones tipo Newman

las nuevas condiciones newmann no son homogéneas. tendremos que homogeneizarlas con el siguiente cambio de variable :

[math] u(x,t)=w(x,t)+2a(t)x^2+b(t)x [/math]

el nuevo problema quedaría de la siguiente manera :

[math] (P)\left\{\begin{matrix}\\w_t-w_{xx}=-3/7 ,x\in\ (0,7), t\gt0\\w_x(0,t)=0, t\gt0\\w_x(7,t)=-3, t\gt0\\w(x,0)=u_0, x\in\ (0,7)\end{matrix}\right. [/math]

Suponiendo que en el instante inicial se verifica :[math] w(x,0)=\left\{\begin{matrix}\\0, x≤5\\3, x\gt5\end{matrix}\right. [/math]

Se comprueba que es un problema bien propuesto analizando su existencia de solución, unicidad de ésta y estabilidad con respecto al dato inicial.

ahora ya se trata de un sistema con condiciones de frontera homogéneas, el problema de autovalores asociado es:

[math] (PA)\left\{\begin{matrix}\\x''+μx=0 ,x\in\ (0,7)\\x'(0)=0\\x'(7)=0\end{matrix}\right. [/math]

Y los autovalores y autofunciones salvo una constante son:

[math]μ_k=\frac{\pi^2k^2}{49}[/math]: [math]φ_k=cos(kπ/7)x[/math] k=0,1,2,3...

Se aproxima el dato inicial y el segundo miembro, se elige Q y se toma:

[math]\displaystyle\sum_{k=0}^Q c_kφ_k(x)=-3/7[/math]:

[math]\displaystyle\sum_{k=0}^Q a_kφ_k(x)=\left\{\begin{matrix}\\0, x≤5\\3, x\gt5\end{matrix}\right.[/math]

donde se deduce que [math]c_k=\frac{\displaystyle\int_{0}^{7} -3/7cos(kπ/7)x\, dx}{\displaystyle\int_{0}^{7} cos^2(kπ/7)x\, dx}[/math]

[math]a_k=\frac{\displaystyle\int_{0}^{7} u_0cos(kπ/7)x\, dx}{\displaystyle\int_{0}^{7} cos^2(kπ/7)x\, dx}[/math]

El problema aproximado queda de la forma:

[math] (P)\left\{\begin{matrix}\\w_t-w_{xx}=\displaystyle\sum_{k=0}^Q c_kφ_k(x) ,x\in\ (0,7), t\gt0\\w_x(0,t)=0, t\gt0\\w_x(7,t)=0, t\gt0\\w(x,0)=w_0=\displaystyle\sum_{k=0}^Q a_kφ_k(x), x\in\ (0,7)\end{matrix}\right. [/math]

La solución aproximada sería:

[math]w(x,t)=\displaystyle\sum_{k=0}^Q T_k(t)φ_k(x)[/math]

Sustituyendo en la ecuación de difusión de una sustancia contaminante e igualando los coeficicientes de Fourier se obtiene el sistema:

[math]\left\{\begin{matrix}\\T'_k(t)+μ_kT_k(t)=c_k\\T_k(0)=a_k\end{matrix}\right.[/math]

Resolviendo el problema de valor inicial, se deduce que [math] T_k(t)=Ce^{-μ_k t}[/math]

La solución aproximada queda de la forma [math]w(x,t)=\displaystyle\sum_{k=0}^Q Ce^{-μ_k t} cos(kπ/7)x[/math]

El valor de la constante C depende de la condición inicial,que es la que condiciona el valor de a_k.

por ultimo solo nos quedaria deshacer el cambio de variable.

Difusión de una sustancia contaminante (Grupo 24C)