Ver el código de «Modelo Térmico de un Edificio.(Grupo 13A)»


{{ TrabajoED | Modelo Térmico de un Edificio.(Grupo 13A) | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED14/15|Curso 2014-15]] | María Aguilera, Paula Martínez, Miguel Sánchez, Laura García, Isabel Roselló, Sarah Boufounas }}
[[Categoría:Ecuaciones Diferenciales]]
[[Categoría:ED14/15]]
[[Categoría:Trabajos 2014-15]]

==Introducción==
Este trabajo consiste en estudiar, mediante un modelo matématico, el comportamiento de la temperatura interior de un edificio en un periodo de 24 horas. Esta temperatura dependerá del calor generado dentro, de la temperatura exterior y del sistema de calefacción o aire acondicionado. Designaremos estos factores como:
:* H(t) → Calor generado por las personas, luces, máquinas...
:* M(t) → Temperatura exterior
:* U(t) → Calentamiento producido por la calefacción o enfriamiento debido al aire acondicionado

==Problema de Cauchy==
Designando la temperatura del edificio como T(t) el problema de valor inicial queda como:
[[Archivo:Jeje2.PNG|centro]]
Donde k es una constante del edificio que depende del número de puertas, ventanas, aislamiento, etc.
===Problema de Cauchy particularizado===
Denominamos t'=1/k al tiempo en que el edificio sufre un cambio de temperatura considerable, este tiempo es denominado como la constante de tiempo del edificio o constante de tiempo de transferencia de calor entre el edificio y el exterior.
Respondiendo a la primera pregunta el tiempo que tarda en cambiar considerablemente la temperatura para un edificio cerrado esta entre 2 y 4 horas.
En este caso particular tomaremos un tiempo t'=3 por lo que la constante k es igual a ⅓, el tiempo a medianoche es t=0 con una To=14, la temperatura exterior M(t) es constante e igual a 8 ºC y el calor generado por las personas H(t) y por la calefaccion U(t) es igual a 0, y nos queda el problema siguiente:
[[Archivo:Jeje3.PNG|centro]]
===Euler Implícito para el problema de valor inicial===
Tomamos una longitud de paso h=0.01 y comparamos el resultado con la ecuación exacta que resolviendo el problema analíticamente da T=6*e^(-1/3*t)+8 y podemos comprobar que es prácticamente idéntica.
[[Archivo:Jojo2.JPG|centro]]
[[Archivo:Jojo1.JPG|centro]]

