MateWiki - Contribuciones del usuario [es]

Modelo Térmico de un Edificio.(Grupo 13A)

2015-03-06T23:01:00Z

P.murios: /* Demostración del enunciado */

{{ 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>

===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.

==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 [T_{D} - 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{KuT_{D}+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[T_{D}-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+KuT_{D} </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{KuT_{D}+kMo+Ho}{k + Ku}\ </math> y <math> T_{p1}= \frac{KuT_{D}+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>

===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]]

Modelo Térmico de un Edificio.(Grupo 13A)

2015-03-06T22:54:44Z

P.murios: /* Demostración del enunciado */

{{ 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]]

Modelo Térmico de un Edificio.(Grupo 13A)

2015-03-06T22:51:39Z

P.murios: /* Demostración del enunciado */

{{ 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(wt)+\frac{w}{k1}\sin(wt)}{1+\frac{w^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(wt)]-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{w}{K_{1}})}^2}cos(wt)-\frac{\frac{KB}{K_{1}^2}w}{1+{(\frac{w}{K_{1}})}^2}sin(wt) </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{w}{K_{1}})}^2}cos(wt)+\frac{\frac{k}{K_{1}^2}w}{1+{(\frac{w}{k_{1}})}^2}sin(wt)]}+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]]

Modelo Térmico de un Edificio.(Grupo 13A)

2015-03-06T22:47:56Z

P.murios:

{{ 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(wt)+\frac{w}{k1}\sin(wt)}{1+\frac{w^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(wt)]-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 ''' T_{p2}'''. La nueva expresión queda:
<math> T_{p2} (t)=-\frac{\frac{kB}{K_{1}}}{1+{(\frac{w}{K_{1}})}^2}cos(wt)-\frac{\frac{KB}{K_{1}^2}w}{1+{(\frac{w}{K_{1}})}^2}sin(wt) </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{w}{K_{1}})}^2}cos(wt)+\frac{\frac{k}{K_{1}^2}w}{1+{(\frac{w}{k_{1}})}^2}sin(wt)]}+ce^{-kt}</math>
Como podemos observar, el factor que multiplica a B1 es ''' F_{1}(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 - Bo + \frac{B}{1+\frac{w^2}{k^2}}\\ </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]]

Modelo Térmico de un Edificio.(Grupo 13A)

2015-03-06T22:47:04Z

P.murios:

{{ 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==
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(wt)+\frac{w}{k1}\sin(wt)}{1+\frac{w^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(wt)]-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 ''' T_{p2}'''. La nueva expresión queda:
<math> T_{p2} (t)=-\frac{\frac{kB}{K_{1}}}{1+{(\frac{w}{K_{1}})}^2}cos(wt)-\frac{\frac{KB}{K_{1}^2}w}{1+{(\frac{w}{K_{1}})}^2}sin(wt) </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{w}{K_{1}})}^2}cos(wt)+\frac{\frac{k}{K_{1}^2}w}{1+{(\frac{w}{k_{1}})}^2}sin(wt)]}+ce^{-kt}</math>
Como podemos observar, el factor que multiplica a B1 es ''' F_{1}(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 - Bo + \frac{B}{1+\frac{w^2}{k^2}}\\ </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]]

Modelo Térmico de un Edificio.(Grupo 13A)

2015-03-06T22:45:56Z

P.murios:

