Ecuación de Calor: "Calentamiento Varilla"

Trabajo realizado por estudiantes
Título	Ecuación del Calor: "Calentamiento de una varilla" (Grupo 5)
Asignatura	Ecuaciones Diferenciales
Curso	Curso 2016-17
Autores	César Blanco Posadas Sara González Bravo Lucía Granados Casado Patricia del Pozo García Marta Nogal Prata Christian Balic Stefanovic
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 Descripción del problema a estudiar

La ecuación del calor fue estudiada y propuesta por ‎Jean-Baptiste Joseph Fourier en 1807, en su memoria sobre la propagación del calor en os cuerpos sólidos. En ella propone ademas el gérmen de lo que pasaría a ser la Teoría de las Series de Fourier. La ecuación del calor es un modelo matemático (quizás el más sencillo) que trata de describir la evolución de la temperatura en un cuerpo sólido. Su expresión matemática aplicada a nuestro objeto (varilla), en este caso un objeto unidimensional, tratando como x la variable de la longitud, sería: \[\frac{\partial U}{\partial t}-\frac{\partial ^2U}{\partial x^2}=0\]

Siendo U(x,t) la temperatura de cada punto de la varilla en cada instante t, consideramos una varilla de longitud L=4 y de un cierto material de espesor constante. La varilla es conductora de calor y los extremos x=0 y x=4 están colocados sobre objetos que mantienen una temperatura constante de 5 y 0 grados respectivamente. En el instante inicial la temperatura de la varilla viene determinada por la siguiente ecuación: \[u(x,0)=e^{-3(x-3)^2}+5+x\]

Las condiciones físicas de la varilla son densidad=1 y conductividad térmica (k)=2. Suponiendo que la temperatura u(x,t) de la varilla satisface la ecuación del calor Ut-Uxx=0, el sistema completo de ecuaciones que satisface u(x,t) será: \[ \left \{\begin{matrix} Ut – 2Uxx = 0 ; t>0, x € (0,4) \\ U(0,t) = 5 ; U(4,t) = 0 ; t>0 \\ U(x,0) = e^{-3(x-3)^2}+5+x ; t>0, x € (0,4) \end{matrix} \right . \]

2 Resolución por distintos métodos

En un primer lugar estudiaremos el problema por el método de diferencias finitas aplicando el método del trapecio en este caso. Con un tamaño de paso h=0.1 para la longitud de la varilla y k=0.1 para la discretización del tiempo construimos el siguiente código:

%concidiones iniciales
a=0;%extremo izquierdo
b=4;%extremo derecho
t0=0;%instante inicial
tM=10;%instante final
q=2;%difusividad térmica
%-------------------------------------------------
ca=@(t)(0*t)+5;%función f(t) de la condición de contorno en x=a
cb=@(t)0*t;%funcion f(t) de la condición de contorno en x=b
u0=@(x)(exp(-3*((x-3).^2)))+5+x;%función f(x) valor inicial
g=@(x,t)0*x;%términos independientes
%-----------------------------------------------
%discretización intervalo[a,b]
h=0.1;
N=round((b-a)/h);%subintervalos en x
x=linspace(a,b,N+1);
x=x'%el vector x tiene N+1 elementos
xx=x(2:end-1);%tomamos los nodos intermedios
%----------------------------------------------------------
%MATRIZ K ( de N-1 x N-1 )
K=(q/h^2)*(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1)-diag(ones(1,N-2),1));
%valor inicial
U0=u0(xx);
%--------------------------------------
%discretización del tiempo
k=0.1;%tamaño de paso del tiempo
M=round((tM-t0)/k);
t=linspace(t0,tM,M+1);%vector con los tiempos
%---------------------------------------------------- 
%MATRIZ CON LAS SOLUCIONES DE LOS DISTINTOS NODOS EN CADA INSTANTE
sol=zeros(N-1,M+1);
sol(:,1)=U0;
%----------------------------------------------------------------
%APLICAMOS MÉTODO DEL TRAPECIO
for i=1:M
G1=g(xx,t(i));%Matriz términos independientes instante i
G2=g(xx,t(i+1));% matriz términos independientes instante i+1
G1(1)=G1(1)+q*ca(t(i))/h^2;
G2(1)=G2(1)+q*ca(t(i+1))/h^2;
G1(end)=G1(end)+q*cb(t(i))/h^2;
G2(end)=G2(end)+q*cb(t(i+1))/h^2;
B= eye(size(K))+((k/2)*K);
B=inv(B);
sol(:,i+1)=B*(sol(:,i)+((k/2)*(G1-K*(sol(:,i))))+(k/2)*G2);
end
%----------------------------------------------------------
Ua=ca(t);
Ub=cb(t);
sol=[Ua;sol;Ub];

