Difusión de un contaminante. Grupo 4

Trabajo realizado por estudiantes
Título	Difusión de una sustancia contaminante. Grupo 4
Asignatura	Ecuaciones Diferenciales
Curso	Curso 2013-14
Autores	Sandra Carrillo del Cura 81, Sergio Castillo Herrero 85, Andrea García Prieto 171, Patricia González Peinado 198, Adrián Salas Calvo 385
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

El proceso irreversible por el cual un grupo de partículas se distribuye de manera uniforme en un medio se denomina difusión, proceso estadísticamente predecible en conjunto siendo totalmente aleatorio el movimiento de cada partícula aislada. La difusión de moléculas de una sustancia disuelta en otra será nuestro objeto de estudio, para lo cual disponemos de un tubo que será detallado posteriormente.

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

1 Modelización del problema

Consideramos un tubo de longitud [math]L[/math] de un cierto material, con sección transversal [math]A[/math] constante. La orientamos en la dirección del eje [math]x[/math], de [math]x = 0[/math] a [math]x = L[/math]. Consideramos una solución compuesta por dos sustancias, de las cuales una de ellas es un contaminante. La sustancia se propaga por el interior de dicho tubo.

Suponemos que el tubo es delgado y que su superficie lateral es aislante, de forma que el contaminante no sale de él, ni puede entrar desde el exterior. Podemos entonces pensar que la concentración es constante a lo largo de cada sección transversal, y ver el tubo como un objeto unidimensional. La concentración va a depender entonces de [math]x[/math] y [math]t[/math]. Designamos por [math]u(x, t)[/math] la concentración de contaminante en la sección de la varilla que dista [math]x ≥ 0[/math] del extremos [math]x = 0[/math] cuando ha pasado un tiempo [math]t ≥ 0[/math].

Tomemos un trozo de tubo entre las secciones [math]x[/math] y [math]x + ∆x[/math], que designaremos por [math][x, x + ∆x][/math]. Suponemos que ∆x es una cantidad muy pequeña.

Designamos por [math]φ(x, t)[/math] el flujo de contaminante, es decir, la cantidad de contaminante que fluye por unidad de tiempo y unidad de área. Como la superficie lateral del tubo está aislada, solamente habría flujo en la dirección del eje [math]X[/math]. Si [math]φ (x, t) \gt 0[/math] pensamos que el flujo va hacia la derecha, si [math]φ (x, t) \lt0[/math] va hacia la izquierda. El flujo de contaminante en [math][x, x + ∆x][/math] está dado por [math]φ (x, t)A − φ (x + ∆x, t)A[/math].

Designamos por [math]f(x, t)[/math] a la concentración de soluto generada por posibles fuentes o sumideros, por unidad de volumen y unidad de tiempo. En el caso de este problema no existen ni fuentes ni sumideros, por lo tanto [math]f(x)=0[/math].

La ley de Fick (análoga a la de Fourier para la temperatura) dice: [math]F=-D \frac{∂u}{∂x}[/math]

Fick demostró empíricamente que: [math] φ(x,t)≅-D \cfrac{(u(x+∆x,t)-u(x,t))}{∆x}[/math]

De aquí se deduce la anteriormente citada ley cuando suponemos que [math]∆x[/math] tiende a [math]0[/math].

Si suponemos que ∆x > 0 y la concentración en un tiempo t es mayor en x + ∆x que en x, entonces u(x + ∆x) − u(x, t) > 0, y si ∆x es pequeño, se tiene [math]ux(x,t)\gt0[/math], y [math]φ (x, t)[/math] será negativo y el flujo irá hacia la izquierda.

Si suponemos que ∆x > 0 y la concentración en un tiempo t es menor en x + ∆x que en x, entonces u(x + ∆x) − u(x, t) < 0, y si ∆x es pequeño, se tiene [math]ux(x,t)\lt0[/math], y [math]φ (x, t)[/math] será positivo y el flujo de calor irá hacia la derecha.

La constante D es el coeficiente de difusión medido en [math]\cfrac{m^2}{s}[/math], que depende de las propiedades físicas de los compuestos. Este coeficiente D puede depender de x si las condiciones físicas varían en las distintas partes de la sección del tubo. También, en principio, podría depender del tiempo. Supondremos que no, pues en caso contrario el problema en términos matemáticos se complica bastante. La supondremos constante e igual a 1.

