Ecuación de Calor en varilla(Grupo 16-C)

Trabajo realizado por estudiantes
Título	Ecuación de Calor en varilla. Grupo 16-C
Asignatura	Ecuaciones Diferenciales
Curso	Curso 2014-15
Autores	Laura Ramos Sangrós Luis Moreno López Álvaro Martínez de Andrés Gonzalo Olmos de la Cruz
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

Se nos presenta una varilla delgada, homogénea, de sección constante y longitud L ocupando el intervalo x € (0,L) con L = 4 m. En los extrema los de la barra están colocados objetos a una temperatura constante de 0 y 5 grados respectivamente, e inicialmente la temperatura de la varilla viene dada por la función: u0 = e^(-6(x-2)^2) + (5-5x/4) Suponiendo que la temperatura de la varilla u(x,t) satisface la ecuación de calor ut-2uxx = 0, El problema (P) con su sistema de ecuaciones que satisface u(x,t) será:

        Ut – 2Uxx = 0   ;                                        t>0, x € (0,4)
        U(0,t) = 5 ; U(4,t) = 0  ;                               t>0
        U(x,0) = e^(-6(x-2)^2) + (5-5x/4)  ;                     t>0, x € (0,4)

En segundo lugar se nos pide resolverlo por el método de diferencias finitas con un tamaño de paso y un intervalo de tiempo de 0.1 por el método del trapecio.

El código en Matlab y su correspondiente gráfico será:

%Trabajo ejercicio 2 
clear all , clc , clf

%Resolverlo por trapecio
%datos iniciales
a=0; b=4; h=0.1;
x=a:h:b;
N=(b-a)/h; %el vector x tiene n+1 elementos
q=2; %término que acompaña a Uxx
%nodos interiores
xx=x(2:N);
xx=xx';
%cond contorno
ua=5; ub=0;
%cond inicial
u0=exp(-6*(xx-2).^2)+5-5*xx/4;
% matriz K
K=q/h^2*(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1) -diag(ones(1,N-2),1));
% término F 
F=0*xx;
F(1)=F(1)+q*ua/h^2;
F(end)=F(end)+q*ub/h^2;
% discretización en t
t0=0;tM=1;
k=0.1; % paso temporal.
t=t0:k:tM;
M=length(t)-1; %número de subintervalos temporales
sol(:,1)=u0;
for i=1:M
%Trapecio
sol(:,1+i)=(eye(size(K))+(k/2)*K)\(sol(:,i)+(k/2)*(-K*sol(:,i)+2*F));
end
%Añado frontera
Ua=q*ua*ones(1,length(t));
Ub=q*ub*ones(1,length(t));
sol=[Ua;sol;Ub];

figure(1)
[Mt,Mx]=meshgrid(t,x); 
surf(Mx,Mt,sol)
xlabel('Longitud');
ylabel('Tiempo');
zlabel('Temperatura');

A continuación se nos pide ver cómo se comporta la temperatura en el punto medio de la barra(X=2). Para ello congelaremos en un gráfico 2D (temperatura/tiempo) en X=2 a lo largo del tiempo.Usando el método de diferencias finitas y trapecio,el código será similar al apartado anterior salvo unos cambios al final:

%Trabajo ejercicio 3 
clear all , clc , clf

%Resolverlo por trapecio

%datos iniciales
a=0; b=4; h=0.1;
x=a:h:b;
N=(b-a)/h; %el vector x tiene n+1 elementos
q=2; %término que acompaña a Uxx
%nodos interiores
xx=x(2:N);
xx=xx';
%cond contorno
ua=5; ub=0;
%cond inicial
u0=exp(-6*(xx-2).^2)+5-5*xx/4;
% matriz K
K=q/h^2*(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1) -diag(ones(1,N-2),1));
% término F 
F=0*xx;
F(1)=F(1)+q*ua/h^2;
F(end)=F(end)+q*ub/h^2;
% discretización en t
t0=0;tM=1;
k=0.1; % paso temporal.
t=t0:k:tM;
M=length(t)-1; %número de subintervalos temporales
sol(:,1)=u0;
for i=1:M
%Trapecio
sol(:,1+i)=(eye(size(K))+(k/2)*K)\(sol(:,i)+(k/2)*(-K*sol(:,i)+2*F));
end
%Añado frontera
Ua=q*ua*ones(1,length(t));
Ub=q*ub*ones(1,length(t));
sol=[Ua;sol;Ub];