%grafico
figure(1)
[Mt,Mx]=meshgrid(t,x)
surf(Mt,Mx,sol);
xlabel('t');
ylabel('x');
figure(2)
xx=2;
fila=round((xx-a)/h+1);
plot(t,sol(fila,:))

Evolución de la temperatura en la varilla frente al tiempo y la posición de los puntos de la misma.

Primeramente lo que obtenemos es una representación en 3D de la evolución de la temperatura en la varilla a lo largo de tiempo. En esta podemos ver como debido a las condiciones iniciales, efectivamente los dos extremos de la varilla mantienen constante su temperatura a lo largo del tiempo (siendo 5 para x=0 y 0 para x=4). Con respecto al resto de la varilla vemos las distintas variaciones que experimentan distintas secciones de ella.

Evolución de la temperatura en el punto medio de la varilla a lo largo del tiempo

En esta gráfica podemos observar como a lo al principio en el punto medio de la varilla se alcanza la temperatura máxima y a medida que discurre el tiempo esta va disminuyendo hasta hacerse nula.

A continuación vamos a estudiar el mismo problema por el método de Euler explícito, Euler implícito y Runge-Kutta de orden 4. Con todo esto compararemos las distintas soluciones y comentaremos las divergencias que podamos encontrar en los resultados.

Primeramente adjuntamos el código correspondiente al método de Euler explícito:

%concidiones iniciales
a=0;%extremo izquierdo
b=4;%extremo derecho
t0=0;%instante inicial
tM=10;%instante final
q=2;%difusividad térmica
%-------------------------------------------------
ca=@(t)(0*t)+5;%función f(t) de la condición de contorno en x=a
cb=@(t)0*t;%función f(t) de la condición de contorno en x=b
u0=@(x)(exp(-3*((x-3).^2)))+5+x;%función f(x) valor inicial
g=@(x,t)0*x;%términos independientes
%-----------------------------------------------
%discretización intervalo[a,b]
h=0.1;
N=round((b-a)/h);%subintervalos en x
x=linspace(a,b,N+1);
x=x'%el vector x tiene N+1 elementos
xx=x(2:end-1);%tomamos los nodos intermedios
%----------------------------------------------------------
%MATRIZ K ( de N-1 x N-1 )
K=(q/h^2)*(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1)-diag(ones(1,N-2),1));
%valor inicial
U0=u0(xx);
%--------------------------------------
%discretización del tiempo
k=0.01;%tamaño de paso del tiempo
M=round((tM-t0)/k);
t=linspace(t0,tM,M+1);%vector con los tiempos
%---------------------------------------------------- 
%MATRIZ CON LAS SOLUCIONES DE LOS DISTINTOS NODOS EN CADA INSTANTE
sol=zeros(N-1,M+1);
sol(:,1)=U0;
%----------------------------------------------------------------
%APLICAMOS MÉTODO DE EULER EXPLÍCITO
for i=1:M
Gi=g(xx,t(i));%Matriz términos independientes instante i

Gi(1)=Gi(1)+q*ca(t(i))/h^2;

Gi(end)=Gi(end)+q*cb(t(i))/h^2;


%k minúscula es tamaño de paso del tiempo
  sol(:,i+1)=sol(:,i)+k*(-K*sol(:,i)+Gi);
%en la ecuación del método pone h
%nuestra h es el tamaño de paso en la x
%el tamaño de paso en el tiempo es k, por eso ponemos k
end
%------------------------------------------------------------
Ua=ca(t);
Ub=cb(t);
sol=[Ua;sol;Ub];
%resultados
xx=2;
fila=round((xx-a)/h+1);
plot(t,sol(fila,:))

Evolución de la temperatura en el punto medio de la varilla frente al tiempo (Inestabilidad)

Evolución de la temperatura en el punto medio de la varilla frente al tiempo

Evolución de la temperatura en la varilla frente al tiempo y la posición de los puntos de la misma (Inestabilidad)

Evolución de la temperatura en la varilla frente al tiempo y la posición de los puntos de la misma

EULER IMPLICITO (cuando comprobemos si esta todo correcto ya lo comentaremos)

