MateWiki - Contribuciones del usuario [es]

Archivo:Placa10.png

2014-05-19T21:55:03Z

Alfonso:

Archivo:Placa02.png

2014-05-19T21:54:43Z

Alfonso:

Archivo:Placa01.png

2014-05-19T21:54:25Z

Alfonso:

Archivo:Placa00.png

2014-05-19T21:54:05Z

Alfonso:

Placa en forma de Anillo

2014-05-19T21:53:43Z

Alfonso:

{{ TrabajoED | Transmisión de calor en placa con forma de anillo. Grupo 11-B | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED13/14|Curso 2013-14]] | Javier Chamizo, Jacobo Campos, Miguel García, Javier Mellado, Alfonso Ascanio }}

Transmisión de calor en una placa con forma de anillo. Es bastante común modelizar el comportamiento del calor en superficies mediante las ecuaciones diferenciales, y además dada la particularidad de la forma de nuestro caso, utilizando coordenadas polares.

==Introducción==

[[Archivo:placadescripcion.png|marco|derecha|Placa en anillo con p ∈ [1, 6]]]

Tenemos la placa que se ve en la imagen superior con las siguientes condiciones. Colocamos en las fronteras interior ρ = 1 y exterior ρ = 6, que sabemos que mantienen una temperatura constante de 0 grados y 10 grados respectivamente como si estuviesen en contacto con un depósito. En nuestro caso vamos a trabajar en coordenadas polares dada la forma de la placa que se presta a ello. Suponemos también que la temperatura u de la placa solo depende de la coordenada radial y el tiempo, es decir u = u(ρ, t). La placa sabemos que satisface la ecuación del calor que se muestra a continuación:

<math> u_t-\Delta u=0</math>

La placa a su vez dispone de las siguientes condiciones iniciales de temperatura u(p,0):