figure(2)
%Comportamiento de la temperatura en el punto medio de la barra
hold on
long=(b-a)/2;
fila=round((long-a)/h+1);
plot(t,sol(fila,:))
xlabel('Tiempo');
ylabel('Temperatura');
hold off

Es el punto donde más variación de temperatura se observa.

Terminamos la primera parte resolviendo este problema (P) con las mismas condiciones que las anteriores y por el método de diferencias finitas, mediante la resolución de los métodos de Euler eplícito, Euler implícito y Runge-Kutta de orden 4, para más tarde comparar sus gráficas como se pide en el enunciado. A primera vista las gráficas obtenidas serán:

%Trabajo ejercicio 4 
clear all , clc , clf

%datos iniciales
a=0; b=4; h=0.1;
x=a:h:b;
N=(b-a)/h; %el vector x tiene n+1 elementos
q=2; %término que acompaña a Uxx
%nodos interiores
xx=x(2:N);
xx=xx';
%cond contorno
ua=5; ub=0;
%cond inicial
u0=exp(-6*(xx-2).^2)+5-5*xx/4;
% matriz K
K=q/h^2*(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1) -diag(ones(1,N-2),1));
% término F
F=0*xx;
F(1)=F(1)+q*ua/h^2;
F(end)=F(end)+q*ub/h^2;
% discretización en t
t0=0;tM=1;
k=0.1 % paso temporal.
t=t0:k:tM;
M=length(t)-1; %número de subintervalos temporales
sol(:,1)=u0;
sol1(:,1)=u0;
sol2(:,1)=u0;
for i=1:M
%Euler explicito
sol(:,1+i)=sol(:,i)+k*(-K*sol(:,i)+F);
%Euler implicito
sol1(:,1+i)=(eye(size(K))+k*K)\(sol1(:,i)+k*F);
%Kutta
k1=-K*sol2(:,i)+F;
k2=-K*(sol2(:,i)+(1/2)*k1*k)+F;
k3=-K*(sol2(:,i)+(1/2)*k2*k)+F;
k4=-K*(sol2(:,i)+k3*k)+F;
sol2(:,i+1)=sol2(:,i)+(k/6)*(k1+k2*2+k3*2+k4);
end
%incluimos cond contorno
Ua=q*ua*ones(1,length(t));
Ub=q*ub*ones(1,length(t));
sol=[Ua;sol;Ub];
sol1=[Ua;sol1;Ub];
sol2=[Ua;sol2;Ub];
%grafico para comparar
figure(1)
[Mt,Mx]=meshgrid(t,x); 
mesh(Mx,Mt,sol) 
title('Euler explícito')
xlabel('Longitud');
ylabel('Tiempo');
zlabel('Temperatura');
figure(2)
mesh(Mx,Mt,sol1)
title(' Euler implícito')
xlabel('Longitud');
ylabel('Tiempo');
zlabel('Temperatura');
figure(3)
mesh(Mx,Mt,sol2)
title('Runge-Kutta')
xlabel('Longitud');
ylabel('Tiempo');
zlabel('Temperatura');

