Revisión actual del 23:11 22 may 2014

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.

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

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

Teniendo en cuenta:

[math]u'_n(t)\simeq\frac{u_{n+1}(t)-u_{n-1}(t)}{h2}[/math]

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=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); 
% 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);
% 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:

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); 
% 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); 
% 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.

Al observar las gráficas vemos que utilizando los distintos métodos obtenemos unas aproximaciones muy parecidas. Todas tienen la misma forma y solo varía la precisión a la hora de ajustarse a la curva, en función del orden del método. El método más preciso es el de Euler Implícito como se demuestra gráficamente y el menos el 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); 
% 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.

Se puede comprobar observando en la siguiente gráfica 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:

[math]u_{xx}(x,t)=0[/math]: [math]u_x(0,t)=0[/math]: [math]u_x(5,t)=0[/math]: [math]u(x,0)=1.2[/math]:

concentración para tiempos grandes.

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

clear all
% 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=10;
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);
% 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

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

De donde se obtiene: \begin{array}{c}dist0=1.4647\\dist1=0.8711\\dist2=0.5854\\dist10=0.0253\end{array}

Para averiguar el tiempo que tarda la concentración en alcanzar el estado estacionario con un error del 5% vamos a calcular la media cuadrática de las distancias desde cada punto de la función al estado estacionario: 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 la fórmula: [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] Haremos esto para cada tiempo t hasta que esta media cuadrática sea menor del 5% del valor estacionario.

clear all
% 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=10;
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); 
% 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 cual es el error del 5%
limite=0.05*1.17;
% Buscamos para que tiempo se alcanza ese valor
for i=1:length(t)
dist=(sqrt(sum((sol(i,:)-1.17).^2)/(N+1)));
if dist<limite
tiempo=(i-1)*dt
break
end
end

Se obtiene un tiempo de 7.8750 segundos.

Al dividir "incremento de x" por 10 (usando N=500) sale un tiempo de 8.1325 segundos.

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);

La gráfica obtenida sería:

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

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 utilizar 2 métodos distintos. Primero vamos a calcular el tiempo que tarda la cantidad de contaminante en reducirse hasta el 5%, y después calculando la media cuadrática igual que hemos hecho en el apartado anterior.

clear all
% Aproximar la ecuacion de difusion
% u_t-u_xx=0, x en (0,L)
% u(0,t)=0 Condicion Dirichlet
% u_x(L,t)=0 Condicion Neumann
% u(x,0)=u0(x)
%%%
% Datos del problema
L=5;
T=40;
% 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
% Calculamos cual es el error del 5%
limite=0.05*5.85;
% Calculamos el area de la funcion en cada tiempo (area = cantidad de contaminante)
for i=1:length(t)
f=sol(i,:);
area=trapz(x,f);
% Buscamos para que tiempo se alcanza ese valor
if area<limite
tiempo=(i-1)*dt
break
end
end

De la primera forma se obtiene un tiempo de 32.225

clear all
% Aproximar la ecuacion de difusion
% u_t-u_xx=0, x en (0,L)
% u(0,t)=0 Condicion Dirichlet
% u_x(L,t)=0 Condicion Neumann
% u(x,0)=u0(x)
%%%
% Datos del problema
L=5;
T=80;
% 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

% Calculamos cual es el error del 5%
limite=0.05*1.17;
% 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 la segunda forma se obtiene un tiempo de 33.275

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}^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π/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.

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

clear all
% Aproximar la ecuacion del difusion con el metodo de Fourier
% u_t-u_xx=0
% u_x(0,t)=0
% u_x(L,t)=0
% u(x,0)=u0(x)
% Datos del problema
L=5;
T=0.5;
% Datos de la aproximacion
Q=1; % Numero de coeficientes de Fourier
N=50; % Numero de intervalos en x
h=L/N;
x=0:h:L;
dt=h/4;
t=0:dt:T;
% Calculamos la aproximacion
sol=zeros(length(t),N+1);
u0=[zeros(1,1+N*3/5),3*ones(1,N*2/5)];
for k=0:1
phi=cos((k*pi/5)*x);
a=(trapz(x,u0.*phi))/(trapz(x,phi.^2));
T=a*exp(-(k^2)*(pi^2)*1/25*t);
sol=sol+T'*phi;
end
[xx,tt]=meshgrid(x,t);
subplot(2,2,1);
surf(xx,tt,sol)