%concidiones iniciales
a=0;%extremo izquierdo
b=4;%extremo derecho
t0=0;%instante inicial
tM=10;%instante final
q=2;%difusividad térmica
%-------------------------------------------------
ca=@(t)(0*t)+5;%función f(t) de la condición de contorno en x=a
cb=@(t)0*t;%función f(t) de la condición de contorno en x=b
u0=@(x)(exp(-3*((x-3).^2)))+5+x;%función f(x) valor inicial
g=@(x,t)0*x;%términos independientes
%-----------------------------------------------
%discretización intervalo[a,b]
h=0.1;
N=round((b-a)/h);%subintervalos en x
x=linspace(a,b,N+1);
x=x'%el vector x tiene N+1 elementos
xx=x(2:end-1);%tomamos los nodos intermedios
%----------------------------------------------------------
%MATRIZ K ( de N-1 x N-1 )
K=(q/h^2)*(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1)-diag(ones(1,N-2),1));
%valor inicial
U0=u0(xx);
%--------------------------------------
%discretización del tiempo
k=0.1;%tamaño de paso del tiempo
M=round((tM-t0)/k);
t=linspace(t0,tM,M+1);%vector con los tiempos
%---------------------------------------------------- 
%MATRIZ CON LAS SOLUCIONES DE LOS DISTINTOS NODOS EN CADA INSTANTE
sol=zeros(N-1,M+1);
sol(:,1)=U0;
%----------------------------------------------------------------
%APLICAMOS METODO DE EULER IMPLÍCITO
for i=1:M

G2=g(xx,t(i+1));% matriz terminos independientes instante i+1
G2(1)=G2(1)+q*ca(t(i+1))/h^2;
G2(end)=G2(end)+q*cb(t(i+1))/h^2;
B= eye(size(K))+((k)*K);%matriz resultante de despejar el metodo implicito
B=inv(B);
%k minuscula es tamaño de paso del tiempo
sol(:,i+1)=B*(sol(:,i)+(k*(G2)));
%en la ecuacion del metodo pone h
%nuestra h es el tamaño de paso en la x
%el tamaño de paso en el tiempo es k, por eso ponemos k
end

%------------------------------------------------------------
Ua=ca(t);
Ub=cb(t);
sol=[Ua;sol;Ub];
%resultados
xx=2;
fila=round((xx-a)/h+1);
plot(t,sol(fila,:))

Evolución de la temperatura en el punto medio de la varilla frente al tiempo

Evolución de la temperatura en la varilla frente al tiempo y la posición de los puntos de la misma.

RUNGE KUTTA

a=0;
b=4;
t0=0;
tM=10;
q=2;
%funciones. Datos de entrada
%ca: Funcion f(t) de la condicion de contorno en x=a (t vector)
ca=@(t)0*t+5;
%cb: Funcion f(t) de la condicion de contorno en x=b (t vector)
cb=@(t)0*t;
%u0: forma analítica f(x) funcion valor inicial (x vector)
u0=@(x)(exp(-3*(x-3).^2))+5+x;
%g: funcion término independiente (x será vector, t escalar)
g=@(x,t)0*x;
Discretización del intervalo [a,b]
h=0.1;
N=round((b-a)/h); %subintervalos en x
x=linspace(a,b,N+1);
x=x';%el vector x tiene N+1 elementos

xx=x(2:end-1);

%Matriz K (de N-1 x N-1)
K=q/h^2*(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1)-diag(ones(1,N-2),1));
 
%valor inicial 
U0=u0(xx);
 
%discretización del tiempo
k=0.01;%paso en t
M=round((tM-t0)/k);%número de subintervalos en t
t=linspace(t0,tM,M+1);%vector con los tiempos 
 
sol=zeros(N-1,M+1);
sol(:,1)=U0;
 
for i=1:M
   
   G1=g(xx,t(i)); 
   G1(1)=G1(1)+q*ca(t(i))/h^2;
   G1(end)=G1(end)+q*cb(t(i))/h^2;
   
   G2=g(xx,t(i)+k/2);
   G2(1)=G2(1)+q*ca(t(i)+k/2)/h^2;
   G2(end)=G2(end)+q*cb(t(i)+k/2)/h^2;
   
   G4=g(xx,t(i+1));
   G4(1)=G4(1)+q*ca(t(i+1))/h^2;
   G4(end)=G4(end)+q*cb(t(i+1))/h^2;
   
   K1=-K*sol(:,i)+G1;
   K2=-K*(sol(:,i)+1/2*K1*k)+G2;
   K3=-K*(sol(:,i)+1/2*K2*k)+G2;
   K4=-K*(sol(:,i)+K3*k)+G4;
   sol(:,i+1)=sol(:,i)+k/6*(K1+2*K2+2*K3+K4);