Como se puede observar, la solución no tiene nada que ver con la realidad,solo concuerda en el caso de Euler Implícito,lo cual induce a pensar que hay un error en el desarrollo de la modelización. Si se revisa el programa, se advierte que la distorsión se debe a que con las condiciones Δx y Δt dadas el problema no cumple la condición de estabilidad necesaria: Δx/(Δt)^2 Y en este caso no se cumple por lo tanto basta con buscar un Δt que haga que cumpla: Δx/(Δt)^2 = k/h^2 = 1/2Δt , con lo que Δt = 0.0025. El nuevo código sería:

%Trabajo ejercicio 4 
clear all , clc , clf

%datos iniciales
a=0; b=4; h=0.1;
x=a:h:b;
N=(b-a)/h; %el vector x tiene n+1 elementos
q=2; %término que acompaña a Uxx
%nodos interiores
xx=x(2:N);
xx=xx';
%cond contorno
ua=5; ub=0;
%cond inicial
u0=exp(-6*(xx-2).^2)+5-5*xx/4;
% matriz K
K=q/h^2*(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1) -diag(ones(1,N-2),1));
% término F
F=0*xx;
F(1)=F(1)+q*ua/h^2;
F(end)=F(end)+q*ub/h^2;
% discretización en t
t0=0;tM=1;
k=0.0025; % paso temporal.
t=t0:k:tM;
M=length(t)-1; %número de subintervalos temporales
sol(:,1)=u0;
sol1(:,1)=u0;
sol2(:,1)=u0;
for i=1:M
%Euler explicito
sol(:,1+i)=sol(:,i)+k*(-K*sol(:,i)+F);
%Euler implicito
sol1(:,1+i)=(eye(size(K))+k*K)\(sol1(:,i)+k*F);
%Kutta
k1=-K*sol2(:,i)+F;
k2=-K*(sol2(:,i)+(1/2)*k1*k)+F;
k3=-K*(sol2(:,i)+(1/2)*k2*k)+F;
k4=-K*(sol2(:,i)+k3*k)+F;
sol2(:,i+1)=sol2(:,i)+(k/6)*(k1+k2*2+k3*2+k4);
end
%incluimos cond contorno
Ua=q*ua*ones(1,length(t));
Ub=q*ub*ones(1,length(t));
sol=[Ua;sol;Ub];
sol1=[Ua;sol1;Ub];
sol2=[Ua;sol2;Ub];
%grafico para comparar
figure(1)
[Mt,Mx]=meshgrid(t,x); 
mesh(Mx,Mt,sol) 
xlabel('Longitud');
ylabel('Tiempo');
zlabel('Temperatura');
figure(2)
mesh(Mx,Mt,sol1)
xlabel('Longitud');
ylabel('Tiempo');
zlabel('Temperatura');
figure(3)
mesh(Mx,Mt,sol2)
title('Kutta');
xlabel('Longitud');
ylabel('Tiempo');
zlabel('Temperatura');

Variación de la temperatura en el tiempo

Otra conclusión interesante que se extrae de este problema es que para tiempos grandes la temperatura apenas cambia en el tiempo. En este caso el término ut se puede despreciar y la temperatura toma un valor estacionario en cada punto de la varilla. Esta ecuación estacionaria será la siguiente al despreciar este término: 2Uxx = 0 Para demostrarlo gráficamente, basta con utilizar el método de Euler implícito, tanteando con distintos tiempos. El código es el mismo que antes pero al variar el valor de t(t=0,1,2,10), se obtienen las siguientes gráficas.

Se puede apreciar que la temperatura apenas varía en el tiempo a partir de t = 10, se mantiene constante y sólo varía en los primeros instantes. La función a la que se parece la solución u(x,t) para tiempos grandes es: ye=5-(5/4)*x, resuelta analíticamente a partir de la ecuación estacionaria antes dicha (2Uxx= 0). La función gráficamente es(es idéntica para cualquier tiempo):

Así la distancia entre la solución transitoria y estacionaria será la resta entre ambas soluciones --> err=abs(ye-sol1) quue en matlab dará una tabla con las distancias en cada instante t por paso(0.1).