Ecuación del calor (Grupo 8)

De MateWiki
Saltar a: navegación, buscar
Trabajo realizado por estudiantes
Título Ecuación del Calor. Grupo 8-B
Asignatura Ecuaciones Diferenciales
Curso Curso 2013-14
Autores Álvaro Conde Gallastegui

Ana Álvarez Cruz

Elena de Frutos Aguilar

Ana Uceda Atán

Diego Pontiveros Bermejo

Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

La ecuación del calor es un modelo matemático que describe la evolución de la temperatura en un cuerpo sólido en función del tiempo y del espacio. Esta ecuación fue propuesta por el físico francés Jean-Baptiste Joseph Fourier en 1807, pero no sería hasta 1822 cuando la academia de ciencias francesas, a la que pertenecía, decidió publicarla bajo el título de “Théorie analytique de la chaleur” (Teoría analítica del calor). Su expresión matemática, simplificada para nuestro caso, en el que tratamos la distribución de temperaturas en una varilla, que es un elemento unidimensional (con x como variable longitud), sería la siguiente: 


[math]\frac{\partial U}{\partial t}-\frac{\partial ^2U}{\partial x^2}=0[/math]


Siendo U (x,t) la temperatura de cada punto de la varilla en cada instante t. Consideramos una varilla de longitud L de un cierto material, de espesor constante. Está orientada en la dirección del eje x, desde x=0 a x=L. La varilla es conductora de calor, por lo que entre dos zonas de ella a diferente temperatura hay un intercambio de energía térmica. En nuestro caso, consideramos una varilla delgada de longitud L=3m, y cuyos extremos (x=0 y x=3) están colocados sobre objetos que mantienen una temperatura constante de 0 y 10º C respectivamente. Suponemos que la varilla es de grosor despreciable y tiene su superficie exterior aislada térmicamente. Podemos entonces pensar que todas las temperaturas son constantes a lo largo de cada sección transversal, y ver la varilla como un objeto unidimensional. Inicialmente la temperatura de la varilla viene dada por U0(x)=10X/3 salvo en su tercio central donde la temperatura ha subido hasta los 100 grados. Designamos por U(x, t) a la temperatura de la sección de la varilla que dista x ≥ 0 del extremo x = 0 cuando ha pasado un tiempo t ≥ 0. La función U (x,t) para t=0 se convierte en una función que llamamos U0(x), y que definimos así:


\begin{cases} u_t-u_{xx}=0, \qquad x\epsilon(0,3), t>0\\ u(0,t)=0, u(3,t)=10, \qquad t>0\\ u(x,0)=u_0(x), \qquad x\epsilon[0,3] \end{cases}

[math] u_0(x)= \begin{cases} 10x/3 & \mbox{si} & x \in \mbox{(0,1)}\cup\mbox{(2,3)} \\ 100 & \mbox{si} & x \in \mbox{(1,2)} \end{cases} [/math]


A continuación aproximaremos la resolución de la ecuación del calor por diferentes métodos y estudiaremos la evolución de las temperaturas en los diferentes puntos de la varilla a lo largo de diferentes intervalos de tiempo. Definamos el Problema propiamente con su ecuación y sus condiciones de frontera:


[math]u_t-u_xx=0 \quad x\in(0,3) \quad t\gt0 [/math]


[math]u(x,0)=\frac{10x}{3} \quad x\in(0,1) U (2,3)[/math]


[math]u(x,0)=100\quad x\in(1,2)[/math]


[math]u(0,t)=0[/math]


[math]u(3,t)=10[/math]


1 Solución mediante el método del Trapecio

%datos del problema
L=3;
h=0.1;
T=2;
%discretizacion espacial
N=L/h;
%vector de puntos en el espacio
x=0:h:L;
xint=h:h:L-h;%nodos interiores
%aproximacion de -Uxx
K=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
K=(1/h^2)*K;
%calculamos F
F=zeros(N-1,1);
F(N-1,1)=10/h^2;
%calculamos U0
U0=(10*xint/3)';
for j=1:length(xint)
    if (xint(j)>1)&(xint(j)<2)
        U0(j)=100;
    end
end

%discretizacion temporal
dt=h; %paso en el espacio
t=0:dt:T;%vector en tiempos