<math>\left\{\begin{matrix}100·(p-1) → p ∈ [1, 2] \\100 → p ∈ [2, 5]\\90·(6-p)+10 → p ∈ [5, 6]\end{matrix}\right.</math>

Además tenemos condiciones tipo Dirichlet en los extremos ya que se mantiene constante la temperatura a lo largo del tiempo en estos al tener objetos que mantienen la temperatura constante en contacto con los mismos:

<math>\left\{\begin{matrix}u(1,t)=0 \\u(6,t)=10\end{matrix}\right.</math>

Debido a que el ángulo teta y la altura Z se mantienen constantes al ser una figura plana, el Laplaciano en coordenadas cilíndricas es el siguiente:

<math>\Delta u=\frac{1}{p} \frac{\partial }{\partial p} (p \frac{\partial u }{\partial p})+\frac{1 }{p^2} \frac{\partial^2 u }{\partial θ^2}+\frac{\partial^2 u }{\partial z^2}= \frac{\partial^2 u }{\partial p^2}+\frac{1 }{p}\frac{\partial u}{\partial p} </math>

Al aplicar el Laplaciano, la ecuación del calor que se va a resolver queda de la forma:
<math> u_t-u_{pp}-\frac{1}{p} u_p=0 </math>
Quedando así similar al sistema que obtendríamos con una barra con longitud entre 1 y 6.

==Resolución del Sistema por diferencias finitas==
Aquí se va a resolver el problema planteado por el método de diferencias finitas con paso de discretización ∆p = 0.1; y para el tiempo usaremos en cada caso los métodos del trapecio, Euler modificado, explícito e implícito. Todos ellos tomando ∆t = ∆p/4 en t ∈ [0, 10]
===Método del trapecio===
{{matlab|codigo=
%Datos del problema
clear all
p0=1;
pf=6;
T0=0;
Tf=10;
%datos discretización
N=50;
h=(pf-p0)/N;
p=p0:h:pf;
pint=p0+h:h:pf-h;
%condiciones iniciales
for i=1:length(pint)
if pint(i)<2
u0(i)=100*(pint(i)-1);
elseif pint(i)<5
u0(i)=100;
else
u0(i)=90*(6-pint(i))+10;
end
end
%Datos discretizacion en espacio
dt=h/4;%paso en tiempo suele ser más pequeño que el paso en espacio
t=T0:dt:Tf;
%diferencias finitas
%Matriz K Aprox. -u_pp
K=diag(2*ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
%K(N-1,N-2)=2;
K=K/h^2;
%Matriz L Aprox. -u_p
L=diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
L=(L*(1/(2*h)));
L=diag(1./pint)*L;
%Matriz F
F=zeros(N-1,1);
F(N-1,1)=-(10/(h^2))+10/(2*h);
%Dato inicial
W0=[u0];
%Metodo del trapecio
WW=W0';
%definimos la matriz sol con la u para pintarla
sol=zeros(length(t),N+1);
sol(1,:)=[0 W0 10];
%iteraciones desde W^j-->W^j+1
for j=1:length(t)-1
WW=(eye(N-1)+dt/2*(K+L))\((eye(N-1)-dt/2*(K+L))*WW-dt*F);
sol(j+1,:)=[0,WW',10];%guardamos solucion
end
[pp,tt]=meshgrid(p,t);
mesh(pp,tt,sol)
}}
[[Archivo:trapeciobueno.png|marco|centro|Distribución de la temperatura en el disco]]
Ahora particularizamos la solución en ρ = 3 y observamos el comportamiento de la temperatura en dichos puntos. La gráfica siguiente muestra una representación 2D temperatura/tiempo
{{matlab|codigo=
%Datos del problema
clear all
p0=1;
pf=6;
T0=0;
Tf=10;
%datos discretización
N=50;
h=(pf-p0)/N;
p=p0:h:pf;
pint=p0+h:h:pf-h;
%condiciones iniciales
for i=1:length(pint)
if pint(i)<2
u0(i)=100*(pint(i)-1);
elseif pint(i)<5
u0(i)=100;
else
u0(i)=90*(6-pint(i))+10;
end
end
%Datos discretizacion en espacio
dt=h/4;%paso en tiempo suele ser más pequeño que el paso en espacio
t=T0:dt:Tf;
%diferencias finitas
%Matriz K Aprox. -u_pp
K=diag(2*ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
%K(N-1,N-2)=2;
K=K/h^2;
%Matriz L Aprox. -u_p
L=diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
L=(L*(1/(2*h)));
L=diag(1./pint)*L;
%Matriz F
F=zeros(N-1,1);
F(N-1,1)=-(10/(h^2))+10/(2*h);
%Dato inicial
W0=[u0];
%Metodo del trapecio
WW=W0';
%definimos la matriz sol con la u para pintarla
sol=zeros(length(t),N+1);
sol(1,:)=[0 W0 10];
%iteraciones desde W^j-->W^j+1
for j=1:length(t)-1
WW=(eye(N-1)+dt/2*(K+L))\((eye(N-1)-dt/2*(K+L))*WW-dt*F);
sol(j+1,:)=[0,WW',10];%guardamos solucion
end
[pp,tt]=meshgrid(p,t);
J=ones(401,51);
plot3(3*J,tt,sol)
}}

[[Archivo:p3.png|marco|centro|Distribución de la temperatura para p=3]]

===Método de Euler explícito===
{{matlab|codigo=
%Datos del problema
clear all
p0=1;
pf=6;
T0=0;
Tf=10;
%datos discretización
N=75;
h=(pf-p0)/N;
p=p0:h:pf;
pint=p0+h:h:pf-h;
%condiciones iniciales
for i=1:length(pint)
if pint(i)<2
u0(i)=100*(pint(i)-1);
elseif pint(i)<5
u0(i)=100;
else
u0(i)=90*(6-pint(i))+10;
end
end
%Datos discretizacion en espacio
dt=h^2/4;%paso en tiempo suele ser más pequeño que el paso en espacio
t=T0:dt:Tf;
%diferencias finitas
%Matriz K Aprox. -u_pp
K=diag(2*ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
%K(N-1,N-2)=2;
K=K/h^2;
%Matriz L Aprox. -u_p
L=diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
L=(L*(1/(2*h)));
L=diag(1./pint)*L;
%Matriz F
F=zeros(N-1,1);
F(N-1,1)=-(10/(h^2))+10/(2*h);
%Dato inicial
W0=[u0];
WW=W0';
%definimos la matriz sol con la u para pintarla
sol=zeros(length(t),N+1);
sol(1,:)=[0 W0 10];
%calculamos W_n+1 a partir de W_n
for j=1:length(t)-1%a continuacion usamos euler explicito
WW=WW-dt*((K+L)*WW+F);
sol(j+1,:)=[0,WW',10]; %guardamos la solucion
end
%dibujamos
[pp,tt]=meshgrid(p,t);
surf(pp,tt,sol)
}}

[[Archivo:explicitobueno.png|marco|centro|Distribución de la temperatura en el disco]]

===Método de Euler implícito===
{{matlab|codigo=
%Datos del problema
clear all
p0=1;
pf=6;
T0=0;
Tf=10;
%datos discretización
N=75;
h=(pf-p0)/N;
p=p0:h:pf;
pint=p0+h:h:pf-h;
%condiciones iniciales
for i=1:length(pint)
if pint(i)<2
u0(i)=100*(pint(i)-1);
elseif pint(i)<5
u0(i)=100;
else
u0(i)=90*(6-pint(i))+10;
end
end
%Datos discretizacion en espacio
dt=h^2/4;%paso en tiempo suele ser más pequeño que el paso en espacio
t=T0:dt:Tf;
%diferencias finitas
%Matriz K Aprox. -u_pp
K=diag(2*ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
%K(N-1,N-2)=2;
K=K/h^2;
%Matriz L Aprox. -u_p
L=diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
L=(L*(1/(2*h)));
L=diag(1./pint)*L;
%Matriz F
F=zeros(N-1,1);
F(N-1,1)=-(10/(h^2))+10/(2*h);
%Dato inicial
W0=[u0];
WW=W0';
%definimos la matriz sol con la u para pintarla
sol=zeros(length(t),N+1);
sol(1,:)=[0 W0 10];
%calculamos W_n+1 a partir de W_n
for j=1:length(t)-1%a continuacion usamos euler implicito
WW=(eye(N-1)+dt*(K+L))\(WW-dt*F);
sol(j+1,:)=[0,WW',10]; %guardamos la solucion
end
%dibujamos
[pp,tt]=meshgrid(p,t);
surf(pp,tt,sol)
}}

[[Archivo:implicitobueno.png|marco|centro|Distribución de la temperatura en el disco]]

===Método de Euler modificado===
{{matlab|codigo=
%Datos del problema
clear all
p0=1;
pf=6;
T0=0;
Tf=10;
%datos discretización
N=75;
h=(pf-p0)/N;
p=p0:h:pf;
pint=p0+h:h:pf-h;
%condiciones iniciales
for i=1:length(pint)
if pint(i)<2
u0(i)=100*(pint(i)-1);
elseif pint(i)<5
u0(i)=100;
else
u0(i)=90*(6-pint(i))+10;
end
end
%Datos discretizacion en espacio
dt=h^2/4;%paso en tiempo suele ser más pequeño que el paso en espacio
t=T0:dt:Tf;
%diferencias finitas
%Matriz K Aprox. -u_pp
K=diag(2*ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
%K(N-1,N-2)=2;
K=K/h^2;
%Matriz L Aprox. -u_p
L=diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
L=(L*(1/(2*h)));
L=diag(1./pint)*L;
%Matriz F
F=zeros(N-1,1);
F(N-1,1)=-(10/(h^2))+10/(2*h);
%Dato inicial
W0=[u0];
WW=W0';
%definimos la matriz sol con la u para pintarla
sol=zeros(length(t),N+1);
sol(1,:)=[0 W0 10];
%calculamos W_n+1 a partir de W_n
for j=1:length(t)-1%a continuacion usamos euler explicito
k1=(K+L)*WW+F;
k2=((dt/2)+(K+L))*WW+F;
WW=WW-(dt/2)*(k1+k2);
sol(j+1,:)=[0,WW',10]; %guardamos la solucion
end
%dibujamos
[pp,tt]=meshgrid(p,t);
surf(pp,tt,sol)
}}

[[Archivo:modificadobueno.png|marco|centro|Distribución de la temperatura en el disco]]

==Evolución de la solución en el tiempo==
Se puede deducir a partir de las siguientes gráficas que para grandes tiempos, la temperatura apenas varía dada la solución caracterizada por una función logarítmica decreciente. En este caso el término <math> u_t</math> se puede despreciar y la temperatura toma un valor estacionario en cada punto de la varilla. Esta ecuación verifica el estado estacionario:
<math> u_{pp}+\frac{1}{p} u_p=0 </math>
Resolviendo la ecuación analíticamente observamos que satisface la función de tipo logarítmica dada por <math> u=\frac{10}{log(6)} log(p) </math>, tal como se observa si representamos la solución del problema añadiendo las siguientes lineas a cualquiera de los métodos realizados anteriormente.
{{matlab|codigo=
plot(p,sol(end,:))
}}

[[Archivo:logaritmica.png|marco|centro|Distribución de la temperatura con el tiempo]]

Este es el código para comparar la diferencia con la solución estacionaria

{{matlab|codigo=
%Datos del problema
clear all
p0=1;
pf=6;
T0=0;
Tf=10;
%datos discretización
h=[0.1,0.01];
color='rg';
hold on
for i=1:2
N=(pf-p0)/h(i);
p=[p0:h(i):pf]';
pint=p0+h(i):h(i):pf-h(i);
%diferencias finitas
%Matriz K Aprox. -u_pp
K=diag(2*ones(N-1,1))-diag(ones(N-2,1),1)-diag(ones(N-2,1),-1);
K=K/h(i)^2;
%Matriz L Aprox. -u_p
L=diag(ones(N-2,1),-1)-diag(ones(N-2,1),1);
L=(L*(1/(2*h(i))));
L=diag(1./pint)*L;
%Matriz F
F=zeros(1,N-1);
F(end)=F(end)+(10/(h(i)^2))+10/(2*5.9*h(i));
u=(K+L)\F';
uu=[T0;u;Tf];
plot(p,uu,color(i));
end
hold off

}}
[[Archivo:1placa.png|marco|centro|Distancia con la solución estacionaria]]

Se observa que la distancia entre la solución estacionaria y la solución u(ρ, t) para t = 0, 1, 2, 10 va aumentando a medida que aumenta el tiempo.

Para t=0 tenemos:

{{matlab|codigo=
%Datos del problema
clear all
p0=1;
pf=6;
T0=0;
Tf=0;
%datos discretización
h=[0.1,0.01];
color='rg';
hold on
for i=1:2
N=(pf-p0)/h(i);
p=[p0:h(i):pf]'/10;
pint=p0+h(i):h(i):pf-h(i);
%diferencias finitas
%Matriz K Aprox. -u_pp
K=diag(2*ones(N-1,1))-diag(ones(N-2,1),1)-diag(ones(N-2,1),-1);
K=K/h(i)^2;
%Matriz L Aprox. -u_p
L=diag(ones(N-2,1),-1)-diag(ones(N-2,1),1);
L=(L*(1/(2*h(i))));
L=diag(1./pint)*L;
%Matriz F
F=zeros(1,N-1);
F(end)=F(end)+(10/(h(i)^2))+10/(2*5.9*h(i));
u=(K+L)\F';
uu=[T0;u;Tf];
plot(p,uu,color(i));
end
hold off
}}

Para t=1:
{{matlab|codigo=
%Datos del problema
clear all
p0=1;
pf=6;
T0=1;
Tf=1;
%datos discretización
h=[0.1,0.01];
color='rg';
hold on
for i=1:2
N=(pf-p0)/h(i);
p=[p0:h(i):pf]'/10;
pint=p0+h(i):h(i):pf-h(i);
%diferencias finitas
%Matriz K Aprox. -u_pp
K=diag(2*ones(N-1,1))-diag(ones(N-2,1),1)-diag(ones(N-2,1),-1);
K=K/h(i)^2;
%Matriz L Aprox. -u_p
L=diag(ones(N-2,1),-1)-diag(ones(N-2,1),1);
L=(L*(1/(2*h(i))));
L=diag(1./pint)*L;
%Matriz F
F=zeros(1,N-1);
F(end)=F(end)+(10/(h(i)^2))+10/(2*5.9*h(i));
u=(K+L)\F';
uu=[T0;u;Tf];
plot(p,uu,color(i));
end
hold off
}}
Para t=2:
{{matlab|codigo=
%Datos del problema
clear all
p0=1;
pf=6;
T0=2;
Tf=2;
%datos discretización
h=[0.1,0.01];
color='rg';
hold on
for i=1:2
N=(pf-p0)/h(i);
p=[p0:h(i):pf]'/10;
pint=p0+h(i):h(i):pf-h(i);
%diferencias finitas
%Matriz K Aprox. -u_pp
K=diag(2*ones(N-1,1))-diag(ones(N-2,1),1)-diag(ones(N-2,1),-1);
K=K/h(i)^2;
%Matriz L Aprox. -u_p
L=diag(ones(N-2,1),-1)-diag(ones(N-2,1),1);
L=(L*(1/(2*h(i))));
L=diag(1./pint)*L;
%Matriz F
F=zeros(1,N-1);
F(end)=F(end)+(10/(h(i)^2))+10/(2*5.9*h(i));
u=(K+L)\F';
uu=[T0;u;Tf];
plot(p,uu,color(i));
end
hold off
}}

Para t=10:
{{matlab|codigo=
%Datos del problema
clear all
p0=1;
pf=6;
T0=10;
Tf=10;
%datos discretización
h=[0.1,0.01];
color='rg';
hold on
for i=1:2
N=(pf-p0)/h(i);
p=[p0:h(i):pf]'/10;
pint=p0+h(i):h(i):pf-h(i);
%diferencias finitas
%Matriz K Aprox. -u_pp
K=diag(2*ones(N-1,1))-diag(ones(N-2,1),1)-diag(ones(N-2,1),-1);
K=K/h(i)^2;
%Matriz L Aprox. -u_p
L=diag(ones(N-2,1),-1)-diag(ones(N-2,1),1);
L=(L*(1/(2*h(i))));
L=diag(1./pint)*L;
%Matriz F
F=zeros(1,N-1);
F(end)=F(end)+(10/(h(i)^2))+10/(2*5.9*h(i));
u=(K+L)\F';
uu=[T0;u;Tf];
plot(p,uu,color(i));
end
hold off
}}

{|
|-
| [[Archivo:placa00.png|thumb|500px|left|t=0]] || [[Archivo:placa01.png|thumb|500px|left|t=1|500px]]
|}

{|
|-
| [[Archivo:placa02.png|thumb|500px|left|t=2]] || [[Archivo:placa10.png|thumb|500px|left|t=10|500px]]
|}

En el caso de dividir la discretización en p por 10 observamos:

{{matlab|codigo=
%Datos del problema
clear all
p0=1;
pf=6;
T0=0;
Tf=10;
%datos discretización
h=[0.1,0.01];
color='rg';
hold on
for i=1:2
N=(pf-p0)/h(i);
p=[p0:h(i):pf]'/10;
pint=p0+h(i):h(i):pf-h(i);
%diferencias finitas
%Matriz K Aprox. -u_pp
K=diag(2*ones(N-1,1))-diag(ones(N-2,1),1)-diag(ones(N-2,1),-1);
K=K/h(i)^2;
%Matriz L Aprox. -u_p
L=diag(ones(N-2,1),-1)-diag(ones(N-2,1),1);
L=(L*(1/(2*h(i))));
L=diag(1./pint)*L;
%Matriz F
F=zeros(1,N-1);
F(end)=F(end)+(10/(h(i)^2))+10/(2*5.9*h(i));
u=(K+L)\F';
uu=[T0;u;Tf];
plot(p,uu,color(i));
end
hold off
}}

[[Archivo:3placa.png|marco|centro|Con la discretización en p una décima parte]]
==Cambio en la condiciones de frontera==

A continuación cambiaremos a unas condiciones tipo Neumann donde en la frontera exterior colocaremos un material aislante. El aislante hace que no haya pérdida de calor en
ese extremo, es decir, el flujo de temperatura en la dirección radial es nulo.

{{matlab|codigo=
L = 5; dr=0.1; N=L/dr; % discretización espacial
r=1+dr:dr:1+L-dr;
r(1)=0; % lo impongo para la primera iteración
i=2:r;
for n=1:length(r)-1 % bucle
r(i+1)= (1/(1+(1/(n*dr))))*(2*r(i)-r(i-1)+(1/(n*dr))*r(i-1));
end
r(length(r))= r(length(r)-1); % condición de contorno tipo newman
U=r; % solución estacionaria
plot(r)
}}

[[Archivo:4placa.png|marco|centro|Solución con las nuevas condiciones de frontera]]
[[Categoría:Ecuaciones Diferenciales]]
[[Categoría:ED13/14]]
[[Categoría:Trabajos 2013-14]]

Placa en forma de Anillo

2014-05-19T21:45:09Z

Alfonso:

{{ TrabajoED | Transmisión de calor en placa con forma de anillo. Grupo 11-B | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED13/14|Curso 2013-14]] | Javier Chamizo, Jacobo Campos, Miguel García, Javier Mellado, Alfonso Ascanio }}

Transmisión de calor en una placa con forma de anillo. Es bastante común modelizar el comportamiento del calor en superficies mediante las ecuaciones diferenciales, y además dada la particularidad de la forma de nuestro caso, utilizando coordenadas polares.

==Introducción==

[[Archivo:placadescripcion.png|marco|derecha|Placa en anillo con p ∈ [1, 6]]]

Tenemos la placa que se ve en la imagen superior con las siguientes condiciones. Colocamos en las fronteras interior ρ = 1 y exterior ρ = 6, que sabemos que mantienen una temperatura constante de 0 grados y 10 grados respectivamente como si estuviesen en contacto con un depósito. En nuestro caso vamos a trabajar en coordenadas polares dada la forma de la placa que se presta a ello. Suponemos también que la temperatura u de la placa solo depende de la coordenada radial y el tiempo, es decir u = u(ρ, t). La placa sabemos que satisface la ecuación del calor que se muestra a continuación:

<math> u_t-\Delta u=0</math>

La placa a su vez dispone de las siguientes condiciones iniciales de temperatura u(p,0):

<math>\left\{\begin{matrix}100·(p-1) → p ∈ [1, 2] \\100 → p ∈ [2, 5]\\90·(6-p)+10 → p ∈ [5, 6]\end{matrix}\right.</math>

Además tenemos condiciones tipo Dirichlet en los extremos ya que se mantiene constante la temperatura a lo largo del tiempo en estos al tener objetos que mantienen la temperatura constante en contacto con los mismos:

<math>\left\{\begin{matrix}u(1,t)=0 \\u(6,t)=10\end{matrix}\right.</math>

Debido a que el ángulo teta y la altura Z se mantienen constantes al ser una figura plana, el Laplaciano en coordenadas cilíndricas es el siguiente:

<math>\Delta u=\frac{1}{p} \frac{\partial }{\partial p} (p \frac{\partial u }{\partial p})+\frac{1 }{p^2} \frac{\partial^2 u }{\partial θ^2}+\frac{\partial^2 u }{\partial z^2}= \frac{\partial^2 u }{\partial p^2}+\frac{1 }{p}\frac{\partial u}{\partial p} </math>

Al aplicar el Laplaciano, la ecuación del calor que se va a resolver queda de la forma:
<math> u_t-u_{pp}-\frac{1}{p} u_p=0 </math>
Quedando así similar al sistema que obtendríamos con una barra con longitud entre 1 y 6.

==Resolución del Sistema por diferencias finitas==
Aquí se va a resolver el problema planteado por el método de diferencias finitas con paso de discretización ∆p = 0.1; y para el tiempo usaremos en cada caso los métodos del trapecio, Euler modificado, explícito e implícito. Todos ellos tomando ∆t = ∆p/4 en t ∈ [0, 10]
===Método del trapecio===
{{matlab|codigo=
%Datos del problema
clear all
p0=1;
pf=6;
T0=0;
Tf=10;
%datos discretización
N=50;
h=(pf-p0)/N;
p=p0:h:pf;
pint=p0+h:h:pf-h;
%condiciones iniciales
for i=1:length(pint)
if pint(i)<2
u0(i)=100*(pint(i)-1);
elseif pint(i)<5
u0(i)=100;
else
u0(i)=90*(6-pint(i))+10;
end
end
%Datos discretizacion en espacio
dt=h/4;%paso en tiempo suele ser más pequeño que el paso en espacio
t=T0:dt:Tf;
%diferencias finitas
%Matriz K Aprox. -u_pp
K=diag(2*ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
%K(N-1,N-2)=2;
K=K/h^2;
%Matriz L Aprox. -u_p
L=diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
L=(L*(1/(2*h)));
L=diag(1./pint)*L;
%Matriz F
F=zeros(N-1,1);
F(N-1,1)=-(10/(h^2))+10/(2*h);
%Dato inicial
W0=[u0];
%Metodo del trapecio
WW=W0';
%definimos la matriz sol con la u para pintarla
sol=zeros(length(t),N+1);
sol(1,:)=[0 W0 10];
%iteraciones desde W^j-->W^j+1
for j=1:length(t)-1
WW=(eye(N-1)+dt/2*(K+L))\((eye(N-1)-dt/2*(K+L))*WW-dt*F);
sol(j+1,:)=[0,WW',10];%guardamos solucion
end
[pp,tt]=meshgrid(p,t);
mesh(pp,tt,sol)
}}
[[Archivo:trapeciobueno.png|marco|centro|Distribución de la temperatura en el disco]]
Ahora particularizamos la solución en ρ = 3 y observamos el comportamiento de la temperatura en dichos puntos. La gráfica siguiente muestra una representación 2D temperatura/tiempo
{{matlab|codigo=
%Datos del problema
clear all
p0=1;
pf=6;
T0=0;
Tf=10;
%datos discretización
N=50;
h=(pf-p0)/N;
p=p0:h:pf;
pint=p0+h:h:pf-h;
%condiciones iniciales
for i=1:length(pint)
if pint(i)<2
u0(i)=100*(pint(i)-1);
elseif pint(i)<5
u0(i)=100;
else
u0(i)=90*(6-pint(i))+10;
end
end
%Datos discretizacion en espacio
dt=h/4;%paso en tiempo suele ser más pequeño que el paso en espacio
t=T0:dt:Tf;
%diferencias finitas
%Matriz K Aprox. -u_pp
K=diag(2*ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
%K(N-1,N-2)=2;
K=K/h^2;
%Matriz L Aprox. -u_p
L=diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
L=(L*(1/(2*h)));
L=diag(1./pint)*L;
%Matriz F
F=zeros(N-1,1);
F(N-1,1)=-(10/(h^2))+10/(2*h);
%Dato inicial
W0=[u0];
%Metodo del trapecio
WW=W0';
%definimos la matriz sol con la u para pintarla
sol=zeros(length(t),N+1);
sol(1,:)=[0 W0 10];
%iteraciones desde W^j-->W^j+1
for j=1:length(t)-1
WW=(eye(N-1)+dt/2*(K+L))\((eye(N-1)-dt/2*(K+L))*WW-dt*F);
sol(j+1,:)=[0,WW',10];%guardamos solucion
end
[pp,tt]=meshgrid(p,t);
J=ones(401,51);
plot3(3*J,tt,sol)
}}

[[Archivo:p3.png|marco|centro|Distribución de la temperatura para p=3]]

===Método de Euler explícito===
{{matlab|codigo=
%Datos del problema
clear all
p0=1;
pf=6;
T0=0;
Tf=10;
%datos discretización
N=75;
h=(pf-p0)/N;
p=p0:h:pf;
pint=p0+h:h:pf-h;
%condiciones iniciales
for i=1:length(pint)
if pint(i)<2
u0(i)=100*(pint(i)-1);
elseif pint(i)<5
u0(i)=100;
else
u0(i)=90*(6-pint(i))+10;
end
end
%Datos discretizacion en espacio
dt=h^2/4;%paso en tiempo suele ser más pequeño que el paso en espacio
t=T0:dt:Tf;
%diferencias finitas
%Matriz K Aprox. -u_pp
K=diag(2*ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
%K(N-1,N-2)=2;
K=K/h^2;
%Matriz L Aprox. -u_p
L=diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
L=(L*(1/(2*h)));
L=diag(1./pint)*L;
%Matriz F
F=zeros(N-1,1);
F(N-1,1)=-(10/(h^2))+10/(2*h);
%Dato inicial
W0=[u0];
WW=W0';
%definimos la matriz sol con la u para pintarla
sol=zeros(length(t),N+1);
sol(1,:)=[0 W0 10];
%calculamos W_n+1 a partir de W_n
for j=1:length(t)-1%a continuacion usamos euler explicito
WW=WW-dt*((K+L)*WW+F);
sol(j+1,:)=[0,WW',10]; %guardamos la solucion
end
%dibujamos
[pp,tt]=meshgrid(p,t);
surf(pp,tt,sol)
}}

[[Archivo:explicitobueno.png|marco|centro|Distribución de la temperatura en el disco]]

===Método de Euler implícito===
{{matlab|codigo=
%Datos del problema
clear all
p0=1;
pf=6;
T0=0;
Tf=10;
%datos discretización
N=75;
h=(pf-p0)/N;
p=p0:h:pf;
pint=p0+h:h:pf-h;
%condiciones iniciales
for i=1:length(pint)
if pint(i)<2
u0(i)=100*(pint(i)-1);
elseif pint(i)<5
u0(i)=100;
else
u0(i)=90*(6-pint(i))+10;
end
end
%Datos discretizacion en espacio
dt=h^2/4;%paso en tiempo suele ser más pequeño que el paso en espacio
t=T0:dt:Tf;
%diferencias finitas
%Matriz K Aprox. -u_pp
K=diag(2*ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
%K(N-1,N-2)=2;
K=K/h^2;
%Matriz L Aprox. -u_p
L=diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
L=(L*(1/(2*h)));
L=diag(1./pint)*L;
%Matriz F
F=zeros(N-1,1);
F(N-1,1)=-(10/(h^2))+10/(2*h);
%Dato inicial
W0=[u0];
WW=W0';
%definimos la matriz sol con la u para pintarla
sol=zeros(length(t),N+1);
sol(1,:)=[0 W0 10];
%calculamos W_n+1 a partir de W_n
for j=1:length(t)-1%a continuacion usamos euler implicito
WW=(eye(N-1)+dt*(K+L))\(WW-dt*F);
sol(j+1,:)=[0,WW',10]; %guardamos la solucion
end
%dibujamos
[pp,tt]=meshgrid(p,t);
surf(pp,tt,sol)
}}

[[Archivo:implicitobueno.png|marco|centro|Distribución de la temperatura en el disco]]

===Método de Euler modificado===
{{matlab|codigo=
%Datos del problema
clear all
p0=1;
pf=6;
T0=0;
Tf=10;
%datos discretización
N=75;
h=(pf-p0)/N;
p=p0:h:pf;
pint=p0+h:h:pf-h;
%condiciones iniciales
for i=1:length(pint)
if pint(i)<2
u0(i)=100*(pint(i)-1);
elseif pint(i)<5
u0(i)=100;
else
u0(i)=90*(6-pint(i))+10;
end
end
%Datos discretizacion en espacio
dt=h^2/4;%paso en tiempo suele ser más pequeño que el paso en espacio
t=T0:dt:Tf;
%diferencias finitas
%Matriz K Aprox. -u_pp
K=diag(2*ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
%K(N-1,N-2)=2;
K=K/h^2;
%Matriz L Aprox. -u_p
L=diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
L=(L*(1/(2*h)));
L=diag(1./pint)*L;
%Matriz F
F=zeros(N-1,1);
F(N-1,1)=-(10/(h^2))+10/(2*h);
%Dato inicial
W0=[u0];
WW=W0';
%definimos la matriz sol con la u para pintarla
sol=zeros(length(t),N+1);
sol(1,:)=[0 W0 10];
%calculamos W_n+1 a partir de W_n
for j=1:length(t)-1%a continuacion usamos euler explicito
k1=(K+L)*WW+F;
k2=((dt/2)+(K+L))*WW+F;
WW=WW-(dt/2)*(k1+k2);
sol(j+1,:)=[0,WW',10]; %guardamos la solucion
end
%dibujamos
[pp,tt]=meshgrid(p,t);
surf(pp,tt,sol)
}}

[[Archivo:modificadobueno.png|marco|centro|Distribución de la temperatura en el disco]]

==Evolución de la solución en el tiempo==
Se puede deducir a partir de las siguientes gráficas que para grandes tiempos, la temperatura apenas varía dada la solución caracterizada por una función logarítmica decreciente. En este caso el término <math> u_t</math> se puede despreciar y la temperatura toma un valor estacionario en cada punto de la varilla. Esta ecuación verifica el estado estacionario:
<math> u_{pp}+\frac{1}{p} u_p=0 </math>
Resolviendo la ecuación analíticamente observamos que satisface la función de tipo logarítmica dada por <math> u=\frac{10}{log(6)} log(p) </math>, tal como se observa si representamos la solución del problema añadiendo las siguientes lineas a cualquiera de los métodos realizados anteriormente.
{{matlab|codigo=
plot(p,sol(end,:))
}}

[[Archivo:logaritmica.png|marco|centro|Distribución de la temperatura con el tiempo]]

Este es el código para comparar la diferencia con la solución estacionaria

{{matlab|codigo=
%Datos del problema
clear all
p0=1;
pf=6;
T0=0;
Tf=10;
%datos discretización
h=[0.1,0.01];
color='rg';
hold on
for i=1:2
N=(pf-p0)/h(i);
p=[p0:h(i):pf]';
pint=p0+h(i):h(i):pf-h(i);
%diferencias finitas
%Matriz K Aprox. -u_pp
K=diag(2*ones(N-1,1))-diag(ones(N-2,1),1)-diag(ones(N-2,1),-1);
K=K/h(i)^2;
%Matriz L Aprox. -u_p
L=diag(ones(N-2,1),-1)-diag(ones(N-2,1),1);
L=(L*(1/(2*h(i))));
L=diag(1./pint)*L;
%Matriz F
F=zeros(1,N-1);
F(end)=F(end)+(10/(h(i)^2))+10/(2*5.9*h(i));
u=(K+L)\F';
uu=[T0;u;Tf];
plot(p,uu,color(i));
end
hold off

}}
[[Archivo:1placa.png|marco|centro|Distancia con la solución estacionaria]]

Se observa que la distancia entre la solución estacionaria y la solución u(ρ, t) para t = 0, 1, 2, 10 va aumentando a medida que aumenta el tiempo.

Para t=0 tenemos:

{{matlab|codigo=
%Datos del problema
clear all
p0=1;
pf=6;
T0=0;
Tf=0;
%datos discretización
h=[0.1,0.01];
color='rg';
hold on
for i=1:2
N=(pf-p0)/h(i);
p=[p0:h(i):pf]'/10;
pint=p0+h(i):h(i):pf-h(i);
%diferencias finitas
%Matriz K Aprox. -u_pp
K=diag(2*ones(N-1,1))-diag(ones(N-2,1),1)-diag(ones(N-2,1),-1);
K=K/h(i)^2;
%Matriz L Aprox. -u_p
L=diag(ones(N-2,1),-1)-diag(ones(N-2,1),1);
L=(L*(1/(2*h(i))));
L=diag(1./pint)*L;
%Matriz F
F=zeros(1,N-1);
F(end)=F(end)+(10/(h(i)^2))+10/(2*5.9*h(i));
u=(K+L)\F';
uu=[T0;u;Tf];
plot(p,uu,color(i));
end
hold off
}}

Para t=1:
{{matlab|codigo=
%Datos del problema
clear all
p0=1;
pf=6;
T0=1;
Tf=1;
%datos discretización
h=[0.1,0.01];
color='rg';
hold on
for i=1:2
N=(pf-p0)/h(i);
p=[p0:h(i):pf]'/10;
pint=p0+h(i):h(i):pf-h(i);
%diferencias finitas
%Matriz K Aprox. -u_pp
K=diag(2*ones(N-1,1))-diag(ones(N-2,1),1)-diag(ones(N-2,1),-1);
K=K/h(i)^2;
%Matriz L Aprox. -u_p
L=diag(ones(N-2,1),-1)-diag(ones(N-2,1),1);
L=(L*(1/(2*h(i))));
L=diag(1./pint)*L;
%Matriz F
F=zeros(1,N-1);
F(end)=F(end)+(10/(h(i)^2))+10/(2*5.9*h(i));
u=(K+L)\F';
uu=[T0;u;Tf];
plot(p,uu,color(i));
end
hold off
}}
Para t=2:
{{matlab|codigo=
%Datos del problema
clear all
p0=1;
pf=6;
T0=2;
Tf=2;
%datos discretización
h=[0.1,0.01];
color='rg';
hold on
for i=1:2
N=(pf-p0)/h(i);
p=[p0:h(i):pf]'/10;
pint=p0+h(i):h(i):pf-h(i);
%diferencias finitas
%Matriz K Aprox. -u_pp
K=diag(2*ones(N-1,1))-diag(ones(N-2,1),1)-diag(ones(N-2,1),-1);
K=K/h(i)^2;
%Matriz L Aprox. -u_p
L=diag(ones(N-2,1),-1)-diag(ones(N-2,1),1);
L=(L*(1/(2*h(i))));
L=diag(1./pint)*L;
%Matriz F
F=zeros(1,N-1);
F(end)=F(end)+(10/(h(i)^2))+10/(2*5.9*h(i));
u=(K+L)\F';
uu=[T0;u;Tf];
plot(p,uu,color(i));
end
hold off
}}

Para t=10:
{{matlab|codigo=
%Datos del problema
clear all
p0=1;
pf=6;
T0=10;
Tf=10;
%datos discretización
h=[0.1,0.01];
color='rg';
hold on
for i=1:2
N=(pf-p0)/h(i);
p=[p0:h(i):pf]'/10;
pint=p0+h(i):h(i):pf-h(i);
%diferencias finitas
%Matriz K Aprox. -u_pp
K=diag(2*ones(N-1,1))-diag(ones(N-2,1),1)-diag(ones(N-2,1),-1);
K=K/h(i)^2;
%Matriz L Aprox. -u_p
L=diag(ones(N-2,1),-1)-diag(ones(N-2,1),1);
L=(L*(1/(2*h(i))));
L=diag(1./pint)*L;
%Matriz F
F=zeros(1,N-1);
F(end)=F(end)+(10/(h(i)^2))+10/(2*5.9*h(i));
u=(K+L)\F';
uu=[T0;u;Tf];
plot(p,uu,color(i));
end
hold off
}}

[[Archivo:2placa.png|marco|centro|Para los diferentes t]]

En el caso de dividir la discretización en p por 10 observamos:

{{matlab|codigo=
%Datos del problema
clear all
p0=1;
pf=6;
T0=0;
Tf=10;
%datos discretización
h=[0.1,0.01];
color='rg';
hold on
for i=1:2
N=(pf-p0)/h(i);
p=[p0:h(i):pf]'/10;
pint=p0+h(i):h(i):pf-h(i);
%diferencias finitas
%Matriz K Aprox. -u_pp
K=diag(2*ones(N-1,1))-diag(ones(N-2,1),1)-diag(ones(N-2,1),-1);
K=K/h(i)^2;
%Matriz L Aprox. -u_p
L=diag(ones(N-2,1),-1)-diag(ones(N-2,1),1);
L=(L*(1/(2*h(i))));
L=diag(1./pint)*L;
%Matriz F
F=zeros(1,N-1);
F(end)=F(end)+(10/(h(i)^2))+10/(2*5.9*h(i));
u=(K+L)\F';
uu=[T0;u;Tf];
plot(p,uu,color(i));
end
hold off
}}

[[Archivo:3placa.png|marco|centro|Con la discretización en p una décima parte]]
==Cambio en la condiciones de frontera==

A continuación cambiaremos a unas condiciones tipo Neumann donde en la frontera exterior colocaremos un material aislante. El aislante hace que no haya pérdida de calor en
ese extremo, es decir, el flujo de temperatura en la dirección radial es nulo.

{{matlab|codigo=
L = 5; dr=0.1; N=L/dr; % discretización espacial
r=1+dr:dr:1+L-dr;
r(1)=0; % lo impongo para la primera iteración
i=2:r;
for n=1:length(r)-1 % bucle
r(i+1)= (1/(1+(1/(n*dr))))*(2*r(i)-r(i-1)+(1/(n*dr))*r(i-1));
end
r(length(r))= r(length(r)-1); % condición de contorno tipo newman
U=r; % solución estacionaria
plot(r)
}}

[[Archivo:4placa.png|marco|centro|Solución con las nuevas condiciones de frontera]]
[[Categoría:Ecuaciones Diferenciales]]
[[Categoría:ED13/14]]
[[Categoría:Trabajos 2013-14]]

Archivo:4placa.png

2014-05-19T21:40:25Z

Alfonso:

Archivo:3placa.png

2014-05-19T21:40:04Z

Alfonso:

Archivo:2placa.png

2014-05-19T21:39:42Z

Alfonso:

Archivo:1placa.png

2014-05-19T21:39:20Z

Alfonso:

Placa en forma de Anillo

2014-05-18T17:26:43Z

Alfonso:

Archivo:Logaritmica.png

2014-05-18T17:25:19Z

Alfonso:

Placa en forma de Anillo

2014-05-18T17:24:58Z

Alfonso:

Archivo:P3.png

2014-05-18T17:12:50Z

Alfonso:

Placa en forma de Anillo

2014-05-18T17:12:31Z

Alfonso:

Archivo:Modificadobueno.png

2014-05-18T17:03:54Z

Alfonso:

Archivo:Implicitobueno.png

2014-05-18T17:03:33Z

Alfonso:

Archivo:Explicitobueno.png

2014-05-18T17:03:02Z

Alfonso:

Archivo:Trapeciobueno.png

2014-05-18T17:02:36Z

Alfonso:

Placa en forma de Anillo

2014-05-18T17:02:00Z

Alfonso:

Placa en forma de Anillo

2014-05-18T16:58:40Z

Alfonso:

Placa en forma de Anillo

2014-05-18T16:53:29Z

Alfonso:

Placa en forma de Anillo

2014-05-14T08:46:20Z

Alfonso:

Placa en forma de Anillo

2014-05-14T08:42:08Z

Alfonso:

Placa en forma de Anillo

2014-05-14T08:40:32Z

Alfonso:

Placa en forma de Anillo

2014-05-14T08:37:02Z

Alfonso:

Placa en forma de Anillo

2014-05-14T08:33:11Z

Alfonso:

Placa en forma de Anillo

2014-05-14T08:28:33Z

Alfonso:

Placa en forma de Anillo

2014-05-14T08:23:35Z

Alfonso:

Placa en forma de Anillo

2014-05-14T08:20:33Z

Alfonso:

Placa en forma de Anillo

2014-05-14T08:18:55Z

Alfonso:

Placa en forma de Anillo

2014-05-14T08:18:01Z

Alfonso:

Placa en forma de Anillo

2014-05-14T08:15:44Z

Alfonso:

Placa en forma de Anillo

2014-05-14T08:15:10Z

Alfonso:

Placa en forma de Anillo

2014-05-14T08:14:42Z

Alfonso:

Placa en forma de Anillo

2014-05-14T08:09:29Z

Alfonso:

Placa en forma de Anillo

2014-05-14T08:07:40Z

Alfonso:

Archivo:Placadescripcion.png

2014-05-14T08:05:47Z

Alfonso:

Placa en forma de Anillo

2014-05-14T08:04:15Z

Alfonso:

Placa en forma de Anillo

2014-05-14T08:01:07Z

Alfonso:

Placa en forma de Anillo

2014-05-14T07:58:19Z

Alfonso:

Placa en forma de Anillo

2014-05-14T07:57:56Z

Alfonso:

Placa en forma de Anillo

2014-05-14T07:57:38Z

Alfonso:

Placa en forma de Anillo

2014-05-14T07:56:38Z

Alfonso:

Placa en forma de Anillo

2014-05-14T07:53:34Z

Alfonso:

Placa en forma de Anillo

2014-05-14T07:49:33Z

Alfonso:

Placa en forma de Anillo

2014-05-14T07:42:00Z

Alfonso:

Placa en forma de Anillo

2014-05-14T06:59:47Z

Alfonso:

Placa en forma de Anillo

2014-05-14T06:58:47Z

Alfonso:

Sistema de masas y muelles

2014-03-05T22:55:48Z

Alfonso:

{{ TrabajoED | Sistema de masas y muelles. Grupo 11-B | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED13/14|Curso 2013-14]] | Javier Mellado, Jacobo Campos, Javier Chamizo, Miguel García, Alfonso Ascanio }}

Sistemas de muelles y masas. Estos sistemas son bastante útiles a la hora de estudiar vibraciones, consiguiendo traducirlo a un lenguaje matemático que sea fácil de estudiar e interpretar.
==Introducción==
En este artículo vamos a estudiar un sistema de 3 masas y 4 muelles en disposición horizontal, entre dos paredes rígidas. Empezaremos estudiándolo primero desde la posición de equilibrio en el caso de no haber rozamiento ni amortiguamiento, por lo que solo actúan las fuerzas restauradoras de los muelles. Tomamos como variables las posiciones de cada masa (x,y,z) en función del tiempo. A partir de ellas escribimos las ecuaciones diferenciales correspondientes a una interpretación según la mecánica de Newton.

[[Archivo:sistema.png|marco|centro|Sistema formado por 3 masas y 4 muelles,anclado a las paredes laterales]]

<math>\left\{\begin{matrix}m_{1}\ddot x=-k_{1}x+k_{2}(y-x)\\m_{2}\ddot y=-k_{2}(y-x)+k_{3}(z-y)\\m_{3}\ddot z=-k_{3}(z-y)-k_{4}(z)\end{matrix}\right.</math>

==Condiciones iniciales y caso particular==

Suponemos que en el instante t=0 las tres masas están desplazadas 0.5, 1 y 0.8 metros hacia la derecha de la posici ́on de equilibrio y se sueltan repentinamente, sin velocidad. Vamos a suponer que k1=4N/m, k2=2N/m, k3=1N/m, k4=3N/m, m1=2kg, m2=1kg, m3=3kg, que la distancia entre las paredes es de 12 metros y que en equilibrio las tres masas están en las posiciones 2.5, 4 y 8.Resolveremos el problema numéricamente utilizando los métodos de Euler modificado y Runge Kutta con pasos de h=0.1m y h=0.025m.

Sustituyendo ahora los valores dados arriba en las ecuaciones obtenemos las siguientes, con las cuales operaremos. Para aplicar el método reduciremos el sistema a uno de primer orden ayudándonos de las siguientes ecuaciones.

<math>\dot x=u\\\dot y=v\\\dot z=w</math>

Y lo juntamos con el ya existente para obtener:

<math>\left\{\begin{matrix}m_{1}\ddot x=-k_{1}x+k_{2}(y-x)\\m_{2}\ddot y=-k_{2}(y-x)+k_{3}(z-y)\\m_{3}\ddot z=-k_{3}(z-y)-k_{4}(z)\\\dot x=u\\\dot y=v\\\dot z=w\end{matrix}\right.</math>

===Método de Euler Modificado===
{{matlab|codigo=
% Euler Modificado h=0.1
% Sistema de 3 masas y 4 muelles
clear all
% Intervalo de tiempo
t0=0;tN=10;
% Datos
k1=4;k2=2;k3=1;k4=3;
m1=2;m2=1;m3=3;
%Condiciones iniciales
x0=0.5;y0=1;z0=0.8;
u0=0;v0=0;w0=0;

N=100; % Pasos para h=0.1
h=(tN-t0)/N; % Intervalo

% Primeros valores
xx(1)=x0;yy(1)=y0;zz(1)=z0;uu(1)=u0;vv(1)=v0;ww(1)=w0;

for n=1:N
% Primera K
k1x=uu(n);
k1y=vv(n);
k1z=ww(n);
k1u=-(k2+k1)/m1*xx(n)+k2/m1*yy(n);
k1v=k2/m2*xx(n)-((k2+k3)/m2)*yy(n)+(k3/m2)*zz(n);
k1w=k3/m3*yy(n)-((k4+k3)/m3)*zz(n);

xa=xx(n)+h/2*k1x;
ya=yy(n)+h/2*k1y;
za=zz(n)+h/2*k1z;
ua=uu(n)+h/2*k1u;
va=vv(n)+h/2*k1v;
wa=ww(n)+h/2*k1w;

% Segunda K
k2x=ua;
k2y=va;
k2z=wa;
k2u=-(k2+k1)/m1*xa+k2/m1*ya;
k2v=k2/m2*xa-((k2+k3)/m2)*ya+(k3/m2)*za;
k2w=k3/m3*ya-((k4+k3)/m3)*za;

xx(n+1)=xx(n)+h*k2x;
yy(n+1)=yy(n)+h*k2y;
zz(n+1)=zz(n)+h*k2z;
uu(n+1)=uu(n)+h*k2u;
vv(n+1)=vv(n)+h*k2v;
ww(n+1)=ww(n)+h*k2w;
end
hold on
t=t0:h:tN;
e=t0:0.01:tN
plot(xx+2.5,t,'-')
plot(yy+4,t,'-r')
plot(zz+8,t,'-k')
plot(2.5,e,'-g')
plot(4,e,'-g')
plot(8,e,'-g')
hold off
axis([0,12,0,10])
hold on
title('Euler Modificado')
legend('Masa m1','Masa m2','Masa m3','Posicion de equilibrio')
hold off
}}

Aqui podemos ver los gráficos obtenidos. El de la izquierda corresponde a un paso de 0.1 y el derecho a uno de 0.025. El código es el mismo para ambos a excepción de N=100 para el primero y N=40 para el segundo.

{|
|-
| [[Archivo:em20.1.png|thumb|600px|left|Posición de cada masa respecto de la posición de equilibrio en intervalo (0,10) con paso de 0.1]] || [[Archivo:em20.25.png|thumb|600px|left|Posición de cada masa respecto de la posición de equilibrio en intervalo (0,10) con paso de 0.025|600px]]
|}

===Método de Runge-Kutta===

El código de Matlab es:
{{matlab|codigo=

% Sistema de 3 masas y 4 muelles
t0=0;tN=10;
% Datos
k1=4;k2=2;k3=1;k4=3;
m1=2;m2=1;m3=3;
%Condiciones iniciales
x0=0.5;y0=1;z0=0.8;u0=0;v0=0;w0=0;

N=100; % Pasos
h=(tN-t0)/N; % Intervalo

% Primeros valores
xx(1)=x0;yy(1)=y0;zz(1)=z0;uu(1)=u0;vv(1)=v0;ww(1)=w0;

% Metodo de Runge-Kutta
for n=1:N
% Primera K
k1x=uu(n);
k1y=vv(n);
k1z=ww(n);
k1u=-(k2+k1)/m1*xx(n)+k2/m1*yy(n);
k1v=k2/m2*xx(n)-((k2+k3)/m2)*yy(n)+(k3/m2)*zz(n);
k1w=k3/m3*yy(n)-((k4+k3)/m3)*zz(n);

% Segunda K
k2x=uu(n)+1/2*h*k1u;
k2y=vv(n)+1/2*h*k1v;
k2z=ww(n)+1/2*h*k1w;
k2u=-(k2+k1)/m1*(xx(n)+1/2*h*k1x)+k2/m1*(yy(n)+1/2*k1y*h);
k2v=k2/m2*(xx(n)+1/2*h*k1x)-((k2+k3)/m2)*(yy(n)+1/2*k1y*h)+(k3/m2)*(zz(n)+1/2*k1z*h);
k2w=k3/m3*(yy(n)+1/2*k1y*h)-((k4+k3)/m3)*(zz(n)+1/2*k1z*h);

% Tercera K
k3x=uu(n)+1/2*h*k2u;
k3y=vv(n)+1/2*h*k2v;
k3z=ww(n)+1/2*h*k2w;
k3u=-(k2+k1)/m1*(xx(n)+1/2*h*k2x)+k2/m1*(yy(n)+1/2*k2y*h);
k3v=k2/m2*(xx(n)+1/2*h*k2x)-((k2+k3)/m2)*(yy(n)+1/2*k2y*h)+(k3/m2)*(zz(n)+1/2*k2z*h);
k3w=k3/m3*(yy(n)+1/2*k2y*h)-((k4+k3)/m3)*(zz(n)+1/2*k2z*h);

% Cuarta K
k4x=uu(n)+h*k3u;
k4y=vv(n)+h*k3v;
k4z=ww(n)+h*k3w;
k4u=-(k2+k1)/m1*(xx(n)+h*k3x)+k2/m1*(yy(n)+k3y*h);
k4v=k2/m2*(xx(n)+h*k3x)-((k2+k3)/m2)*(yy(n)+k3y*h)+(k3/m2)*(zz(n)+k3z*h);
k4w=k3/m3*(yy(n)+h*k3y)-((k4+k3)/m3)*(zz(n)+k3z*h);

xx(n+1)=xx(n)+(h/6)*(k1x+2*k2x+2*k3x+k4x);
yy(n+1)=yy(n)+(h/6)*(k1y+2*k2y+2*k3y+k4y);
zz(n+1)=zz(n)+(h/6)*(k1z+2*k2z+2*k3z+k4z);
uu(n+1)=uu(n)+(h/6)*(k1u+2*k2u+2*k3u+k4u);
vv(n+1)=vv(n)+(h/6)*(k1v+2*k2v+2*k3v+k4v);
ww(n+1)=ww(n)+(h/6)*(k1w+2*k2w+2*k3w+k4w);
end

e=t0:0.01:tN
t=t0:h:tN;
axis([0,12,0,10])
hold on
title('Runge-Kutta')
plot(xx+2.5,t,'-')
plot(yy+4,t,'-r')
plot(zz+8,t,'-k')
plot(2.5,e,'-g')
plot(4,e,'-g')
plot(8,e,'-g')
legend('Masa m1','Masa m2','Masa m3','Posicion de equilibrio')
hold off
}}
Aqui podemos ver los gráficos obtenidos. El superior corresponde a un paso de 0.1 y el inferior a uno de 0.025. El código es el mismo para ambos a excepción de N=100 para el primero y N=40 para el segundo.

{|
|-
| [[Archivo:rk20.1.png|thumb|600px|left|Posición de cada masa respecto de la posición de equilibrio en intervalo (0,10) con paso de 0.1]] || [[Archivo:rk20.25.png|thumb|600px|left|Posición de cada masa respecto de la posición de equilibrio en intervalo (0,10) con paso de 0.025|600px]]
|}

==Otras condiciones iniciales==

Ahora vamos a imponer unas nuevas CI para representar dos nuevas situaciones. Dichas condiciones son:

{|
|-
| A<math>\left\{\begin{matrix}x(0)=0.5 \\\dot x(0)=1\\y(0)=1 \\\dot y(0)=1\\z(0)=-0.5 \\\dot z=-1\\\end{matrix}\right.</math> || B<math>\left\{\begin{matrix}x(0)=0.5 \\\dot x(0)=-1\\y(0)=1 \\\dot y(0)=0\\z(0)=-0.5 \\\dot z=0.5\\\end{matrix}\right.</math>
|}

===Método del trapecio===
Utilizando las condiciones iniciales A y con un h=0.2 metros durante los primeros 20 segundos obtenemos:

{{matlab|codigo=
%Sistema EDOS lineal
%Metodo de Trapecio Apartado A
%%datos del problema
t0=0;
tf=20;
y0=[-0.5 1 1 1 -0.5 1]';
k1=4;
k2=2;
k3=1;
k4=3;
m1=2;
m2=1;
m3=3;
u=-1;
v=1;
w=1;

A=[0 1 0 0 0 0; -(k2+k1)/m1 0 k2/m1 0 0 0; 0 0 0 1 0 0; k2/m2 0 -((k2+k3)/m2) 0 (k3/m2) 0; 0 0 0 0 0 1 ; 0 0 k3/m3 0 -((k4+k3)/m3) 0];
%%discretizacion
N=100;
h=(tf-t0)/N;
%%Vectores de tiempo y solucion approx.
t=t0:h:tf;
y=zeros(6,N+1);%%matriz de ceros con 3 filas (por no. de ecus)y N+1 columnas
%%Inicializacion
y(:,1)=y0;
yy=y0;
%%necesito gn y gn+1
for n=1:N
gn=[0 0 0 0 0 0]';%vector con terminos libres de y
gnp1=[0 0 0 0 0 0]';
yy=(eye(6)-h/2*A)\(((eye(6)+h/2*A)*yy)+h/2*(gn+gnp1));
y(:,n+1)=yy;
end
clf

e=t0:0.05:tf
hold on
title('Trapecio A')
plot(y(1,:)+2.5,t,'-')
plot(y(1,:)+4,t,'-r')
plot(y(1,:)+8,t,'-k')
plot(2.5,e,'-g')
plot(4,e,'-g')
plot(8,e,'-g')
legend('Masa m1','Masa m2','Masa m3','Posicion de equilibrio')
hold off
}}

Si bien usando los mismos parámetros de h y tiempo con las condiciones B tenemos:

{{matlab|codigo=

%Sistema EDOS lineal
%Metodo de Trapecio Apartado B
%%datos del problema
t0=0;
tf=20;
y0=[0.5 -1 1 0 -0.5 0.5]';
k1=4;
k2=2;
k3=1;
k4=3;
m1=2;
m2=1;
m3=3;
u=-1;
v=0;
w=0.5;

A=[0 1 0 0 0 0; -(k2+k1)/m1 0 k2/m1 0 0 0; 0 0 0 1 0 0; k2/m2 0 -((k2+k3)/m2) 0 (k3/m2) 0; 0 0 0 0 0 1 ; 0 0 k3/m3 0 -((k4+k3)/m3) 0];
%%discretizacion
N=100;
h=(tf-t0)/N;
%%Vectores de tiempo y solucion approx.
t=t0:h:tf;
y=zeros(6,N+1);%%matriz de ceros con 3 filas (por no. de ecus)y N+1 columnas
%%Inicializacion
y(:,1)=y0;
yy=y0;
%%necesito gn y gn+1
for n=1:N
gn=[0 0 0 0 0 0]';%vector con terminos libres de y
gnp1=[0 0 0 0 0 0]';
yy=(eye(6)-h/2*A)\(((eye(6)+h/2*A)*yy)+h/2*(gn+gnp1));
y(:,n+1)=yy;
end
clf

e=t0:0.05:tf
hold on
title('Trapecio B')
plot(y(1,:)+2.5,t,'-')
plot(y(1,:)+4,t,'-r')
plot(y(1,:)+8,t,'-k')
plot(2.5,e,'-g')
plot(4,e,'-g')
plot(8,e,'-g')
legend('Masa m1','Masa m2','Masa m3','Posicion de equilibrio')
hold off
}}

{|
|-
| [[Archivo:tresatrap.png|thumb|600px|left|Posición de cada masa respecto de la posición de equilibrio en intervalo (0,20) condiciones A]] || [[Archivo:tresbtrap.png|thumb|600px|left|Posición de cada masa respecto de la posición de equilibrio en intervalo (0,20) condiciones B|600px]]
|}

===Método de Euler Modificado===
Ahora utilizaremos el método de Euler Modificado junto con las condiciones A con un h=0.05 metros y un tiempo de 10 segundos:

-->Para la primera y tercera masas tenemos:

{{matlab|codigo=

% Euler Modificado Apartado A con Masas 1 y 3
% Sistema de 3 masas y 4 muelles
clear all
% Intervalo de tiempo
t0=0;
tN=10;
% Datos
k1=4;
k2=2;
k3=1;
k4=3;
m1=2;
m2=1;
m3=3;
%Condiciones iniciales
x0=-0.5;
y0=1;
z0=-0.5;
u0=1;
v0=1;
w0=1;

N=50; % Pasos
h=(tN-t0)/N; % Intervalo

% Primeros valores
xx(1)=x0;
yy(1)=y0;
zz(1)=z0;
uu(1)=u0;
vv(1)=v0;
ww(1)=w0;

for n=1:N
% Primera K
k1x=uu(n);
k1y=vv(n);
k1z=ww(n);
k1u=-(k2+k1)/m1*xx(n)+k2/m1*yy(n);
k1v=k2/m2*xx(n)-((k2+k3)/m2)*yy(n)+(k3/m2)*zz(n);
k1w=k3/m3*yy(n)-((k4+k3)/m3)*zz(n);

xa=xx(n)+h/2*k1x;
ya=yy(n)+h/2*k1y;
za=zz(n)+h/2*k1z;
ua=uu(n)+h/2*k1u;
va=vv(n)+h/2*k1v;
wa=ww(n)+h/2*k1w;

% Segunda K
k2x=ua;
k2y=va;
k2z=wa;
k2u=-(k2+k1)/m1*xa+k2/m1*ya;
k2v=k2/m2*xa-((k2+k3)/m2)*ya+(k3/m2)*za;
k2w=k3/m3*ya-((k4+k3)/m3)*za;

xx(n+1)=xx(n)+h*k2x;
yy(n+1)=yy(n)+h*k2y;
zz(n+1)=zz(n)+h*k2z;
uu(n+1)=uu(n)+h*k2u;
vv(n+1)=vv(n)+h*k2v;
ww(n+1)=ww(n)+h*k2w;
end
hold on
t=t0:h:tN;
plot(zz+8-(xx+2),t,'-k')
hold off
axis([0,12,0,10])
hold on
title('Euler Modificado Apartado A con Masas 1 y 3')

legend('Dif entre m3,m1')
hold off
}}

-->Para todas las masas tenemos:

{{matlab|codigo=

% Euler Modificado Apartado A con Todas las Masas
% Sistema de 3 masas y 4 muelles
clear all
% Intervalo de tiempo
t0=0;
tN=10;
% Datos
k1=4;
k2=2;
k3=1;
k4=3;
m1=2;
m2=1;
m3=3;
%Condiciones iniciales
x0=-0.5;
y0=1;
z0=-0.5;
u0=1;
v0=1;
w0=1;

N=50; % Pasos
h=(tN-t0)/N; % Intervalo

% Primeros valores
xx(1)=x0;
yy(1)=y0;
zz(1)=z0;
uu(1)=u0;
vv(1)=v0;
ww(1)=w0;

for n=1:N
% Primera K
k1x=uu(n);
k1y=vv(n);
k1z=ww(n);
k1u=-(k2+k1)/m1*xx(n)+k2/m1*yy(n);
k1v=k2/m2*xx(n)-((k2+k3)/m2)*yy(n)+(k3/m2)*zz(n);
k1w=k3/m3*yy(n)-((k4+k3)/m3)*zz(n);

xa=xx(n)+h/2*k1x;
ya=yy(n)+h/2*k1y;
za=zz(n)+h/2*k1z;
ua=uu(n)+h/2*k1u;
va=vv(n)+h/2*k1v;
wa=ww(n)+h/2*k1w;

% Segunda K
k2x=ua;
k2y=va;
k2z=wa;
k2u=-(k2+k1)/m1*xa+k2/m1*ya;
k2v=k2/m2*xa-((k2+k3)/m2)*ya+(k3/m2)*za;
k2w=k3/m3*ya-((k4+k3)/m3)*za;

xx(n+1)=xx(n)+h*k2x;
yy(n+1)=yy(n)+h*k2y;
zz(n+1)=zz(n)+h*k2z;
uu(n+1)=uu(n)+h*k2u;
vv(n+1)=vv(n)+h*k2v;
ww(n+1)=ww(n)+h*k2w;
end
hold on
t=t0:h:tN;
plot((yy+4)-(xx+2.5),t,'-k')
plot((zz+8)-(yy+4)+((yy+4)-(xx+2.5)),t,'-g')
hold off
axis([0,12,0,10])
hold on
title('Euler Modificado')

legend('Dif entre m1,m2','Dif entre m2,m3')
hold off
}}

{|
|-
| [[Archivo:tresa13.png|thumb|600px|left|Diferencia relativa entre las posiciones m1,m3]] || [[Archivo:tresa123.png|thumb|600px|left|Diferencia relativa entre las posiciones m1,m2,m3|600px]]
|}

Variando a continuación solo las condiciones iniciales utilizando las condiciones de B:

-->Para la primera y tercera masas tenemos:

{{matlab|codigo=
% Euler Modificado Apartado B con Masas 1 y 3
% Sistema de 3 masas y 4 muelles
clear all
% Intervalo de tiempo
t0=0;
tN=10;
% Datos
k1=4;
k2=2;
k3=1;
k4=3;
m1=2;
m2=1;
m3=3;
%Condiciones iniciales
x0=0.5;
y0=1;
z0=-0.5;
u0=-1;
v0=0;
w0=0.5;

N=50; % Pasos
h=(tN-t0)/N; % Intervalo

% Primeros valores
xx(1)=x0;
yy(1)=y0;
zz(1)=z0;
uu(1)=u0;
vv(1)=v0;
ww(1)=w0;

for n=1:N
% Primera K
k1x=uu(n);
k1y=vv(n);
k1z=ww(n);
k1u=-(k2+k1)/m1*xx(n)+k2/m1*yy(n);
k1v=k2/m2*xx(n)-((k2+k3)/m2)*yy(n)+(k3/m2)*zz(n);
k1w=k3/m3*yy(n)-((k4+k3)/m3)*zz(n);

xa=xx(n)+h/2*k1x;
ya=yy(n)+h/2*k1y;
za=zz(n)+h/2*k1z;
ua=uu(n)+h/2*k1u;
va=vv(n)+h/2*k1v;
wa=ww(n)+h/2*k1w;

% Segunda K
k2x=ua;
k2y=va;
k2z=wa;
k2u=-(k2+k1)/m1*xa+k2/m1*ya;
k2v=k2/m2*xa-((k2+k3)/m2)*ya+(k3/m2)*za;
k2w=k3/m3*ya-((k4+k3)/m3)*za;

xx(n+1)=xx(n)+h*k2x;
yy(n+1)=yy(n)+h*k2y;
zz(n+1)=zz(n)+h*k2z;
uu(n+1)=uu(n)+h*k2u;
vv(n+1)=vv(n)+h*k2v;
ww(n+1)=ww(n)+h*k2w;
end
hold on
t=t0:h:tN;
plot(zz+8-(xx+2),t,'-k')
hold off
axis([0,12,0,10])
hold on
title('Euler Modificado Apartado B con Masas 1 y 3')

legend('Dif entre m3,m1')
hold off
}}
-->Para todas las masas tenemos:

{{matlab|codigo=

% Euler Modificado Apartado B con Todas las Masas
% Sistema de 3 masas y 4 muelles
clear all
% Intervalo de tiempo
t0=0;
tN=10;
% Datos
k1=4;
k2=2;
k3=1;
k4=3;
m1=2;
m2=1;
m3=3;
%Condiciones iniciales
x0=0.5;
y0=1;
z0=-0.5;
u0=-1;
v0=0;
w0=0.5;

N=50; % Pasos
h=(tN-t0)/N; % Intervalo

% Primeros valores
xx(1)=x0;
yy(1)=y0;
zz(1)=z0;
uu(1)=u0;
vv(1)=v0;
ww(1)=w0;

for n=1:N
% Primera K
k1x=uu(n);
k1y=vv(n);
k1z=ww(n);
k1u=-(k2+k1)/m1*xx(n)+k2/m1*yy(n);
k1v=k2/m2*xx(n)-((k2+k3)/m2)*yy(n)+(k3/m2)*zz(n);
k1w=k3/m3*yy(n)-((k4+k3)/m3)*zz(n);

xa=xx(n)+h/2*k1x;
ya=yy(n)+h/2*k1y;
za=zz(n)+h/2*k1z;
ua=uu(n)+h/2*k1u;
va=vv(n)+h/2*k1v;
wa=ww(n)+h/2*k1w;

% Segunda K
k2x=ua;
k2y=va;
k2z=wa;
k2u=-(k2+k1)/m1*xa+k2/m1*ya;
k2v=k2/m2*xa-((k2+k3)/m2)*ya+(k3/m2)*za;
k2w=k3/m3*ya-((k4+k3)/m3)*za;

xx(n+1)=xx(n)+h*k2x;
yy(n+1)=yy(n)+h*k2y;
zz(n+1)=zz(n)+h*k2z;
uu(n+1)=uu(n)+h*k2u;
vv(n+1)=vv(n)+h*k2v;
ww(n+1)=ww(n)+h*k2w;
end
hold on
t=t0:h:tN;
plot((yy+4)-(xx+2.5),t,'-k')
plot((zz+8)-(yy+4)+((yy+4)-(xx+2.5)),t,'-g')
hold off
axis([0,12,0,10])
hold on
title('Euler Modificado')

legend('Dif entre m1,m2','Dif entre m2,m3')
hold off
}}

{|
|-
| [[Archivo:tresb13.png|thumb|600px|left|Diferencia relativa entre las posiciones m1,m3]] || [[Archivo:tresb123.png|thumb|600px|left|Diferencia relativa entre las posiciones m1,m2,m3|600px]]
|}

==Estudio de un sistema en medio viscoso y con fuerza exterior==
En muchos de los casos que estudian vibraciones es necesario tener en cuenta un amortiguador que va frenando las vibraciones. Vamos a estudiar el caso de una única masa (m=1kg), <math>\mu=c=1</math> , k=4 y una fuerza exterior <math>f(t)=2·exp(-0.01t)·sin(2t)</math>. Las condiciones iniciales son x(0)=0.3 v(0)=0.3
Aplicando las ecuaciones de Newton obtenemos la siguiente ecuación diferencial:
:<math>\ddot x+\dot x +kx -2e^{(-0.01t)·sin(2t)}=0</math>

Para resolver numericamente este problema aplicamos el método del trapecio. Necesitamos pasar de la ecuación de segundo orden a un sistema de dos ecuaciones de primer orden, el cual queda así:

<math>\left\{\begin{matrix}\dot v+v+4x-2e^{(-0.01t)·sin(2t)}=0\\\dot x=v\end{matrix}\right.</math>

Si operamos el sistema aplicando la fórmula del trapecio y despejamos dejando como incógnita los términos en (n+1) resolvemos con el siguiente código de paso h=0.2:

{{matlab|codigo=
a=0
b=60
N=300
h=(b-a)/N
t=a:h:b;
t0=a
x0=0.3
v0=0.3
k=4;
c=1;
x(1)=x0;v(1)=v0;E(1)=2*(x0^2)+0.5*(v0^2)

for n=1:N
A=[x(n)-c*h/2*v(n)+h/2*k*x(n)+h/2*(2*exp(-0.01*t(n))*sin(2*t(n)))+2*exp(-0.01*t(n+1)*sin(2*t(n+1)));
c*(1+h/2)*v(n)]
B=[1-h/2*k,c*h/2;0,c]
rest=A\B;
x(n+1)=rest(1);
v(n+1)=rest(2);
end

E=2*(x.^2)+0.5*(v.^2);
semilogx(t,E,'g')
}}

Ahora vemos la ganancia de energía del sistema a medida que transcurre el tiempo, en escala logarítmica. Primero vemos para c=1 y a continuación para c=3. Es bastante razonable que se disipe más energía con c más alta porque la viscosidad del medio es mayor y sería necesaria una mayor energía para mantener el movimiento. Sin embargo en nuestro caso se aprecia que aumenta y se puede deber a que la fuerza externa aporta más energía al sistema de la que disipa el medio viscoso.

[[Archivo:energiac1.png|marco|centro|Energía del sistema en el primer minuto para c=1]] || [[Archivo:energiac3.png|marco|centro|Energía del sistema en el primer minuto para c=3]]

[[Categoría:Ecuaciones Diferenciales]]
[[Categoría:ED13/14]]
[[Categoría:Trabajos 2013-14]]