Utilizando la ley de Fick, la concentración en la sección transversal del tubo que dista x del extremo [math]x = 0[/math] cuando ha pasado un tiempo [math]t[/math], [math]u(x, t)[/math], satisface: [math]u_t (x,t) = \cfrac{∂}{∂x} Du_x+f(x,t)[/math]: [math]x ∈(0,L),t\gt0[/math]

es decir, la ecuación en derivadas parciales: [math]u_t= \cfrac{∂}{∂x} Du_x+f(x,t)[/math]: [math]u_t-Du_{xx} = f(x,t)[/math]

Como hemos señalado anteriormente [math]f(x,t)=0[/math] y [math]D=1[/math]: [math]u_t-u_{xx}=0[/math]

Obtenemos de esta forma una ecuación de difusión similar a la del calor.

2 Planteamiento del sistema

El sistema de ecuaciones que debe satisfacer u(x,t) para que el problema esté bien propuesto es el siguiente:

[math] \left\{\begin{matrix}\\u_t-Du_{xx}=0\\u_x(0,t)=0\\u_x(L,t)=0\\u(x,0)=u_0\end{matrix}\right. [/math]

Supondremos D constante e igual a 1 y L con valor 5.

3 Resolución del sistema

Método de diferencias finitas con ∆x=0.1 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]

3.1 Método del trapecio

% Difusión de una sustancia contaminante
% u_t-u_xx=0, x en (0,L)
% u_x(0,t)=0
% u_x(L,t)=0
% u(0,t)=u0(x)
% ______________________________________
%
clear all
% Datos del problema
L=5;
T=5;
% Discretización espacial
h=0.1;
N=L/h;
% Vector de puntos en el espacio
x=0:h:L;
% Aproximación de -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=(1/h^2)*K;
% Vector F
F=zeros(N+1,1);
% Vector u0
u0=[zeros(31,1);3*ones(20,1)];
% Discretización temporal
dt=h/4;
t=0:dt:T; %Vector de tiempos
% Método del trapecio
A=eye(N+1)+dt/2*K;
U(1,:)=u0';
uu=u0;
for n=1:length(t)-1
    uu=A\((eye(N+1)-dt/2*K)*uu+dt*F);
    U(n+1,:)=uu';
end
% Dibujamos
[xx,tt]= meshgrid(x,t);
surf(xx,tt,U)

3.2 Método de Euler explícito

% Difusión de una sustancia contaminante
% u_t-u_xx=0, x en (0,L)
% u_x(0,t)=0
% u_x(L,t)=0
% u(0,t)=u0(x)
% ______________________________________
%
clear all
% Datos del problema
L=5;
T=5;
% Discretización espacial
h=0.1;
N=L/h;
% Vector de puntos en el espacio
x=0:h:L;
% Aproximación de -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=(1/h^2)*K;
% Vector F
F=zeros(N+1,1);
% Vector u0
u0=[zeros(31,1);3*ones(20,1)];
% Discretización temporal
dt=h^2/2;
t=0:dt:T; %Vector de tiempos
% Método de Euler explícito
U(1,:)=u0';
uu=u0;
for n=1:length(t)-1
    uu=uu+dt*(-K*uu+F);
    U(n+1,:)=uu';
end
% Dibujamos
[xx,tt]= meshgrid(x,t);
surf(xx,tt,U)

3.3 Método de Euler implícito

% Difusión de una sustancia contaminante
% u_t-u_xx=0, x en (0,L)
% u_x(0,t)=0
% u_x(L,t)=0
% u(0,t)=u0(x)
% ______________________________________
%
clear all
% Datos del problema
L=5;
T=5;
% Discretización espacial
h=0.1;
N=L/h;
% Vector de puntos en el espacio
x=0:h:L;
% Aproximación de -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=(1/h^2)*K;
% Vector F
F=zeros(N+1,1);
% Vector u0
u0=[zeros(31,1);3*ones(20,1)];
% Discretización temporal
dt=h^2/2;
t=0:dt:T; %Vector de tiempos
% Método de Euler implícito
A=eye(N+1)+dt*K;
U(1,:)=u0';
uu=u0;
for n=1:length(t)-1
    uu=A\(uu+dt*F);
    U(n+1,:)=uu';
end
% Dibujamos
[xx,tt]= meshgrid(x,t);
surf(xx,tt,U)

3.4 Método de Euler modificado