%    
   
end 
 
%añadimos las condiciones Dirichlet
Ua=ca(t);
Ub=cb(t);
sol=[Ua;sol;Ub]%las añadimos

%gráfico 3D
[Mt,Mx]=meshgrid(t,x);
surf(Mt,Mx,sol);
xlabel('t');
ylabel('x');

%resultados
xx=2;
fila=round((xx-a)/h+1);
plot(t,sol(fila,:))

Evolución de la temperatura en el punto medio de la varilla frente al tiempo (Inestabilidad)

Evolución de la temperatura en el punto medio de la varilla frente al tiempo

Evolución de la temperatura en la varilla frente al tiempo y la posición de los puntos de la misma (Inestabilidad)

Evolución de la temperatura en la varilla frente al tiempo y la posición de los puntos de la misma

Compobacion Grafica (APARTADO 5)

a=0;%extremo izquierdo
b=4;%extremo derecho
t0=0;%instante inicial
tM=10;%instante final
q=2;%difusividad termica
%-------------------------------------------------
ca=@(t)(0*t)+5;%función f(t) de la condición de contorno en x=a
cb=@(t)0*t;%función f(t) de la condición de contorno en x=b
u0=@(x)(exp(-3*((x-3).^2)))+5+x;%función f(x) valor inicial
g=@(x,t)0*x;%términos indepencientes
%-----------------------------------------------
%discretización intervalo[a,b]
h=0.1;
N=round((b-a)/h);%subintervalos en x
x=linspace(a,b,N+1);
x=x';
xx=x(2:end-1);%tomamos los nodos intermedios
%----------------------------------------------------------
%MATRIZ K ( de N-1 x N-1 )
K=(q/h^2)*(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1)-diag(ones(1,N-2),1));
%valor inicial
U0=u0(xx);
%--------------------------------------
%discretización del tiempo
k=0.1;%tamaño de paso del tiempo
M=round((tM-t0)/k);
t=linspace(t0,tM,M+1);%vector con los tiempos
%---------------------------------------------------- 
%MATRIZ CON LAS SOLUCIONES DE LOS DISTINTOS NODOS EN CADA INSTANTE
sol=zeros(N-1,M+1);
sol(:,1)=U0;
%----------------------------------------------------------------
%APLICAMOS METODO DEL TRAPECIO
for i=1:M

G1=g(xx,t(i));%Matriz terminos independientes instante i
G2=g(xx,t(i+1));% matriz terminos independientes instante i+1
G1(1)=G1(1)+q*ca(t(i))/h^2;
G2(1)=G2(1)+q*ca(t(i+1))/h^2;
G1(end)=G1(end)+q*cb(t(i))/h^2;
G2(end)=G2(end)+q*cb(t(i+1))/h^2;
B= eye(size(K))+((k/2)*K);%matriz resultante de despejar el metodo implicito
B=inv(B);

sol(:,i+1)=B*(sol(:,i)+((k/2)*(G1-K*(sol(:,i))))+(k/2)*G2);

