Modelo térmico aplicado a un edificio. Grupo 8C

Trabajo realizado por estudiantes
Título	Modelo térmico aplicado a un edificio. Grupo 8-C
Asignatura	Ecuaciones Diferenciales
Curso	Curso 2014-15
Autores	Miguel Arbaizar, Cristino Pérez , Javier Mellado , Alejandro López
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 Introducción

A lo largo de este trabajo, formularemos un modelo matemático para estudiar el comportamiento de la temperatura en el interior de un edificio durante un periodo de 24 horas. La temperatura dependerá de las siguientes variables:

• Temperatura exterior [M(t)]
• Calor producido en el interior por las personas, luces, máquinas, etc. [H(t)]
• Calentamiento y enfriamiento por la calefacción y/o aire acondicionado [U(t)]

Utilizaremos un análisis compartimental según el cual consideraremos a nuestro edificio como un solo compartimento en el que su temperatura estará designada por T(t).

2 Problema de valor inicial (P.V.I)

La razón de cambio de la temperatura quedará determinada según el problema de valor inicial o de Cauchy:

[math] \left \{ \begin{matrix} T'=k(M(t)-T)+H(t)+U(t), t\gt0 \\ T(to)=T_0 \end{matrix} \right . [/math]

Donde k es una constante real que depende de las condiciones del edificio tales como el número de ventanas, aislamiento, etc.

3 Estudio de la temperatura interior

3.1 En función exclusivamente de la temperatura exterior [M(t)]

En este caso vamos a considerar que al final del día (to=0), la temperatura exterior permanece constante (M=8), y tanto el calor producido en el interior como el aportado por la calefacción/aire acondicionado es nulo (H=U=0). La temperatura para to=0 (medianoche) será de 14 grados ºC.

La constante de tiempo del edificio nos da una medida de cuánto tarda en cambiar considerablemente la temperatura del edificio. Esta constante también es conocida como la constante de tiempo de transferencia de calor entre el edificio y el exterior. La expresión será [math]k=\frac{1}{t}[/math], que en esta primera hipótesis será [math]\frac{1}{3}[/math].

El problema de valor inicial (P.V.I) en este caso será:

[math] \left \{ \begin{matrix} T'=\frac{1}{3}*(8-T), t\gt0 \\ T(0)=14 \end{matrix} \right . [/math]

3.1.1 Método de Euler Implícito

Con la ayuda de Matlab, vamos a resolver el P.V.I. con un tamaño de paso h=0.01

%Euler implícito
clear all
%Datos del programa
t0 = 0;
tN = 24;
y0 = 14;
h = 0.01;
%Definimos la variable independiente
t = t0:h:tN;
%Guardamos los valores de la solución aproximada en el siguiente vector
y = zeros(1,length(t));
y(1) = y0;
%Bucle Euler implícito
for i = 1:length(t)-1
    y(i+1)=(1/(1+1/3*h)*(y(i)+8/3*h));
end
%grafica
figure(1)
hold on
plot(t,y)
legend('Tinterior Euler ímplícito','Location','best');
xlabel('Tiempo en horas');
ylabel('Temperatura en ºC');
axis([0,24,6,16])
hold off

centro

Como podemos observar, la temperatura del edificio parte de 14º, para después bajar rapidamente (en las 5 primeras horas del día aproximadamente) hasta estabilizarse en 8º. Es decir, nuestro edificio pierde a lo largo de la noche su calor interior hasta igualarse a la temperatura exterior.

Tambien es fácil comprobar que a mayor diferencia de temperatura entre el interior y exterior del edificio, más rápido cambiará la función T. Cuando las dos temperaturas comienzan a ser similares, alrededor de las 5 de la mañana, la temperatura interior tiende a 8º y la velocidad de enfriamiento se reduce considerablemente hasta que se iguala la temperatura exterior e interior a lo largo del resto del día.

3.2 En función de la temperatura exterior [M(t)] y el calor producido en el interior [H(t)]

3.2.1 Solución teórico-analítica

3.2.1.1 Demostración de la solución del P.V.I.

En este caso vamos a resolver analíticamente el P.V.I. siguiente: [math] \left \{ \begin{matrix} T'=k(Mo-B*cos(wt)-T)+Ho, t\gt0 \\ T(to)=T_0 \end{matrix} \right . [/math]

Operando la ecuación diferencial obtenemos: [math]T'+k*T=k*Mo-K*B*cos(wt)+Ho[/math]Se trata de una ecuación lineal de primer orden en T. La solución de este problema de valor inicial será del tipo: [math]T(t)=T_{H}(t)+T_{P}(t)[/math] En un primer lugar calculamos la homogénea: [math]T'+k*T=0 \rightarrow \frac{dT}{T}=-k*dt \rightarrow T_{H}=c*e^{-k*t}[/math] En segundo lugar calculamos la solución particular, que se realizará aplicando el principio de superposición de las soluciones: [math]f_{1}(t)=k*Mo+Ho; f_{2}(t)=-k*B*cos(wt) \rightarrow T_{P}(t)=T_{P_1}(t)+T_{P_2}(t)[/math] La solución para la primera particular será, sustituyendo en la completa: [math]T_{P_1}=A; T'_{P_1}=0 \rightarrow T'_{P_1}+k*T_{P_1}=k*Mo+Ho \rightarrow k*A=k*Mo+Ho \rightarrow A=Mo+\frac{Ho}{k} \rightarrow T_{P_1}(t)=Mo+\frac{Ho}{k}[/math] La solución para la segunda particular se obtendrá tanteando con: [math]T_{P_2}(t)=a*cos(wt)+b*sen(wt); T'_{P_2}(t)=-a*w*sen(wt)+b*w*cos(wt)[/math]Sustituyendo en la completa: [math]-a*w*sen(wt)+b*w*cos(wt)+k*a*cos(wt)+k*b*sen(wt)=-k*B*cos(wt)[/math]Identificando coeficientes obtenemos el sistema: [math] \left \{ \begin{matrix} b*w+k*a=-k*B \\ k*b-a*w=0 \end{matrix} \right . [/math] Del sistema, aplicando Cramer, obtenemos: [math]b=\frac{\begin{vmatrix} -k*B & k \\ 0 & -w \end{vmatrix}}{\begin{vmatrix} w & k \\ k & -w \end{vmatrix}}=\frac{w*k*B}{-w^2-k^2}=\frac{\frac{w}{k}*B}{1+{(\frac{w}{k})}^2};a=\frac{\begin{vmatrix} w & -k*B \\ k & 0 \end{vmatrix}}{\begin{vmatrix} w & k \\ k & -w \end{vmatrix}}=\frac{k^2*B}{-w^2-k^2}=\frac{-B}{1+{(\frac{w}{k})}^2}[/math] La segunda solucion particular será: [math]T_{P_2}(t)=\frac{-B}{1+{(\frac{w}{k})}^2}*cos(wt)-\frac{\frac{w}{k}*B}{1+{(\frac{w}{k})}^2}*sin(wt)[/math] Finalmente la solución completa será: [math]T(t)=Mo+\frac{Ho}{k}-B*{[\frac{1}{1+{(\frac{w}{k})}^2}*cos(wt)+\frac{\frac{w}{k}}{1+{(\frac{w}{k})}^2}*sin(wt)]}+c*e^{-k*t}[/math] Identificando los coeficientes obtenidos con los de la solución propuesta obtenemos los valores que pretendíamos demostrar:[math]Bo=Mo+\frac{Ho}{k};F(t)=\frac{1}{1+{(\frac{w}{k})}^2}*cos(wt)+\frac{\frac{w}{k}}{1+{(\frac{w}{k})}^2}*sin(wt)[/math] Sustituyendo la condición inicial en la completa ([math]T(0)=T_0[/math]) se obtiene el último coeficiente pedido:[math]T_0=Mo+\frac{Ho}{k}-B*{[\frac{1}{1+{(\frac{w}{k})}^2}*cos(w*0)+\frac{\frac{w}{k}}{1+{(\frac{w}{k})}^2}*sin(w*0)]}+c*e^{-k*0} \rightarrow c=T_0-Bo+\frac{B}{1+{(\frac{w}{k})}^2}[/math]

3.2.1.2 Demostración de Bo es aproximadamente igual a Mo

La temperatura media del edificio queda definida como: [math]Tmedia(t)=\frac{\int_{0}^{24} T(t)\, dt}{24}=\frac{\int_{0}^{24} (Bo-B*F(t)+c*e^{-k*t})\, dt}{24}\simeq\frac{\int_{0}^{24} Bo\, dt}{24}\simeq Bo[/math] Siendo:[math]\int_{0}^{24} F(t)\, dt=0;:\int_{0}^{24} c*e^{-k*t}\, dt\simeq 0[/math] Ya que la primera de ellas está formada por funciones periódicas y la segunda indica el enunciado que es nula.

3.2.1.3 Demostración de la variación senoidal del edificio tiene un retraso igual a [math]\boldsymbol{\alpha}[/math]

Si [math]F(t)={(1+{(\frac{w}{k})}^{2})}^{-\frac{1}{2}}*{(cos(wt)+\frac{w}{k}*sin(wt))}={(1+{(\frac{w}{k})}^{2})}^{-\frac{1}{2}}*cos{(wt-\boldsymbol{\alpha})}[/math]

Haciendo el desarrollo de [math]cos(wt)+\frac{w}{k}*sin(wt)=cos(wt)*cos(\boldsymbol{\alpha})+sin(wt)*sin(\boldsymbol{\alpha})[/math]

Llegamos a la igualdad[math] \left \{ \begin{matrix} cos(\boldsymbol{\alpha})=1 (1) \\ sen(\boldsymbol{\alpha})=\frac{w}{k} (2) \end{matrix} \right . [/math]

De donde finalmente se sigue, si dividimos [math]\frac{(2)}{(1)}[/math] obtenemos que[math]\frac{sen(\boldsymbol{\alpha})}{cos(\boldsymbol{\alpha})}=tan(\boldsymbol{\alpha})=\frac{w}{k}[/math]

3.2.2 Solución numérica

En este ocasión vamos a considerar que sí existe una razón de calentamiento adicional [H(t)], que supondremos constante a lo largo del día, a H(t)=Ho=3. De la misma manera que en el apartado anterior, el calor aportado por la calefacción/aire acondicionado es nulo (U=0), es decir, nos vamos a adentrar en el estudio de la temperatura del edificio durante las estaciones de primavera y otoño, cuando no se precisan calefactores o aparatos de aire acondicionado para alcanzar la temperatura deseada. La constante de tiempo será [math]k=\frac{1}{3}[/math] y la temperatura para to=0 (medianoche) será de 13 grados ºC.

A diferencia del caso anterior, la temperatura exterior [M(t)] no permanece constante, sino que varía en forma de onda senoidal durante un periodo de 24 horas, con su mínimo en t=0 (medianoche) y su máximo en t=12 (mediodía), de la siguiente manera: [math]M(t)=Mo-B*cos(wt)[/math]En el caso que nos ocupa, Mo recibirá el valor de 7ºC, B el valor de 5ºC y w el valor de [math]\frac{2*π}{24}=\frac{π}{12}[/math]radianes/hora. La expresión de la temperatura exterior quedará definida de la siguiente forma: [math]M(t)=7 - 5*cos(\frac{π}{12}*t)[/math]

A partir de la situación anteriormente descrita, planteamos el problema de valor inicial (P.V.I):

[math] \left \{ \begin{matrix} T'=\frac{1}{3}*(7-5*cos(\frac{π}{12}*t)-T)+3, t\gt0\\ T(0)=13 \end{matrix} \right . [/math]

3.2.2.1 Método de Euler

Con la ayuda de Matlab, vamos a resolver el P.V.I. con un tamaño de paso h1=0.1, h2=0.01 y h3=0.001. Sin embargo, por la baja resolución de la gráfica adjunta y las pequeñas diferencias milimétricas entre las funciones obtenidas para cada tamaño de paso, apenas son apreciables las diferencias entre ellas.

%Euler
clear all
%Datos del problema
t0 = 0;
tN = 24;
y0 = 13;
h1 = 0.1;
h2 = 0.01;
h3 = 0.001;
%Definimos las variables independientes
t1 = t0:h1:tN;
t2 = t0:h2:tN;
t3 = t0:h3:tN;
%Guardamos los valores de la solución aproximada en los siguientes vectores
%Euler 0.1
a = zeros(1,length(t1));
a(1) = y0;
%Euler 0.01
b = zeros(1,length(t2));
b(1) = y0;
%Euler 0.001
c = zeros(1,length(t3));
c(1) = y0;
for i = 1:length(t1)-1
    %Euler 0.1
    a(i+1) = a(i)+h1*(7/3-5/3*cos(pi/12*t1(i))-a(i)/3+3);
end
for j = 1:length(t2)-1
    %Euler 0.01
    b(j+1) = b(j)+h2*(7/3-5/3*cos(pi/12*t2(j))-b(j)/3+3);
end
for m = 1:length(t3)-1
    %Euler 0.001
    c(m+1) = c(m)+h3*(7/3-5/3*cos(pi/12*t3(m))-c(m)/3+3);
end
%Temperatura exterior
M = 7-5*cos((pi/12)*t3);
%Gráfica
figure(1)
hold on
plot(t1,a,t2,b,t3,c,t3,M)
legend('Tinterior Euler 0.1','Tinterior Euler 0.01','Tinterior Euler 0.001','Texterior','Location','best');
xlabel('Tiempo en horas');
ylabel('Temperatura en ºC');
axis([0,24,0,21])
hold off

centro

Mediante esta gráfica, podemos extraer varios aspectos interesantes de este estudio. Primeramente, podríamos decir que la temperatura exterior esta actuando sobre el edificio como modulador de su temperatura, aumentando y disminuyendo ésta de forma similar, con variaciones muy parecidas. Es importante señalar que la expresión del P.V.I. obtenida confiere a la función de la temperatura interior T un comportamiento más "equilibrado" que la temperatura exterior, pues ambas fluctúan alrededor de los 16ºC y los 7ºC respectivamente, pero abarcando un rango de temperaturas menor : Tε[12,20] y Mε[2,12].

Sin embargo, los periodos de variación de ambas funciones son iguales, de 12 horas exactamente. Tambien, debemos mencionar el retraso de α horas de la temperatura interior respecto a la exterior, que conllevan pequeños periodos durante los que la temperatura esta disminuyendo dentro del edificio y aumentando al de fuera. (Este curioso fenómeno producido se trata con mayor profundidad en la parte analítica de este trabajo "3.2.1.3 Demostración de la variación senoidal del edificio tiene un retraso igual a α").

3.2.2.2 Método de Runge-Kutta de orden 4

Con la ayuda de Matlab,también vamos a resolver el P.V.I con un tamaño de paso h1=0.1, h2=0.01 y h3=0.001. Como el P.V.I. es el mismo que el resuelto en el punto "3.2.2.1 Método de Euler", las gráficas serán casi idénticas salvo por el hecho de que Runge-Kutta es un método numérico de orden 4 y aproxima la solución mucho mejor que el método de Euler.

%RK4
clear all
%Datos del problema
t0 = 0;
tN = 24;
y0 = 13;
h1 = 0.1;
h2 = 0.01;
h3 = 0.001;
%Definimos las variables independientes
t1 = t0:h1:tN;
t2 = t0:h2:tN;
t3 = t0:h3:tN;
%Guardamos los valores de la solución aproximada en los siguientes vectores
%RK4 0.1
r = zeros(1,length(t1));
r(1) = y0;
%RK4 0.01
r1 = zeros(1,length(t2));
r1(1) = y0;
%RK4 0.001
r2 = zeros(1,length(t3));
r2(1) = y0;
for i = 1:length(t1)-1
    %RK4 0.1
    k1 = (1/3)*(7-5*cos((pi/12)*t1(i))-r(i))+3;
    k2 = (1/3)*(7-5*cos((pi/12)*(t1(i)+1/2*h1))-r(i)-1/2*k1*h1)+3;
    k3 = (1/3)*(7-5*cos((pi/12)*(t1(i)+1/2*h1))-r(i)-1/2*k2*h1)+3;
    k4 = (1/3)*(7-5*cos((pi/12)*(t1(i)+h1))-r(i)-k3*h1)+3;
    r(i+1) = r(i)+h1/6*(k1+2*k2+2*k3+k4);
end
for j = 1:length(t2)-1
    %RK4 0.01
    k1 = (1/3)*(7-5*cos((pi/12)*t2(j))-r1(j))+3;
    k2 = (1/3)*(7-5*cos((pi/12)*(t2(j)+1/2*h2))-r1(j)-1/2*k1*h2)+3;
    k3 = (1/3)*(7-5*cos((pi/12)*(t2(j)+1/2*h2))-r1(j)-1/2*k2*h2)+3;
    k4 = (1/3)*(7-5*cos((pi/12)*(t2(j)+h2))-r1(j)-k3*h2)+3;
    r1(j+1) = r1(j)+h2/6*(k1+2*k2+2*k3+k4);
end
for m = 1:length(t3)-1
    %RK4 0.001
    k1 = (1/3)*(7-5*cos((pi/12)*t3(m))-r2(m))+3; 
    k2 = (1/3)*(7-5*cos((pi/12)*(t3(m)+1/2*h3))-r2(m)-1/2*k1*h3)+3;
    k3 = (1/3)*(7-5*cos((pi/12)*(t3(m)+1/2*h3))-r2(m)-1/2*k2*h3)+3;
    k4 = (1/3)*(7-5*cos((pi/12)*(t3(m)+h3))-r2(m)-k3*h3)+3;
    r2(m+1) = r2(m)+h3/6*(k1+2*k2+2*k3+k4);
end
%Temperatura exterior
M = 7-5*cos((pi/12)*t3);
%Gráfica
figure (2)
hold on 
plot(t1,r,t2,r1,t3,r2,t3,M)
legend('Tinterior RK4 0.1','Tinterior RK4 0.01','Tinterior RK4 0.001','Texterior','Location','best');
xlabel('Tiempo en horas');
ylabel('Temperatura en ºC');
axis([0,24,0,21])
hold off

centro

Como ya hemos comentado antes de resolver el problema por este método, las gráficas obtenidas son prácticamente idénticas a las obtenidas por el método de Euler y recibirán la misma intepretación.

3.3 En función de la temperatura exterior [M(t)], el calor producido en el interior [H(t)] y el efecto de la calefacción y/o aire acondicionado [U(t)]

3.3.1 Solución teórico-analítica

3.3.1.1 Demostración de la solución del P.V.I.

En este caso vamos a resolver analíticamente el P.V.I. siguiente: [math] \left \{ \begin{matrix} T'=k(Mo-B*cos(wt)-T)+Ho+Ku*{[T_{D}-T]}, t\gt0 \\ T(to)=T_0 \end{matrix} \right . [/math]

Operando la ecuación diferencial obtenemos: [math]T'+k*T+Ku*T=Ho+k*Mo+Ku*T_{D}-K*B*cos(wt)[/math]Igual qu en el caso anterior, se trata de una ecuación lineal de primer orden en T. La solución de este problema de valor inicial será del tipo: [math]T(t)=T_{H}(t)+T_{P}(t)[/math] En un primer lugar calculamos la homogénea: [math]T'+k*T+Ku*T=0 \rightarrow \frac{dT}{T}=-{(k+Ku)}*dt \rightarrow T_{H}=c*e^{-{(k+Ku)}*t}[/math]Donde [math]K_{1}=(k+Ku)[/math]. En segundo lugar calculamos la solución particular, que se realizará aplicando el principio de superposición de las soluciones: [math]f_{1}(t)=k*Mo+Ho+Ku*T_{D}; f_{2}(t)=-k*B*cos(wt) \rightarrow T_{P}(t)=T_{P_1}(t)+T_{P_2}(t)[/math] La solución para la primera particular será, sustituyendo en la completa: [math]T_{P_1}=A; T'_{P_1}=0 \rightarrow T'_{P_1}+K_{1}*T_{P_1}=k*Mo+Ho+Ku*T_{D} \rightarrow K_{1}*A=k*Mo+Ho+Ku*T_{D} \rightarrow A=\frac{k*Mo+Ho+Ku*T_{D}}{K_{1}}[/math]La primera solución particular será: [math]T_{P_1}(t)=\frac{k*Mo+Ho+Ku*T_{D}}{K_{1}}[/math] La solución para la segunda particular se obtendrá tanteando con: [math]T_{P_2}(t)=a*cos(wt)+b*sen(wt); T'_{P_2}(t)=-a*w*sen(wt)+b*w*cos(wt)[/math]Sustituyendo en la completa: [math]-a*w*sen(wt)+b*w*cos(wt)+K_{1}*a*cos(wt)+K_{1}*b*sen(wt)=-k*B*cos(wt)[/math]Identificando coeficientes obtenemos el sistema: [math] \left \{ \begin{matrix} b*w+K_{1}*a=-k*B \\ K_{1}*b-a*w=0 \end{matrix} \right . [/math] Del sistema, aplicando Cramer, obtenemos: [math]b=\frac{\begin{vmatrix} -k*B & K_{1} \\ 0 & -w \end{vmatrix}}{\begin{vmatrix} w & K_{1} \\ K_{1} & -w \end{vmatrix}}=\frac{w*k*B}{-w^2-{K_{1}}^2}=-\frac{\frac{K*B}{K_{1}^2}*w}{1+{(\frac{w}{K_{1}})}^2};a=\frac{\begin{vmatrix} w & -k*B \\ K_{1} & 0 \end{vmatrix}}{\begin{vmatrix} w & K_{1} \\ K_{1} & -w \end{vmatrix}}=\frac{w*K_{1}*B}{-w^2-{K_{1}}^2}=\frac{\frac{k*B}{K_{1}}}{1+{(\frac{w}{K_{1}})}^2}[/math] La segunda solución particular será: [math]T_{P_2}(t)=-\frac{\frac{k*B}{K_{1}}}{1+{(\frac{w}{K_{1}})}^2}*cos(wt)-\frac{\frac{K*B}{K_{1}^2}*w}{1+{(\frac{w}{K_{1}})}^2}*sin(wt)[/math] Finalmente la solución completa será: [math]T(t)=\frac{k*Mo+Ho+Ku*T_{D}}{K_{1}}-B*{[\frac{\frac{k}{K_{1}}}{1+{(\frac{w}{K_{1}})}^2}*cos(wt)+\frac{\frac{K}{K_{1}^2}*w}{1+{(\frac{w}{K_{1}})}^2}*sin(wt)]}+c*e^{-k*t}[/math] Identificando los coeficientes obtenidos con los de la solución propuesta obtenemos los valores que pretendíamos demostrar:[math]B_{2}=\frac{k*Mo+Ho+Ku*T_{D}}{K_{1}};F_{1}(t)=\frac{cos(wt)+\frac{w}{K_{1}}*sin(wt)}{1+{(\frac{w}{K_{1}})}^2};B_{1}=\frac{k*B}{K_{1}}[/math] Sustituyendo la condición inicial en la completa ([math]T(0)=T_0[/math]) se obtiene el último coeficiente pedido:[math]T_0=\frac{k*Mo+Ho+Ku*T_{D}}{K_{1}}-\frac{k*B}{K_{1}}*{[\frac{cos(w*0)+\frac{w}{K_{1}}*sin(w*0)}{1+{(\frac{w}{K_{1}})}^2}]}+c*e^{-k*0} \rightarrow c=T_0-B_{2}+B_{1}*F_{1}(0)[/math]

3.3.1.2 Demostración de Bo es aproximadamente igual a Mo

La temperatura media del edificio queda definida como: [math]Tmedia(t)=\frac{\int_{0}^{24} T(t)\, dt}{24}=\frac{\int_{0}^{24} (B_{2}-B_{1}*F_{1}(t)+c*e^{-K_{1}*t})\, dt}{24}\simeq\frac{\int_{0}^{24} B_{2}\, dt}{24}\simeq B_{2}[/math] Siendo:[math]\int_{0}^{24} B_{2}*F_{1}(t)\, dt=0;:\int_{0}^{24} c*e^{-K_{1}*t}\, dt\simeq 0[/math] Ya que la primera de ellas está formada por funciones periódicas y la segunda indica el enunciado que es nula.

[math]Temperatura media\simeq B_{2}\simeq \frac{k*Mo+Ho+Ku*T_{D}}{K_{1}}\simeq T_{D} \leftrightarrow \left \{ \begin{matrix} \frac{Ku}{K_{1}}=\frac{Ku}{k+Ku}\simeq 1 & \mbox{si }k\mbox{\lt\lt\ltKu} \\ \frac{k*Mo+Ho}{K_{1}}\simeq 0 \end{matrix} \right . [/math]

3.3.2 Solución numérica

A partir de este momento el estudio de la variación de la temperatura en el interior del edificio lo realizaremos suponiendo que está afectado por todas las variables. Existe una razón de calentamiento adicional [H(t)], que seguimos suponiendo constante (H=3). La temperatura para t0=0 (medianoche) será de 24 grados ºC.

La temperatura exterior [M(t)] sigue variando de forma de senoidal durante un periodo de 24 horas, con su mínimo en t=0 (medianoche) y su máximo en t=12 (mediodía), de la siguiente manera: [math]M(t)=Mo-B*cos(wt)[/math]En esta situación, Mo recibirá el valor de 4ºC, B el valor de 5ºC y w el valor de [math]\frac{2*π}{24}=\frac{π}{12}[/math]radianes/hora. La expresión de la temperatura exterior quedará definida de la siguiente forma: [math]M(t)=4 - 5*cos(\frac{π}{12}*t)[/math]

Suponemos que en el edificio del apartado anterior está instalado un termostato sencillo que se utiliza para comparar la temperatura real interior del edificio con una temperatura deseada [math]T_{D}[/math]. Si la temperatura es menor que la deseada, el calefactor proporciona calentamiento; si es mayor se mantiene apagado. Si la temperatura es mayor que la deseada, el aire acondicionado proporciona enfriamiento; en caso contrario se encuentra apagado. Suponemos que la cantidad de calor o enfriamiento suministrado es proporcional a la diferencia de temperatura [math]U(t)=Ku*[T_{D}−T(t)][/math], donde Ku es una constante positiva. En esta situación la constante Ku tiene un valor de [math]k=\frac{7}{8}[/math] y la temperatura deseada es de 22 grados ºC.

Por lo tanto el P.V.I. que modeliza el esta situación será:

[math] \left \{ \begin{matrix} T'=\frac{1}{3}*(4-5*cos(\frac{π}{12}*t)-T)+3+\frac{7}{8}*(22-T), t\gt0\\ T(0)=24 \end{matrix} \right . [/math]

3.3.2.1 Método de Euler

Con la ayuda de Matlab, vamos a resolver el P.V.I con un tamaño de paso h1=0.1, h2=0.01 y h3=0.001. Sin embargo, por la baja resolución de la gráfica adjunta y las pequeñas diferencias milimétricas entre las funciones obtenidas para cada tamaño de paso, apenas son apreciables las diferencias entre ellas.

%Euler
clear all
%Datos del problema
t0 = 0;
tN = 24;
y0 = 24;
h1 = 0.1;
h2 = 0.01;
h3 = 0.001;
%Definimos las variables independientes
t1 = t0:h1:tN;
t2 = t0:h2:tN;
t3 = t0:h3:tN;
%Guardamos los valores de la solución aproximada en los siguientes vectores
%Euler 0.1
a = zeros(1,length(t1));
a(1) = y0;
%Euler 0.01
b = zeros(1,length(t2));
b(1) = y0;
%Euler 0.001
c = zeros(1,length(t3));
c(1) = y0;
for i = 1:length(t1)-1
    %Euler 0.1
    a(i+1) = a(i)+h1*(4/3-5/3*cos(pi/12*t1(i))-a(i)/3+3+7/8*(22-a(i)));
end
for j = 1:length(t2)-1
    %Euler 0.01
    b(j+1) = b(j)+h2*(4/3-5/3*cos(pi/12*t2(j))-b(j)/3+3+7/8*(22-b(j)));
end
for m = 1:length(t3)-1
    %Euler 0.001
    c(m+1) = c(m)+h3*(4/3-5/3*cos(pi/12*t3(m))-c(m)/3+3+7/8*(22-c(m)));
end
%Temperatura exterior
M = 7-5*cos((pi/12)*t3);
%Grafica
figure(1)
hold on
plot(t1,a,t2,b,t3,c,t3,M)
legend('Tinterior Euler 0.1','Tinterior Euler 0.01','Tinterior Euler 0.001','Texterior','Location','best');
xlabel('Tiempo en horas');
ylabel('Temperatura en ºC');
axis([0,24,0,25])
hold off

centro

En esta gráfica se observa como la temperatura exterior sigue actuando sobre el edificio como modulador de su temperatura, aumentando y disminuyendo aunque de forma mucho menos acusada que en el caso anterior. Esto es debido al efecto tanto de la calefacción como del aire acondicionado, entre ambos mantienen la temperatura próxima a la temperatura deseada [math]T_{D}[/math]. La temperatura media del interior del edificio está en torno a unos 20ºC, cerca de la [math]T_{D}[/math] pero que no se llega a alcanzar ya que la temperatura media exterior es de 6ºC lo que hace difícil calentar el edificio hasta la [math]T_{D}[/math]. Durante las horas centrales del día, cuando la temperatura exterior se maximiza y alcanza unos 13ºC, la temperatura interior del edificio alcanza la [math]T_{D}[/math].

Sin embargo, los periodos de variación de ambas funciones son iguales, de 12 horas exactamente. Tambien, debemos mencionar el retraso de [math]\alpha[/math] horas de la temperatura interior respecto a la exterior, que conllevan pequeños periodos durante los que la temperatura está aumentando dentro del edificio y disminuyendo la exterior, aunque del mismo modo este fenómeno es mucho menor que en el anterior caso ya que está atenuado por la calefacción y/o aire acondicionado y prácticamente es imperceptible.

3.3.2.2 Método de Runge-Kutta de orden 4

Con la ayuda de Matlab, también vamos a resolver el P.V.I con un tamaño de paso h1=0.1, h2=0.01 y h3=0.001. Como el P.V.I. es el mismo que el resuelto en el punto "3.3.2.1 Método de Euler", las gráficas serán casi idénticas salvo por el hecho de que Runge-Kutta es un método numérico de orden 4 y aproxima la solución mucho mejor que el método de Euler.

%RK4
clear all
%Datos del problema
t0 = 0;
tN = 24;
y0 = 24;
h1 = 0.1;
h2 = 0.01;
h3 = 0.001;
%Definimos las variables independientes
t1 = t0:h1:tN;
t2 = t0:h2:tN;
t3 = t0:h3:tN;
%Guardamos los valores de la solución aproximada en los siguientes vectores
%RK4 0.1
r = zeros(1,length(t1));
r(1) = y0;
%RK4 0.01
r1 = zeros(1,length(t2));
r1(1) = y0;
%RK4 0.001
r2 = zeros(1,length(t3));
r2(1) = y0;
for i = 1:length(t1)-1
    %RK4 0.1
    k1 = (1/3)*(4-5*cos((pi/12)*t1(i))-r(i))+3+7/8*(22-r(i)); 
    k2 = (1/3)*(4-5*cos((pi/12)*(t1(i)+1/2*h1))-r(i)-1/2*k1*h1)+3+7/8*(22-r(i)-1/2*k1*h1);
    k3 = (1/3)*(4-5*cos((pi/12)*(t1(i)+1/2*h1))-r(i)-1/2*k2*h1)+3+7/8*(22-r(i)-1/2*k2*h1);
    k4 = (1/3)*(4-5*cos((pi/12)*(t1(i)+h1))-r(i)-k3*h1)+3+7/8*(22-r(i)-k3*h1);
    r(i+1) = r(i)+h1/6*(k1+2*k2+2*k3+k4);
end
for j = 1:length(t2)-1
    %RK4 0.01
    k1 = (1/3)*(4-5*cos((pi/12)*t2(j))-r1(j))+3+7/8*(22-r1(j)); 
    k2 = (1/3)*(4-5*cos((pi/12)*(t2(j)+1/2*h2))-r1(j)-1/2*k1*h2)+3+7/8*(22-r1(j)-1/2*k1*h2);
    k3 = (1/3)*(4-5*cos((pi/12)*(t2(j)+1/2*h2))-r1(j)-1/2*k2*h2)+3+7/8*(22-r1(j)-1/2*k1*h2);
    k4 = (1/3)*(4-5*cos((pi/12)*(t2(j)+h2))-r1(j)-k3*h2)+3+7/8*(22-r1(j)-k3*h2);
    r1(j+1) = r1(j)+h2/6*(k1+2*k2+2*k3+k4);
end
for m = 1:length(t3)-1
    %RK4 0.001
    k1 = (1/3)*(4-5*cos((pi/12)*t3(m))-r2(m))+3+7/8*(22-r2(m)); 
    k2 = (1/3)*(4-5*cos((pi/12)*(t3(m)+1/2*h3))-r2(m)-1/2*k1*h3)+3+7/8*(22-r2(m)-1/2*k1*h3);
    k3 = (1/3)*(4-5*cos((pi/12)*(t3(m)+1/2*h3))-r2(m)-1/2*k2*h3)+3+7/8*(22-r2(m)-1/2*k1*h3);
    k4 = (1/3)*(4-5*cos((pi/12)*(t3(m)+h3))-r2(m)-k3*h3)+3+7/8*(22-r2(m)-k3*h3);
    r2(m+1) = r2(m)+h3/6*(k1+2*k2+2*k3+k4);
end
%Temperatura exterior
M = 7-5*cos((pi/12)*t3);
%Grafica
figure (2)
hold on 
plot(t1,r,t2,r1,t3,r2,t3,M)
legend('Tinterior RK4 0.1','Tinterior RK4 0.01','Tinterior RK4 0.001','Texterior','Location','best');
xlabel('Tiempo en horas');
ylabel('Temperatura en ºC');
axis([0,24,0,25])
hold off

centro

Como ya hemos comentado antes de resolver el problema por este método, las gráficas obtenidas son prácticamente idénticas a las obtenidas por el método de Euler y recibirán la misma intepretación.

4 Estudio de la temperatura interior del edificio separado en dos zonas, A y B

Ahora vamos a suponer que el edificio está dividido en dos zonas, zona A y zona B. La temperatura en la zona A estará designada por [Ta(t)] en tanto que la temperatura de la zona B estará designada por [Tb(t)]. La temperatura inicial en la zona A será de 20 grados ºC mientras que la temperatura inicial de la zona B será de 18 grados ºC

En la zona A la razón de calentamiento adicional será [Ha(t)] y el efecto de la calefacción y/o aire acondicionado [Ua(t)]. En la zona B la razón de calentamiento adicional será [Hb(t)] y el efecto de la calefacción y/o aire acondicionado [Ub(t)]. En este caso, para las zonas A y B respectivamente: [math] \left \{ \begin{matrix} Ha(t)=10*\frac{t}{1+t} \\ Ua(t)=\frac{5}{3}*(22-Ta(t)) \end{matrix} \right . [/math]

[math] \left \{ \begin{matrix} Hb(t)=4*cos(\frac{π}{6}*t)\\ Ub(t)=\frac{13}{7}*(23-Tb(t)) \end{matrix} \right . [/math]

La temperatura exterior variará en forma de onda senoidal durante un periodo de 24 horas, con su mínimo en t=0 (medianoche) y su máximo en t=12 (mediodía), de la siguiente manera tanto en la zona A como B: [math]M(t)=Mo-B*cos(wt)[/math]En esta situación será igual tanto para la zona A como para la zona B y su valor será: [math]M(t)=2 - 7*cos(\frac{π}{12}*t)[/math]

La constante de tiempo de cada zona nos da una medida de cuánto tarda en cambiar considerablemente la temperatura de dicha zona. Esta constante también es conocida como la constante de tiempo de transferencia de calor entre el edificio y el exterior. Esta constante tendrá un valor de [math]\frac{1}{4}[/math] entre la zona A y el exterior, un valor de [math]\frac{1}{5}[/math] entre la zona B y el exterior y un valor de [math]\frac{1}{2}[/math] entre la zona A y la zona B.

Para la resolución de este caso, tendremos que resolver un sistema de ecuaciones diferenciales, que será el siguiente:

[math] \left \{ \begin{matrix} Ta'= \frac{1}{4}*(2-7*cos(\frac{π}{12}*t)-Ta) + 10*\frac{t}{1+t} + \frac{5}{3}*(22-Ta) + \frac{1}{2}*(Tb-Ta) \\ Tb'= \frac{1}{5}*(2-7*cos(\frac{π}{12}*t)-Tb) + 4*cos(\frac{π}{6}*t) + \frac{13}{7}*(23-Tb) - \frac{1}{2}*(Tb-Ta) \end{matrix} \right . [/math]

4.1 Método de Euler

Con la ayuda de Matlab, vamos a resolver el sistema de ecuaciones diferenciales con un tamaño de paso h1=0.1, h2=0.01 y h3=0.001. Sin embargo, por la baja resolución de la gráfica adjunta y las pequeñas diferencias milimétricas entre las funciones obtenidas para cada tamaño de paso, apenas son apreciables las diferencias entre ellas.

centro

%Euler
clear all
%Datos del problema
t0 = 0;
tN = 24;
y0 = [20,18]';
h1 = 0.1;
h2 = 0.01;
h3 = 0.001;
%Definimos las variables independientes
t1 = t0:h1:tN;
t2 = t0:h2:tN;
t3 = t0:h3:tN;
%Guardamos los valores de la solución aproximada en los siguientes vectores
%Euler 0.1
y = zeros(2,length(t1));
y(:,1) = y0;
%Euler 0.01
y1 = zeros(2,length(t2));
y1(:,1) = y0;
%Euler 0.001
y2 = zeros(2,length(t3));
y2(:,1) = y0;
%Definimos la matriz de coeficientes de Ta y Tb
A = [-29/12,1/2;1/2,-179/70];
%Bucles
for i = 1:length(t1)-1
    %Euler 0.1
    y(:,i+1) = y(:,i)+h1*(A*y(:,i)+[(223/6)+(10*t1(i)/(1+t1(i)))-(7/4)*cos((pi/12)*t1(i));(1509/35)-(7/5)*cos((pi/12)*t1(i))+4*cos((pi/6)*t1(i))]);
end
for j = 1:length(t2)-1
    %Euler 0.01
    y1(:,j+1) = y1(:,j)+h2*(A*y1(:,j)+[(223/6)+(10*t2(j)/(1+t2(j)))-(7/4)*cos((pi/12)*t2(j)),(1509/35)-(7/5)*cos((pi/12)*t2(j))+4*cos((pi/6)*t2(j))]');
end
for m = 1:length(t3)-1
    %Euler 0.001
    y2(:,m+1) = y2(:,m)+h3*(A*y2(:,m)+[(223/6)+(10*t3(m)/(1+t3(m)))-(7/4)*cos((pi/12)*t3(m)),(1509/35)-(7/5)*cos((pi/12)*t3(m))+4*cos((pi/6)*t3(m))]');
end
%Temperatura exterior
M = 2-7*cos((pi/12)*t3);
figure(1)
hold on
plot(t1,y(1,:),t1,y(2,:),t2,y1(1,:),t2,y1(2,:),t3,y2(1,:),t3,y2(2,:),t3,M)
legend('Ta Euler 0.1','Tb Euler 0.1','Ta Euler 0.01','Tb Euler 0.01','Ta Euler 0.001','Tb Euler 0.001','Texterior','Location','best');
xlabel('Tiempo en horas');
ylabel('Temperatura en ºC');
axis([0,24,-7,27])
hold off

Inicialmente la temperatura de los dos edificios varía de forma similar hasta que en un cierto instante de tiempo se igualan. Durante el transcurso del tiempo, la temperatura de la zona A (gráfica verde) tiende a estabilizarse, aunque sigue apreciándose una leve variación periódica que es influida por la variación en forma de onda senoidal de la temperatura exterior. Sin embargo, la temperatura de la zona B (gráfica azul) sufre variaciones periódicas de su temperatura mucho más acusadas ya que en este lado del edificio, a demás de la variación de la temperatura exterior, el calor producido en el interior también tiene una variación en forma de onda senoidal. Finalmente, podemos observar que la temperatura máxima de ambas zonas se produce cuando la temperatura exterior alcanza también su máximo.

4.2 Método de Runge-Kutta de orden 4

Con la ayuda de Matlab, también vamos a resolver el sistema de ecuaciones diferenciales con un tamaño de paso h1=0.1, h2=0.01 y h3=0.001. Como el sistema es el mismo que el resuelto en el punto "4.1 Método de Euler", las gráficas serán casi idénticas salvo por el hecho de que Runge-Kutta es un método numérico de orden 4 y aproxima la solución mucho mejor que el método de Euler.

%RK4
clear all
%Datos del problema
t0 = 0;
tN = 24;
y0 = [20,18]';
h1 = 0.1;
h2 = 0.01;
h3 = 0.001;
%Definimos las variables independientes
t1 = t0:h1:tN;
t2 = t0:h2:tN;
t3 = t0:h3:tN;
%Guardamos los valores de la solución aproximada en los siguientes vectores
%RK4 0.1
z = zeros(2,length(t1));
z(:,1) = y0;
%RK4 0.01
z1 = zeros(2,length(t2));
z1(:,1) = y0;
%RK4 0.001
z2 = zeros(2,length(t3));
z2(:,1) = y0;
%Definimos la matriz de coeficientes de Ta y Tb
A = [-29/12,1/2;1/2,-179/70];
%Bucles
for i = 1:length(t1)-1
    %RK4 0.1
    k1 = A*z(:,i)+[(223/6)+(10*t1(i)/(1+t1(i)))-(7/4)*cos((pi/12)*t1(i));(1509/35)-(7/5)*cos((pi/12)*t1(i))+4*cos((pi/6)*t1(i))];
    k2 = A*(z(:,i)+h1/2*k1)+[(223/6)+(10*(t1(i)+h1/2)/(1+(t1(i)+h1/2))-(7/4)*cos((pi/12)*(t1(i)+h1/2)));(1509/35)-(7/5)*cos((pi/12)*(t1(i)+h1/2))+4*cos((pi/6)*(t1(i)+h1/2))];
    k3 = A*(z(:,i)+h1/2*k2)+[(223/6)+(10*(t1(i)+h1/2)/(1+(t1(i)+h1/2))-(7/4)*cos((pi/12)*(t1(i)+h1/2)));(1509/35)-(7/5)*cos((pi/12)*(t1(i)+h1/2))+4*cos((pi/6)*(t1(i)+h1/2))];
    k4 = A*(z(:,i)+h1*k3)+[(223/6)+(10*(t1(i)+h1)/(1+(t1(i)+h1))-(7/4)*cos((pi/12)*(t1(i)+h1)));(1509/35)-(7/5)*cos((pi/12)*(t1(i)+h1))+4*cos((pi/6)*(t1(i)+h1))];
    z(:,i+1) = z(:,i)+h1/6*(k1+2*k2+2*k3+k4);
end
for j = 1:length(t2)-1
    %RK4 0.01
    k1 = A*z1(:,j)+[(223/6)+(10*t2(j)/(1+t2(j)))-(7/4)*cos((pi/12)*t2(j));(1509/35)-(7/5)*cos((pi/12)*t2(j))+4*cos((pi/6)*t2(j))]; 
    k2 = A*(z1(:,j)+h2/2*k1)+[(223/6)+(10*(t2(j)+h2/2)/(1+(t2(j)+h2/2))-(7/4)*cos((pi/12)*(t2(j)+h2/2)));(1509/35)-(7/5)*cos((pi/12)*(t2(j)+h2/2))+4*cos((pi/6)*(t2(j)+h2/2))];
    k3 = A*(z1(:,j)+h2/2*k2)+[(223/6)+(10*(t2(j)+h2/2)/(1+(t2(j)+h2/2))-(7/4)*cos((pi/12)*(t2(j)+h2/2)));(1509/35)-(7/5)*cos((pi/12)*(t2(j)+h2/2))+4*cos((pi/6)*(t2(j)+h2/2))];
    k4 = A*(z1(:,j)+h2*k3)+[(223/6)+(10*(t2(j)+h2)/(1+(t2(j)+h2))-(7/4)*cos((pi/12)*(t2(j)+h2)));(1509/35)-(7/5)*cos((pi/12)*(t2(j)+h2))+4*cos((pi/6)*(t2(j)+h2))];
    z1(:,j+1) = z1(:,j)+h2/6*(k1+2*k2+2*k3+k4);
end
for m = 1:length(t3)-1
    %RK4 0.001
    k1 = A*z2(:,m)+[(223/6)+(10*t3(m)/(1+t3(m)))-(7/4)*cos((pi/12)*t3(m));(1509/35)-(7/5)*cos((pi/12)*t3(m))+4*cos((pi/6)*t3(m))]; 
    k2 = A*(z2(:,m)+h3/2*k1)+[(223/6)+(10*(t3(m)+h3/2)/(1+(t3(m)+h3/2))-(7/4)*cos((pi/12)*(t3(m)+h3/2)));(1509/35)-(7/5)*cos((pi/12)*(t3(m)+h3/2))+4*cos((pi/6)*(t3(m)+h3/2))];
    k3 = A*(z2(:,m)+h3/2*k2)+[(223/6)+(10*(t3(m)+h3/2)/(1+(t3(m)+h3/2))-(7/4)*cos((pi/12)*(t3(m)+h3/2)));(1509/35)-(7/5)*cos((pi/12)*(t3(m)+h3/2))+4*cos((pi/6)*(t3(m)+h3/2))];
    k4 = A*(z2(:,m)+h3*k3)+[(223/6)+(10*(t3(m)+h3)/(1+(t3(m)+h3))-(7/4)*cos((pi/12)*(t3(m)+h3)));(1509/35)-(7/5)*cos((pi/12)*(t3(m)+h3))+4*cos((pi/6)*(t3(m)+h3))];
    z2(:,m+1) = z2(:,m)+h3/6*(k1+2*k2+2*k3+k4);
end
%Temperatura exterior
M = 2-7*cos((pi/12)*t3);
%Gráfica
figure(2)
hold on
plot(t1,z(1,:),t1,z(2,:),t2,z1(1,:),t2,z1(2,:),t3,z2(1,:),t3,z2(2,:),t3,M)
legend('Ta RK4 0.1','Tb RK4 0.1','Ta RK4 0.01','Tb RK4 0.01','Ta RK4 0.001','Tb RK4 0.001','Texterior','Location','best');
xlabel('Tiempo en horas');
ylabel('Temperatura en ºC');
axis([0,24,-7,27])
hold off

centro

Como ya hemos comentado antes de resolver el problema por este método, las gráficas obtenidas son prácticamente idénticas a las obtenidas por el método de Euler y recibirán la misma intepretación.

4.3 Método de Euler Implícito

Con la ayuda de Matlab, también vamos a resolver el sistema de ecuaciones diferenciales con un tamaño de paso h1=0.1, h2=0.01 y h3=0.001. Como el sistema es el mismo que el resuelto en el punto "4.1 Método de Euler", las gráficas serán casi idénticas salvo por el hecho de que Euler Implícito es un método numérico de orden 2 y aproxima la solución mejor que el método de Euler.

centro

%Euler implícito
clear all
%Datos del problema
t0 = 0;
tN = 24;
y0 = [20,18]';
h1 = 0.1;
h2 = 0.01;
h3 = 0.001;
%Definimos las variables independientes
t1 = t0:h1:tN;
t2 = t0:h2:tN;
t3 = t0:h3:tN;
%Guardamos los valores de la solución aproximada en los siguientes vectores
%Euler Implícito 0.1
e = zeros(2,length(t1));
e(:,1) = y0;
%Euler 0.01
e1 = zeros(2,length(t2));
e1(:,1) = y0;
%Euler 0.001
e2 = zeros(2,length(t3));
e2(:,1) = y0;
%Definimos la matriz de coeficientes de Ta y Tb
A = [-29/12,1/2;1/2,-179/70];
%Bucles
for i = 1:length(t1)-1
    %Euler implícito 0.1
    b = e(:,i)+h1*[(223/6)+(10*t1(i+1)/(1+t1(i+1)))-(7/4)*cos((pi/12)*t1(i+1));(1509/35)-(7/5)*cos((pi/12)*t1(i+1))+4*cos((pi/6)*t1(i+1))];
    e(:,i+1) = (eye(2)-h1*A)\b;
end
for j = 1:length(t2)-1
    %Euler implícito 0.01
    b1 = e1(:,j)+h2*[(223/6)+(10*t2(j+1)/(1+t2(j+1)))-(7/4)*cos((pi/12)*t2(j+1));(1509/35)-(7/5)*cos((pi/12)*t2(j+1))+4*cos((pi/6)*t2(j+1))];
    e1(:,j+1) = (eye(2)-h2*A)\b1;
end
for m = 1:length(t3)-1
    %Euler implícito 0.001
    b2 = e2(:,m)+h3*[(223/6)+(10*t3(m+1)/(1+t3(m+1)))-(7/4)*cos((pi/12)*t3(m+1));(1509/35)-(7/5)*cos((pi/12)*t3(m+1))+4*cos((pi/6)*t3(m+1))];
    e2(:,m+1) = (eye(2)-h3*A)\b2;
end
%Temperatura exterior
M = 2-7*cos((pi/12)*t3);
%Gráfica
figure(3)
hold on
plot(t1,e(1,:),t1,e(2,:),t2,e1(1,:),t2,e1(2,:),t3,e2(1,:),t3,e2(2,:),t3,M)
legend('Ta Euler implícito 0.1','Tb Euler implícito 0.1','Ta Euler implícito 0.01','Tb Euler implícito 0.01','Ta Euler implícito 0.001','Tb Euler implícito 0.001','Texterior','Location','best');
xlabel('Tiempo en horas');
ylabel('Temperatura en ºC');
axis([0,24,-7,27])
hold off

Como ya hemos comentado antes de resolver el problema por este método, las gráficas obtenidas son prácticamente idénticas a las obtenidas por el método de Euler y recibirán la misma intepretación.