%metodo del trapecio
I=eye(N-1);
M=I+(dt./2)*K;
U(1,:)=[0   U0' 10];
uu=U0;
for n=1:length(t)-1
    uu=M\(uu+dt/2*(F)+dt/2*(-K*uu+F));
    U(n+1,:)=[0 uu' 10];
end
%dibujamos la solucion (apartado 3)



[xx,tt]=meshgrid(x,t);
 
figure(1)
surf(xx,tt,U)
xlabel('espacio')
ylabel('tiempo')
zlabel('Temperatura')
figure(2)
plot(t,U(:,15))
Variación de la temperatura en funcion del tiempo y la posición de la barra
Variación de la temperatura respecto al tiempo en el punto medio de la varilla

En la primera gráfica se observa la diferencia de temperatura en la barra debido a las condiciones iniciales y cómo va disminuyendo su temperatura con el paso del tiempo. Como se aprecia en la gráfica de la derecha, el punto medio de la barra sufre una disminución en la temperatura a lo largo del tiempo. Hay un descenso de la temperatura mas acentuada en los primeros instantes, que se debe al contraste existente entre la barra y el entorno, el cual actúa como sumidero de calor. En ambas representaciones se ven unos picos o asientos térmicos que pueden ser debidos a la inexactitud del método empleado. Con un paso en el tiempo menor y un mayor intervalo, estos errores son menos apreciables.

2 Solución mediante el método de Euler implicito

%datos del problema
L=3;
h=0.1;
T=2;
%discretizacion espacial
N=L/h;
%vector de puntos en el espacio
x=0:h:L;
xint=h:h:L-h;%nodos interiores
%aproximacion de -Uxx
K=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
K=(1/h^2)*K;
%calculamos F
F=zeros(N-1,1);
F(N-1,1)=10/h^2;
%calculamos U0
U0=(10*xint/3)';
for j=1:length(xint)
    if (xint(j)>1)&(xint(j)<2)
        U0(j)=100;
    end
end

%discretizacion temporal
dt=h; %paso en el espacio
t=0:dt:T;%vector en tiempos



%metodo de euler implicito
I=eye(N-1);
M=I+dt.*K;
U(1,:)=[0   U0' 10];
uu=U0;
for n=1:length(t)-1

    uu=M\(uu+dt*(F));
    U(n+1,:)=[0 uu' 10];
end




%dibujamos la solucion (apartado 3)

[xx,tt]=meshgrid(x,t);
 
figure(1)
surf(xx,tt,U)
xlabel('espacio')
ylabel('tiempo')
zlabel('Temperatura')
figure(2)
plot(t,U(:,15))
Variación de la temperatura en funcion del tiempo y la posición de la barra
Variación de la temperatura respecto al tiempo en el punto medio de la varilla

3 Solución mediante el método de Euler explícito

%datos del problema
L=3;
h=0.1;
T=2;
%discretizacion espacial
N=L/h;
%vector de puntos en el espacio
x=0:h:L;
xint=h:h:L-h;%nodos interiores
%aproximacion de -Uxx
K=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
K=(1/h^2)*K;
%calculamos F
F=zeros(N-1,1);
F(N-1,1)=10/h^2;
%calculamos U0
U0=(10*xint/3)';
for j=1:length(xint)
    if (xint(j)>1)&(xint(j)<2)
        U0(j)=100;
    end
end

%discretizacion temporal
dt=h^2/2; %paso en el espacio 
t=0:dt:T;%vector en tiempos



%metodo de euler explicito
I=eye(N-1);
M=I-dt.*K;
U(1,:)=[0   U0' 10];
uu=U0;
for n=1:length(t)-1

    uu=M*uu+dt*F;
    U(n+1,:)=[0 uu' 10];
end

%dibujamos la solucion (apartado 3)
[xx,tt]=meshgrid(x,t);
 
figure(1)
surf(xx,tt,U)
xlabel('espacio')
ylabel('tiempo')
zlabel('Temperatura')
figure(2)
plot(t,U(:,15))
Si dt=h, el paso en el tiempo es demasiado grande y el programa no responde. Por ello, es necesario reducir su valor para que se ejecute el programa correctamente. En nuestro caso hemos tomado dt=h^2/2.
Variación de la temperatura en funcion del tiempo y la posición de la barra
Variación de la temperatura respecto al tiempo en el punto medio de la varilla

4 Solución mediante el método de Runge-Kutta

%datos del problema
L=3;
h=0.1;
T=2;
%discretizacion espacial
N=L/h;
%vector de puntos en el espacio
x=0:h:L;
xint=h:h:L-h;%nodos interiores
%aproximacion de -Uxx
K=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
K=(1/h^2)*K;
%calculamos F
F=zeros(N-1,1);
F(N-1,1)=10/h^2;
%calculamos U0
U0=(10*xint/3)';
for j=1:length(xint)
    if (xint(j)>1)&(xint(j)<2)
        U0(j)=100;
    end
end

%discretizacion temporal
dt=h^2/2; %paso en el espacio
t=0:dt:T;%vector en tiempos



%metodo de runge-kuta
I=eye(N-1);
M=I-dt.*K;
U(1,:)=[0   U0' 10];
uu=U0;
for n=1:length(t)-1
 k1=-K*uu+F;
 k2=-K*(uu+k1*dt/2)+F;
 k3=-K*(uu+k2*dt/2)+F;
 k4=-K*(uu+k3*dt)+F;
 uu=uu+(dt/6)*(k1+2*k2+2*k3+k4);
 U(n+1,:)=[0 uu' 10];

end




%dibujamos la solucion (apartado 3)
[xx,tt]=meshgrid(x,t);
 
figure(1)
surf(xx,tt,U)
xlabel('espacio')
ylabel('tiempo')
zlabel('Temperatura')
figure(2)
plot(t,U(:,15))
Variación de la temperatura en función del tiempo y la posición de la barra
Variación de la temperatura respecto al tiempo en el punto medio de la varilla
Debido a la elección de un paso extremadamente pequeño, las líneas de nivel de la superficie se encuentran tan próximas que hacen imposible la visualización de la superficie. Por ello, empleamos el comando "shading flat" para eliminar de la gráfica estas lineas y que se distinga mejor.
Variación de la temperatura en función del tiempo y la posición de la barra

Una vez resuelta la ecuación por los cuatro métodos anteriores podemos deducir que el método de Euler implícito es el que mejor aproxima,ya que tanto Runge-kutta como euler explícito necesitan un paso menor en la discretización para mejorar la aproximación.


5 Estado estacionario para tiempos grandes

%datos del problema
L=3;
h=0.1;
T=18;
%discretizacion espacial
N=L/h;
%vector de puntos en el espacio
x=0:h:L;
xint=h:h:L-h;%nodos interiores
%aproximacion de -Uxx
K=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
K=(1/h^2)*K;
%calculamos F
F=zeros(N-1,1);
F(N-1,1)=10/h^2;
%calculamos U0
U0=(10*xint/3)';
for j=1:length(xint)
    if (xint(j)>1)&(xint(j)<2)
        U0(j)=100;
    end
end

%discretizacion temporal
dt=h^2/2; %paso en el espacio
t=0:dt:T;%vector en tiempos



%metodo de runge-kuta
I=eye(N-1);
M=I-dt.*K;
U(1,:)=[0   U0' 10];
uu=U0;
for n=1:length(t)-1
 k1=-K*uu+F;
 k2=-K*(uu+k1*dt/2)+F;
 k3=-K*(uu+k2*dt/2)+F;
 k4=-K*(uu+k3*dt)+F;
 uu=uu+(dt/6)*(k1+2*k2+2*k3+k4);
 U(n+1,:)=[0 uu' 10];

end
 ut1=U(2001,:);
 ut2=U(201,:);
 ut3=U(401,:);
 ut4=U(1,:);
 ut=10.*x/3;
 ur1=ut1-ut;
ur2=ut2-ut;
ur3=ut3-ut;
ur4=ut4-ut;

for i=1:length(t)
    urt=U(i,:);
    if max(urt-ut)<=0.05
    inst=(i-1)/200;
    break
    end
end
inst
%nos dará el instante exacto en el que la temperatura se regulariza hasta ser independiente del tiempo, concretamente 6,395 segundos
%dibujamos la solucion (apartado 3)
[xx,tt]=meshgrid(x,t);

figure(1)
surf(xx,tt,U)
shading flat
xlabel('espacio')
ylabel('tiempo')
zlabel('Temperatura')
figure(2)
plot(x,ut1)
hold on
plot(x,ut,'k')
figure(3)
plot(x,ut2)
hold on
plot(x,ut,'k')
figure(4)
plot(x,ut3)
hold on
plot(x,ut,'k')
figure(5)
plot(x,ut4)
hold on
plot(x,ut,'k')
figure(6)
plot(x,ur1)
figure(7)
plot(x,ur2)
figure(8)
plot(x,ur3)
figure(9)
plot(x,ur4)
Variación de la temperatura en funcion del tiempo y la posición de la barra
Distribución de temperaturas en t=0 comparada con la solución estacionaria
Diferencia entre la distribución de temperaturas en t=0 y la solución estacionaria
Distribución de temperaturas en t=1 comparada con la solución estacionaria
Diferencia entre la distribución de temperaturas en t=1 y la solución estacionaria
Distribución de temperaturas en t=2 comparada con la solución estacionaria
Diferencia entre la distribución de temperaturas en t=2 y la solución estacionaria
Distribución de temperaturas en t=10 comparada con la solución estacionaria
Diferencia entre la distribución de temperaturas en t=10 y la solución estacionaria
Distribución de temperaturas en t=6,395 comparada con la solución estacionaria
Diferencia entre la distribución de temperaturas en t=6,395 y la solución estacionaria

Cómo comparar las distribuciones de temperatura con la solución estacionaria.


Definimos un vector (en este programa lo hemos llamado ut) que contenga una fila de la matriz U, es decir, la distribución de temperaturas en todos los puntos de la varilla para un instante determinado. Con T=10 y dt=0,005, sabemos que la matriz tiene 2001 filas. Tomamos los instantes 0, 1, 2 y 10, que se corresponden con las filas 1, 201, 401 y 2001 y comparamos la distribución de temperaturas en ese instante con la solución estacionaria que hemos calculado. Sabemos que la ecuación del calor es regularizante, por lo que la distribución de temperaturas para tiempos muy grandes tiende a ser independiente del instante, dependiendo únicamente de la posición. En este caso, dadas las condiciones iniciales, un extremo se encuentra a 0ºC y el opuesto a 10ºC, y dado que la varilla es homogénea y mide 3m, es fácilmente deducible que si no existen más intercambios de calor, la distribución final en función de la distancia al extremo será: U(x)=10X/3. Podemos ver en las gráficas que la diferencia entre esta ecuación y la distribución para t=10 es prácticamente inapreciable, por lo que, debido a que se trata de una aproximación, podemos asegurar que la distribución estacionaria de temperaturas se alcanza antes del instante t=10. Pero, ¿en que instante se alcanza exactamente? Mediante un sencillo bucle (líneas 54 a 60) comparamos todas las filas de la matriz U con la solución estacionaria hasta que la máxima diferencia entre ellas sea menor al 0,05%. Una vez se alcance esta diferencia, el bucle se detiene y nos da el numero de fila que corresponde. Obtenemos así la fila 1280, que se corresponde con el instante t=6,395, instante en el que podemos afirmar que la distribución de temperaturas empieza a ser independiente del tiempo con un error del 0,05%.

6 Varilla con extremo derecho aislado térmicamente

[math] \begin{cases} u_t-u_{xx}=0, \qquad x\epsilon(0,3), t\gt0\\ u(0,t)=0, u_x(3,t)=0, \qquad t\gt0\\ u(x,0)=u_0(x), \qquad x\epsilon[0,3] \end{cases} [/math]

En este apartado tenemos una condición tipo Newmann (u_x(3,t)=0) lo que indica que el extremo derecho está térmicamente aislado por lo que el flujo de calor es nulo. La diferencia de este programa radica en la variación en la condiciones del nuevo sistema, en este caso la condición tipo Newmann no nos da información sobre la temperatura de la varilla en su extremo derecho (L=3) provocando que la dimensión de la matriz solución aumente pasando de N-1 a N. 

%Datos del problema
L=3;
T=10;
 
%Datos de la discretización
h=0.1;
N=L/h;
x=0:h:L;
xint=h:h:L; %excluyo un extremo porque tengo la condicion(??)
dt=h/2;
t=0:dt/10:T;
%Discretización espacial
%aprox 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=(10*xint/3)';
for j=1:length(xint)
    if (xint(j)>1)&(xint(j)<2)
        U0(j)=100;
    end
end
 
%Definmos F
F=zeros(N,1);
%tenemos U'+KU=F ahora tenemos que aproximar temporalmete
%Discretización temporal
sol(1,:)=[0,U0'];%guardamos la solucion
U=U0; %inicialización
%calculamos U_n ---> U_n+1
 
%Euler implícito 
for n=1:length(t)-1
    U=(eye(N)+dt*K)\(U+dt*F);
    sol(n+1,:)=[0,U'];%guardamos la solucion
end
for i=1:length(t)
    urt=sol(i,:);
    a=max(urt);
    if a<=0.05
    b=i;
    inst=(b-1)/200;
    break
    end
end
inst
b
%dibujamos
[xx,tt]=meshgrid(x,t);
figure(1)
surf(xx,tt,sol)
shading flat
y=zeros(1,length(x));
figure(2)
plot(x,urt)
hold on 
plot(x,y,'r')
Variación de la temperatura en funcion del tiempo y la posición de la barra
Diferencia en la distribución de temperaturas entre la solución estacionaria y la varilla en t=2,525

De la misma forma que en el apartado anterior, obtenemos el instante en el cual la distribución de temperaturas se vuelve estacionaria con un bucle que compare cada fila de la matriz solución de la ecuación con la solución de la ecuación que resulta al eliminar el término ut.

En este caso sabemos que dicha solución estacionaria es además independiente de la posición del punto en la varilla. Concretamente, dadas las condiciones iniciales, es una función constante de valor 0 en todos los puntos, lo que simplifica considerablemente nuestros cálculos.

En la segunda gráfica se representan los valores que toma la temperatura de todos los puntos de la varilla en el instante obtenido como comienzo del estado estacionario (t=2,525). Se puede apreciar que en ningún momento superan el valor 0,05, lo que nos permite afirmar que a partir de ese momento comienza el estado estacionario con un error del 5%.

7 Aproximación por Serie de Fourier

Para resolver el problema por el método de Fourier hallamos los autovalores y autofunciones correspondientes:

[math] \begin{cases} \varphi’’(x)+\lambda\varphi(x)=0\\ \varphi(0)=0, \varphi’(3)=0 \end{cases} [/math]


Los autovalores que cumplen la ecuación anterior son [math]\mu_k[/math] y sus respectivas autofunciones [math]\varphi_k(x)[/math]

[math]{\mu_k=(k-{1\over2})^2{\pi^2\over9}}[/math];
[math]{\varphi_k(x)=\sin(k-{1\over2}){\pi\over3}x}, \qquad k=1,2,3…N[/math]

Debido a la unicidad de los coeficientes de Fourier obtenemos:

[math]{a_k={\int_{0}^{3}u(x,0)\varphi_k(x)dx\over\int_{0}^{3}\varphi_k^2(x)dx}}[/math]

La resolución numérica consiste en dar como solución aproximada una función del tipo:


[math]{u(x,t)= \sum_{k=1}^Q a_ke^{-(k-{1\over2})^2{\pi^2\over9}t}\sin(k-{1\over2}){\pi\over3}x}[/math]

Q representa el número de elementos de la serie de Fourier, cuya aproximación es mayor cuanto mayor sea, esto puede apreciarse en las siguientes gráficas con Q= 1,3,5,10 y 20. Hemos añadido una gráfica para Q=100 en la que se observa con mayor claridad. 


%u_t-u_xx= 0
%u(0,t)=0
%u_x(3,t)=0
%u(x,0)=u0(x)

%Datos del problema
L=3;
T=10;
%Datos de la aproximación
Q=1; %numero de coef de Fourier
h=0.01; %numero de int en x
N=L/h;
x=0:h:L;
dt=h; %mismo paso en timpo y espacio
t=0:dt:T;
f=((10*x)/3).*(x<=1)+100*(1<x).*(x<2)+((10*x)/3).*(x>=2);
%calculamos la aproximacion
sol=zeros(length(t),N+1);


for k=1:Q
    phi=sin(((k-1/2)*pi*x)/3); %autofuncio hallada antes a mano
    c=0; %coeficiente
    a=(trapz(x,f.*phi))/trapz(x,phi.^2);
    T=a*exp(-((k-1/2)^2*pi^2/9)*t);
    sol=sol+T'*phi;
end

[xx,tt]=meshgrid(x,t);
surf(xx,tt,sol)
shading flat
figure(2)
plot(x, sol(1,:),'r')
hold on 
plot(x,f)


Q=1

Aproximación de la distribución de la temperatura inicial con el coeficiente Q=1
Aproximación de la temperatura en funcion del tiempo y la posición de la barra para t=10
Aproximación de la temperatura respecto al tiempo en el punto medio de la varilla para t=0.5

Q=3

Aproximación de la distribución de la temperatura inicial con el coeficiente Q=3
Aproximación de la temperatura respecto al tiempo en el punto medio de la varilla para t=10
Aproximación de la temperatura respecto al tiempo en el punto medio de la varilla para t=0.5

Q=5

Aproximación de la distribución de la temperatura inicial con el coeficiente Q=5
Aproximación de la temperatura respecto al tiempo en el punto medio de la varilla para t=10
Aproximación de la temperatura respecto al tiempo en el punto medio de la varilla para t=0.5

Q=10

Aproximación de la distribución de la temperatura inicial con el coeficiente Q=10
Aproximación de la temperatura respecto al tiempo en el punto medio de la varilla para t=10
Aproximación de la temperatura respecto al tiempo en el punto medio de la varilla para t=0.5

Q=20

Aproximación de la distribución de la temperatura inicial con el coeficiente Q=20
Aproximación de la temperatura respecto al tiempo en el punto medio de la varilla para t=10
Aproximación de la temperatura respecto al tiempo en el punto medio de la varilla para t=0.5

Q=100

Aproximación de la distribución de la temperatura inicial con el coeficiente Q=100
Aproximación de la temperatura respecto al tiempo en el punto medio de la varilla para t=10
Aproximación de la temperatura respecto al tiempo en el punto medio de la varilla para t=0.5

8 Pérdida de calor de la varilla a través del aire a Tª=16

[math] \begin{cases} u_t-u_{xx}+u-16=0,\ x \in \mbox{[0,3]} & t\gt0 \\ u(0,t)=0; u(3,t)=10; \end{cases} [/math] 

%datos del problema
 
L=3;
 
h=0.1;
 
T=10;
 
%discretizacion espacial
 
N=L/h;
 
%vector de puntos en el espacio
 
x=0:h:L;
 
xint=h:h:L-h;%nodos interiores
 
%aproximacion de -Uxx
 
K=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
 
K=(1/h^2)*K;
K=eye(N-1)+K;
 
%calculamos F
 
F=16*ones(N-1,1);
 
F(N-1,1)=16+10/h^2;
 
%calculamos U0
 
U0=(10*xint/3)';
 
for j=1:length(xint)
 
    if (xint(j)>1)&(xint(j)<2)
 
        U0(j)=100;
 
    end
 
end
 
 
 
%discretizacion temporal
 
dt=h^2/2; %paso en el espacio
 
t=0:dt:T;%vector en tiempos
 
%metodo de runge-kuta
 
I=eye(N-1);
 
M=I-dt.*K;
 
U(1,:)=[0   U0' 10];
 
uu=U0;
 
for n=1:length(t)-1
 
 k1=-K*uu+F;
 
 k2=-K*(uu+k1*dt/2)+F;
 
 k3=-K*(uu+k2*dt/2)+F;
 
 k4=-K*(uu+k3*dt)+F;
 
 uu=uu+(dt/6)*(k1+2*k2+2*k3+k4);
 
 U(n+1,:)=[0 uu' 10];
 
end
urt=U(2001,:);
 
[xx,tt]=meshgrid(x,t);

yy=zeros(length(t),length(x));
y=(2*(-3*exp(3-x)+8*exp(6-x)-8*exp(x)+3*exp(x+3)+8-8*exp(6))/(1-exp(6)));

for j=1:length(t)
    yy(j,:)=y;
end
for i=1:length(t)
    urt=U(i,:);
    a=max(abs(urt-y));
    if a<=0.001
    b=i;
    inst=((b-1)/200);
    break
    end
end
inst
b

figure(1)
surf(xx,tt,U)
shading flat
 
xlabel('espacio')
 
ylabel('tiempo')
 
zlabel('Temperatura')
hold on
surf(xx,tt,yy)
Apartado8fig2.jpg
Apartado8fig1.jpg

En este caso nos encontramos con unas condiciones iniciales diferentes, que consisten en una temperatura constante del aire de 16ºC. Esto quiere decir que tras un intervalo de tiempo suficientemente largo como para que se produzca un gran intercambio de calor, la temperatura de cada punto de la varilla no dependerá del instante en el que se mida. Sabemos de inicio que el extremo x=0 se encuentra a 0ºC y el x=3 se encuentra a 10ºC. Además, estando el entorno a 16ºC, la temperatura en cada punto de la varilla a partir del instante en que se vuelve estacionaria seguirá una función que hemos calculado resolviendo la ecuación diferencial resultante de igualar el término u_t a 0. Así, obtenemos la ecuación siguiente:

[math]u(x)=2\frac{\ -3e^(3-x)+8e^(6-x)-8e^x+3e^(x+3)+8-8e^6}{\ 1-e^6}[/math]

Mediante un proceso análogo al del apartado 5, comparamos la distribución de temperaturas en cada instante con la solución estacionaria obtenida hasta que sea menor a 10-3. El resultado obtenido nos dice que a partir de un instante t=4,54s, se puede afirmar que la distribución de temperaturas es independiente del tiempo con un error inferior al 0,1%. En la gráfica hemos representado la distribución de temperaturas junto con esa distribución estacionaria (en negro). Gracias a una vista cenital se puede apreciar perfectamente el instante en que la distribución se vuelve estacionaria.


8.1 Diferentes condiones de contorno

Nuestro problema inicial tenía unas condiciones de contorno con valores constantes, en el extremo izquierdo (x=0)se mantiene a 0 y en el derecho (x=3) a 10. En cambio, en este último apartado el extremo izquierdo está en contacto con un depósito o material cuya temperatura en el instante t es 10sen(t). En el extremo derecho, cuya condición es tipo Neumann, está entrando una cantidad de calor por unidad de tiempo igual a 1. 

clear all
%datos del problema

L=3;
h=0.1;
T=10;
%discretizacion espacial
N=L/h;
%vector de puntos en el espacio
x=0:h:L; 
xint=h:h:L;%nodos interiores
 
%discretizacion temporal
 
dt=h^2/2; %paso en el espacio
 
t=0:dt:T;%vector en tiempos
%aproximacion de -Uxx
K=2*diag(ones(1,N))-diag(ones(1,N-1),1)-diag(ones(1,N-1),-1);
K(N,N-1)=-2;
K=(1/h^2)*K;
K=eye(N)+K;
%calculamos F
g=10*sin(t);

F=16*ones(N,1);
 
F(N,1)=16+2/h;

 
%calculamos U0
 
U0=(10*xint/3);
 
for j=1:length(xint)
 
    if (xint(j)>1)&(xint(j)<2)
 
        U0(j)=100;
 
    end
 
end
 
 


U(1,:)=[g(1) U0];

%metodo de runge-kuta
 
I=eye(N);
 
M=I-dt.*K;
 

 
uu=U0';
 
for n=1:length(t)-1
 F(1,1)=16+g(n)./(h^2);
 
 k1=-K*uu+F;
 
 k2=-K*(uu+k1*dt/2)+F;
 
 k3=-K*(uu+k2*dt/2)+F;
 
 k4=-K*(uu+k3*dt)+F;
 
 uu=uu+(dt/6)*(k1+2*k2+2*k3+k4);
 
 U(n+1,:)=[g(n+1) uu'];
 
end
 

[xx,tt]=meshgrid(x,t);
 figure(1)
surf(xx,tt,U)
shading flat
 
xlabel('espacio')
 
ylabel('tiempo')
 
zlabel('Temperatura')
Anauceda.jpg


Como se observa en la gráfica, la temperatura con respecto al tiempo en el extremo x=0 varía de forma sinusoidal con un periodo aproximado de unos 5 segundos. Sin embargo, en el otro extremo y debido a la condición Neumann establecida en el enunciado, la temperatura apenas varía.

Obtenido de «https://mat.caminos.upm.es/w/index.php?title=Ecuación_del_calor_(Grupo_8)&oldid=13598»