for k=0:3
phi=cos((k*pi/5)*x);
a=(trapz(x,u0.*phi))/(trapz(x,phi.^2));
T=a*exp(-(k^2)*(pi^2)*1/25*t);
sol=sol+T'*phi;
end
[xx,tt]=meshgrid(x,t);
subplot(2,2,2);
surf(xx,tt,sol)


for k=0:5
phi=cos((k*pi/5)*x);
a=(trapz(x,u0.*phi))/(trapz(x,phi.^2));
T=a*exp(-(k^2)*(pi^2)*1/25*t);
sol=sol+T'*phi;
end
[xx,tt]=meshgrid(x,t);
subplot(2,2,3);
surf(xx,tt,sol)

for k=0:10
phi=cos((k*pi/5)*x);
a=(trapz(x,u0.*phi))/(trapz(x,phi.^2));
T=a*exp(-(k^2)*(pi^2)*1/25*t);
sol=sol+T'*phi;
end
[xx,tt]=meshgrid(x,t);
subplot(2,2,4);
surf(xx,tt,sol)

Las gráficas obtenidas serían:

Solución por Fourier para 1,2,3,5 y 10 términos.

El código en MATLAB del método de Fourier con 20 términos de la serie sería:

clear all
% Aproximar la ecuacion del difusion con el metodo de Fourier
% u_t-u_xx=0
% u_x(0,t)=0
% u_x(L,t)=0
% u(x,0)=u0(x)
% Datos del problema
L=5;
T=0.5;
% Datos de la aproximacion
Q=20; % Numero de coeficientes de Fourier
N=50; % Numero de intervalos en x
h=L/N;
x=0:h:L;
dt=h/4;
t=0:dt:T;
% Calculamos la aproximacion
sol=zeros(length(t),N+1);
u0=[zeros(1,1+N*3/5),3*ones(1,N*2/5)];
for k=0:Q
phi=cos((k*pi/5)*x);
a=(trapz(x,u0.*phi))/(trapz(x,phi.^2));
T=a*exp(-(k^2)*(pi^2)*1/25*t);
sol=sol+T'*phi;
end
[xx,tt]=meshgrid(x,t);
surf(xx,tt,sol)

Y su gráfica correspondiente:

Solución por Fourier para 20 términos.

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

Solución por Fourier para T=5.

2.3 Variación de las condiciones de frontera

Dado el problema

[math] (P)\left\{\begin{matrix}\\u_t-u_{xx}=0 ,x\in\ (0,5), t\gt0\\u_x(0,t)=3, t\gt0\\u(5,t)=10sent, 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 es un problema bien propuesto analizando su existencia de solución, unicidad de ésta y estabilidad con respecto al dato inicial.

Para obtener una aproximación de la solución del problema planteado, se resuelve numéricamente dicho problema empleando el método de diferencias finitas (método de Euler implícito).

El código en MATLAB sería:

clear all
% Aproximar la ecuacion de difusion de una sustancia contaminante
% u_t-u_xx=0, x en (0,L)
% u_x(0,t)=3 Condicion Neumann
% u(L,t)=g(t) Condicion Dirichlet
% u(x,0)=u0(x)
%%%
% Datos del problema
L=5;
T=10;
% Datos de la discretizacion
N=50;
h=L/N;
x=0:h:L; % Vector de nodos
xint=h:h:L-h; % Vector de nodos interiores 
% Discretizacion temporal. 
dt=h^2/2; 
t=0:dt:T;
% Funcion g(t)
g=10*sin(t);
% Aproximacion de -*u_xx
K=2*diag(ones(1,N))-diag(ones(1,N-1),1)-diag(ones(1,N-1),-1);
K(1,2)=-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);
F(1)=-3*2/h;
% Guardamos la solucion
sol(1,:)=[U0' g(1)]; 
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
F(N)=g(j)/h^2;
U=(eye(N)+dt*K)\(U+dt*F);
sol(j+1,:)=[U' g(j+1)]; % Guardamos la solucion
end
% Dibujamos
[xx,tt]=meshgrid(x,t);
surf(xx,tt,sol);

Y la gráfica obtenida es la siguiente:

Solución con las nuevas condiciones de frontera.