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');