% Difusión de una sustancia contaminante
% u_t-u_xx=0, x en (0,L)
% u_x(0,t)=0
% u_x(L,t)=0
% u(0,t)=u0(x)
% ______________________________________
%
clear all
% Datos del problema
L=5;
T=5;
% Discretización espacial
h=0.1;
N=L/h;
% Vector de puntos en el espacio
x=0:h:L;
% Aproximación de -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=(1/h^2)*K;
% Vector F
F=zeros(N+1,1);
% Vector u0
u0=[zeros(31,1);3*ones(20,1)];
% Discretización temporal
dt=h^2/2;
t=0:dt:T; %Vector de tiempos
% Método de Euler modificado
U(1,:)=u0';
uu=u0;
for n=1:length(t)-1
    k1=-K*uu;
    k2=-K*(uu+dt*k1);
    uu=uu+dt/2*(k1+k2);
    U(n+1,:)=uu';
end
% Dibujamos
[xx,tt]= meshgrid(x,t);
surf(xx,tt,U)

3.5 Comparativa de métodos

En términos de precisión cabe destacar que es más adecuado el método de Euler implítico frente a Euler Modificado, o Runge-Kutta de orden 2. El método que presenta peores resultados es Euler explítico dado que es muy inestable para sistemas.

Se muestra a continuación gráficas explicativas entre los tres métodos:

Error entre Euler implícito y Euler modificado

Error entre Euler explícito y Euler implícito

Error entre Euler explícito y Euler modificado

4 Conservación de la masa total de contaminante

Para comprobar que la masa se conserva a lo largo del tiempo, dado que estamos ante la ecuación [math]u_t-u_{xx}=0[/math] , basta con resolver la siguiente integral:

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

que tendrá solución trivial debido a las condiciones frontera.

La primera integral,que hace referencia a la concentración [math]u(x, t)[/math] , da como resultado una constante, que ha sido resuelta numéricamente por el método del trapecio.

A=U';
integral=trapz(A);
concentracion=U(:,26);
subplot(1,2,1)
plot(t,integral)
xlabel('Tiempo')
ylabel('Masa total de contaminante')
title('Evolución de la cantidad de contaminante en el sistema')
subplot(1,2,2)
plot(t,concentracion)
xlabel('Tiempo')
ylabel('Concentración en el punto medio del tubo')
title('Evolución de la concentración en el punto medio del tubo')

Queda ilustrado en las gráficas perfectamente la evolución de la cantidad de contaminante de sistema, frente al comportamiento de la concentración en el punto medio del tubo a lo largo del tiempo.

5 Resultados para tiempos grandes

Despreciamos [math]u_t[/math]:

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

Usamos una T=100 y hacemos los cálculos con el método del trapecio:

clear all
% Datos del problema
L=5;
T=100;
% Discretización espacial
h=0.1;
N=L/h;
% Vector de puntos en el espacio
x=0:h:L;
% Aproximación de -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=(1/h^2)*K;
% Vector F
F=zeros(N+1,1);
% Vector u0
u0=[zeros(31,1);3*ones(20,1)];
% Discretización temporal
dt=h/4;
t=0:dt:T; %Vector de tiempos
% Método del trapecio
A=eye(N+1)+dt/2*K;
U(1,:)=u0';
uu=u0;
for n=1:length(t)-1
    uu=A\((eye(N+1)-dt/2*K)*uu+dt*F);
    U(n+1,:)=uu';
end
%Solución estacionaria
a=U(length(t),:);
b=U(1,:);
c=U((length(t)-1)/100+1,:);
d=U((length(t)-1)*2/100+1,:);
e=U((length(t)-1)*10/100+1,:);
hold on
plot(x,a,'b','linewidth',2);
plot(x,b,'-y','linewidth',1);
plot(x,c,'-g','linewidth',1);
plot(x,d,'-k','linewidth',1);
plot(x,e,'-r','linewidth',1);
legend('100','0','1','2','10');
xlabel('Posición')
ylabel('Concentracion')
hold off

Y en la gráfica comprobamos las distintas distancias.

6 Cálculo del tiempo necesario para alcanzar el estado estacionario con un error del 5%

