|
|
| Línea 5: |
Línea 5: |
| | | | |
| | ∂U/∂t-(∂^2 U)/(∂x^2 )=0 | | ∂U/∂t-(∂^2 U)/(∂x^2 )=0 |
| − | <math>\vec{u}=(u_1,u_2)=(\frac{\partial \psi}{\partial y},-\frac{\partial \psi}{\partial x})</math> | + | <math>\(\frac{\partial \psi}{\partial y},-\frac{\partial \psi}{\partial x})</math> |
| | | | |
| | | | |
Revisión del 13:18 9 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:
∂U/∂t-(∂^2 U)/(∂x^2 )=0
[math]\(\frac{\partial \psi}{\partial y},-\frac{\partial \psi}{\partial x})[/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í:
{█(U0(x)=10x/3 ∀ x∈ [├ 0,1)┤U(2,3]@U0(x)=100 ∀x∈[2,3]@-)┤
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.
'
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];%VER
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))
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];%VER
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];%VER
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.
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];%VER
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 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) |
|
|
|
