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