Ecuación del calor (Grupo 8)
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]
Contenido
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))
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 foco térmico. 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))
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))
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))
5 apartado 5
%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)
6 Solución apartado 6
%Datos del problema
L=3;
T=3;
%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: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
%dibujamos
[xx,tt]=meshgrid(x,t);
surf(xx,tt,sol)
7 apartado 7. Serie de Fourier
%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
Q=3
Q=5
Q=10
Q=20
Q=100
8 apartado 8
%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;
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
[xx,tt]=meshgrid(x,t);
surf(xx,tt,U)
shading flat
xlabel('espacio')
ylabel('tiempo')
zlabel('Temperatura')
9 apartado 9
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);%primA O No prima
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')