clear all
% Datos del problema
L=5;
T=50;
% Discretización espacial
h=0.1;
N=L/h;
% Vector de puntos en el espacio
x=0:h:L;
% Aproximación de -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=(1/h^2)*K;
% Vector F
F=zeros(N+1,1);
% Vector u0
u0=[zeros(31,1);3*ones(20,1)];
% Discretización temporal
dt=h^2/2;
t=0:dt:T; %Vector de tiempos
% Método del trapecio
A=eye(N+1)+dt/2*K;
U(1,:)=u0';
uu=u0;
for n=1:length(t)-1;
    uu=A\((eye(N+1)-dt/2*K)*uu+dt*F);
    U(n+1,:)=uu';
end
% Estado estacionario de la concentración
est=U(length(t),:);
clear A F K L N T U dt h n t u0 uu x
% Repetimos el cálculo y lo paramos al alcanzar un error máximo del 5%
% Datos del problema
L=5;
T=50;
% Discretización espacial
h=0.1;
N=L/h;
% Vector de puntos en el espacio
x=0:h:L;
% Aproximación de -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=(1/h^2)*K;
% Vector F
F=zeros(N+1,1);
% Vector u0
u0=[zeros(31,1);3*ones(20,1)];
% Discretización temporal
dt=h^2/2;
t=0:dt:T; %Vector de tiempos
% Método del trapecio
A=eye(N+1)+dt/2*K;
U(1,:)=u0';
uu=u0;
n=1;
while max((abs(uu'-est))./est)>0.05
    uu=A\((eye(N+1)-dt/2*K)*uu+dt*F);
    U(n+1,:)=uu';
    n=n+1;
end
% Tiempo invertido
a=size(U);
b=a(1);
tiempo=(b-1)*dt
clear A F K L N T U a b dt h n t u0 uu x
% Cálculo del tiempo necesario para alcanzar el estado estacionario con un
% error del 5% para un intervalo espacial 10 veces más pequeño
% Datos del problema
L=5;
T=50;
% Discretización espacial
h=0.1/10;
N=L/h;
% Vector de puntos en el espacio
x=0:h:L;
% Aproximación de -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=(1/h^2)*K;
% Vector F
F=zeros(N+1,1);
% Vector u0
u0=[zeros(301,1);3*ones(200,1)];
% Discretización temporal
dt=(h*10)^2/2; % Mantenemos el mismo intervalo temporal
t=0:dt:T; %Vector de tiempos
% Método del trapecio
A=eye(N+1)+dt/2*K;
U(1,:)=u0';
uu=u0;
for n=1:length(t)-1;
    uu=A\((eye(N+1)-dt/2*K)*uu+dt*F);
    U(n+1,:)=uu';
end
% Estado estacionario de la concentración
est=U(length(t),:);
clear A F K L N T U a b dt h n t u0 uu x
% Comprobamos si varía el tiempo al hacer h 10 veces más pequeño
% Datos del problema
L=5;
T=50;
% Discretización espacial
h=0.1/10;
N=L/h;
% Vector de puntos en el espacio
x=0:h:L;
% Aproximación de -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=(1/h^2)*K;
% Vector F
F=zeros(N+1,1);
% Vector u0
u0=[zeros(301,1);3*ones(200,1)];
% Discretización temporal
dt=(h*10)^2/2; % Mantenemos el mismo intervalo temporal
t=0:dt:T; %Vector de tiempos
% Método del trapecio
A=eye(N+1)+dt/2*K;
U(1,:)=u0';
uu=u0;
n=1;
while max((abs(uu'-est))./est)>0.05
    uu=A\((eye(N+1)-dt/2*K)*uu+dt*F);
    U(n+1,:)=uu';
    n=n+1;
end
% Tiempo invertido
a=size(U);
b=a(1);
nuevo_tiempo=(b-1)*dt

El tiempo que tarda en alcanzar la concentración el estado estacionario con un error del 5% es 8.68 segundos con h=0.1. Si cambiamos la longitud de paso a 0.01, el nuevo tiempo es practicamente invariable, 8.645 segundos.

7 Caso de colocar un limpiador en un extremo

Colocar un limpiador en el extremo izquierdo se traduce físicamente en que la concentración en x=0 es ahora nula. Por lo tanto el problema a resolver ahora es el siguiente:

[math] \left\{\begin{matrix}\\u_t-Du_{xx}=0\\u(0,t)=0 & t≥0\\u_x(L,t)=0 & t≥0\\u(x,0)=u_0 & xє(0,L)\end{matrix}\right. [/math]

7.1 Método de diferencias finitas

[math] \left\{\begin{matrix}\\u'_n(t) + \frac {-u_n + 2u_n(t) + u_{n+1}(t)}{h^2}\\\frac {u_{n+1}(t)-u_{n+1}(t)}{2h}\\u_0(t)=0\\u_n(0)=u_0\end{matrix}\right. [/math]

7.2 Resolución numérica

Al cambiar las condiciones de contorno, el valor estacionario de la concentración [math] u(x,t)[/math] es 0 es todos los puntos del tubo, lo que se traduce en que el contaminante desaparece por completo cuando el tiempo es muy grande.

Será preciso calcular el error relativo en cada punto del tubo, obtener el valor absoluto y escoger el máximo valor.

[math]error \ relativo=\cfrac{valor \ real-valor \ estacionario}{valor \ estacionario}=5 [/math]%

El tiempo que se tarda en alcanzar el valor del 5% es dependiente del valor que consideramos como estacionario, que en nuestro caso será la concentración correspondiente a T=100 segundos desde que se inicia el ensayo (elegimos este tiempo ya que es suficientemente pequeño el valor de la concentración [math] u(x,t)[/math]. El resultado que se obtiene de tiempo para alcanzar dicho error del 5% es 99.49 segundos, quedando reflejado a continuación:

% Difusión de una sustancia contaminante con limpiador en el extremo de la
% izquierda
% u_t-u_xx=0, x en (0,L)
% u(0,t)=0
% u_x(L,t)=0
% u(0,t)=u0(x)
% ______________________________________
%
clear all
% Datos del problema
L=5;
T=100;
% Discretización espacial
h=0.1;
N=L/h;
% Vector de puntos en el espacio
x=0:h:L;
% Aproximación de -u_xx
K=2*diag(ones(1,N+1))-diag(ones(1,N),1)-diag(ones(1,N),-1);
K(N+1,N)=-2;
K=(1/h^2)*K;
% Vector F
F=zeros(N+1,1);
% Vector u0
u0=[zeros(31,1);3*ones(20,1)];
% Discretización temporal
dt=h^2/2;
t=0:dt:T; %Vector de tiempos
% Método del trapecio
A=eye(N+1)+dt/2*K;
U(1,:)=u0';
uu=u0;
for n=1:length(t)-1
    uu=A\((eye(N+1)-dt/2*K)*uu+dt*F);
    U(n+1,:)=uu';
end
% Dibujamos
[xx,tt]= meshgrid(x,t);
surf(xx,tt,U)
% Valor estacionario de la concentración en el tubo
est=U(length(t),:);
clear A F K L N T U dt h n t u0 uu x xx tt
% Repetimos el cálculo y lo paramos al alcanzar un error máximo del 5%
% Datos del problema
L=5;
T=100;
% Discretización espacial
h=0.1;
N=L/h;
% Vector de puntos en el espacio
x=0:h:L;
% Aproximación de -u_xx
K=2*diag(ones(1,N+1))-diag(ones(1,N),1)-diag(ones(1,N),-1);
K(N+1,N)=-2;
K=(1/h^2)*K;
% Vector F
F=zeros(N+1,1);
% Vector u0
u0=[zeros(31,1);3*ones(20,1)];
% Discretización temporal
dt=h^2/2;
t=0:dt:T; %Vector de tiempos
% Método del trapecio
A=eye(N+1)+dt/2*K;
U(1,:)=u0';
uu=u0;
n=1;
while max((abs(uu'-est))./est)>0.05
    uu=A\((eye(N+1)-dt/2*K)*uu+dt*F);
    U(n+1,:)=uu';
    n=n+1;
end
% Tiempo invertido
a=size(U);
b=a(1);
tiempo=(b-1)*dt

8 Método de Fourier

Para la resolución del sistema mediante el método de Fourier:

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

Buscaremos soluciones del tipo [math]u(x,t)=X(x)T(t)[/math]:

[math]XT^{'}-X^{''} T=0[/math]: [math]\cfrac{T'}{T}=\cfrac{X^{''}}{X}=-μ[/math]

El problema de autovalores asociado, será:

[math] \left\{\begin{matrix}\\X^{''}-μX=0 \\X'(0)=0 \\X^{'} (L)=0\end{matrix}\right. [/math]:

[math]μ_k=(\cfrac{kπ}{L})^2[/math]:

[math]φ_k=cos⁡(\cfrac{kπ}{L} x)[/math]

Por otro lado, la ecuación diferencial que debe cumplir [math]T(t)[/math] es:

[math] T^{'}_k (t)-μ_k T_k (t)=0[/math]

Resolviéndola obtenemos:

[math] T_k (t)=C_k e^{-μ_k t}[/math]

Entonces, la solución aproximada mediante el método de Fourier es la siguiente:

[math]u(x,t)= \sum_{k=0}^K C_k e^{-μt} φ_k [/math]

A continuación se muestra el código de matlab que resuelve el problema mediante Fourier usando 1, 3, 5, 10 y 20 términos de la serie.

% Difusión de una sustancia contaminante
% u_t-u_xx=0, x en (0,L)
% u_x(0,t)=0
% u_x(L,t)=0
% u(0,t)=u0(x)
% Método de Fourier
% ______________________________________
%
clear all
% Datos del problema
L=5; % Longitud del tubo
a=0; b=5; % Extremos
T=5; % Tiempo final
% Numero de autofunciones de la serie = M
M=20;
% Mallado
% En x --------> Paso espacial = h
h=(b-a)/50;
x=a:h:b;
% En t --------> Paso temporal = j
j=h/4;
t=0:j:T;
[Mx,Mt]=meshgrid(x,t);
f=[zeros(1,31),3*ones(1,20)]; % Valor inicial
U=0;
for k=0:M
    p=cos(k*pi*x/L); % Autofunciones
    Ck=trapz(x,f.*p)/trapz(x,p.^2);
    U=U+Ck*(exp(-(k^2*pi^2/L^2)*Mt).*cos(k*pi*Mx/L));
end
surf(Mx,Mt,U)
xlabel('Espacio');
ylabel('Tiempo');

Para obtener los diferentes resultados basta con introducir en M el número de términos de la serie de Fourier que vayamos a usar. Para M=20 se obtiene un resultado similar a los obtenidos anteriormente con otros métodos

y si lo calculamos para M=50 se aprecia que la gráfica se suaviza y se asemeja aun más.

Las aproximaciones en t=0.5 para 1, 3, 5, 10 y 20 términos de Fourier serán

Se aprecia que para M=1 (en verde) se da la peor aproximación.

9 Cambio en las condiciones de frontera

El problema pasa a ser

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

Este sistema tiene una condición de contorno tipo Dirichlet en el extremo de la varilla [math]x=L[/math]. Esto quiere decir que este extremo está en contacto con una fuente de contaminante siendo la concentración en el instante t en dicho extremo [math]u=10sent[/math]

La otra condición es de tipo Neumann en cero [math]-Du_x(0,t)=φ(t)[/math]

En este caso [math]φ(t)=-3\lt0[/math], además [math]D=1\gt0[/math], entonces debe verificarse [math]u_x(0,t)\gt0[/math]. Cogiendo [math]∆x\gt0[/math] y pequeño obtenemos

[math]u(0+∆x,t)-u(0,t)\gt0[/math]

[math]u(0+∆x,t)\gtu(0,t)[/math]

lo que indica que el flujo de contaminante irá hacia la izquierda,interpretado de forma que por el extremo x=0 está saliendo una cantidad de contaminante en el instante t igual a φ=3.

% Difusión de una sustancia contaminante
% u_t-u_xx=0, x en (0,L)
% u_x(0,t)=3
% u(L,t)=f(t)
% u(0,t)=u0(x)
% ______________________________________
%
clear all
% Datos del problema
L=5;
T=10;
% Discretización espacial
h=0.1;
N=L/h;
% Vector de puntos en el espacio
x=0:h:L;
xint=0:h:L-h; % Nodos interiores
% Aproximación 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=(1/h^2)*K;
% Vector F
F=zeros(N,1);
F(1)=-3*2/h;
% Vector u0
u0=[zeros(31,1);3*ones(19,1)];
% Discretización temporal
dt=h^2/2;
t=0:dt:T; %Vector de tiempos
% Función f(t)
f=10*sin(t);
% Método de Euler explícito
U(1,:)=[u0' f(1)];
uu=u0;
for n=1:length(t)-1
    F(N)=f(n)/h^2;
    uu=uu+dt*(-K*uu+F);
    U(n+1,:)=[uu' f(n+1)];
end
% Dibujamos
[xx,tt]= meshgrid(x,t);
surf(xx,tt,U)