Diferencia entre revisiones de «Ecuación del calor (Grupo 8)»
| Línea 1: | Línea 1: | ||
| − | + | 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: | 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: | ||
Revisión del 18:25 15 may 2014
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
- 2 Solución mediante el método de Euler implicito
- 3 Solución mediante el método de Euler explícito
- 4 Solución mediante el método de Runge-Kutta
- 5 Estado estacionario para tiempos grandes
- 6 Solución apartado 6
- 7 apartado 7. Serie de Fourier
- 8 apartado 8
- 9 apartado 9
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 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))
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))
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)
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 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
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 Newmann, 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);%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')