{{ 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 termostato==
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(wt)+\frac{w}{k1}\sin(wt)}{1+\frac{w^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(wt)]-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 ''' T_{p2}'''. La nueva expresión queda:
<math> T_{p2} (t)=-\frac{\frac{kB}{K_{1}}}{1+{(\frac{w}{K_{1}})}^2}cos(wt)-\frac{\frac{KB}{K_{1}^2}w}{1+{(\frac{w}{K_{1}})}^2}sin(wt) </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{w}{K_{1}})}^2}cos(wt)+\frac{\frac{k}{K_{1}^2}w}{1+{(\frac{w}{k_{1}})}^2}sin(wt)]}+ce^{-kt}</math>
Como podemos observar, el factor que multiplica a B1 es ''' F_{1}(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 - Bo + \frac{B}{1+\frac{w^2}{k^2}}\\ </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.

==Modelo con una variación senoidal con termostato==
===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]]

Estudio del movimiento de un sistema de partículas. Grupo 8-A

2014-12-04T19:29:13Z

P.murios:

{{ TrabajoED | Estudio del movimiento de un sistema de partículas. Grupo 8-A | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC14/15|2014-15]] | Laura García Pedraz, Paula Martínez Urios, Sarah Boufounas, Isabel Roselló Colom }}

==Introducción==

El presente trabajo pretende analizar el movimiento de un sistema de partículas así como varios aspectos físicos tales como la velocidad y momento angular, la energía cinética y otros que veremos a continuación.

En todo el estudio que realizamos utilizaremos el programa informático Matlab para realizar los cálculos necesarios y dibujar figuras que ilustren todos los pasos del análisis que realizamos.

Las aclaraciones de algunos pasos más simples aparecen como comentarios en Matlab por si pudiesen ayudar al lector en la comprensión del proceso seguido.

==Sistema de partículas inicial==

Comenzamos dibujando el sistema de partículas inicial. Se trata de diez partículas con coordenadas en un sistema de tres dimensiones que forman una nube de puntos en este espacio. En todo momento trabajaremos en la base ortonormal {e1,e2,e3}

{{matlab|codigo=
%Introducimos las coordenadas de los puntos
P=[0 1/5 2/5 3/5 4/5 1 6/5 7/5 8/5 9/5;0 sin(pi/4) sin(pi/2) sin(3*pi/4) sin(pi) sin(5*pi/4) sin(3*pi/2) sin(7*pi/4) sin(2*pi) sin(9*pi/4);0 pi/30 pi/15 pi/10 4*pi/30 pi/6 pi/5 7*pi/30 4*pi/15 3*pi/10];
%Los puntos son las columnas de la matriz P. P es una matriz de 3x10
figure(1)
% Dibujamos los puntos
plot3(P(1,:),P(2,:),P(3,:),'o-','MarkerFaceColor','b')
% Dibujamos los ejes de la región de trabajo
axis([-2,2,-2,2,-2,2])
}}

==Cálculo del centro de masas del sistema==

[[Archivo:Centromasas.jpg|200px|thumb|left|Sistema de partículas inicial (puntos azules) y centro de masas del sistema (punto verde).]]
A continuación calcularemos el centro de masas de la nube de puntos. Nuestro sistema de masas es discreto y está formado por un conjunto de masas puntuales, todas ellas de valor igual a 10. Para el cálculo del centro de masas utilizaremos la fórmula que conocemos de física:

∑110(rimi)/M
donde M es la masa total (de valor 100) y ri el vector posición de cada partícula i.
{{matlab|codigo=
%Con una iteración en Matlab podemos crear un programa que calcule el centro de masas.
sumR=P(:,1);
for n=2:10
sumR=sumR+P(:,n);
end
sumR
cm=sumR/10
hold on
%Dibujamos el centro de masas calculado
plot3(cm(1),cm(2),cm(3),'g.')
}}

==Rotaciones del sistema de puntos alrededor de un eje y de valor el de un ángulo determinado==

[[Archivo:TransformadoT1.jpg|200px|thumb|left|Sistema de partículas inicial (puntos azules) y transformados de los puntos según el eje ω=e3 y ángulo θ= π/16 (puntos amarillos).]]
Lo siguiente que queremos hacer es aplicar a nuestro sistema una rotación de θ= π/16 alrededor del eje ω=e3. Para ello calcularemos previamente la matriz de componentes asociada a esta rotación.

{{matlab|codigo=
%Introducimos las matrices que intervienen en el cálculo de la matriz de rotación según nuestros conocimientos de Teoría de Campos.
I=[1 0 0;0 1 0;0 0 1];
EE1=[0 0 0;0 0 0;0 0 1];
EX1=[0 -1 0;1 0 0;0 0 0];
%Calculamos la matriz de componentes asociada a la rotación
R1=cos(pi/16)*I+(1-cos(pi/16))*EE1+sin(pi/16)*EX1;
%La matriz T1, estará formada por los vectores transformados de los iniciales según el eje e3 y el ángulo pi/16 que rotan los puntos iniciales alrededor del eje
%T1 será una matriz de 10x3
T1=zeros(10,3) ;
for n=1:10
T1(n,:)=R1*P(:,n);
end
T1
% Igual que hemos hecho en el primer apartado, dibujamos los vectores tras la rotación
plot3(T1(:,1),T1(:,2),T1(:,3),'o-','MarkerFaceColor','y') % Dibujar los transformados con eje e3
}}

[[Archivo:Figura_trabajo.jpg|200px|thumb|left|Sistema inicial de partículas (azul), centro de masas (verde), simulación 1 (amarillo), simulación 2 (rojo), simulación 3 (cián), simulación 4 (magenta).]]
En vez del eje ω=e3 que hemos decidido tomar, podemos tomar otros ejes como por ejemplo el ω=e1 (simulación 2), el ω=e2 (simulación 3), o el ω=e1+e2+e3 (simulación 4). En todas estas simulaciones hemos supuesto que el ángulo girado era el mismo (θ= π/16), pero se puede cambiar con un sencillo cambio en los comandos de Matlab.
El proceso es el mismo en cada caso, pero cambian las matrices que intervienen en el cálculo de la matriz de componentes asociada a cada rotación.

{{matlab|codigo=
%Repetimos el proceso
%Introducimos las nuevas matrices
EE2=[1 0 0;0 0 0;0 0 0];
EX2=[0 0 0;0 0 -1;0 1 0];
EE3=[0 0 0;0 1 0;0 0 0];
EX3=[0 0 1;0 0 0;-1 0 0];
EE4=[1 1 1;1 1 1;1 1 1];
EX4=[0 -1 1;1 0 -1;-1 1 0];
%Calculamos la matriz de componentes asociada a cada rotación
R2=cos(pi/16)*I+(1-cos(pi/16))*EE2+sin(pi/16)*EX2;
R3=cos(pi/16)*I+(1-cos(pi/16))*EE3+sin(pi/16)*EX3;
R4=cos(pi/16)*I+(1-cos(pi/16))*EE4+sin(pi/16)*EX4;
%Calculamos los transformados en las tres simulaciones
T2=zeros(10,3) ;
for n=1:10
T2(n,:)=R2*P(:,n);
end
T2
T3=zeros(10,3) ;
for n=1:10
T3(n,:)=R3*P(:,n);
end
T3
T4=zeros(10,3) ;
for n=1:10
T4(n,:)=R4*P(:,n);
end
T4
% Dibujamos los vectores transformados por las distintas rotaciones
plot3(T2(:,1),T2(:,2),T2(:,3),'o-','MarkerFaceColor','r') % Dibujar los transformados con eje e1
plot3(T3(:,1),T3(:,2),T3(:,3),'o-','MarkerFaceColor','c') % Dibujar los transformados con eje e2
plot3(T4(:,1),T4(:,2),T4(:,3),'o-','MarkerFaceColor','m') % Dibujar los transformados con eje e1+e2+e3
}}

[[Categoría:Teoría de Campos]]
[[Categoría:TC14/15]]

Estudio del movimiento de un sistema de partículas. Grupo 8-A

2014-12-04T19:26:30Z

P.murios: /* Cálculo del centro de masas del sistema */

{{ TrabajoED | Estudio del movimiento de un sistema de partículas. Grupo 8-A | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC14/15|2014-15]] | Laura García Pedraz, Paula Martínez Urios, Sarah Boufounas, Isabel Roselló Colom }}

==Introducción==

El presente trabajo pretende analizar el movimiento de un sistema de partículas así como varios aspectos físicos tales como la velocidad y momento angular, la energía cinética y otros que veremos a continuación.

En todo el estudio que realizamos utilizaremos el programa informático Matlab para realizar los cálculos necesarios y dibujar figuras que ilustren todos los pasos del análisis que realizamos.

Las aclaraciones de algunos pasos más simples aparecen como comentarios en Matlab por si pudiesen ayudar al lector en la comprensión del proceso seguido.

==Sistema de partículas inicial==

Comenzamos dibujando el sistema de partículas inicial. Se trata de diez partículas con coordenadas en un sistema de tres dimensiones que forman una nube de puntos en este espacio. En todo momento trabajaremos en la base ortonormal {e1,e2,e3}

{{matlab|codigo=
%Introducimos las coordenadas de los puntos
P=[0 1/5 2/5 3/5 4/5 1 6/5 7/5 8/5 9/5;0 sin(pi/4) sin(pi/2) sin(3*pi/4) sin(pi) sin(5*pi/4) sin(3*pi/2) sin(7*pi/4) sin(2*pi) sin(9*pi/4);0 pi/30 pi/15 pi/10 4*pi/30 pi/6 pi/5 7*pi/30 4*pi/15 3*pi/10];
%Los puntos son las columnas de la matriz P. P es una matriz de 3x10
figure(1)
% Dibujamos los puntos
plot3(P(1,:),P(2,:),P(3,:),'o-','MarkerFaceColor','b')
% Dibujamos los ejes de la región de trabajo
axis([-2,2,-2,2,-2,2])
}}

==Cálculo del centro de masas del sistema==

[[Archivo:Centromasas.jpg|200px|thumb|left|Sistema de partículas inicial (puntos azules) y centro de masas del sistema (punto verde).]]
A continuación calcularemos el centro de masas de la nube de puntos. Nuestro sistema de masas es discreto y está formado por un conjunto de masas puntuales, todas ellas de valor igual a 10. Para el cálculo del centro de masas utilizaremos la fórmula que conocemos de física:

∑110(rimi)/M
donde M es la masa total (de valor 100) y ri el vector posición de cada partícula i.
{{matlab|codigo=
%Con una iteración en Matlab podemos crear un programa que calcule el centro de masas.
sumR=P(:,1);
for n=2:10
sumR=sumR+P(:,n);
end
sumR
cm=sumR/10
hold on
%Dibujamos el centro de masas calculado
plot3(cm(1),cm(2),cm(3),'g.')
}}

==Rotaciones del sistema de puntos alrededor de un eje y de valor el de un ángulo determinado==

[[Archivo:TransformadoT1.jpg|200px|thumb|left|Sistema de partículas inicial (puntos azules) y transformados de los puntos según el eje ω=e3 y angulo θ= π/16 (puntos amarillos).]]
Lo siguiente que queremos hacer es aplicar a nuestro sistema una rotación de pi/16 alrededor del eje ω=e3. Para ello calcularemos previamente la matriz de componentes asociada a esta rotación.

{{matlab|codigo=
%Introducimos las matrices que intervienen en el cálculo de la matriz de rotación según nuestros conocimientos de Teoría de Campos.
I=[1 0 0;0 1 0;0 0 1];
EE1=[0 0 0;0 0 0;0 0 1];
EX1=[0 -1 0;1 0 0;0 0 0];
%Calculamos la matriz de componentes asociada a la rotación
R1=cos(pi/16)*I+(1-cos(pi/16))*EE1+sin(pi/16)*EX1;
%La matriz T1, estará formada por los vectores transformados de los iniciales según el eje e3 y el ángulo pi/16 que rotan los puntos iniciales alrededor del eje
%T1 será una matriz de 10x3
T1=zeros(10,3) ;
for n=1:10
T1(n,:)=R1*P(:,n);
end
T1
% Igual que hemos hecho en el primer apartado, dibujamos los vectores tras la rotación
plot3(T1(:,1),T1(:,2),T1(:,3),'o-','MarkerFaceColor','y') % Dibujar los transformados con eje e3
}}

[[Archivo:Figura_trabajo.jpg|200px|thumb|left|Sistema inicial de partículas (azul), centro de masas (verde), simulación 1 (amarillo), simulación 2 (rojo), simulación 3 (cián), simulación 4 (magenta).]]
En vez del eje ω=e3 que hemos decidido tomar, podemos tomar otros ejes como por ejemplo el ω=e1 (simulación 2), el ω=e2 (simulación 3), o el ω=e1+e2+e3 (simulación 4). En todas estas simulaciones hemos supuesto que el ángulo girado era el mismo, pero se puede cambiar con un sencillo cambio en los comandos de Matlab.
El proceso es el mismo en cada caso, pero cambian las matrices que intervienen en el cálculo de la matriz de componentes asociada a cada rotación.

{{matlab|codigo=
%Repetimos el proceso
%Introducimos las nuevas matrices
EE2=[1 0 0;0 0 0;0 0 0];
EX2=[0 0 0;0 0 -1;0 1 0];
EE3=[0 0 0;0 1 0;0 0 0];
EX3=[0 0 1;0 0 0;-1 0 0];
EE4=[1 1 1;1 1 1;1 1 1];
EX4=[0 -1 1;1 0 -1;-1 1 0];
%Calculamos la matriz de componentes asociada a cada rotación
R2=cos(pi/16)*I+(1-cos(pi/16))*EE2+sin(pi/16)*EX2;
R3=cos(pi/16)*I+(1-cos(pi/16))*EE3+sin(pi/16)*EX3;
R4=cos(pi/16)*I+(1-cos(pi/16))*EE4+sin(pi/16)*EX4;
%Calculamos los transformados en las tres simulaciones
T2=zeros(10,3) ;
for n=1:10
T2(n,:)=R2*P(:,n);
end
T2
T3=zeros(10,3) ;
for n=1:10
T3(n,:)=R3*P(:,n);
end
T3
T4=zeros(10,3) ;
for n=1:10
T4(n,:)=R4*P(:,n);
end
T4
% Dibujamos los vectores transformados por las distintas rotaciones
plot3(T2(:,1),T2(:,2),T2(:,3),'o-','MarkerFaceColor','r') % Dibujar los transformados con eje e1
plot3(T3(:,1),T3(:,2),T3(:,3),'o-','MarkerFaceColor','c') % Dibujar los transformados con eje e2
plot3(T4(:,1),T4(:,2),T4(:,3),'o-','MarkerFaceColor','m') % Dibujar los transformados con eje e1+e2+e3
}}

[[Categoría:Teoría de Campos]]
[[Categoría:TC14/15]]

Estudio del movimiento de un sistema de partículas. Grupo 8-A

2014-12-04T19:25:44Z

P.murios: /* Cálculo del centro de masas del sistema */

{{ TrabajoED | Estudio del movimiento de un sistema de partículas. Grupo 8-A | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC14/15|2014-15]] | Laura García Pedraz, Paula Martínez Urios, Sarah Boufounas, Isabel Roselló Colom }}

==Introducción==

El presente trabajo pretende analizar el movimiento de un sistema de partículas así como varios aspectos físicos tales como la velocidad y momento angular, la energía cinética y otros que veremos a continuación.

En todo el estudio que realizamos utilizaremos el programa informático Matlab para realizar los cálculos necesarios y dibujar figuras que ilustren todos los pasos del análisis que realizamos.

Las aclaraciones de algunos pasos más simples aparecen como comentarios en Matlab por si pudiesen ayudar al lector en la comprensión del proceso seguido.

==Sistema de partículas inicial==

Comenzamos dibujando el sistema de partículas inicial. Se trata de diez partículas con coordenadas en un sistema de tres dimensiones que forman una nube de puntos en este espacio. En todo momento trabajaremos en la base ortonormal {e1,e2,e3}

{{matlab|codigo=
%Introducimos las coordenadas de los puntos
P=[0 1/5 2/5 3/5 4/5 1 6/5 7/5 8/5 9/5;0 sin(pi/4) sin(pi/2) sin(3*pi/4) sin(pi) sin(5*pi/4) sin(3*pi/2) sin(7*pi/4) sin(2*pi) sin(9*pi/4);0 pi/30 pi/15 pi/10 4*pi/30 pi/6 pi/5 7*pi/30 4*pi/15 3*pi/10];
%Los puntos son las columnas de la matriz P. P es una matriz de 3x10
figure(1)
% Dibujamos los puntos
plot3(P(1,:),P(2,:),P(3,:),'o-','MarkerFaceColor','b')
% Dibujamos los ejes de la región de trabajo
axis([-2,2,-2,2,-2,2])
}}

==Cálculo del centro de masas del sistema==

[[Archivo:Centromasas.jpg|200px|thumb|left|Sistema de partículas inicial (puntos azules) y centro de masas del sistema (punto verde).]]
A continuación calcularemos el centro de masas de la nube de puntos. Nuestro sistema de masas es discreto y está formado por un conjunto de masas puntuales, todas ellas de valor igual a 10. Para el cálculo del centro de masas utilizaremos la fórmula que conocemos de física:

∑110<(/sup>rimi)/M
donde M es la masa total (de valor 100) y ri el vector posición de cada partícula i.
{{matlab|codigo=
%Con una iteración en Matlab podemos crear un programa que calcule el centro de masas.
sumR=P(:,1);
for n=2:10
sumR=sumR+P(:,n);
end
sumR
cm=sumR/10
hold on
%Dibujamos el centro de masas calculado
plot3(cm(1),cm(2),cm(3),'g.')
}}

==Rotaciones del sistema de puntos alrededor de un eje y de valor el de un ángulo determinado==

[[Archivo:TransformadoT1.jpg|200px|thumb|left|Sistema de partículas inicial (puntos azules) y transformados de los puntos según el eje ω=e3 y angulo θ= π/16 (puntos amarillos).]]
Lo siguiente que queremos hacer es aplicar a nuestro sistema una rotación de pi/16 alrededor del eje ω=e3. Para ello calcularemos previamente la matriz de componentes asociada a esta rotación.

{{matlab|codigo=
%Introducimos las matrices que intervienen en el cálculo de la matriz de rotación según nuestros conocimientos de Teoría de Campos.
I=[1 0 0;0 1 0;0 0 1];
EE1=[0 0 0;0 0 0;0 0 1];
EX1=[0 -1 0;1 0 0;0 0 0];
%Calculamos la matriz de componentes asociada a la rotación
R1=cos(pi/16)*I+(1-cos(pi/16))*EE1+sin(pi/16)*EX1;
%La matriz T1, estará formada por los vectores transformados de los iniciales según el eje e3 y el ángulo pi/16 que rotan los puntos iniciales alrededor del eje
%T1 será una matriz de 10x3
T1=zeros(10,3) ;
for n=1:10
T1(n,:)=R1*P(:,n);
end
T1
% Igual que hemos hecho en el primer apartado, dibujamos los vectores tras la rotación
plot3(T1(:,1),T1(:,2),T1(:,3),'o-','MarkerFaceColor','y') % Dibujar los transformados con eje e3
}}

[[Archivo:Figura_trabajo.jpg|200px|thumb|left|Sistema inicial de partículas (azul), centro de masas (verde), simulación 1 (amarillo), simulación 2 (rojo), simulación 3 (cián), simulación 4 (magenta).]]
En vez del eje ω=e3 que hemos decidido tomar, podemos tomar otros ejes como por ejemplo el ω=e1 (simulación 2), el ω=e2 (simulación 3), o el ω=e1+e2+e3 (simulación 4). En todas estas simulaciones hemos supuesto que el ángulo girado era el mismo, pero se puede cambiar con un sencillo cambio en los comandos de Matlab.
El proceso es el mismo en cada caso, pero cambian las matrices que intervienen en el cálculo de la matriz de componentes asociada a cada rotación.

{{matlab|codigo=
%Repetimos el proceso
%Introducimos las nuevas matrices
EE2=[1 0 0;0 0 0;0 0 0];
EX2=[0 0 0;0 0 -1;0 1 0];
EE3=[0 0 0;0 1 0;0 0 0];
EX3=[0 0 1;0 0 0;-1 0 0];
EE4=[1 1 1;1 1 1;1 1 1];
EX4=[0 -1 1;1 0 -1;-1 1 0];
%Calculamos la matriz de componentes asociada a cada rotación
R2=cos(pi/16)*I+(1-cos(pi/16))*EE2+sin(pi/16)*EX2;
R3=cos(pi/16)*I+(1-cos(pi/16))*EE3+sin(pi/16)*EX3;
R4=cos(pi/16)*I+(1-cos(pi/16))*EE4+sin(pi/16)*EX4;
%Calculamos los transformados en las tres simulaciones
T2=zeros(10,3) ;
for n=1:10
T2(n,:)=R2*P(:,n);
end
T2
T3=zeros(10,3) ;
for n=1:10
T3(n,:)=R3*P(:,n);
end
T3
T4=zeros(10,3) ;
for n=1:10
T4(n,:)=R4*P(:,n);
end
T4
% Dibujamos los vectores transformados por las distintas rotaciones
plot3(T2(:,1),T2(:,2),T2(:,3),'o-','MarkerFaceColor','r') % Dibujar los transformados con eje e1
plot3(T3(:,1),T3(:,2),T3(:,3),'o-','MarkerFaceColor','c') % Dibujar los transformados con eje e2
plot3(T4(:,1),T4(:,2),T4(:,3),'o-','MarkerFaceColor','m') % Dibujar los transformados con eje e1+e2+e3
}}

[[Categoría:Teoría de Campos]]
[[Categoría:TC14/15]]

Estudio del movimiento de un sistema de partículas. Grupo 8-A

2014-12-04T19:25:13Z

P.murios: /* Cálculo del centro de masas del sistema */

{{ TrabajoED | Estudio del movimiento de un sistema de partículas. Grupo 8-A | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC14/15|2014-15]] | Laura García Pedraz, Paula Martínez Urios, Sarah Boufounas, Isabel Roselló Colom }}

==Introducción==

El presente trabajo pretende analizar el movimiento de un sistema de partículas así como varios aspectos físicos tales como la velocidad y momento angular, la energía cinética y otros que veremos a continuación.

En todo el estudio que realizamos utilizaremos el programa informático Matlab para realizar los cálculos necesarios y dibujar figuras que ilustren todos los pasos del análisis que realizamos.

Las aclaraciones de algunos pasos más simples aparecen como comentarios en Matlab por si pudiesen ayudar al lector en la comprensión del proceso seguido.

==Sistema de partículas inicial==

Comenzamos dibujando el sistema de partículas inicial. Se trata de diez partículas con coordenadas en un sistema de tres dimensiones que forman una nube de puntos en este espacio. En todo momento trabajaremos en la base ortonormal {e1,e2,e3}

{{matlab|codigo=
%Introducimos las coordenadas de los puntos
P=[0 1/5 2/5 3/5 4/5 1 6/5 7/5 8/5 9/5;0 sin(pi/4) sin(pi/2) sin(3*pi/4) sin(pi) sin(5*pi/4) sin(3*pi/2) sin(7*pi/4) sin(2*pi) sin(9*pi/4);0 pi/30 pi/15 pi/10 4*pi/30 pi/6 pi/5 7*pi/30 4*pi/15 3*pi/10];
%Los puntos son las columnas de la matriz P. P es una matriz de 3x10
figure(1)
% Dibujamos los puntos
plot3(P(1,:),P(2,:),P(3,:),'o-','MarkerFaceColor','b')
% Dibujamos los ejes de la región de trabajo
axis([-2,2,-2,2,-2,2])
}}

==Cálculo del centro de masas del sistema==

[[Archivo:Centromasas.jpg|200px|thumb|left|Sistema de partículas inicial (puntos azules) y centro de masas del sistema (punto verde).]]
A continuación calcularemos el centro de masas de la nube de puntos. Nuestro sistema de masas es discreto y está formado por un conjunto de masas puntuales, todas ellas de valor igual a 10. Para el cálculo del centro de masas utilizaremos la fórmula que conocemos de física:

∑110(rimi)/M
donde M es la masa total (de valor 100) y ri el vector posición de cada partícula i.
{{matlab|codigo=
%Con una iteración en Matlab podemos crear un programa que calcule el centro de masas.
sumR=P(:,1);
for n=2:10
sumR=sumR+P(:,n);
end
sumR
cm=sumR/10
hold on
%Dibujamos el centro de masas calculado
plot3(cm(1),cm(2),cm(3),'g.')
}}

==Rotaciones del sistema de puntos alrededor de un eje y de valor el de un ángulo determinado==

[[Archivo:TransformadoT1.jpg|200px|thumb|left|Sistema de partículas inicial (puntos azules) y transformados de los puntos según el eje ω=e3 y angulo θ= π/16 (puntos amarillos).]]
Lo siguiente que queremos hacer es aplicar a nuestro sistema una rotación de pi/16 alrededor del eje ω=e3. Para ello calcularemos previamente la matriz de componentes asociada a esta rotación.

{{matlab|codigo=
%Introducimos las matrices que intervienen en el cálculo de la matriz de rotación según nuestros conocimientos de Teoría de Campos.
I=[1 0 0;0 1 0;0 0 1];
EE1=[0 0 0;0 0 0;0 0 1];
EX1=[0 -1 0;1 0 0;0 0 0];
%Calculamos la matriz de componentes asociada a la rotación
R1=cos(pi/16)*I+(1-cos(pi/16))*EE1+sin(pi/16)*EX1;
%La matriz T1, estará formada por los vectores transformados de los iniciales según el eje e3 y el ángulo pi/16 que rotan los puntos iniciales alrededor del eje
%T1 será una matriz de 10x3
T1=zeros(10,3) ;
for n=1:10
T1(n,:)=R1*P(:,n);
end
T1
% Igual que hemos hecho en el primer apartado, dibujamos los vectores tras la rotación
plot3(T1(:,1),T1(:,2),T1(:,3),'o-','MarkerFaceColor','y') % Dibujar los transformados con eje e3
}}

[[Archivo:Figura_trabajo.jpg|200px|thumb|left|Sistema inicial de partículas (azul), centro de masas (verde), simulación 1 (amarillo), simulación 2 (rojo), simulación 3 (cián), simulación 4 (magenta).]]
En vez del eje ω=e3 que hemos decidido tomar, podemos tomar otros ejes como por ejemplo el ω=e1 (simulación 2), el ω=e2 (simulación 3), o el ω=e1+e2+e3 (simulación 4). En todas estas simulaciones hemos supuesto que el ángulo girado era el mismo, pero se puede cambiar con un sencillo cambio en los comandos de Matlab.
El proceso es el mismo en cada caso, pero cambian las matrices que intervienen en el cálculo de la matriz de componentes asociada a cada rotación.

{{matlab|codigo=
%Repetimos el proceso
%Introducimos las nuevas matrices
EE2=[1 0 0;0 0 0;0 0 0];
EX2=[0 0 0;0 0 -1;0 1 0];
EE3=[0 0 0;0 1 0;0 0 0];
EX3=[0 0 1;0 0 0;-1 0 0];
EE4=[1 1 1;1 1 1;1 1 1];
EX4=[0 -1 1;1 0 -1;-1 1 0];
%Calculamos la matriz de componentes asociada a cada rotación
R2=cos(pi/16)*I+(1-cos(pi/16))*EE2+sin(pi/16)*EX2;
R3=cos(pi/16)*I+(1-cos(pi/16))*EE3+sin(pi/16)*EX3;
R4=cos(pi/16)*I+(1-cos(pi/16))*EE4+sin(pi/16)*EX4;
%Calculamos los transformados en las tres simulaciones
T2=zeros(10,3) ;
for n=1:10
T2(n,:)=R2*P(:,n);
end
T2
T3=zeros(10,3) ;
for n=1:10
T3(n,:)=R3*P(:,n);
end
T3
T4=zeros(10,3) ;
for n=1:10
T4(n,:)=R4*P(:,n);
end
T4
% Dibujamos los vectores transformados por las distintas rotaciones
plot3(T2(:,1),T2(:,2),T2(:,3),'o-','MarkerFaceColor','r') % Dibujar los transformados con eje e1
plot3(T3(:,1),T3(:,2),T3(:,3),'o-','MarkerFaceColor','c') % Dibujar los transformados con eje e2
plot3(T4(:,1),T4(:,2),T4(:,3),'o-','MarkerFaceColor','m') % Dibujar los transformados con eje e1+e2+e3
}}

[[Categoría:Teoría de Campos]]
[[Categoría:TC14/15]]

Estudio del movimiento de un sistema de partículas. Grupo 8-A

2014-12-04T19:24:35Z

P.murios:

{{ TrabajoED | Estudio del movimiento de un sistema de partículas. Grupo 8-A | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC14/15|2014-15]] | Laura García Pedraz, Paula Martínez Urios, Sarah Boufounas, Isabel Roselló Colom }}

==Introducción==

El presente trabajo pretende analizar el movimiento de un sistema de partículas así como varios aspectos físicos tales como la velocidad y momento angular, la energía cinética y otros que veremos a continuación.

En todo el estudio que realizamos utilizaremos el programa informático Matlab para realizar los cálculos necesarios y dibujar figuras que ilustren todos los pasos del análisis que realizamos.

Las aclaraciones de algunos pasos más simples aparecen como comentarios en Matlab por si pudiesen ayudar al lector en la comprensión del proceso seguido.

==Sistema de partículas inicial==

Comenzamos dibujando el sistema de partículas inicial. Se trata de diez partículas con coordenadas en un sistema de tres dimensiones que forman una nube de puntos en este espacio. En todo momento trabajaremos en la base ortonormal {e1,e2,e3}

{{matlab|codigo=
%Introducimos las coordenadas de los puntos
P=[0 1/5 2/5 3/5 4/5 1 6/5 7/5 8/5 9/5;0 sin(pi/4) sin(pi/2) sin(3*pi/4) sin(pi) sin(5*pi/4) sin(3*pi/2) sin(7*pi/4) sin(2*pi) sin(9*pi/4);0 pi/30 pi/15 pi/10 4*pi/30 pi/6 pi/5 7*pi/30 4*pi/15 3*pi/10];
%Los puntos son las columnas de la matriz P. P es una matriz de 3x10
figure(1)
% Dibujamos los puntos
plot3(P(1,:),P(2,:),P(3,:),'o-','MarkerFaceColor','b')
% Dibujamos los ejes de la región de trabajo
axis([-2,2,-2,2,-2,2])
}}

==Cálculo del centro de masas del sistema==

[[Archivo:Centromasas.jpg|200px|thumb|left|Sistema de partículas inicial (puntos azules) y centro de masas del sistema (punto verde).]]
A continuación calcularemos el centro de masas de la nube de puntos. Nuestro sistema de masas es discreto y está formado por un conjunto de masas puntuales, todas ellas de valor igual a 10. Para el cálculo del centro de masas utilizaremos la fórmula que conocemos de física:

∑110rimi)/M
donde M es la masa total (de valor 100) y ri el vector posición de cada partícula i.
{{matlab|codigo=
%Con una iteración en Matlab podemos crear un programa que calcule el centro de masas.
sumR=P(:,1);
for n=2:10
sumR=sumR+P(:,n);
end
sumR
cm=sumR/10
hold on
%Dibujamos el centro de masas calculado
plot3(cm(1),cm(2),cm(3),'g.')
}}

==Rotaciones del sistema de puntos alrededor de un eje y de valor el de un ángulo determinado==

[[Archivo:TransformadoT1.jpg|200px|thumb|left|Sistema de partículas inicial (puntos azules) y transformados de los puntos según el eje ω=e3 y angulo θ= π/16 (puntos amarillos).]]
Lo siguiente que queremos hacer es aplicar a nuestro sistema una rotación de pi/16 alrededor del eje ω=e3. Para ello calcularemos previamente la matriz de componentes asociada a esta rotación.

{{matlab|codigo=
%Introducimos las matrices que intervienen en el cálculo de la matriz de rotación según nuestros conocimientos de Teoría de Campos.
I=[1 0 0;0 1 0;0 0 1];
EE1=[0 0 0;0 0 0;0 0 1];
EX1=[0 -1 0;1 0 0;0 0 0];
%Calculamos la matriz de componentes asociada a la rotación
R1=cos(pi/16)*I+(1-cos(pi/16))*EE1+sin(pi/16)*EX1;
%La matriz T1, estará formada por los vectores transformados de los iniciales según el eje e3 y el ángulo pi/16 que rotan los puntos iniciales alrededor del eje
%T1 será una matriz de 10x3
T1=zeros(10,3) ;
for n=1:10
T1(n,:)=R1*P(:,n);
end
T1
% Igual que hemos hecho en el primer apartado, dibujamos los vectores tras la rotación
plot3(T1(:,1),T1(:,2),T1(:,3),'o-','MarkerFaceColor','y') % Dibujar los transformados con eje e3
}}

[[Archivo:Figura_trabajo.jpg|200px|thumb|left|Sistema inicial de partículas (azul), centro de masas (verde), simulación 1 (amarillo), simulación 2 (rojo), simulación 3 (cián), simulación 4 (magenta).]]
En vez del eje ω=e3 que hemos decidido tomar, podemos tomar otros ejes como por ejemplo el ω=e1 (simulación 2), el ω=e2 (simulación 3), o el ω=e1+e2+e3 (simulación 4). En todas estas simulaciones hemos supuesto que el ángulo girado era el mismo, pero se puede cambiar con un sencillo cambio en los comandos de Matlab.
El proceso es el mismo en cada caso, pero cambian las matrices que intervienen en el cálculo de la matriz de componentes asociada a cada rotación.

{{matlab|codigo=
%Repetimos el proceso
%Introducimos las nuevas matrices
EE2=[1 0 0;0 0 0;0 0 0];
EX2=[0 0 0;0 0 -1;0 1 0];
EE3=[0 0 0;0 1 0;0 0 0];
EX3=[0 0 1;0 0 0;-1 0 0];
EE4=[1 1 1;1 1 1;1 1 1];
EX4=[0 -1 1;1 0 -1;-1 1 0];
%Calculamos la matriz de componentes asociada a cada rotación
R2=cos(pi/16)*I+(1-cos(pi/16))*EE2+sin(pi/16)*EX2;
R3=cos(pi/16)*I+(1-cos(pi/16))*EE3+sin(pi/16)*EX3;
R4=cos(pi/16)*I+(1-cos(pi/16))*EE4+sin(pi/16)*EX4;
%Calculamos los transformados en las tres simulaciones
T2=zeros(10,3) ;
for n=1:10
T2(n,:)=R2*P(:,n);
end
T2
T3=zeros(10,3) ;
for n=1:10
T3(n,:)=R3*P(:,n);
end
T3
T4=zeros(10,3) ;
for n=1:10
T4(n,:)=R4*P(:,n);
end
T4
% Dibujamos los vectores transformados por las distintas rotaciones
plot3(T2(:,1),T2(:,2),T2(:,3),'o-','MarkerFaceColor','r') % Dibujar los transformados con eje e1
plot3(T3(:,1),T3(:,2),T3(:,3),'o-','MarkerFaceColor','c') % Dibujar los transformados con eje e2
plot3(T4(:,1),T4(:,2),T4(:,3),'o-','MarkerFaceColor','m') % Dibujar los transformados con eje e1+e2+e3
}}

[[Categoría:Teoría de Campos]]
[[Categoría:TC14/15]]

Estudio del movimiento de un sistema de partículas. Grupo 8-A

2014-12-04T19:19:19Z

P.murios:

{{ TrabajoED | Estudio del movimiento de un sistema de partículas. Grupo 8-A | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC14/15|2014-15]] | Laura García Pedraz, Paula Martínez Urios, Sarah Boufounas, Isabel Roselló Colom }}

==Introducción==

El presente trabajo pretende analizar el movimiento de un sistema de partículas así como varios aspectos físicos tales como la velocidad y momento angular, la energía cinética y otros que veremos a continuación.

En todo el estudio que realizamos utilizaremos el programa informático Matlab para realizar los cálculos necesarios y dibujar figuras que ilustren todos los pasos del análisis que realizamos.

Las aclaraciones de algunos pasos más simples aparecen como comentarios en Matlab por si pudiesen ayudar al lector en la comprensión del proceso seguido.

==Sistema de partículas inicial==

Comenzamos dibujando el sistema de partículas inicial. Se trata de diez partículas con coordenadas en un sistema de tres dimensiones que forman una nube de puntos en este espacio. En todo momento trabajaremos en la base ortonormal {e1,e2,e3}

{{matlab|codigo=
%Introducimos las coordenadas de los puntos
P=[0 1/5 2/5 3/5 4/5 1 6/5 7/5 8/5 9/5;0 sin(pi/4) sin(pi/2) sin(3*pi/4) sin(pi) sin(5*pi/4) sin(3*pi/2) sin(7*pi/4) sin(2*pi) sin(9*pi/4);0 pi/30 pi/15 pi/10 4*pi/30 pi/6 pi/5 7*pi/30 4*pi/15 3*pi/10];
%Los puntos son las columnas de la matriz P. P es una matriz de 3x10
figure(1)
% Dibujamos los puntos
plot3(P(1,:),P(2,:),P(3,:),'o-','MarkerFaceColor','b')
% Dibujamos los ejes de la región de trabajo
axis([-2,2,-2,2,-2,2])
}}

==Cálculo del centro de masas del sistema==

[[Archivo:Centromasas.jpg|200px|thumb|left|Sistema de partículas inicial (puntos azules) y centro de masas del sistema (punto verde).]]
A continuación calcularemos el centro de masas de la nube de puntos. Nuestro sistema de masas es discreto y está formado por un conjunto de masas puntuales, todas ellas de valor igual a 10. Para el cálculo del centro de masas utilizaremos la fórmula que conocemos de física.

{{matlab|codigo=
%Con una iteración en Matlab podemos crear un programa que calcule el centro de masas.
sumR=P(:,1);
for n=2:10
sumR=sumR+P(:,n);
end
sumR
cm=sumR/10
hold on
%Dibujamos el centro de masas calculado
plot3(cm(1),cm(2),cm(3),'g.')
}}

==Rotaciones del sistema de puntos alrededor de un eje y de valor el de un ángulo determinado==

[[Archivo:TransformadoT1.jpg|200px|thumb|left|Sistema de partículas inicial (puntos azules) y transformados de los puntos según el eje ω=e3 y angulo θ= π/16 (puntos amarillos).]]
Lo siguiente que queremos hacer es aplicar a nuestro sistema una rotación de pi/16 alrededor del eje ω=e3. Para ello calcularemos previamente la matriz de componentes asociada a esta rotación.

{{matlab|codigo=
%Introducimos las matrices que intervienen en el cálculo de la matriz de rotación según nuestros conocimientos de Teoría de Campos.
I=[1 0 0;0 1 0;0 0 1];
EE1=[0 0 0;0 0 0;0 0 1];
EX1=[0 -1 0;1 0 0;0 0 0];
%Calculamos la matriz de componentes asociada a la rotación
R1=cos(pi/16)*I+(1-cos(pi/16))*EE1+sin(pi/16)*EX1;
%La matriz T1, estará formada por los vectores transformados de los iniciales según el eje e3 y el ángulo pi/16 que rotan los puntos iniciales alrededor del eje
%T1 será una matriz de 10x3
T1=zeros(10,3) ;
for n=1:10
T1(n,:)=R1*P(:,n);
end
T1
% Igual que hemos hecho en el primer apartado, dibujamos los vectores tras la rotación
plot3(T1(:,1),T1(:,2),T1(:,3),'o-','MarkerFaceColor','y') % Dibujar los transformados con eje e3
}}

[[Archivo:Figura_trabajo.jpg|200px|thumb|left|Sistema inicial de partículas (azul), centro de masas (verde), simulación 1 (amarillo), simulación 2 (rojo), simulación 3 (cián), simulación 4 (magenta).]]
En vez del eje ω=e3 que hemos decidido tomar, podemos tomar otros ejes como por ejemplo el ω=e1 (simulación 2), el ω=e2 (simulación 3), o el ω=e1+e2+e3 (simulación 4). En todas estas simulaciones hemos supuesto que el ángulo girado era el mismo, pero se puede cambiar con un sencillo cambio en los comandos de Matlab.
El proceso es el mismo en cada caso, pero cambian las matrices que intervienen en el cálculo de la matriz de componentes asociada a cada rotación.

{{matlab|codigo=
%Repetimos el proceso
%Introducimos las nuevas matrices
EE2=[1 0 0;0 0 0;0 0 0];
EX2=[0 0 0;0 0 -1;0 1 0];
EE3=[0 0 0;0 1 0;0 0 0];
EX3=[0 0 1;0 0 0;-1 0 0];
EE4=[1 1 1;1 1 1;1 1 1];
EX4=[0 -1 1;1 0 -1;-1 1 0];
%Calculamos la matriz de componentes asociada a cada rotación
R2=cos(pi/16)*I+(1-cos(pi/16))*EE2+sin(pi/16)*EX2;
R3=cos(pi/16)*I+(1-cos(pi/16))*EE3+sin(pi/16)*EX3;
R4=cos(pi/16)*I+(1-cos(pi/16))*EE4+sin(pi/16)*EX4;
%Calculamos los transformados en las tres simulaciones
T2=zeros(10,3) ;
for n=1:10
T2(n,:)=R2*P(:,n);
end
T2
T3=zeros(10,3) ;
for n=1:10
T3(n,:)=R3*P(:,n);
end
T3
T4=zeros(10,3) ;
for n=1:10
T4(n,:)=R4*P(:,n);
end
T4
% Dibujamos los vectores transformados por las distintas rotaciones
plot3(T2(:,1),T2(:,2),T2(:,3),'o-','MarkerFaceColor','r') % Dibujar los transformados con eje e1
plot3(T3(:,1),T3(:,2),T3(:,3),'o-','MarkerFaceColor','c') % Dibujar los transformados con eje e2
plot3(T4(:,1),T4(:,2),T4(:,3),'o-','MarkerFaceColor','m') % Dibujar los transformados con eje e1+e2+e3
}}

[[Categoría:Teoría de Campos]]
[[Categoría:TC14/15]]

Estudio del movimiento de un sistema de partículas. Grupo 8-A

2014-12-04T16:25:37Z

P.murios:

{{ TrabajoED | Estudio del movimiento de un sistema de partículas. Grupo 8-A | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC14/15|2014-15]] | Laura García Pedraz, Paula Martínez Urios, Sarah Boufounas, Isabel Roselló Colom }}

==Introducción==

El presente trabajo pretende analizar el movimiento de un sistema de partículas así como varios aspectos físicos tales como la velocidad y momento angular, la energía cinética y otros que veremos a continuación.

En todo el estudio que realizamos utilizaremos el programa informático Matlab para realizar los cálculos necesarios y dibujar figuras que ilustren todos los pasos del análisis que realizamos.

Las aclaraciones de algunos pasos más simples aparecen como comentarios en Matlab por si pudiesen ayudar al lector en la comprensión del proceso seguido.

==Sistema de partículas inicial==

Comenzamos dibujando el sistema de partículas inicial. Se trata de diez partículas con coordenadas en un sistema de tres dimensiones que forman una nube de puntos en este espacio.

{{matlab|codigo=
%Introducimos las coordenadas de los puntos
P=[0 1/5 2/5 3/5 4/5 1 6/5 7/5 8/5 9/5;0 sin(pi/4) sin(pi/2) sin(3*pi/4) sin(pi) sin(5*pi/4) sin(3*pi/2) sin(7*pi/4) sin(2*pi) sin(9*pi/4);0 pi/30 pi/15 pi/10 4*pi/30 pi/6 pi/5 7*pi/30 4*pi/15 3*pi/10];
%Los puntos son las columnas de la matriz P. P es una matriz de 3x10
figure(1)
% Dibujamos los puntos
plot3(P(1,:),P(2,:),P(3,:),'o-','MarkerFaceColor','b')
% Dibujamos los ejes de la región de trabajo
axis([-2,2,-2,2,-2,2])
}}

==Cálculo del centro de masas del sistema==

[[Archivo:Centromasas.jpg|200px|thumb|left|Sistema de partículas inicial (puntos azules) y centro de masas del sistema (punto verde).]]
A continuación calcularemos el centro de masas de la nube de puntos. Nuestro sistema de masas es discreto y está formado por un conjunto de masas puntuales, todas ellas de valor igual a 10. Para el cálculo del centro de masas utilizaremos la fórmula que conocemos de física.

{{matlab|codigo=
%Con una iteración en Matlab podemos crear un programa que calcule el centro de masas.
sumR=P(:,1);
for n=2:10
sumR=sumR+P(:,n);
end
sumR
cm=sumR/10
hold on
%Dibujamos el centro de masas calculado
plot3(cm(1),cm(2),cm(3),'g.')
}}

==Rotaciones del sistema de puntos alrededor de un eje y de valor el de un ángulo determinado==

[[Archivo:TransformadoT1.jpg|200px|thumb|left|Sistema de partículas inicial (puntos azules) y transformados de los puntos según el eje e3 y el ángulo pi/16 (puntos amarillos).]]
Lo siguiente que queremos hacer es aplicar a nuestro sistema una rotación de pi/16 alrededor del eje e3. Para ello calcularemos previamente la matriz de componentes asociada a esta rotación.

{{matlab|codigo=
%Introducimos las matrices que intervienen en el cálculo de la matriz de rotación según nuestros conocimientos de Teoría de Campos.
I=[1 0 0;0 1 0;0 0 1];
EE1=[0 0 0;0 0 0;0 0 1];
EX1=[0 -1 0;1 0 0;0 0 0];
%Calculamos la matriz de componentes asociada a la rotación
R1=cos(pi/16)*I+(1-cos(pi/16))*EE1+sin(pi/16)*EX1;
%La matriz T1, estará formada por los vectores transformados de los iniciales según el eje e3 y el ángulo pi/16 que rotan los puntos iniciales alrededor del eje
%T1 será una matriz de 10x3
T1=zeros(10,3) ;
for n=1:10
T1(n,:)=R1*P(:,n);
end
T1
% Igual que hemos hecho en el primer apartado, dibujamos los vectores tras la rotación
plot3(T1(:,1),T1(:,2),T1(:,3),'o-','MarkerFaceColor','y') % Dibujar los transformados con eje e3
}}

[[Archivo:Figura_trabajo.jpg|200px|thumb|left|Sistema inicial de partículas (azul), centro de masas (verde), simulación 1 (amarillo), simulación 2 (rojo), simulación 3 (cián), simulación 4 (magenta).]]
En vez del eje e3 que hemos decidido tomar, podemos tomar otros ejes como por ejemplo el e1 (simulación 2), el e2 (simulación 3), o el e1+e2+e3 (simulación 4). En todas estas simulaciones hemos supuesto que el ángulo girado era el mismo, pero se puede cambiar con un sencillo cambio en los comandos de Matlab.
El proceso es el mismo en cada caso, pero cambian las matrices que intervienen en el cálculo de la matriz de componentes asociada a cada rotación.

{{matlab|codigo=
%Repetimos el proceso
%Introducimos las nuevas matrices
EE2=[1 0 0;0 0 0;0 0 0];
EX2=[0 0 0;0 0 -1;0 1 0];
EE3=[0 0 0;0 1 0;0 0 0];
EX3=[0 0 1;0 0 0;-1 0 0];
EE4=[1 1 1;1 1 1;1 1 1];
EX4=[0 -1 1;1 0 -1;-1 1 0];
%Calculamos la matriz de componentes asociada a cada rotación
R2=cos(pi/16)*I+(1-cos(pi/16))*EE2+sin(pi/16)*EX2;
R3=cos(pi/16)*I+(1-cos(pi/16))*EE3+sin(pi/16)*EX3;
R4=cos(pi/16)*I+(1-cos(pi/16))*EE4+sin(pi/16)*EX4;
%Calculamos los transformados en las tres simulaciones
T2=zeros(10,3) ;
for n=1:10
T2(n,:)=R2*P(:,n);
end
T2
T3=zeros(10,3) ;
for n=1:10
T3(n,:)=R3*P(:,n);
end
T3
T4=zeros(10,3) ;
for n=1:10
T4(n,:)=R4*P(:,n);
end
T4
% Dibujamos los vectores transformados por las distintas rotaciones
plot3(T2(:,1),T2(:,2),T2(:,3),'o-','MarkerFaceColor','r') % Dibujar los transformados con eje e1
plot3(T3(:,1),T3(:,2),T3(:,3),'o-','MarkerFaceColor','c') % Dibujar los transformados con eje e2
plot3(T4(:,1),T4(:,2),T4(:,3),'o-','MarkerFaceColor','m') % Dibujar los transformados con eje e1+e2+e3
}}

[[Categoría:Teoría de Campos]]
[[Categoría:TC14/15]]

Estudio del movimiento de un sistema de partículas. Grupo 8-A

2014-12-04T16:23:14Z

P.murios:

{{ TrabajoED | Estudio del movimiento de un sistema de partículas. Grupo 8-A | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC14/15|2014-15]] | Laura García Pedraz, Paula Martínez Urios, Sarah Boufounas, Isabel Roselló Colom }}

==Introducción==

El presente trabajo pretende analizar el movimiento de un sistema de partículas así como varios aspectos físicos tales como la velocidad y momento angular, la energía cinética y otros que veremos a continuación.

En todo el estudio que realizamos utilizaremos el programa informático Matlab para realizar los cálculos necesarios y dibujar figuras que ilustren todos los pasos del análisis que realizamos.

Las aclaraciones de algunos pasos más simples aparecen como comentarios en Matlab por si pudiesen ayudar al lector en la comprensión del proceso seguido.

==Sistema de partículas inicial==

Comenzamos dibujando el sistema de partículas inicial. Se trata de diez partículas con coordenadas en un sistema de tres dimensiones que forman una nube de puntos en este espacio.

{{matlab|codigo=
%Introducimos las coordenadas de los puntos
P=[0 1/5 2/5 3/5 4/5 1 6/5 7/5 8/5 9/5;0 sin(pi/4) sin(pi/2) sin(3*pi/4) sin(pi) sin(5*pi/4) sin(3*pi/2) sin(7*pi/4) sin(2*pi) sin(9*pi/4);0 pi/30 pi/15 pi/10 4*pi/30 pi/6 pi/5 7*pi/30 4*pi/15 3*pi/10];
%Los puntos son las columnas de la matriz P. P es una matriz de 3x10
figure(1)
% Dibujamos los puntos
plot3(P(1,:),P(2,:),P(3,:),'o-','MarkerFaceColor','b')
% Dibujamos los ejes de la región de trabajo
axis([-2,2,-2,2,-2,2])
}}

==Cálculo del centro de masas del sistema==

[[Archivo:Centromasas.jpg|200px|thumb|left|Sistema de partículas inicial (puntos azules) y centro de masas del sistema (punto verde).]]
A continuación calcularemos el centro de masas de la nube de puntos. Nuestro sistema de masas es discreto y está formado por un conjunto de masas puntuales, todas ellas de valor igual a 10. Para el cálculo del centro de masas utilizaremos la fórmula que conocemos de física:

FÓRMULA

{{matlab|codigo=
%Con una iteración en Matlab podemos crear un programa que calcule el centro de masas.
sumR=P(:,1);
for n=2:10
sumR=sumR+P(:,n);
end
sumR
cm=sumR/10
hold on
%Dibujamos el centro de masas calculado
plot3(cm(1),cm(2),cm(3),'g.')
}}

==Rotaciones del sistema de puntos alrededor de un eje y de valor el de un ángulo determinado==

[[Archivo:TransformadoT1.jpg|200px|thumb|left|texto descriptivo]]
Lo siguiente que queremos hacer es aplicar a nuestro sistema una rotación de pi/16 alrededor del eje e3. Para ello calcularemos previamente la matriz de componentes asociada a esta rotación.

{{matlab|codigo=
%Introducimos las matrices que intervienen en el cálculo de la matriz de rotación según nuestros conocimientos de Teoría de Campos.
I=[1 0 0;0 1 0;0 0 1];
EE1=[0 0 0;0 0 0;0 0 1];
EX1=[0 -1 0;1 0 0;0 0 0];
%Calculamos la matriz de componentes asociada a la rotación
R1=cos(pi/16)*I+(1-cos(pi/16))*EE1+sin(pi/16)*EX1;
%La matriz T1, estará formada por los vectores transformados de los iniciales según el eje e3 y el ángulo pi/16 que rotan los puntos iniciales alrededor del eje
%T1 será una matriz de 10x3
T1=zeros(10,3) ;
for n=1:10
T1(n,:)=R1*P(:,n);
end
T1
% Igual que hemos hecho en el primer apartado, dibujamos los vectores tras la rotación
plot3(T1(:,1),T1(:,2),T1(:,3),'o-','MarkerFaceColor','y') % Dibujar los transformados con eje e3
}}

[[Archivo:Figura_trabajo.jpg|200px|thumb|left|Sistema inicial de partículas (azul), centro de masas (verde), simulación 1 (amarillo), simulación 2 (rojo), simulación 3 (cián), simulación 4 (magenta).]]
En vez del eje e3 que hemos decidido tomar, podemos tomar otros ejes como por ejemplo el e1 (simulación 2), el e2 (simulación 3), o el e1+e2+e3 (simulación 4). En todas estas simulaciones hemos supuesto que el ángulo girado era el mismo, pero se puede cambiar con un sencillo cambio en los comandos de Matlab.
El proceso es el mismo en cada caso, pero cambian las matrices que intervienen en el cálculo de la matriz de componentes asociada a cada rotación.

{{matlab|codigo=
%Repetimos el proceso
%Introducimos las nuevas matrices
EE2=[1 0 0;0 0 0;0 0 0];
EX2=[0 0 0;0 0 -1;0 1 0];
EE3=[0 0 0;0 1 0;0 0 0];
EX3=[0 0 1;0 0 0;-1 0 0];
EE4=[1 1 1;1 1 1;1 1 1];
EX4=[0 -1 1;1 0 -1;-1 1 0];
%Calculamos la matriz de componentes asociada a cada rotación
R2=cos(pi/16)*I+(1-cos(pi/16))*EE2+sin(pi/16)*EX2;
R3=cos(pi/16)*I+(1-cos(pi/16))*EE3+sin(pi/16)*EX3;
R4=cos(pi/16)*I+(1-cos(pi/16))*EE4+sin(pi/16)*EX4;
%Calculamos los transformados en las tres simulaciones
T2=zeros(10,3) ;
for n=1:10
T2(n,:)=R2*P(:,n);
end
T2
T3=zeros(10,3) ;
for n=1:10
T3(n,:)=R3*P(:,n);
end
T3
T4=zeros(10,3) ;
for n=1:10
T4(n,:)=R4*P(:,n);
end
T4
% Dibujamos los vectores transformados por las distintas rotaciones
plot3(T2(:,1),T2(:,2),T2(:,3),'o-','MarkerFaceColor','r') % Dibujar los transformados con eje e1
plot3(T3(:,1),T3(:,2),T3(:,3),'o-','MarkerFaceColor','c') % Dibujar los transformados con eje e2
plot3(T4(:,1),T4(:,2),T4(:,3),'o-','MarkerFaceColor','m') % Dibujar los transformados con eje e1+e2+e3
}}

[[Categoría:Teoría de Campos]]
[[Categoría:TC14/15]]

Archivo:Figura trabajo.jpg

2014-12-04T16:20:53Z

P.murios: Sistema inicial de partículas (azul), centro de masas (verde), simulación 1 (amarillo), simulación 2 (rojo), simulación 3 (cián), simulación 4 (magenta).

Sistema inicial de partículas (azul), centro de masas (verde), simulación 1 (amarillo), simulación 2 (rojo), simulación 3 (cián), simulación 4 (magenta).

Archivo:TransformadoT1.jpg

2014-12-04T16:17:29Z

P.murios: Sistema de partículas inicial y transformados de los puntos según el eje e3 y el ángulo pi/16

Sistema de partículas inicial y transformados de los puntos según el eje e3 y el ángulo pi/16

Estudio del movimiento de un sistema de partículas. Grupo 8-A

2014-12-04T16:14:49Z

P.murios:

{{ TrabajoED | Estudio del movimiento de un sistema de partículas. Grupo 8-A | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC14/15|2014-15]] | Laura García Pedraz, Paula Martínez Urios, Sarah Boufounas, Isabel Roselló Colom }}

==Introducción==

El presente trabajo pretende analizar el movimiento de un sistema de partículas así como varios aspectos físicos tales como la velocidad y momento angular, la energía cinética y otros que veremos a continuación.

En todo el estudio que realizamos utilizaremos el programa informático Matlab para realizar los cálculos necesarios y dibujar figuras que ilustren todos los pasos del análisis que realizamos.

Las aclaraciones de algunos pasos más simples aparecen como comentarios en Matlab por si pudiesen ayudar al lector en la comprensión del proceso seguido.

==Sistema de partículas inicial==

Comenzamos dibujando el sistema de partículas inicial. Se trata de diez partículas con coordenadas en un sistema de tres dimensiones que forman una nube de puntos en este espacio.

{{matlab|codigo=
%Introducimos las coordenadas de los puntos
P=[0 1/5 2/5 3/5 4/5 1 6/5 7/5 8/5 9/5;0 sin(pi/4) sin(pi/2) sin(3*pi/4) sin(pi) sin(5*pi/4) sin(3*pi/2) sin(7*pi/4) sin(2*pi) sin(9*pi/4);0 pi/30 pi/15 pi/10 4*pi/30 pi/6 pi/5 7*pi/30 4*pi/15 3*pi/10];
%Los puntos son las columnas de la matriz P. P es una matriz de 3x10
figure(1)
% Dibujamos los puntos
plot3(P(1,:),P(2,:),P(3,:),'o-','MarkerFaceColor','b')
% Dibujamos los ejes de la región de trabajo
axis([-2,2,-2,2,-2,2])
}}

==Cálculo del centro de masas del sistema==

[[Archivo:Centromasas.jpg|200px|thumb|left|Sistema de partículas inicial (puntos azules) y centro de masas del sistema (punto verde).]]
A continuación calcularemos el centro de masas de la nube de puntos. Nuestro sistema de masas es discreto y está formado por un conjunto de masas puntuales, todas ellas de valor igual a 10. Para el cálculo del centro de masas utilizaremos la fórmula que conocemos de física:

FÓRMULA

{{matlab|codigo=
%Con una iteración en Matlab podemos crear un programa que calcule el centro de masas.
sumR=P(:,1);
for n=2:10
sumR=sumR+P(:,n);
end
sumR
cm=sumR/10
hold on
%Dibujamos el centro de masas calculado
plot3(cm(1),cm(2),cm(3),'g.')
}}

==Rotaciones del sistema de puntos alrededor de un eje y de valor el de un ángulo determinado==

Lo siguiente que queremos hacer es aplicar a nuestro sistema una rotación de pi/16 alrededor del eje e3. Para ello calcularemos previamente la matriz de componentes asociada a esta rotación.

{{matlab|codigo=
%Introducimos las matrices que intervienen en el cálculo de la matriz de rotación según nuestros conocimientos de Teoría de Campos.
I=[1 0 0;0 1 0;0 0 1];
EE1=[0 0 0;0 0 0;0 0 1];
EX1=[0 -1 0;1 0 0;0 0 0];
%Calculamos la matriz de componentes asociada a la rotación
R1=cos(pi/16)*I+(1-cos(pi/16))*EE1+sin(pi/16)*EX1;
%La matriz T1, estará formada por los vectores transformados de los iniciales según el eje e3 y el ángulo pi/16 que rotan los puntos iniciales alrededor del eje
%T1 será una matriz de 10x3
T1=zeros(10,3) ;
for n=1:10
T1(n,:)=R1*P(:,n);
end
T1
% Igual que hemos hecho en el primer apartado, dibujamos los vectores tras la rotación
plot3(T1(:,1),T1(:,2),T1(:,3),'o-','MarkerFaceColor','y') % Dibujar los transformados con eje e3
}}

En vez del eje e3 que hemos decidido tomar, podemos tomar otros ejes como por ejemplo el e1 (simulación 2), el e2 (simulación 3), o el e1+e2+e3 (simulación 4). En todas estas simulaciones hemos supuesto que el ángulo girado era el mismo, pero se puede cambiar con un sencillo cambio en los comandos de Matlab.
El proceso es el mismo en cada caso, pero cambian las matrices que intervienen en el cálculo de la matriz de componentes asociada a cada rotación.

{{matlab|codigo=
%Repetimos el proceso
%Introducimos las nuevas matrices
EE2=[1 0 0;0 0 0;0 0 0];
EX2=[0 0 0;0 0 -1;0 1 0];
EE3=[0 0 0;0 1 0;0 0 0];
EX3=[0 0 1;0 0 0;-1 0 0];
EE4=[1 1 1;1 1 1;1 1 1];
EX4=[0 -1 1;1 0 -1;-1 1 0];
%Calculamos la matriz de componentes asociada a cada rotación
R2=cos(pi/16)*I+(1-cos(pi/16))*EE2+sin(pi/16)*EX2;
R3=cos(pi/16)*I+(1-cos(pi/16))*EE3+sin(pi/16)*EX3;
R4=cos(pi/16)*I+(1-cos(pi/16))*EE4+sin(pi/16)*EX4;
%Calculamos los transformados en las tres simulaciones
T2=zeros(10,3) ;
for n=1:10
T2(n,:)=R2*P(:,n);
end
T2
T3=zeros(10,3) ;
for n=1:10
T3(n,:)=R3*P(:,n);
end
T3
T4=zeros(10,3) ;
for n=1:10
T4(n,:)=R4*P(:,n);
end
T4
% Dibujamos los vectores transformados por las distintas rotaciones
plot3(T2(:,1),T2(:,2),T2(:,3),'o-','MarkerFaceColor','r') % Dibujar los transformados con eje e1
plot3(T3(:,1),T3(:,2),T3(:,3),'o-','MarkerFaceColor','c') % Dibujar los transformados con eje e2
plot3(T4(:,1),T4(:,2),T4(:,3),'o-','MarkerFaceColor','m') % Dibujar los transformados con eje e1+e2+e3
}}

[[Categoría:Teoría de Campos]]
[[Categoría:TC14/15]]

Estudio del movimiento de un sistema de partículas. Grupo 8-A

2014-12-04T16:13:20Z

P.murios:

{{ TrabajoED | Estudio del movimiento de un sistema de partículas. Grupo 8-A | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC14/15|2014-15]] | Laura García Pedraz, Paula Martínez Urios, Sarah Boufounas, Isabel Roselló Colom }}

==Introducción==

El presente trabajo pretende analizar el movimiento de un sistema de partículas así como varios aspectos físicos tales como la velocidad y momento angular, la energía cinética y otros que veremos a continuación.

En todo el estudio que realizamos utilizaremos el programa informático Matlab para realizar los cálculos necesarios y dibujar figuras que ilustren todos los pasos del análisis que realizamos.

Las aclaraciones de algunos pasos más simples aparecen como comentarios en Matlab por si pudiesen ayudar al lector en la comprensión del proceso seguido.

==Sistema de partículas inicial==

Comenzamos dibujando el sistema de partículas inicial. Se trata de diez partículas con coordenadas en un sistema de tres dimensiones que forman una nube de puntos en este espacio.

{{matlab|codigo=
%Introducimos las coordenadas de los puntos
P=[0 1/5 2/5 3/5 4/5 1 6/5 7/5 8/5 9/5;0 sin(pi/4) sin(pi/2) sin(3*pi/4) sin(pi) sin(5*pi/4) sin(3*pi/2) sin(7*pi/4) sin(2*pi) sin(9*pi/4);0 pi/30 pi/15 pi/10 4*pi/30 pi/6 pi/5 7*pi/30 4*pi/15 3*pi/10];
%Los puntos son las columnas de la matriz P. P es una matriz de 3x10
figure(1)
% Dibujamos los puntos
plot3(P(1,:),P(2,:),P(3,:),'o-','MarkerFaceColor','b')
% Dibujamos los ejes de la región de trabajo
axis([-2,2,-2,2,-2,2])
}}

==Cálculo del centro de masas del sistema==

[[Archivo:Centromasas.jpg|200px|thumb|left|Sistema de partículas inicial (puntos azules) y centro de masas del sistema (punto verde).]]
A continuación calcularemos el centro de masas de la nube de puntos. Nuestro sistema de masas es discreto y está formado por un conjunto de masas puntuales, todas ellas de valor igual a 10. Para el cálculo del centro de masas utilizaremos la fórmula que conocemos de física:

FÓRMULA

{{matlab|codigo=
%Con una iteración en Matlab podemos crear un programa que calcule el centro de masas.
sumR=P(:,1);
for n=2:10
sumR=sumR+P(:,n);
end
sumR
cm=sumR/10
hold on
%Dibujamos el centro de masas calculado
plot3(cm(1),cm(2),cm(3),'g.')
}}

==Rotaciones del sistema de puntos alrededor de un eje y de valor el de un ángulo determinado==

Lo siguiente que queremos hacer es aplicar a nuestro sistema una rotación de pi/16 alrededor del eje e3. Para ello calcularemos previamente la matriz de componentes asociada a esta rotación.

{{matlab|codigo=
%Introducimos las matrices que intervienen en el cálculo de la matriz de rotación según nuestros conocimientos de Teoría de Campos.
I=[1 0 0;0 1 0;0 0 1];
EE1=[0 0 0;0 0 0;0 0 1];
EX1=[0 -1 0;1 0 0;0 0 0];
%Calculamos la matriz de componentes asociada a la rotación
R1=cos(pi/16)*I+(1-cos(pi/16))*EE1+sin(pi/16)*EX1;
%La matriz T1, estará formada por los vectores transformados de los iniciales según el eje e3 y el ángulo pi/16 que rotan los puntos iniciales alrededor del eje
%T1 será una matriz de 10x3
T1=zeros(10,3) ;
for n=1:10
T1(n,:)=R1*P(:,n);
end
T1
% Igual que hemos hecho en el primer apartado, dibujamos los vectores tras la rotación
plot3(T1(:,1),T1(:,2),T1(:,3),'o-','MarkerFaceColor','y') % Dibujar los transformados con eje e3
}}

En vez del eje e3 que hemos decidido tomar, podemos tomar otros ejes como por ejemplo el e1 (simulación 2), el e2 (simulación 3), o el e1+e2+e3 (simulación 4). En todas estas simulaciones hemos supuesto que el ángulo girado era el mismo, pero se puede cambiar con un sencillo cambio en los comandos de Matlab.
El proceso es el mismo en cada caso, pero cambian las matrices que intervienen en el cálculo de la matriz de componentes asociada a cada rotación.

{{matlab|codigo=
%Repetimos el proceso
%Introducimos las nuevas matrices
EE2=[1 0 0;0 0 0;0 0 0];
EX2=[0 0 0;0 0 -1;0 1 0];
EE3=[0 0 0;0 1 0;0 0 0];
EX3=[0 0 1;0 0 0;-1 0 0];
EE4=[1 1 1;1 1 1;1 1 1];
EX4=[0 -1 1;1 0 -1;-1 1 0];
%Calculamos la matriz de componentes asociada a cada rotación
R2=cos(pi/16)*I+(1-cos(pi/16))*EE2+sin(pi/16)*EX2;
R3=cos(pi/16)*I+(1-cos(pi/16))*EE3+sin(pi/16)*EX3;
R4=cos(pi/16)*I+(1-cos(pi/16))*EE4+sin(pi/16)*EX4;
%Calculamos los transformados en las tres simulaciones
T2=zeros(10,3) ;
for n=1:10
T2(n,:)=R2*P(:,n);
end
T2

T3=zeros(10,3) ;
for n=1:10
T3(n,:)=R3*P(:,n);
end
T3

T4=zeros(10,3) ;
for n=1:10
T4(n,:)=R4*P(:,n);
end
T4

% Dibujamos los vectores transformados por las distintas rotaciones
plot3(T2(:,1),T2(:,2),T2(:,3),'o-','MarkerFaceColor','r') % Dibujar los transformados con eje e1
plot3(T3(:,1),T3(:,2),T3(:,3),'o-','MarkerFaceColor','c') % Dibujar los transformados con eje e2
plot3(T4(:,1),T4(:,2),T4(:,3),'o-','MarkerFaceColor','m') % Dibujar los transformados con eje e1+e2+e3
}}

[[Categoría:Teoría de Campos]]
[[Categoría:TC14/15]]

Estudio del movimiento de un sistema de partículas. Grupo 8-A

2014-12-04T16:12:06Z

P.murios:

{{ TrabajoED | Estudio del movimiento de un sistema de partículas. Grupo 8-A | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC14/15|2014-15]] | Laura García Pedraz, Paula Martínez Urios, Sarah Boufounas, Isabel Roselló Colom }}

==Introducción==

El presente trabajo pretende analizar el movimiento de un sistema de partículas así como varios aspectos físicos tales como la velocidad y momento angular, la energía cinética y otros que veremos a continuación.

En todo el estudio que realizamos utilizaremos el programa informático Matlab para realizar los cálculos necesarios y dibujar figuras que ilustren todos los pasos del análisis que realizamos.

Las aclaraciones de algunos pasos más simples aparecen como comentarios en Matlab por si pudiesen ayudar al lector en la comprensión del proceso seguido.

==Sistema de partículas inicial==

Comenzamos dibujando el sistema de partículas inicial. Se trata de diez partículas con coordenadas en un sistema de tres dimensiones que forman una nube de puntos en este espacio.

{{matlab|codigo=
%Introducimos las coordenadas de los puntos
P=[0 1/5 2/5 3/5 4/5 1 6/5 7/5 8/5 9/5;0 sin(pi/4) sin(pi/2) sin(3*pi/4) sin(pi) sin(5*pi/4) sin(3*pi/2) sin(7*pi/4) sin(2*pi) sin(9*pi/4);0 pi/30 pi/15 pi/10 4*pi/30 pi/6 pi/5 7*pi/30 4*pi/15 3*pi/10];
%Los puntos son las columnas de la matriz P. P es una matriz de 3x10
figure(1)
% Dibujamos los puntos
plot3(P(1,:),P(2,:),P(3,:),'o-','MarkerFaceColor','b')
% Dibujamos los ejes de la región de trabajo
axis([-2,2,-2,2,-2,2])
}}

==Cálculo del centro de masas del sistema==

[[Archivo:Centromasas.jpg|200px|thumb|left|texto descriptivo]]
A continuación calcularemos el centro de masas de la nube de puntos. Nuestro sistema de masas es discreto y está formado por un conjunto de masas puntuales, todas ellas de valor igual a 10. Para el cálculo del centro de masas utilizaremos la fórmula que conocemos de física:

FÓRMULA

{{matlab|codigo=
%Con una iteración en Matlab podemos crear un programa que calcule el centro de masas.
sumR=P(:,1);
for n=2:10
sumR=sumR+P(:,n);
end
sumR
cm=sumR/10
hold on
%Dibujamos el centro de masas calculado
plot3(cm(1),cm(2),cm(3),'g.')
}}

==Rotaciones del sistema de puntos alrededor de un eje y de valor el de un ángulo determinado==

Lo siguiente que queremos hacer es aplicar a nuestro sistema una rotación de pi/16 alrededor del eje e3. Para ello calcularemos previamente la matriz de componentes asociada a esta rotación.

{{matlab|codigo=
%Introducimos las matrices que intervienen en el cálculo de la matriz de rotación según nuestros conocimientos de Teoría de Campos.
I=[1 0 0;0 1 0;0 0 1];
EE1=[0 0 0;0 0 0;0 0 1];
EX1=[0 -1 0;1 0 0;0 0 0];
%Calculamos la matriz de componentes asociada a la rotación
R1=cos(pi/16)*I+(1-cos(pi/16))*EE1+sin(pi/16)*EX1;
%La matriz T1, estará formada por los vectores transformados de los iniciales según el eje e3 y el ángulo pi/16 que rotan los puntos iniciales alrededor del eje
%T1 será una matriz de 10x3
T1=zeros(10,3) ;
for n=1:10
T1(n,:)=R1*P(:,n);
end
T1
% Igual que hemos hecho en el primer apartado, dibujamos los vectores tras la rotación
plot3(T1(:,1),T1(:,2),T1(:,3),'o-','MarkerFaceColor','y') % Dibujar los transformados con eje e3
}}

En vez del eje e3 que hemos decidido tomar, podemos tomar otros ejes como por ejemplo el e1 (simulación 2), el e2 (simulación 3), o el e1+e2+e3 (simulación 4). En todas estas simulaciones hemos supuesto que el ángulo girado era el mismo, pero se puede cambiar con un sencillo cambio en los comandos de Matlab.
El proceso es el mismo en cada caso, pero cambian las matrices que intervienen en el cálculo de la matriz de componentes asociada a cada rotación.

{{matlab|codigo=
%Repetimos el proceso
%Introducimos las nuevas matrices
EE2=[1 0 0;0 0 0;0 0 0];
EX2=[0 0 0;0 0 -1;0 1 0];
EE3=[0 0 0;0 1 0;0 0 0];
EX3=[0 0 1;0 0 0;-1 0 0];
EE4=[1 1 1;1 1 1;1 1 1];
EX4=[0 -1 1;1 0 -1;-1 1 0];
%Calculamos la matriz de componentes asociada a cada rotación
R2=cos(pi/16)*I+(1-cos(pi/16))*EE2+sin(pi/16)*EX2;
R3=cos(pi/16)*I+(1-cos(pi/16))*EE3+sin(pi/16)*EX3;
R4=cos(pi/16)*I+(1-cos(pi/16))*EE4+sin(pi/16)*EX4;
%Calculamos los transformados en las tres simulaciones
T2=zeros(10,3) ;
for n=1:10
T2(n,:)=R2*P(:,n);
end
T2

T3=zeros(10,3) ;
for n=1:10
T3(n,:)=R3*P(:,n);
end
T3

T4=zeros(10,3) ;
for n=1:10
T4(n,:)=R4*P(:,n);
end
T4

% Dibujamos los vectores transformados por las distintas rotaciones
plot3(T2(:,1),T2(:,2),T2(:,3),'o-','MarkerFaceColor','r') % Dibujar los transformados con eje e1
plot3(T3(:,1),T3(:,2),T3(:,3),'o-','MarkerFaceColor','c') % Dibujar los transformados con eje e2
plot3(T4(:,1),T4(:,2),T4(:,3),'o-','MarkerFaceColor','m') % Dibujar los transformados con eje e1+e2+e3
}}

[[Categoría:Teoría de Campos]]
[[Categoría:TC14/15]]

Archivo:Centromasas.jpg

2014-12-04T16:07:58Z

P.murios: Sistema de partículas inicial (puntos azules) y centro de masas del sistema (punto verde).

Sistema de partículas inicial (puntos azules) y centro de masas del sistema (punto verde).

Estudio del movimiento de un sistema de partículas. Grupo 8-A

2014-12-04T15:50:57Z

P.murios:

{{ TrabajoED | Estudio del movimiento de un sistema de partículas. Grupo 8-A | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC14/15|2014-15]] | Laura García Pedraz, Paula Martínez Urios, Sarah Boufounas, Isabel Roselló Colom }}

==Introducción==

El presente trabajo pretende analizar el movimiento de un sistema de partículas así como varios aspectos físicos tales como la velocidad y momento angular, la energía cinética y otros que veremos a continuación.

En todo el estudio que realizamos utilizaremos el programa informático Matlab para realizar los cálculos necesarios y dibujar figuras que ilustren todos los pasos del análisis que realizamos.

Las aclaraciones de algunos pasos más simples aparecen como comentarios en Matlab por si pudiesen ayudar al lector en la comprensión del proceso seguido.

==Sistema de partículas inicial==

Comenzamos dibujando el sistema de partículas inicial. Se trata de diez partículas con coordenadas en un sistema de tres dimensiones que forman una nube de puntos en este espacio.

{{matlab|codigo=
%Introducimos las coordenadas de los puntos
P=[0 1/5 2/5 3/5 4/5 1 6/5 7/5 8/5 9/5;0 sin(pi/4) sin(pi/2) sin(3*pi/4) sin(pi) sin(5*pi/4) sin(3*pi/2) sin(7*pi/4) sin(2*pi) sin(9*pi/4);0 pi/30 pi/15 pi/10 4*pi/30 pi/6 pi/5 7*pi/30 4*pi/15 3*pi/10];
%Los puntos son las columnas de la matriz P. P es una matriz de 3x10
figure(1)
% Dibujamos los puntos
plot3(P(1,:),P(2,:),P(3,:),'o-','MarkerFaceColor','b')
% Dibujamos los ejes de la región de trabajo
axis([-2,2,-2,2,-2,2])
}}

==Cálculo del centro de masas del sistema==

A continuación calcularemos el centro de masas de la nube de puntos. Nuestro sistema de masas es discreto y está formado por un conjunto de masas puntuales, todas ellas de valor igual a 10. Para el cálculo del centro de masas utilizaremos la fórmula que conocemos de física:

FÓRMULA

{{matlab|codigo=
%Con una iteración en Matlab podemos crear un programa que calcule el centro de masas.
sumR=P(:,1);
for n=2:10
sumR=sumR+P(:,n);
end
sumR
cm=sumR/10
hold on
%Dibujamos el centro de masas calculado
plot3(cm(1),cm(2),cm(3),'g.')
}}

==Rotaciones del sistema de puntos alrededor de un eje y de valor el de un ángulo determinado==

Lo siguiente que queremos hacer es aplicar a nuestro sistema una rotación de pi/16 alrededor del eje e3. Para ello calcularemos previamente la matriz de componentes asociada a esta rotación.

{{matlab|codigo=
%Introducimos las matrices que intervienen en el cálculo de la matriz de rotación según nuestros conocimientos de Teoría de Campos.
I=[1 0 0;0 1 0;0 0 1];
EE1=[0 0 0;0 0 0;0 0 1];
EX1=[0 -1 0;1 0 0;0 0 0];
%Calculamos la matriz de componentes asociada a la rotación
R1=cos(pi/16)*I+(1-cos(pi/16))*EE1+sin(pi/16)*EX1;
%La matriz T1, estará formada por los vectores transformados de los iniciales según el eje e3 y el ángulo pi/16 que rotan los puntos iniciales alrededor del eje
%T1 será una matriz de 10x3
T1=zeros(10,3) ;
for n=1:10
T1(n,:)=R1*P(:,n);
end
T1
% Igual que hemos hecho en el primer apartado, dibujamos los vectores tras la rotación
plot3(T1(:,1),T1(:,2),T1(:,3),'o-','MarkerFaceColor','y') % Dibujar los transformados con eje e3
}}

En vez del eje e3 que hemos decidido tomar, podemos tomar otros ejes como por ejemplo el e1 (simulación 2), el e2 (simulación 3), o el e1+e2+e3 (simulación 4). En todas estas simulaciones hemos supuesto que el ángulo girado era el mismo, pero se puede cambiar con un sencillo cambio en los comandos de Matlab.
El proceso es el mismo en cada caso, pero cambian las matrices que intervienen en el cálculo de la matriz de componentes asociada a cada rotación.

{{matlab|codigo=
%Repetimos el proceso
%Introducimos las nuevas matrices
EE2=[1 0 0;0 0 0;0 0 0];
EX2=[0 0 0;0 0 -1;0 1 0];
EE3=[0 0 0;0 1 0;0 0 0];
EX3=[0 0 1;0 0 0;-1 0 0];
EE4=[1 1 1;1 1 1;1 1 1];
EX4=[0 -1 1;1 0 -1;-1 1 0];
%Calculamos la matriz de componentes asociada a cada rotación
R2=cos(pi/16)*I+(1-cos(pi/16))*EE2+sin(pi/16)*EX2;
R3=cos(pi/16)*I+(1-cos(pi/16))*EE3+sin(pi/16)*EX3;
R4=cos(pi/16)*I+(1-cos(pi/16))*EE4+sin(pi/16)*EX4;
%Calculamos los transformados en las tres simulaciones
T2=zeros(10,3) ;
for n=1:10
T2(n,:)=R2*P(:,n);
end
T2

T3=zeros(10,3) ;
for n=1:10
T3(n,:)=R3*P(:,n);
end
T3

T4=zeros(10,3) ;
for n=1:10
T4(n,:)=R4*P(:,n);
end
T4

% Dibujamos los vectores transformados por las distintas rotaciones
plot3(T2(:,1),T2(:,2),T2(:,3),'o-','MarkerFaceColor','r') % Dibujar los transformados con eje e1
plot3(T3(:,1),T3(:,2),T3(:,3),'o-','MarkerFaceColor','c') % Dibujar los transformados con eje e2
plot3(T4(:,1),T4(:,2),T4(:,3),'o-','MarkerFaceColor','m') % Dibujar los transformados con eje e1+e2+e3
}}

[[Categoría:Teoría de Campos]]
[[Categoría:TC14/15]]

Estudio del movimiento de un sistema de partículas. Grupo 8-A

2014-12-04T15:23:46Z

P.murios:

Estudio del movimiento de un sistema de partículas. Grupo 8-A

2014-12-04T15:22:14Z

P.murios:

Estudio del movimiento de un sistema de partículas. Grupo 8-A

2014-12-04T15:20:37Z

P.murios:

Estudio del movimiento de un sistema de partículas. Grupo 8-A

2014-12-04T15:08:32Z

P.murios:

Estudio del movimiento de un sistema de partículas. Grupo 8-A

2014-12-04T14:49:36Z

P.murios:

Estudio del movimiento de un sistema de partículas. Grupo 8-A

2014-12-04T14:46:16Z

P.murios:

<pre><nowiki>{{ TrabajoED | Estudio del movimiento de un sistema de partículas. Grupo 8-A | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC14/15|2013-14]] | Laura García Pedraz, Paula Martínez Urios, Sarah Boufounas, Isabel Roselló Colom }}</nowiki></pre>

1. Sistema de partículas inicial.

<pre><nowiki>[[Categoría:Teoría de Campos]]</nowiki></pre>
<pre><nowiki>[[Categoría:TC14/15]]</nowiki></pre>

Estudio del movimiento de un sistema de partículas. Grupo 8-A

2014-12-04T14:43:21Z

P.murios: Página creada con «<pre><nowiki>{{ TrabajoED | Estudio del movimiento de un sistema de partículas. Grupo 8-A | Teoría de Campos|:Categoría:TC14/15|2013-...»

<pre><nowiki>{{ TrabajoED | Estudio del movimiento de un sistema de partículas. Grupo 8-A | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC14/15|2013-14]] | Nuestros nombres }}</nowiki></pre>
1. Sistema de partículas inicial.

<pre><nowiki>[[Categoría:Teoría de Campos]]</nowiki></pre>
<pre><nowiki>[[Categoría:TC14/15]]</nowiki></pre>

Visualización de campos escalares y vectoriales en elasticidad grupo 14A

2013-12-10T18:39:46Z

P.murios:

[[Archivo:Gsubtheta.jpg|miniaturadeimagen|default]]
{{beta}}

'''Visualización de campos escalares y vectoriales en elasticidad grupo 14A'''

=='''1. Mallado del sólido'''==

{{matlab|codigo=

h=0.1; %Establecemos paso de muestra

u=1:h:2; %Valor de u según el muestreo del intervalo

v=0:h:2*pi+h; %Valor de v según el muestreo del intervalo

[uu,vv]=meshgrid(u,v); %Matrices de la placa en coordenadas polares

figure(1) %Abrimos la ventana donde dibujaremos el mallado

xx=uu.*cos(vv); %Parametrización

yy=uu.*sin(vv); %Parametrización

mesh(xx,yy,0*xx) %Dibujamos el mallado

axis([-3,3,-3,3]) %Establecemos los límites de los ejes x e y en la representación de la placa

view(2) %visualización

}}
[[Archivo:Malladog14.jpg|miniaturadeimagen|Mallado del Sólido]]

=='''2. Campo de temperaturas'''==

{{matlab|codigo=

h=0.1; %Establecemos paso de muestra

u=1:h:2; %Valor de u (ro) según el muestreo del intervalo

v=0:h:2*pi+h; %Valor de v (theta) según el muestreo del interval

[uu,vv]=meshgrid(u,v); %Matrices de la placa en coordenadas polares

figure(1) %Abrimos la ventana donde dibujaremos el mallado

xx=uu.*cos(vv); %Parametrización

yy=uu.*sin(vv); %Parametrización

f=-log(uu+0.1); %Campo escalar

surf(xx,yy,f) %Dibujamos la superficie

axis([-3,3,-3,3]) %Establecemos los límites de los ejes x e y en la representación de la placa

colorbar %Hacemos aparecer en la figura la leyenda de los colores

view(2) %Visualización
}}
[[Archivo:Temperaturasg14.jpg|miniaturadeimagen|Campo de temperaturas]]

=='''5. Desplazamiento del sólido'''==

La posición de nuestra placa circular sólida viene dada en cada punto por la ecuación:

r(x,y,t)= r0(x,y,t) + <math> \vec u</math>. Siendo <math> \vec u</math>. nuestro campo <math> \vec u(\rho,\theta)=\frac{\sin(\pi \theta/2)}{30\rho}\vec g_{\rho}.
</math>

{{matlab|codigo=

h=0.1; % establecemos un paso de muestra

u=1:h:2; % valor u según el paso de muestreo en el intervalo [1,2]

v=0:h:2*pi+h; % valor v según el muestreo en el intervalo [0,2*pi]

[rho,theta]=meshgrid(u,v); % matrices en componentes polares de la placa sobre la que se trabaja

xx=rho.*cos(theta); % parametrizamos

yy=rho.*sin(theta);

fx=(sin((pi*theta)/2))./(30*rho); % componente x(rho,theta) de mi campo vectorial

fy=0*yy; % componente y(rho,theta) del campo vectorial

quiver(xx,yy,fx,fy) % representamos el campo vectorial sobre la placa

axis([-3,3,-3,3]) % establecemos los límites de los ejes x e y en la representación de la placa

view % visualización

}}

La figura muestra la visualización tras aplicar el campo a los puntos. En ella se aprecia la dirección que cada uno de los puntos de nuestra placa sufre frente a su posición inicial.
Una vez aplicado el desplazamiento podremos comparar posiciones iniciales y finales.
[[Archivo:Desplazamientoanillo.jpg|miniaturadeimagen|derecha|La figura muestra el desplazamiento sobre los puntos]]

{{matlab|codigo=

h=0.1 % paso de muestreo

u=1:h:2; % valor u según el muestreo en el intervalo [1,2]

v=0:h:2*pi+h; % valor v según el muestreo en el intervalo [0,2*pi]

[rho,theta]=meshgrid(u,v); % matrices de las componentes polares de la placa

xx=rho.*cos(theta); % parametrización

yy=rho.*sin(theta);

subplot(1,2,1)

mesh(xx,yy,0*yy)

axis([-4,4,-4,4])

subplot(1,2,2)

fx=(sin((pi*theta)/2))./(30*rho);

fy=0*yy;

mesh(xx+fx,yy+fy,yy*0)

axis([-4,4,-4,4])

}}

Como puede apreciarse, la variación es mínima. Ésto es debido a la ecuación del campo aplicada. En el punto de mayor deformación es a penas 1/30 para el caso del círculo interior de radio 1, y 1/60 en el caso del exterior por ser de radio 2. (Siendo por tanto evidente que el anillo se deforma más cuanto más cerca nos encontremos de su centro).

Sólo con una vista en detalle podremos observar que el desplazamiento sufrido rompe el anillo ya que desplaza la posición de la variable theta=0º, frente a theta= 360º.
[[Archivo:Comparacionanillos.jpg|miniaturadeimagen|derecha|La figura muestra la comparación tras sufrir el desplazamiento]]

[[Archivo:Detallecomp.jpg|miniaturadeimagen|derecha|Detalle del desplazamiento sufrido por los puntos]]

=='''6. Divergencia del Campo'''==

La divergencia es una medida del cambio de volumen local debido al desplazamiento:

<math>\bigtriangledown{}\cdot{}\overrightarrow{u}
=\frac{1}{\rho}[(\frac{\partial}{\partial\rho})(\rho\cdot{}\frac{\sin(\frac{\pi\theta}{2})}{30\rho})] =0</math>

El volumen de control se comprime, expande o queda igual. Si <math>\bigtriangledown{}\cdot{}\overrightarrow{u}
=0</math> para todo punto del dominio el campo se llama incompresible.

Al ser la divergencia 0 , todos los puntos del solido tienen la misma y el cambio de volumen experimentado por el sólido es nulo, ya que experimenta cambios transversales que no alteran su volumen.

=='''7. Rotacional del Campo'''==

El rotacional es un operador vectorial que muestra la tendencia de un campo vectorial a inducir rotación alrededor de un punto.

En nuestro caso el rotacional es el siguiente:

<math>\overrightarrow{\bigtriangledown{}}\times{}\overrightarrow{u}
=\frac{1}{\rho}\begin{vmatrix} \overrightarrow{g_\rho} & \overrightarrow{g_\theta} & \overrightarrow{g_z} \\ \frac{\partial}{\partial\rho} & \frac{\partial}{\partial\theta} & \frac{\partial}{\partial z} \\ \frac{sin\frac{\pi\theta}{2}}{30\rho} & 0 & 0 \end{vmatrix}=-(\frac{\pi}{60\rho^2})(cos\frac{\pi\theta}{2})\overrightarrow{g_z}</math>

Y el módulo del rotacional es:

<math>|\overrightarrow{\bigtriangledown{}}\times{}\overrightarrow{u}|
=(\frac{\pi}{60\rho^2})(cos\frac{\pi\theta}{2})</math>

{{matlab|codigo=
h=0.1; % paso de muestreo

u=1:h:2; % valor u según el muestreo en el intervalo [1,2]

v=0:h:2*pi+h; % valor v según el muestreo en el intervalo [0,2*pi]

[rho,theta]=meshgrid(u,v); % matrices de las componentes polares de la placa

xx=rho.*cos(theta); % parametrización

yy=rho.*sin(theta);

f=(pi./((rho.^2)*60)).*(cos((pi.*theta)./2)); % campo escalar del rotacional

subplot(1,2,1)

surf(xx,yy,f) % representación del campo escalar sobre la placa

axis([-3,3,-3,3]) % límites de los ejes x e y en la representación de la placa

view % visualización del rotacional sobre la placa

subplot(1,2,2)

surf(xx,yy,f) % representación del campo escalar sobre la placa

axis([-3,3,-3,3]) % límites de los ejes x e y en la representación de la placa

view % visualización del rotacional sobre la placa
}}
[[Archivo:Rotacional plano y 3D G14.jpeg|miniaturadeimagen|Vista en planta y en perspectiva del rotacional]]

=='''8.Tensiones'''==

En un medio elástico lineal, isótropo y homogéneo los desplazamientos permiten escribir el tensor de tensiones
σij a través de la fórmula:
σij=λ<math>\bigtriangledown{}\cdot{}\overrightarrow{u}</math>+2μεij
Donde εij es la parte simétrica del vector <math>\bar{u}</math>. Tomamos λ=μ=1. Y los vectores <math>\bar{i}</math>= <math>\bar{gρ}</math> y <math>\bar{j}</math>=<math>\bar{gθ/ρ}</math>.
En este apartado se nos pide dibujar las tensiones normales en la dirección que marca el eje <math>\bar{i}</math>, es decir <math>\bar{<math>gρ</math>}</math>·σij·<math>\bar{gρ}</math>= .
y las tensiones normales en la dirección que marca el eje <math>\bar{j}</math>, es decir <math>\bar{gθ/ρ}</math>·σij·<math>\bar{gθ/ρ}</math>= ((sin(pi*vv./2))./(30*(xx.^2+yy.^2))).
.
{{matlab|codigo=

% Tensiones normales en la dirección que marca el eje <math>\bar{i}</math>.
h=0.1; %Establecemos paso de muestra

u=1:h:2; %Valor de u según el muestreo del intervalo

v=0:h:2*pi+h; %Valor de v según el muestreo del intervalo

[uu,vv]=meshgrid(u,v); %Matrices de la placa en coordenadas polares

figure(1) %Abrimos la ventana donde dibujaremos el mallado

xx=uu.*cos(vv); %Parametrización

yy=uu.*sin(vv); %Parametrización

f=((-2.*sin(pi*vv./2))./(30*(xx.^2+yy.^2))); % Campo gρ

surf(xx,yy,f) % Dibujamos la malla

axis([-2,2,-2,2]) % seleccionamos región

axis equal

colorbar

title('Tensión normal en la dirección gρ')

view(1)

subplot(1,2,2)

f=((sin(pi*vv./2))./(30*(xx.^2+yy.^2))); % Campo gθ

surf(xx,yy,f) % Dibujamos la malla

axis([-2,2,-2,2]) % Seleccionamos la región

axis equal

colorbar

title('Tensión normal en la dirección gθ')

v
iew(2)
}}
[[Archivo:Gsubtheta.jpg|miniaturadeimagen|default]]
[[Archivo:Tensionesnormalesagsubrho.jpg|miniaturadeimagen|default]]

[[Categoría:Teoría de Campos]]