Difusión de una sustancia contaminante (Grupo 3B)

Trabajo realizado por estudiantes
Título	Difusión de una sustancia contaminante. Grupo 3B
Asignatura	Ecuaciones Diferenciales
Curso	Curso 2013-14
Autores	María Bartol Calderón Rodrigo Bellot Rodríguez Margarita Santiago Ruiz Rocío Santos Rodrigo
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 PLANTEAMIENTO DEL PROBLEMA

1.1 Modelización

Se considera una solución en un tubo largo compuesta por dos sustancias de las cuales una de ellas es un contaminante. Se orienta el tubo en la dirección del eje x, desde x=0 hasta x=L, siendo L la longitud del tubo, cuyo valor es igual a 5 metros. Dicho tubo presenta sección transversal constante.

Se denota por u la concentración de contaminante en cada posición del tubo, medida en mol/m²s. La concentración es la misma en cualquier punto de la sección transversal del tubo. Por tanto, va a depender únicamente de dos variables u=u(x,t).

Se designa por F(x,t) el flujo de difusión del contaminante, siendo éste la cantidad de sustancia contaminante que fluye por unidad de tiempo y unidad de área. La superficie del tubo está aislada, por lo que solamente hay flujo de contaminante en la dirección del eje x. En los extremos se coloca un aislante que hace que F(x,t) al exterior del tubo sea nulo.

Dado que la sección del tubo es constante, la concentración del contaminante va a depender únicamente del tiempo y el número de moles, siendo éste último proporcional a la masa de la sustancia.

Teniendo en cuenta el principio de conservación de masa, se deduce, que 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. Al ser nulo el último sumando, la variación de u va a estar influenciada únicamente por el flujo de contaminante.

La ley de Flick establece 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]

donde D es el coeficiente de difusión (medido en m²/s), que depende de las propiedades químicas de los compuestos, y que se supone constante D=1.

1.2 Ecuación de difusión

Utilizando la ley de Flick y aplicando el principio de conservación de masa, se deduce que la concentración de la sustancia contaminante en una sección del tubo que dista x del extremo x=0, cuando pasa un tiempo t debe satisfacer:

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

Imponiendo que D es el coeficiente de difusión y su valor es la unidad, la ecuación de difusión que describe cómo la sustancia contaminada se dispersa, queda de la forma:

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

Para la deducción de dicha fórmula, se tienen en cuenta las siguientes consideraciones.

Se define [math] A [/math] como la variable de superficie, y se toma 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]

Y al derivar en el tiempo:

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

De esta forma se obtiene el primer término de la ecuación de difusión.

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.

Suponiendo 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 flujo 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:

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

1.3 Problema bien propuesto

Dado el problema

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

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

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

1.3.1 Existencia de solución

Se debe asegurar que el problema admite una solución u(x,t) mediante el método de separación de variables. Para ello se realizan los siguientes pasos:

1. Reducción a un problema con condiciones de frontera homogéneas

Para poder aplicar el método de separación de variables las condiciones de frontera del problema deben ser homogéneas. En este caso no es necesario realizar ningún cambio de variable para conseguirlo.

2. Resolución del problema de autovalores asociado a (P)

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

Los autovalores y autofunciones salvo una constante son:

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

3. Cálculo de la serie de Fourier de f(x,t)=0 en términos de las autofunciones del problema de autovalores asociado

[math]S f(x,t)= \displaystyle\sum_{k=0}^∞ c_kφ_k(x)=0\cdotφ_1(x)+0\cdotφ_2(x)+0\cdotφ_3(x)+...[/math]

4. Ensayo con soluciones de la forma [math]u(x,t)=\displaystyle\sum_{k=0}^∞ T_k(t)φ_k(x)[/math]

Ya que [math]φ'_k(0)=0[/math] y [math]φ'_k(5)=0[/math], cualesquiera que sean las funciones [math]T_k(t)[/math], u(x,t) satisface las condiciones de frontera.

Impongo ahora que satisfaga la ecuación en derivadas parciales [math]u_t(x,t)- u_{xx}(x,t)= 0 [/math]:

[math]\displaystyle\sum_{k=0}^∞ T'_k(t)φ_k(x)+\displaystyle\sum_{k=0}^∞ T_k(t)μ_kφ_k(x)=0[/math]: [math]\displaystyle\sum_{k=0}^∞ [T'_k(t)+μ_kT_k(t)]φ_k(x)=0[/math]

Es decir, [math]\displaystyle\sum_{k=0}^∞ [T'_k(t)+μ_kT_k(t)]φ_k(x)=\displaystyle\sum_{k=0}^∞ c_kφ_k(x)[/math]:

Por la unicidad de los coeficientes de Fourier se cumple [math]T'_k(t)+μ_kT_k(t)=0[/math], t>0, k=0,1,2,3...

Por lo tanto si [math]T_k(t)[/math] es una solución cualquiera de la ecuación anterior, la función [math]u(x,t)=\displaystyle\sum_{k=0}^∞ T_k(t)cos(kπ/5)x[/math] satisface la ecuación en derivadas parciales y las condiciones de frontera de (P).

5. La solución [math]u(x,t)=\displaystyle\sum_{k=0}^∞ T_k(t)φ_k(x)[/math] debe satisfacer también la condición inicial

[math]u(x,0)=\displaystyle\sum_{k=0}^∞ T_k(0)φ_k(x)=\left\{\begin{matrix}\\0, x≤3\\3, x\gt3\end{matrix}\right.[/math]

Calculando la serie de Fourier de la condición inicial y aplicando la unicidad de los coeficientes de Fourier se obtiene el valor de [math]T_k(0)[/math].

Si [math]T_k(t) k=0,1,2... [/math] es solución del problema de valor incial, formado por [math]T'_k(t)+μ_kT_k(t)=0[/math] y [math]T_k(0)[/math], entonces:

[math]u(x,t)=\displaystyle\sum_{k=0}^∞ T_k(t)cos(kπ/5)x,   xɶ(0,5), t\gt0[/math] es solución de (P).

1.3.2 Unicidad de solución

Se utiliza el método de la energía o método de los multiplicadores para demostrar que si existe una solución del problema debe ser única.

Por la linealidad de (P) se verifica que (P) admite una única solución si y solamente si (PH) admite como única solución la trivial, siendo

[math] (PH)\left\{\begin{matrix}\\u_t-u_{xx}=0 ,xɶ(0,5), t\gt0\\u_x(0,t)=0, t\gt0\\u_x(L,t)=0, t\gt0\\u(x,0)=0, xɶ(0,5)\end{matrix}\right. [/math]

Sea w(x,t) una solución de (PH), se demuestra la unicidad de solución de (P) comprobando que w(x,t)=0.

Si w(x,t) satisface el problema :[math] (PH)\left\{\begin{matrix}\\w_t-w_{xx}=0 ,xɶ(0,5), t\gt0\\w_x(0,t)=0, t\gt0\\w_x(L,t)=0, t\gt0\\w(x,0)=0, xɶ(0,5)\end{matrix}\right. [/math]

entonces [math]w_t(x,t)-w_{xx}(x,t)=0 ,xɶ(0,5), t\gt0[/math]

Multiplicamos la ecuación en derivadas parciales por w(x,t) e integramos en la variable x en (0,5):

[math]\displaystyle\int_{0}^{5} w_t(x,t)w(x,t)\, dx - \displaystyle\int_{0}^{5} w_{xx}(x,t)w(x,t)\, dx =0 t\gt0 [/math]

Haciendo integración por partes en la segunda integral anterior, y aplicando que [math]w(x,t)[/math] satisface [math]w(0,t)=w(5,t)=0[/math], se obtiene:

[math]\displaystyle\int_{0}^{5} w_t(x,t)w(x,t)\, dx + \displaystyle\int_{0}^{5} w_x^2(x,t)\, dx =0 t\gt0 [/math]

Que implica:

[math]\displaystyle\int_{0}^{5} w_t(x,t)w(x,t)\, dx= - \displaystyle\int_{0}^{5} w_x^2(x,t)\, dx =0 t\gt0 [/math]

Recordando que:

[math]d/\partial t (w_x^2(x,t))=2w_t(x,t)w(x,t)[/math]

y que:

[math] - \displaystyle\int_{0}^{5} w_x^2(x,t)\, dx≤0,  t\gt0[/math]

se halla la función:

[math] E(t)= \displaystyle\int_{0}^{5} w^2(x,t)\, dx[/math]

Derivando la función comprobamos que se verifica que [math]E'(t)≤0[/math], es decir, que [math] E(t)[/math] es decreciente. Como [math]w_x^2(x,t)≥0[/math], entonces [math]E(t)≥0, t≥0 [/math]

Y ya que [math]w(x,0)=0[/math],:

[math] E(0)= \displaystyle\int_{0}^{5} w^2(x,0)\, dx = \displaystyle\int_{0}^{5} 0\, dx[/math]

Por ser [math] E(t)[/math] decreciente, [math]0 ≤ E(t) ≤ E(0) =0[/math]

Lo anterior implica que:

[math] E(t)= \displaystyle\int_{0}^{5} w^2(x,t)\, dx = 0[/math]

Ya que [math]w^2(x,t)≥0[/math] la igualdad anterior sólo se puede dar si [math]w^2(x,t)=0[/math] y por lo tanto [math]w(x,t)≥0[/math]

1.3.3 Estabilidad

Un problema estable facilita su tratamiento numérico, por ello, se debe comprobar que sea estable con respecto al dato inicial.

Si se considera el problema

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

Se debe demostrar que u₀ y h(x) son "próximas", entonces las soluciones de (P) y de (P') serán "próximas". Si u(x,t) y v(x,t) son las soluciones de (P) y (P'), entonces se verifica:

[math] sup_{t≥0}\displaystyle\int_{0}^{5} \left | u(x,t)-v(x,t) \right |\ ^2 dx ≤ C sup_{x\in\ [0,5]} \left | u_0 - h(x) \right | \ ^2[/math]

Se demuestra la inecuación utilizando el método de la energía.

[math] w(x,t)= u(x,t)-v(x,t)[/math] es solución del problema

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

y por tanto satisface

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

Multiplicamos la ecuación en derivadas parciales por w(x,t) e integramos en la variable x en (0,5):

[math]\displaystyle\int_{0}^{5} w_t(x,t)w(x,t)\, dx - \displaystyle\int_{0}^{5} w_{xx}(x,t)w(x,t)\, dx =0 t\gt0 [/math]

Haciendo integración por partes en la segunda integral anterior, y aplicando que [math]w(x,t)[/math] satisface [math]w(0,t)=w(5,t)=0[/math], se obtiene:

[math]\displaystyle\int_{0}^{5} w_t(x,t)w(x,t)\, dx + \displaystyle\int_{0}^{5} w_x^2(x,t)\, dx =0 t\gt0 [/math]

Que implica:

[math]\displaystyle\int_{0}^{5} w_t(x,t)w(x,t)\, dx= - \displaystyle\int_{0}^{5} w_x^2(x,t)\, dx =0 t\gt0 [/math]

Recordando que:

[math]d/\partial t (w_x^2(x,t))=2w_t(x,t)w(x,t)[/math]

y que:

[math] - \displaystyle\int_{0}^{5} w_x^2(x,t)\, dx≤0,  t\gt0[/math]

se halla la función:

[math] E(t)= \displaystyle\int_{0}^{5} w^2(x,t)\, dx[/math]

Derivando la función comprobamos que se verifica que [math]E'(t)≤0[/math], es decir, que [math] E(t)[/math] es decreciente. Como [math]w_x^2(x,t)≥0[/math], entonces [math]E(t)≥0, t≥0 [/math]

Y ya que [math]w(x,0)=0[/math],:

[math] E(0)= \displaystyle\int_{0}^{5} w^2(x,0)\, dx = \displaystyle\int_{0}^{5} 0\, dx[/math]

Por ser [math] E(t)[/math] decreciente, [math]0 ≤ E(t) ≤ E(0) =0[/math]

Lo anterior implica que:

[math] E(t)= \displaystyle\int_{0}^{5} w^2(x,t)\, dx = 0[/math]

o lo que es lo mismo:

[math] \displaystyle\int_{0}^{5} \left | u(x,t)-v(x,t) \right |\ ^2 dx ≤ \displaystyle\int_{0}^{5} \left | u_0 - h(x) \right |\ ^2 dx ≤ sup_{x\in\ [0,5]} \left | u_0 - h(x) \right | \ ^2[/math]

Dado que el lado derecho de la expresión no depende de t, tomando supremos en t se obtiene:

[math] sup_{t≥0}\displaystyle\int_{0}^{5} \left | u(x,t)-v(x,t) \right |\ ^2 dx ≤ sup_{x\in\ [0,5]} \left | u_0 - h(x) \right | \ ^2[/math]

que es lo que inicialmente se pide verificar.

Por último, se debe tener en cuenta que el resultado de estabilidad es más fuerte que el de unicidad, en el sentido de que aquel implica este.

2 RESOLUCIÓN DEL PROBLEMA

Para obtener una aproximación de la solución del problema planteado, se resuelve numéricamente dicho problema empleando dos métodos:

• Método de diferencias finitas.
  
• Método de Fourier.
  

2.1 Método de diferencias finitas

El método de diferencias finitas proporciona una aproximación de la solución numérica de un sistema continuo.

Para aplicar este método en primer lugar se realiza una partición del intervalo [0,5] definiendo h como la anchura de paso.

En segundo lugar, se aplica la ecuación diferencial a un nodo interior Xn, sustituyendo la derivada u_xx por esta aproximación de segundo orden::

[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]

Teniendo en cuenta esta aproximación, y la notación [math] u_n(t) = u(x_n,t)[/math] , resulta el siguiente sistema de N+1 ecuaciones, obtenido al aplicar la ecuación diferencial para cada n:

\begin{array}{c}u'_n(t)+\frac{-u_{n-1}(t)+2u_n(t)+u_{n+1}(t)}{h^2}=0\end{array} n=1,2,3...,N-1

Al aplicar las condiciones de contorno se reducen a N-1 ecuaciones:

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

\begin{array}{c}U'+KU=F=0\\U(0)=U^0\end{array}

Siendo [math]U^0[/math] las condiciones iniciales dadas.

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

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,5][/math], con [math]\triangle x=0.1[/math]

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

% Aproximar la ecuacion de difusion de una sustancia contaminante
% 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 del problema
L=5;
T=5;
q=1;
% Datos de la discretizacion espacial
N=50;
h=L/N;
% Vector de nodos en en espacio
x=0:h:L; 
% Aproximacion de -q*u_xx
K=2*diag(ones(1,N+1))-diag(ones(1,N),1)-diag(ones(1,N),-1);
 K(1,2)=-2;
 K(N+1,N)=-2; %?¿?
 K=q*K/h^2;
 % Calculamos u0
 u0=[zeros(1,1+N*3/5),3*ones(1,N*2/5)]';
 % Discretizacion temporal
dt=h/4; % Paso
t=0:dt:T; % Vector de tiempos
% Definimos F
F=zeros(N+1,1); %sin nodos interiores
% Guardamos la solucion
sol(1,:)=u0'; 
U=u0; % Inicializacion
% Calculamos U_n ----> U_n+1
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.

En segundo lugar, hallamos la solución numérica del problema utilizando el método de Euler explícito, el método de Euler implícito y por último el método de Euler modificado, y comparamos los diferentes resultados obtenidos.

El código en MATLAB de los tres métodos sería:

MÉTODO DE EULER EXPLÍCITO

% Aproximar la ecuacion de la difusión de una sustancia contaminante
% 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 del problema
L=5;
T=5;
q=1;
% Datos de la discretizacion espacial
N=50;
h=L/N;
% Vector de nodos en en espacio
x=0:h:L; 
% Aproximacion de -q*u_xx
K=2*diag(ones(1,N+1))-diag(ones(1,N),1)-diag(ones(1,N),-1);
 K(1,2)=-2;
 K(N+1,N)=-2; %?¿?
 K=q*K/h^2;
 % Calculamos u0
 u0=[zeros(1,1+N*3/5),3*ones(1,N*2/5)]';
 % Discretizacion temporal
dt=h/20; % Paso para euler explicito
t=0:dt:T; % Vector de tiempos
% Definimos F
F=zeros(N+1,1); %sin nodos interiores
% Guardamos la solucion
sol(1,:)=u0'; 
U=u0; % Inicializacion
% Calculamos U_n ----> U_n+1
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:

Archivo:Método Euler Explícito.png

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

MÉTODO DE EULER IMPLÍCITO

% Aproximar la ecuacion de difusion de una sustancia contaminante
% 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 del problema
L=5;
T=5;
q=1;
% Datos de la discretizacion espacial
N=50;
h=L/N;
% Vector de nodos en en espacio
x=0:h:L; 
% Aproximacion de -q*u_xx
K=2*diag(ones(1,N+1))-diag(ones(1,N),1)-diag(ones(1,N),-1);
 K(1,2)=-2;
 K(N+1,N)=-2; %?¿?
 K=q*K/h^2;
 % Calculamos u0
 u0=[zeros(1,1+N*3/5),3*ones(1,N*2/5)]';
 % Discretizacion temporal
dt=h/4; % Paso
t=0:dt:T; % Vector de tiempos
% Definimos F
F=zeros(N+1,1); %sin nodos interiores
% Guardamos la solucion
sol(1,:)=u0'; 
U=u0; % Inicializacion
% Calculamos U_n ----> U_n+1
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)\(U+dt*F); % metodo de Euler implicito
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 de Euler Implícito.

MÉTODO DE EULER MODIFICADO

% Aproximar la ecuacion de difusion
% 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 del problema
L=5;
T=5;
q=1;
% Datos de la discretizacion espacial
N=50;
h=L/N;
% Vector de nodos en en espacio
x=0:h:L; 
% Aproximacion de -q*u_xx
K=2*diag(ones(1,N+1))-diag(ones(1,N),1)-diag(ones(1,N),-1);
 K(1,2)=-2;
 K(N+1,N)=-2; %?¿?
 K=q*K/h^2;
 % Calculamos u0
 u0=[zeros(1,1+N*3/5),3*ones(1,N*2/5)]';
 % Discretizacion temporal
dt=h/20; % Paso
t=0:dt:T; % Vector de tiempos
% Definimos F
F=zeros(N+1,1); %sin nodos interiores
% Guardamos la solucion
sol(1,:)=u0'; 
U=u0; % Inicializacion
% Calculamos U_n ----> U_n+1
for j=1:length(t)-1 % Voy hasta la longitud de t menos 1 porque ya conozco un valor
% Calculamos k1=h*f(tn,yn)
k1=-K*U;
% Calculamos k2=h*f(tn+h/2,yn+1/2*k1)
k2=-K*(U+dt*k1);
% Meto los datos en el vector
U=U+(dt/2)*(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 Euler Modificado.

2.1.1 Conservación de la masa total de contaminante

Se deduce a partir de la ecuación u(x,t) 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
% Aproximar la ecuacion de la difusion de una sustancia contaminante
% 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 del problema
L=5;
T=5;
q=1;
% Datos de la discretizacion espacial
N=500;
h=L/N;
% Vector de nodos en en espacio
x=0:h:L; 
% Aproximacion de -q*u_xx
K=2*diag(ones(1,N+1))-diag(ones(1,N),1)-diag(ones(1,N),-1);
 K(1,2)=-2;
 K(N+1,N)=-2; %?¿?
 K=q*K/h^2;
 % Calculamos u0
 u0=[zeros(1,1+N*3/5),3*ones(1,N*2/5)]';
 % Discretizacion temporal
dt=h/4; % Paso
t=0:dt:T; % Vector de tiempos
% Definimos F
F=zeros(N+1,1); %sin nodos interiores
% Guardamos la solucion
sol(1,:)=u0'; 
U=u0; % Inicializacion
% Calculamos U_n ----> U_n+1
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)\(U+dt*F); % metodo implicito
sol(j+1,:)= U'; % Guardamos solucion
end
% Calculamos el area de la funcion en cada tiempo (area = cantidad de contaminante)
areas=zeros(1,length(t));
for i=1:length(t)
f=sol(i,:);
aa=trapz(x,f);
area(i)=aa;
end
% Dibujamos la cantidad de contaminante a lo largo del tiempo
plot(t,area)
xlabel('Tiempo')
ylabel('Cantidad de contaminante')

En la gráfica obtenida se observa cómo la masa total de contaminante se conserva a lo largo del tiempo.

Conservación masa total de contaminante.

Podemos ver el comportamiento de la concentración en el punto medio del tubo en la siguiente gráfica:

Concentración de contaminante en el punto medio del tubo.

En la siguiente gráfica se observa que para tiempos grandes el contaminante se va extendiendo de manera uniforme hasta tener concentración constante. Además, la solución del problema debe parecerse a la función u(x,t)=1.2 y cumplir que en el estado estacionario u_x(x,t)=0.

concentración para tiempos grandes.

2.1.2 Variación de las condiciones de frontera

Se coloca un limpiador que elimina el contaminante en el extremo izquierdo del tubo, lo que nos produce un cambio en las condiciones de frontera del problema, siendo necesario modelizar la nueva situación. Esta variación provoca 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,5), t\gt0\\u(0,t)=0, t\gt0\\u_x(5,t)=0, t\gt0\\u(x,0)=u_0, x\in\ (0,5)\end{matrix}\right. [/math]

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

Se comprueba que está bien propuesto analizando la existencia de solución, la unicidad de ésta y su estabilidad, de la misma forma que se hizo anteriormente.

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
% Aproximar la ecuacion de difusion de una sustancia contaminante
% u_t-u_xx=0, x en (0,L)
% u(0,t)=0 Condicion Dirichlet
% u_x(L,t)=0 Condicion Neumann
% u(0,t)=u0(x)
%%%
% Datos del problema
L=5;
T=5;
% Datos de la discretizacion
N=50;
h=L/N;
x=0:h:L; % Vector de nodos
xint=h:h:L; % Vector de nodos interiores 
% Discretizacion temporal. 
dt=h/4; 
t=0:dt:T;
% Aproximacion de -*u_xx
K=2*diag(ones(1,N))-diag(ones(1,N-1),1)-diag(ones(1,N-1),-1);
K(N,N-1)=-2;
K=K/h^2;
% Calculamos U0 (condicion en el instante 0)
U0=[zeros(1,N*3/5),3*ones(1,N*2/5)]';
% Definimos F
F=zeros(N,1);
% Guardamos la solucion
sol(1,:)=[0,U0']; 
U=U0; % Inicializacion
% Calculamos U_n ----> U_n+1
for j=1:length(t)-1 % voy hasta la longitud de t menos 1 porque ya conozco un valor
U=(eye(N)+dt*K)\(U+dt*F);
sol(j+1,:)=[0,U']; % Guardamos la solucion
end
% Dibujamos
[xx,tt]=meshgrid(x,t);
surf(xx,tt,sol);

Solución con la nueva condición de frontera.

2.2 Método de Fourier

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,5), t\gt0\\u_x(0,t)=0, t\gt0\\u_x(5,t)=0, t\gt0\\u(x,0)=u_0, x\in\ (0,5)\end{matrix}\right. [/math]

Cuyas condiciones iniciales son: [math] u(x,0)=\left\{\begin{matrix}\\0, x≤3\\3, x\gt3\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}\\y''+μy=0 ,x\in\ (0,5)\\y'(0)=0\\y'(5)=0\end{matrix}\right. [/math]

Y los autovalores y autofunciones salvo una constante son:

[math]μ_k=\frac{\pi^2k^2}{25}[/math]: [math]φ_k=cos(kπ/5)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≤3\\3, x\gt3\end{matrix}\right.[/math]

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

El problema aproximado queda de la forma:

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

La solución aproximada sería:

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

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

[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 obtiene que [math] T_k(t)=Ce^{-μ_k t}[/math]

La solución aproximada tendría la forma [math]u(x,t)=\displaystyle\sum_{k=0}^∞ Ce^{-μ_k t} cos(kπ/5)x[/math]

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