==Modelo con una variación senoidal==
En este apartado consideramos la razon de calentamiento  <math>\ H(t)=Ho=cte </math>, no hay calefaccion o enfriamiento U(t)=0, y la temperatura exterior M(t) varia de forma senoidal en un perido de 24 horas con un mínimo en t=0 horas y un máximo en t=12 horas. <math>\ M(t)=Mo - Bcos(ωt) </math> con Mo y B constantes >0 y <math>\ ω=\frac{π}{12}\ </math> radianes/hora y un <math>\ T(0)=To </math>.
Para demostrar que:  <math> T(t)=Bo - BF(t) + Ce^{-kt} </math> siendo \[\left\{\begin{matrix}\ Bo = Mo + \frac{Ho}{k}\ , & \\ F(t)= \frac{\cos(ωt)+\frac{ω}{k}\sin(ωt)}{1+\frac{ω^2}{k^2}}\ \, \\ C= To - Bo + \frac{B}{1+\frac{ω^2}{k^2}}\\ & \end{matrix}\right.\] 
Tenemos que resolver el problema de valor inicial: \[\left\{\begin{matrix}\ T'=k[Mo-Bcos(ωt)-T]+Ho, \\ T(0)=To & \end{matrix}\right.\]
Es una ecuación lineal de primer orden y la resolvemos como <math> Tgeneral=Thomogenea+Tparticular </math>.
:Para la homogénea <math> T'+kT=0 </math> la solución es <math> T_{h}=Ce^{-kt} </math>
:Para resolver la solución particular aplicamos el principio de la superposición y obtenemos dos soluciones particulares una para <math> y_{1}(t)=kMo+Ho </math> y otra para <math> y_{2}(t)=-kBcos(ωt) </math>
Tanteamos una solucion <math> T_{p1}=Α \ y \ T'_{p1}=0 </math> y resolviendo nos da <math> Α=Mo+\frac{Ho}{k}\ </math> y <math> T_{p1}=Mo+\frac{Ho}{k}\ </math>
Y seguidamente para la segunda solucion particular tanteamos <math> T_{p2}=fcos(ωt)+gsin(ωt) </math> que metiendola en la ecuación <math> T'_{p2}+kT_{p2}=-kBcos(ωt)</math> y resolviendo por Cramer obtenemos <math> f=\frac{ωkB}{-{ω}^2-{k}^2}\ </math> <math> g=\frac{{k}^2B}{-{ω}^2-{k}^2}\ </math> y por tanto <math> T_{p2}=\frac{-B}{1+{(\frac{ω}{k})}^2}cos(ωt)-\frac{\frac{ω}{k}B}{1+{(\frac{ω}{k})}^2}sin(ωt)</math>


La solución <math> T(t)=T_{h}+T_{p1}+T_{p2} </math> nos queda <math> T(t)=Mo+\frac{Ho}{k}\ -B[\frac{1}{1+{(\frac{ω}{k})}^2}cos(ωt)+\frac{\frac{ω}{k}}{1+{(\frac{ω}{k})}^2}sin(ωt)] +Ce^{-kt} </math> e imponiendo la condición inicial <math> T(0)=To </math> se obtiene un <math> C=To -Bo+\frac{B}{1+{(\frac{ω}{k})}^2}\ </math> exactamente como dice el enunciado.

El siguiente apartado pide demostrar que <math> Bo </math> es aproximadamente la temperatura media diaria dentro del edificio, y que con <math> Ho=0 </math> entonces <math> Bo=Mo </math>

Así que planteamos la siguiente ecuación → <math> T_{media}(t) = \frac{\int_0^{24}{T(t)dt}}{24}\ = \frac{\int_0^{24}{Bo-BF(t)+Ce^{-kt})dt}}{24}\ ≈ \frac{\int_0^{24}{Bo dt}}{24}\ ≈ Bo </math>   al ser  <math> \int_0^{24}{F(t)dt} = 0 </math> y <math>  \int_0^{24}{Ce^{-kt}dt} ≈ 0 </math>

En el tercer apartado tenemos que demostrar que la variación senoidal del edificio tiene un retraso de '''α''' horas respecto a la variación exterior con un <math> F(t) = (1+\frac{ω^2}{k^2})^{-1/2} cos(ωt-α), t>0 </math> y que <math> tan(α)=\frac{ω}{k}\ </math>

De modo que nos queda <math> F(t) = (1+\frac{ω^2}{k^2})^{-1/2}(cos(ωt) + \frac{ω}{k}\sin(ωt)) = (1+\frac{ω^2}{k^2})^{-1/2}cos(ωt-α)  </math> → <math> cos(ωt) + \frac{ω}{k}\sin(ωt) = cos(ωt)cos(α) + sin(ωt)sin(α) </math>

\[\left\{\begin{matrix}\ cos(ωt)=cos(ωt)cos(α) & \frac{w}{k}\sin(ωt)=sin(ωt)sin(α) & \end{matrix}\right.\]   \[\left\{\begin{matrix}\ 1=cos(α) & \frac{w}{k}=sin(α) & \end{matrix}\right.\]

Y por lo tanto <math> tan(α)=\frac{ω}{k}\ </math>

==Modelo con una variación senoidal con termostato==
===Demostración del enunciado===
Analizamos a continuación un segundo caso. Supongamos ahora que el edificio dispone de un termostato que adecúa la temperatura interior del edificio en función de la temperatura deseada (TD). Para ello el sistema mide la temperatura real y aporta calor o frío al interior, o, en caso de corresponder con la temperatura deseada, permanece apagado.  La cantidad de calor o enfriamiento suministrado vendrá dado por la función U(t), proporcional a esta diferencia de temperatura. 
 <math> U(t)=Ku [Td - T(t)] </math>
Demostraremos a continuación que la solución del nuevo problema de valor inicial (PVI) es:
<math> T(t) = B2-B1F1(t)+Ce^{-k1t} </math>
Donde:
\[\left\{\begin{matrix}\ K1 = k+Ku\ , & \\ B2 = \frac{KuTd+kMo+Ho}{K1}\ \, \\ B1 = \frac{kB}{K1}\ \,  \\ F1(t)= \frac{\cos(ωt)+\frac{ω}{k1}\sin(ωt)}{1+\frac{ω^2}{k^2}}\ \, \\ C= To - B2 + F1(0) \ & \end{matrix}\right.\]   
Para demostrarlo paso a paso comenzamos resolviendo el problema de valor inicial: <math> T' = k[(Mo-Bcos(ωt)]-T]+Ho+Ku[Td-T] </math>
Es una ecuación lineal de primer orden y la resolvemos como <math> Tgeneral=Thomogenea+Tparticular </math>.
:La solución de la ecuación homogénea <math> T'+KT + KuT = 0 </math>, empleando el método de separación de variables es <math> T = Ce^{-(k+Ku)} </math> . Como podemos observar, el término <math> k+Ku=K1 </math> . Queda así resuelta la primera parte de esta cuestión. 
:Para resolver la solución particular aplicamos el principio de la superposición y obtenemos dos soluciones particulares una para <math> y_{1}(t)=kMo+Ho+KuTd </math> y otra para <math> y_{2}(t)=-kBcos(ωt) </math>
Empleando el método de resolución por tanteo en <math> T_{p1}=Α \ y \ T'_{p1}=0 </math> obtenemos <math> A = \frac{KuTd+kMo+Ho}{k + Ku}\ </math> y <math> T_{p1}= \frac{KuTd+kMo+Ho}{K1}\ </math>.
Como podemos ver, la expresión B2 coincide con esta primera solución particular hallada.
Del mismo modo, para la segunda solución particular tanteamos <math> T_{p2}=fcos(ωt)+gsin(ωt) </math> . Derivando esta expresión  math> T_{p2}=-fsen(ωt)+gωcos (ωt) </math> y sustituyéndola en la EDOL llegamos a un sistema de ecuaciones que es fácilmente resoluble por Cramer, obteniendo los valores de f y g y sustituyéndolos en la anterior expresión de ''' Tp2'''. La nueva expresión queda:
<math> T_{p2} (t)=-\frac{\frac{kB}{K_{1}}}{1+{(\frac{ω}{K_{1}})}^2}cos(ωt)-\frac{\frac{KB}{K_{1}^2}ω}{1+{(\frac{ω}{K_{1}})}^2}sin(ωt) </math>

La solución <math> T(t)=T_{h}+T_{p1}+T_{p2} </math> nos queda <math>T(t)=\frac{kMo+Ho+kUT_{D}}{K_{1}}-B1{[\frac{\frac{k}{K_{1}}}{1+{(\frac{ω}{K_{1}})}^2}cos(ωt)+\frac{\frac{k}{K_{1}^2}ω}{1+{(\frac{ω}{k_{1}})}^2}sin(ωt)]}+Ce^{-kt}</math>
Como podemos observar, el factor que multiplica a B1 es  ''' F1(t)'''.
Finalmente, imponiendo la condición inicial <math> T(0)=To </math> obtenemos el término correspondiente al término C dado por el enunciado:
<math> C = To – B2 + +F1(0) </math>

===Representación de la temperatura interior y exterior del edificio===

Mediante el método de Euler y de Runge-kutta 4 hacemos un programa en Matlab que nos representa la solución para los datos:
H0=3, M0=7, k=1/3, B=5, T0=13

{{matlab|codigo=
%DATOS DEL PROBLEMA
t0=0;
tN=24;
y0=13;
H0=3;
M0=7;
k=1/3;
B=5;
w=pi/12;
B0=M0+(H0/k);
C=y0-B0+(B/(1+(w/k)^2));
h=input('Introduce tamaño de paso');
N=(tN-t0)/h;
% Definimos la variable independiente
t=t0:h:tN;

F=(1/(1+(w/k)^2))*(cos(w*t)+(sin(w*t)*(w/k)));
%solucion exacta
x=B0-B*F+C*(exp(-k*t));
%Ahora vamos a guardar los valores de las soluciones aproximadas en varios vectores (diferentes métodos)
y=zeros(1,N+1); %euler
y(1)=y0;

u=zeros(1,N+1); %runge-kutta
u(1)=y0;


for i=1:N
 y(i+1)=y(i)+h*((k*((M0-(B*cos(w*t(i))))-y(i)))+3); %euler

 U1=(k*((M0-(B*cos(w*t(i))))-y(i)))+3;
 U2=(k*((M0-(B*cos(w*(t(i)+(0.5*h))))-(y(i)+(0.5*U1*h))))+3);
 U3=k*((M0-(B*cos(w*(t(i)+(0.5*h))))-(y(i)+(0.5*U2*h))))+3;
 U4=k*((M0-(B*cos(w*(t(i)+(h))))-(y(i)+(U3*h))))+3;
 u(i+1)=u(i)+h/6*(U1+2*U2+2*U3+U4); %runge-kutta4
end
M=M0-(B*cos((w)*t));
figure(1)
hold on
plot(t,y,'b')
plot(t,u,'g')
plot(t,M,'r')
plot(t,x,'y')
legend('Temperatura interior, Euler','Temperatura interior, Runge-kutta','Temperatura exterior','Location','best')
xlabel('Tiempo en horas');
ylabel('Temperatura en ºC');
hold off
Mediay=mean(y)
Mediau=mean(u)
MediaM=mean(M)
Mediax=mean(x)
Maxy=max(y)
Maxu=max(u)
MaxM=max(M)
Maxx=max(x)
Miny=min(y)
Minu=min(u)
MinM=min(M)
Minx=min(x)
}}

'''Grafico Euler y Runge-kutta 4 h=0.1 :'''
[[Archivo:2h01.png|centro]]

[[Archivo:T2h0.png|centro]]

La variación de la temperatura interior del edificio en este gráfico la hemos obtenido mediante los métodos de Euler y Runge-Kutta 4 con un h=0.1. 
Se observa en la gráfica que ambos métodos nos proporcionan una función muy similar a la solución exacta. La temperatura máxima interior obtenida con Euler es 19.96 º C y con Runge-kutta es 19.99 º C,  alcanzadas a las 14.5 horas y 14.6 horas respectivamente. Este dato se aproxima mucho a la temperatura máxima de la solución exacta, que es 19.93  Así mismo la temperatura mínima obtenida por Euler es 12.08 ºC y por Runge-Kutta es 12.12ºC, alcanzadas ambas a las 2.6 horas, también muy próximas a los 12.11 ºC de la solución exacta.
El periodo de la función es de 24 horas y el rango en el que varían las temperaturas es de 7.874. La temperatura media obtenida es   15.9996 y  16.0322 para cada uno de los métodos y 16 en el caso de la exacta. 
En cuanto a la variación de la temperatura exterior, tiene su mínimo, 2ºC, a medianoche (0 horas y 24 horas) y su máximo, 12 grados a las 12 horas (mediodía) y también tiene un periodo de 24 horas. El rango en el que varían las temperaturas en este caso es 10.
La temperatura exterior influye en la temperatura interior del edificio de manera que a medida que se produce una variación en la exterior, se produce una variación similar pero atenuada  y con cierto retraso en la temperatura interior del edificio.

'''Grafico Euler y Runge-Kutta 4 h=0.01 :'''
[[Archivo:2h001.png|centro]]

[[Archivo:T2h00.png|centro]]

Al emplear un tamaño de paso menor que en el anterior caso, observamos que la diferencia entre la función representada por el método de Euler y la representada por el método de Runge-kutta 4 es inferior que en el caso anterior, y también su variación con respecto a la solución exacta como se puede comprobar en los datos concretos como la temperatura máxima, la mínima y la media.

'''Grafico Euler y Runge-Kutta 4 h=0.001 :'''
[[Archivo:2h0001.png|centro]]

[[Archivo:T2h000.png|centro]]

Del mismo modo que ocurrió en el caso anterior, al disminuir aún más el tamaño de paso las soluciones obtenidas se acercan más a la realidad, y por tanto los datos como las temperaturas mínima, máxima y media coinciden en ambos métodos y en este caso también coinciden con los de la exacta.


===Aproximación de la temperatura media a B2 y a la temperatura deseada===
A continuación vamos a demostrar que la temperatura media es igual a B2

Suponemos un valor medio de 24 horas:
[[Archivo:Susu1.PNG|centro]]
[[Archivo:Susu3.PNG|centro]]

Seguidamente el ejercio nos pide que demostremos que la temperatura media está muy próxima a la temperatura deseada.
[[Archivo:Susu4.PNG|centro]]

===Representación de la temperatura interior y exterior del edificio===
Mediante el método de Euler y de Runge-kutta 4 hacemos un programa en Matlab que nos representa la solución para los datos:
H0=3, M0=4, k=1/3, B=5, T0=24,KU=7/8, TD=22

{{matlab|codigo=
%DATOS DEL PROBLEMA
t0=0;
tN=24;
y0=24;
H0=3;
M0=4;
k=1/3;
B=5;
KU=7/8;
TD=22;
w=pi/12;
B0=M0+(H0/k);

K1=k+KU;
B2=(1/K1)*((KU*TD)+(k*M0)+H0);
B1=(k*B)/K1;
h=input('Introduce tamaño de paso');
N=(tN-t0)/h;

% Definimos la variable independiente
t=t0:h:tN;

F1=(1/(1+(w/K1)^2))*(cos(w*t)+(sin(w*t)*(w/K1)));
C=y0-B2+(1/(1+(w/K1)^2))*(cos(w*0)+(sin(w*0)*(w/K1)));
%solucion exacta
x=B2-B1*F1+C*(exp(-K1*t));
%Ahora vamos a guardar los valores de las soluciones aproximadas en varios vectores (diferentes métodos)
y=zeros(1,N+1); %euler
y(1)=y0;

u=zeros(1,N+1); %runge-kutta
u(1)=y0;
F1=zeros(1,N+1);
for i=1:N
y(i+1)=y(i)+h*((k*((M0-(B*cos(w*t(i))))-y(i)))+3+(KU*(TD-y(i)))); %euler

 U1=(k*((M0-(B*cos(w*t(i))))-y(i)))+3+(KU*(TD-y(i)));
 U2=(k*((M0-(B*cos(w*(t(i)+(0.5*h))))-(y(i)+(0.5*U1*h))))+3+(KU*(TD-(y(i)+(0.5*U1*h)))));
 U3=k*((M0-(B*cos(w*(t(i)+(0.5*h))))-(y(i)+(0.5*U2*h))))+3+(KU*(TD-(y(i)+(0.5*U2*h))));
 U4=k*((M0-(B*cos(w*(t(i)+(h))))-(y(i)+(U3*h))))+3+(KU*(TD-(y(i)+(U3*h))));
 u(i+1)=u(i)+h/6*(U1+2*U2+2*U3+U4); %runge-kutta4
end
M=M0-(B*cos((w)*t));
figure(1)
hold on
plot(t,y,'b')
plot(t,u,'g')
plot(t,M,'r')
plot(t,x,'y')
legend('Temperatura interior, Euler','Temperatura interior, Runge-kutta','Temperatura exterior','Location','best')
xlabel('Tiempo en horas');
ylabel('Temperatura en ºC');
Mediay=mean(y)
Mediau=mean(u)
MediaM=mean(M)
Mediax=mean(x)
Maxy=max(y)
Maxu=max(u)
MaxM=max(M)
Maxx=max(x)
Miny=min(y)
Minu=min(u)
MinM=min(M)
Minx=min(x)
}}

'''Grafico Euler y Runge-Kutta 4 h=0.1 :'''
[[Archivo:3h01.png|centro]]

[[Archivo:T3h0.png|centro]]

La variación de la temperatura interior del edificio en este gráfico la hemos obtenido mediante los métodos de Euler y Runge-Kutta 4 con un h=0.1. 
Se observa en la gráfica que ambos métodos nos proporcionan una función muy similar a la solución exacta. La temperatura máxima interior obtenida con Euler es 24 º C y con Runge-kutta es también 24 º C. Este dato es bastante póximo a la temperatura máxima de la solución exacta, que es 23.64ºC.  Así mismo la temperatura mínima obtenida por Euler es 18.19 ºC y por Runge-Kutta es 18.53ºC, también muy próximas a los 18.2 ºC de la solución exacta.
El periodo de la función es de 24 horas y el rango en el que varían las temperaturas es de 7.874. La temperatura media obtenida es   19.71 y  20.04  para cada uno de los métodos y 19.7 en el caso de la exacta. 
Por tanto observamos que para este tamaño de paso la solución más aproximada nos la proporciona el método de Euler.
En cuanto a la variación de la temperatura exterior, tiene su mínimo, 2ºC, a medianoche (0 horas y 24 horas) y su máximo, 12 grados a las 12 horas (mediodía) y también tiene un periodo de 24 horas. El rango en el que varían las temperaturas en este caso es 10.
Podemos observar que la temperatura exterior influye en la temperatura interior del edificio de manera que a medida que se produce una variación en la exterior, se produce una variación similar pero atenuada  y con cierto retraso en la temperatura interior del edificio, pero en este caso, al haber un aparato que regula la temperatura de forma que se aproxime a una temperatura deseada observamos que el rango de temperaturas en el que varía la temperatura interior del edificio es de 5.43, mucho menor que en el del apartado anterior y que la temperatura media es de 19.7 ºC valor bastante cercano a los 22ºC de la temperatura deseada.

'''Grafico Euler y Runge-Kutta 4 h=0.01 :'''
[[Archivo:3h001.png|centro]]

[[Archivo:T3h00.png|centro]]

Al emplear un tamaño de paso menor que en el anterior caso, observamos que la diferencia entre la función representada por el método de Euler y la representada por el método de Runge-kutta 4 es inferior que en el caso anterior, y también su variación con respecto a la solución exacta como se puede comprobar en los datos concretos como la temperatura máxima, la mínima y la media.
En este caso también obtenemos resultados más aproximados con el método de Euler.

'''Grafico Euler y Runge-Kutta 4 h=0.001 :'''
[[Archivo:3h0001.png|centro]]

[[Archivo:T3h000.png|centro]]

Del mismo modo que ocurrió en el caso anterior, al disminuir aún más el tamaño de paso las soluciones obtenidas se acercan más a la realidad, y por tanto los datos como las temperaturas mínima, máxima y media coinciden en ambos métodos y son prácticamente iguales a los de la exacta.

==Estudio de la temperatura del edificio dividido en dos zonas==

El edificio consta ahora de dos zonas A y B. Designaremos a la temperatura en la zona A por Ta(t), y la temperatura en la zona B por Tb(t). Como dato del enunciado sabemos que la temperatura inicial en la zona A es de 20 grados y  la de la zona B es de 18 grados.
Ha(t) designara a la razón de calentamiento adicional en la zona A y  Ua(t) el efecto de la calefacción y/o aire acondicionado.
Asi mismo, Hb(t) designara a la razón de calentamiento adicional en la zona B  y  Ub(t) el efecto de la calefacción y/o aire acondicionado.

Sabemos que las constantes de tiempo de transferencia de calor son de 4 h entre A y el exterior, 5 h entre B y el exterior y 2 h entre A y B. 
Despejando de t=1/k obtenemos : 1/4 , 1/5 y 1/2 respectivamente.

Recordamos que el problema de valor inicial (PVI) o problema de Cauchy es el siguiente sistema de dos ecuaciones:

[[Archivo:Mariaecus2.PNG]]
 
La temperatura exterior  variara en forma de onda senoidal durante un t[0,24] de este modo :
M(t) - Mo - B * cos ( wt ) y nos queda que : M(t) - 2- 7 * cos (( pi/12)*t)

[[Archivo:Mariaecus3.PNG]]

En nuestro caso, el problema de Cauchy que satisface las temperaturas Ta(t) y Tb(t) es:

[[Archivo:Mariaecus1.PNG]]

===Solución numérica===

Para resolver el sistema de ecuaciones resultante nos ayudaremos de matlab. Se ha creado un programa que nos dará las soluciones a este sistema y nos dibujará las gráficas de las temperaturas (temperatura interior en zona A, temperatura interior en zona B y temperatura exterior) en función del tiempo durante 24 horas. Para ello deberemos dar un número de paso (h). En nuestro caso se ha utilizado h=0.1;h=0.01;h=0.001. Cuanto mayor sea el número de paso, mayor será la exactitud de la solución pues más se acercará esta a la solución exacta, la cual se lograría integrando. Se utilizan tres métodos distintos: Euler, Euler implícito y Runge-Kutta de orden 4. Cada uno de ellos nos dará una solución distinta, pues todos son métodos aproximados, aunque todos ellos nos darán soluciones muy similares, que serán también similares a la solución exacta.
Por último cabe señalar que al ser un sistema en lugar de una sola ecuación diferencial estos programas deberán de operan con matrices.

{{matlab|codigo=
%Datos del problema
t0=0; tN=24; h=0.001; y0=[20;18]; %Variaremos h dándole como valores 0.1;0.01;0.001
%Definimos la variable independiente
t=t0:h:tN ;
%Preparamos la matriz de la solución
x=zeros(2,length(t)); %Primera fila Temperatura A, segunda fila Temperatura B, Euler
%Inicializamos
x(:,1)=y0;
y=zeros(2,length(t)); %Primera fila Temperatura A, segunda fila Temperatura B,  Euler modificado
%Inicializamos
y(:,1)=y0;
z=zeros(2,length(t)); %Primera fila Temperatura A, segunda fila Temperatura B, Runge Kutta
%Inicializamos
z(:,1)=y0;
%Aplicamos el método
for i =1:length(t)-1
%Comenzacmos con Euler
 x(:,i+1)=x(:,i)+h*[(1/4)*(2-7*cos((pi/12)*t(i))-x(1,i))+((10*t(i))/(1+t(i)))+(5/3)*(22-x(1,i))+1/2*(x(2,i)-x(1,i)) ; 1/5*(2-7*cos(pi/12*t(i))-x(2,i))+4*cos(pi/6*t(i))+13/7*(22-x(2,i))-1/2*(x(2,i)-x(1,i))];
%Ahora Euler implícito
K=[(1/4)*(2-7*cos((pi/12)*t(i))-y(1,i))+((10*t(i))/(1+t(i)))+(5/3)*(22-y(1,i))+1/2*(y(2,i)-y(1,i)) ; 1/5*(2-7*cos(pi/12*t(i))-y(2,i))+4*cos(pi/6*t(i))+13/7*(22-y(2,i))-1/2*(y(2,i)-y(1,i))];
y(:,i+1)=y(:,i)+h*[(1/4)*(2-7*cos((pi/12)*t(i))-(y(1,i))+1/2*h*K(1,1))+((10*t(i))/(1+t(i)))+(5/3)*(22-(y(1,i)+1/2*h*K(1,1)))+1/2*((y(2,i)+1/2*h*K(2,1))-(y(1,i)+1/2*h*K(1,1))) ; 1/5*(2-7*cos(pi/12*t(i))-(y(2,i)+1/2*h*K(2,1)))+4*cos(pi/6*t(i))+13/7*(22-(y(2,i)+1/2*h*K(2,1)))-1/2*((y(2,i)+1/2*h*K(2,1))-(y(1,i)+1/2*h*K(1,1)))];
%Ahora Runge-Kutta
K1=[1/4*(2-7*cos((pi/12)*t(i))-z(1,i))+(10*t(i)/(1+t(i)))+5/3*(22-z(1,i))+1/2*(z(2,i)-z(1,i)) ; 1/5*(2-7*cos(pi/12*t(i))-z(2,i))+4*cos(pi/6*t(i))+13/7*(22-z(2,i))-1/2*(z(2,i)-z(1,i))];
K2=[1/4*(2-7*cos((pi/12)*t(i))-(z(1,i)+1/2*h*K1(1)))+(10*t(i)/(1+t(i)))+5/3*(22-(z(1,i)+1/2*h*K1(1)))+1/2*((z(2,i)+1/2*h*K1(2))-(z(1,i)+1/2*h*K1(1))) ; 1/5*(2-7*cos(pi/12*t(i))-(z(2,i)+1/2*h*K1(2)))+4*cos(pi/6*t(i))+13/7*(22-(z(2,i)+1/2*h*K1(2)))-1/2*((z(2,i)+1/2*h*K1(2))-(z(1,i)+1/2*h*K1(1)))];
K3=[1/4*(2-7*cos((pi/12)*t(i))-(z(1,i)+1/2*h*K2(1)))+(10*t(i)/(1+t(i)))+5/3*(22-(z(1,i)+1/2*h*K2(1)))+1/2*((z(2,i)+1/2*h*K2(2))-(z(1,i)+1/2*h*K2(1))) ; 1/5*(2-7*cos(pi/12*t(i))-(z(2,i)+1/2*h*K2(2)))+4*cos(pi/6*t(i))+13/7*(22-(z(2,i)+1/2*h*K2(2)))-1/2*((z(2,i)+1/2*h*K2(2))-(z(1,i)+1/2*h*K2(1)))];
K4=[1/4*(2-7*cos((pi/12)*t(i))-(z(1,i)+h*K3(1)))+(10*t(i)/(1+t(i)))+5/3*(22-(z(1,i)+h*K3(1)))+1/2*((z(2,i)+h*K3(2))-(z(1,i)+h*K3(1))) ; 1/5*(2-7*cos(pi/12*t(i))-(z(2,i)+h*K3(2)))+4*cos(pi/6*t(i))+13/7*(22-(z(2,i)+h*K3(2)))-1/2*((z(2,i)+h*K3(2))-(z(1,i)+h*K3(1)))];
z(:,i+1)=z(:,i)+h/6*(K1+2*K2+2*K3+K4);
end
%Temperatura exterior
M=2-7*cos(pi*t/12);
%Dibujamos las gráficas
figure (1)
hold on
plot(t,M,'c-*')
plot(t,x(1,:),'y')
plot(t,x(2,:),'m')
plot(t,y(1,:),'r')
plot(t,y(2,:),'g')
plot(t,z(1,:),'b')
plot(t,z(2,:),'k')
title('Paso h=0.001')
legend('Tª Exterior','TªA Euler','TªB Euler','TªA Euler implícito','TªB Euler implícito','TªA Runge-Kutta','TªB Runge-Kutta','Location','best')
xlabel('Tiempo en horas')
ylabel('Temperatura A y Temperatura B')
hold off
}}
====Euler, Euler implícito y Runge-Kutta de orden 4 con h=0.1 ====
En la siguiente gráfica se muestran las temperaturas interiores en cada zona calculadas con h=0.1 por los tres métodos citados, así como la temperatura del exterior. Se puede observar que con este tamaño de paso, que es relativamente pequeño las gráficas de cada temperatura interior obtenidas con los distintos métodos difieren entre ellas. Aunque la diferencia no es muy grande si es apreciable por lo que no sería un tamaño de paso demasiado conveniente. 
[[Archivo:H01.jpg|marco|centro]]
La siguiente imagen es un zoom de anterior gráfica para que se pueda apreciar con mayor detalle la diferencia de las soluciones halladas con cada método.
[[Archivo:H01zoom.jpg|marco|centro]]

==== Euler, Euler implícito y Runge-Kutta de orden 4 con h=0.01  ====
En la siguiente gráfica se muestran las temperaturas interiores en cada zona calculadas con h=0.01 por los tres métodos citados, así como la temperatura del exterior. Se puede observar que con este tamaño de paso, que ya comienza a ser más considerable las gráficas de cada temperatura interior obtenidas con los distintos métodos prácticamente no difieren entre ellas. A primera vista y sin realizar un zoom se podría decir que las curvas están superpuestas. Este tamaño ya puede resultarnos lo suficientemente aproximado.
[[Archivo:H001.jpg|marco|centro]]
En la siguiente imagen se puede observar que para apreciar la separación de las curvas el zoom que se debe realizar es mucho mayor que el anterior.
[[Archivo:H001zoom.jpg|marco|centro]]

==== Euler, Euler implícito y Runge-Kutta de orden 4 con h=0.001 ====
En la siguiente gráfica se muestran las temperaturas interiores en cada zona calculadas con h=0.001 por los tres métodos citados, así como la temperatura del exterior. Se puede observar que con este tamaño de paso las curvas se acercan aún más que en el caso anterior, por lo cual este tamaño de paso nos daría un salución aún más aproximada a la realidad.
[[Archivo:H0001.jpg|marco|centro]]
En la imagen se observa que para lograr ver que las curvas no están realmente superpuestas, el zoom que hay que realizar es enorme y por supuesto incluso mayor que el anterior.
[[Archivo:H0001zoom.jpg|marco|centro]]