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 | |
Contenido
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);%como cuando hacemos N=round((b-a)/h)
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
[Mt,Mx]=meshgrid(t,x)
surf(Mt,Mx,sol);
xlabel('t');
ylabel('x');
xx=2;
fila=round((xx-a)/h+1);
plot(t,sol(fila,:))