end
%------------------------------------------------------------
Ua=ca(t);
Ub=cb(t);
sol=[Ua;sol;Ub];
%------------------------------------------------------------------------
U=@(x,t)((-5/4)*x)+5;% f aproximada para valores estacionarios
%------------------------------------------------------------
%hallamos instante en que la T pasa a estacionaria
for i=1:length(t)
urt=sol(i,:);
ut=U(x,t(i));
if max(urt-ut)<=0.05
inst=t(i);
break
end
end
inst %instante en el que la temperatura permanece estacionaria
%CALCULAMOS VARIACIONES ENTRE LA SOLUCIÓN POR TRAPECIO Y LA APROXIMADA
%instantes
ta=0;
tb=1;
tc=2;
td=10;
%números de las columnas correspondientes a nuestros instantes de tiempo
cola=round((ta-t0)/k+1);
colb=round((tb-t0)/k+1);
colc=round((tc-t0)/k+1);
cold=round((td-t0)/k+1);
ut=U(x);%vector soluciones estacionarias
solu=[sol(:,cola),sol(:,colb),sol(:,colc),sol(:,cold),ut];
%  sol2 = [Tªs en t=0, Tªs en t=1,........., Tªs ESTACIONARIAS]
%DISTANCIA  para t=0
dista=(solu(:,1)-solu(:,5))';
%DISTANCIA  para t=1
distb=(solu(:,2)-solu(:,5))';
%DISTANCIA  para t=2
distc=(solu(:,3)-solu(:,5))';
%DISTANCIA  para t=10
distd=(solu(:,4)-solu(:,5))';
%diagramas----------------------------------------------------------------
subplot(2,2,1);
hold on
plot(x,ut,'r')
plot(x,solu(:,1)','b')
xlabel('longitud de la barra')
ylabel('temperatura')
title('distancia U(x,t) y solución estacionaria en t=0')
legend('solución estacionaria','solución numérica)');
hold off
%%--------------------------------------------------------------
subplot(2,2,2);
hold on
plot(x,ut,'r')
plot(x,solu(:,2)','b')
xlabel('longitud de la barra')
ylabel('temperatura')
title('distancia U(x,t) y solución estacionaria en t=1')
legend('solución estacionaria','solución numérica)');
hold off
%%--------------------------------------------------------------
subplot(2,2,3);
hold on
plot(x,ut,'r')
plot(x,solu(:,3)','b')
xlabel('longitud de la barra')
ylabel('temperatura')
title('distancia U(x,t) y solución estacionaria en t=2')
legend('solución estacionaria','solución numérica)');
hold off
%%--------------------------------------------------------------
subplot(2,2,4);
hold on
plot(x,ut,'r')
plot(x,solu(:,4)','b')
xlabel('longitud de la barra')
ylabel('temperatura')
title('distancia U(x,t) y solución estacionaria en t=10')
legend('solución estacionaria','solución numérica)');
hold off
%%-----------------------------------------------------------

FOURIER (apartado 7)

a=0;%extremo izquierdo
b=4;%extremo derecho
t0=0;%instante inicial
tM=10;%instante final
h=0.1;k=0.1;
N=round((b-a)/h);
M=round((tM-t0)/k);


x=linspace(a,b,N+1);
t=linspace(t0,tM,M+1);

u0=(exp(-3*((x-3).^2)))+5+x;

[xx,tt]=meshgrid(x,t);

%-------------------------------------------------------
U=0;
Q=1
for j=1:Q
    
    p=sin(j*pi/2*x)
    aj=trapz(x,u0.*p)/trapz(x,p.*p);
   
    U=U+aj*exp(-(j^2)*(pi^2)/4*tt).*sin(j*pi/2*xx);
end

U1=U+5-(5/4)*xx;
plot(x,U1(round(0.5/h+1),:));

hold on

Q=3;
for j=1:Q
    p=sin(j*pi/2*x)
    aj=trapz(x,u0.*p)/trapz(x,p.*p);
    U=U+aj*exp(-(j^2)*(pi^2)/4*tt).*sin(j*pi/2*xx);
end

U1=U+5-(5/4)*xx;
plot(x,U1(round(0.5/h+1),:));

Q=5;
for j=1:Q
    p=sin(j*pi/2*x)
    aj=trapz(x,u0.*p)/trapz(x,p.*p);
    U=U+aj*exp(-(j^2)*(pi^2)/4*tt).*sin(j*pi/2*xx);
end

U1=U+5-(5/4)*xx;
plot(x,U1(round(0.5/h+1),:));


Q=10;
for j=1:Q
    p=sin(j*pi/2*x)
    aj=trapz(x,u0.*p)/trapz(x,p.*p);
    U=U+aj*exp(-(j^2)*(pi^2)/4*tt).*sin(j*pi/2*xx);
end

U1=U+5-(5/4)*xx;
plot(x,U1(round(0.5/h+1),:));

Q=20;
for j=1:Q
    p=sin(j*pi/2*x)
    aj=trapz(x,u0.*p)/trapz(x,p.*p);
    U=U+aj*exp(-(j^2)*(pi^2)/4*tt).*sin(j*pi/2*xx);
end

U1=U+5-(5/4)*xx;
plot(x,U1(round(0.5/h+1),:));

legend('Q=1','Q=3','Q=5','Q=10','Q=20')
title('Serie de Fourier con 1,3,5,10,20 términos de la serie')
hold off

Con este gráfico queremos representar una vez pasadas 0.5 unidades de tiempo la temperatura a lo largo de la barra. Hemos hecho la representación usando Q= 1,3,5,10 y 20 como términos de la Serie de Fourier. La aproximación más exacta la obtenemos con Q=20 y todas son proporcionales decreciendo cumpliendo así las condiciones de contorno impuestas en el enunciado u(0,t)=5 y u(4,t)=0.