Ver el código de «Ecuacion de calor (Grupo 25C)»

{{ TrabajoED | ECUACIÓN DE CALOR. GRUPO 25C | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED14/15|Curso 2014-15]]|Gálvez Aparici, Antonio<br />Megino León, Guillermo<br />Popa, Silviu<br />Sistac Ara, Alejandro<br />Veiga López, Roberto}}

== Introducción ==

Se considera una placa en forma de anillo comprendida entre los radios ρ=1 y ρ=3 y unos objetos colocados en estas fronteras que mantienen la placa a una temperatura constante de 3 y 0 grados respectivamente. 


[[Archivo:Placa.jpg|marco|centro|Representación gráfica de la placa]]


Consideramos también que su temperatura inicial sigue la siguiente función:

[[Archivo:Condicioninicial.jpg|marco|centro]]

Todas estas características llevan al desarrollo de un problema de calor cuyo sistema se va a plantear a continuación.

==Sistema de ecuaciones en coordenadas polares==

En este caso se va a suponer que la temperatura "U" de la placa solamente depende de la coordenada radial y el tiempo, descartándose así la coordenada cíclica theta. Además se sabe que "U" satisface la siguiente ecuación:

[[Archivo:EcuacionU.jpg|marco|centro]]

Sabiendo esto se puede desarrollar un sistema de ecuaciones en coordenadas polares que representa el problema planteado.

Al no tener dependencia de la variable “θ” el problema que nos queda es de contorno en “ρ”:


<math> \Delta U = U_{ρρ} +\frac{1 }{ρ} U_ρ;    (1) </math>

Aplicando la ecuación que nuestro problema tiene que satisfacer:


<math> U_t - \Delta U = 0;    (2)</math>

Sustituyendo (1) en (2):

<math> U_t - U_{ρρ} -\frac{1}{ρ} U_ρ = 0 </math>


Añadiendo las condiciones de contorno e inicial nos queda el siguiente planteamiento:

<math> \left\{\begin{matrix} U_t - U_{ρρ} -\frac{1 }{ρ} U_ρ = 0  \\ U(1,t)=3 ;  U(3,t)=0 \\U(ρ,0)=e^{-8(ρ-2)^2}+3-ρ \\ ρ∈(1,3);t>0 \end{matrix}\right.</math>

==Resolución del problema==

Se resuelve por el método de diferencias finitas discretizando la variable espacial "ρ" y la variable temporal "t". En este caso tenemo dos matrices K1 y K2 que acompañan el vector "U" ya que el planteamiento lo exige por las apariciones de las derivadas primera y segunda de U con respecto de "ρ". 

===Método del trapecio===

En este caso se pide realizar el metodo con los siguientes pasos:

[[Archivo:Pl4.jpg|centro]]

{{matlab|codigo=
clear all
clc
%Discretización espacial
a=1;b=3;
h=0.1;
N=(b-a)/h;
r=a:h:b;
rr=r(2:N); %Elementos interiores del vector 
%Definicion de las variables K1 K2
K1=(1/(h^2))*(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1)-diag(ones(1,N-2),1));
K2=(1/(2*h))*(-diag(ones(1,N-2),1)+diag(ones(1,N-2),-1));
G=diag(1./rr);
K2=G*K2;
%Definicion de F y valores de contorno
F=zeros(1,N-1);Ua=3;Ub=0;
F(1)=F(1)+Ua*((1/h^2)-(1/(2*h*rr(1))));
F(end)=F(end)+Ub*((1/h^2)+(1/(2*h*rr(end))));
%Discretizacion temporal
t0=0;tM=10;
k=h/4;
t=t0:k:tM;
M=length(t)-1;
%Definicion del valor inicial
U0=exp(-8.*((rr-2).^2))+3-rr;
sol(:,1)=U0;
%Solucion
for i=1:M
    sol(:,i+1)=(eye(size(K1))+(k/2)*(K1+K2))\((eye(size(K1))-(k/2)*(K1+K2))*sol(:,i)+k*F');
end
%Grafico
UA=Ua*ones(1,length(t));
UB=Ub*ones(1,length(t));
sol=[UA;sol;UB];
[Mr,Mt]=meshgrid(r,t);
s=surf(Mr,Mt,sol');
set(s,'edgecolor','flat');
colormap summer
}}

[[Archivo:Graficatrapecio.jpg|marco|centro|Gráfico del problema mediante el método del trapecio]]

===Temperatura para ρ=2===

Para comprobar como cambia la temperatura en función del tiempo para ρ=2 hay que implementar el siguiente codigo en Matlab, a continuación del anterior:
{{matlab|codigo=
d=(2-a)/h+1;
plot(t,sol(d,:));
}}

Como resultado nos sale el siguiente gráfico donde se observa que para un radio 2, la temperatura baja desde los dos grados hasta los 1.1 para luego mantenerse constante:

[[Archivo:Graficaro2.jpg|marco|centro|Gráfico para ρ=2]]

===Método de Euler explícito===
Para obtener tanto este método como el de Euler implícito y el de Heun, lo unico que hay que modificar es el apartado "%Solucion" del código del método del trapecio.

{{matlab|codigo=
%Solucion
for i=1:M
    sol(:,i+1)=sol(:,i)+k*(-(K1+K2)*sol(:,i)+F');
end
}}

En este caso se han tenido problemas con la obtención total de los datos y del gráfico para un tiempo de 10 segundos ya que el método es muy impreciso y da valores anormalmente altos. Para un t=1 se ha podido obtener el siguiente gráfico:

[[Archivo:Graficoeulerexplicito.jpg|marco|centro|Gráfico con el método de Euler explítcito]]

===Método de Euler implícito===

{{matlab|codigo=
%Solucion
for i=1:M
    sol(:,i+1)=(eye(size(K1))+k*(K1+K2))\(sol(:,i)+k*F');
end
}}

La gráfica obtenida es muy parecida a la del método del trapecio:

[[Archivo:GraficaEulerimplicito.jpg|marco|centro|Gráfico meiante el método de Euler implícito]]

===Método de Heun===

El código empleado ha sido el siguiente:
{{matlab|codigo=
%Solucion
for i=1:M
    Q1=(-(K1+K2))*sol(:,1)+F';
    Q2=(-(K1+K2))*(sol(:,1)+k*Q1)+F';
    sol(:,i+1)=sol(:,1)+(k/2)*(Q1+Q2);
end
}}

[[Archivo:GraficoHeun.jpg|marco|centro|Gráfico mediante el método de Heun]]

==Planteamiento del problema para tiempos grandes==

En este caso es el mismo planteamiento que el anterior pero sin la variable dependiente del tiempo. El problema se reduce a un problema de contorno que se resuelve mediante el siguiente código:

{{matlab|codigo=
clear all
clc
%Discretización espacial
a=1;b=3;
h=0.1;
N=(b-a)/h;
r=a:h:b;
rr=r(2:N); %Elementos interiores del vector 
%Definicion de las variables K1 K2
K1=(1/(h^2))*(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1)-diag(ones(1,N-2),1));
K2=(1/(2*h))*(-diag(ones(1,N-2),1)+diag(ones(1,N-2),-1));
G=diag(1./rr);
K2=G*K2;
%Definicion de F y valores de contorno
F=zeros(1,N-1);
Ua=3;Ub=0;
F(1)=F(1)+Ua*((1/h^2)-(1/(2*h*rr(1))));
F(end)=F(end)+Ub*((1/h^2)+(1/(2*h*rr(end))));
%Discretizacion temporal
Y=(K1+K2)\F';
Y=[Ua;Y;Ub];
plot(r,Y);
}}

La gráfica resultante nos muestra a qué debería parecerse la función de temperatura para tiempos grandes:

[[Archivo:Asgag.jpg|marco|centro]]

Volviendo a plantear el sistema de antes con Euler implícito y añadiendo los valores temporales para los cuales queremos estudiar el error de los resultados

{{matlab|codigo=
clear all
clc
%Discretización espacial
a=1;b=3;
h=0.1;
N=(b-a)/h;
r=a:h:b;
rr=r(2:N); %Elementos interiores del vector 
%Definicion de las variables K1 K2
K1=(1/(h^2))*(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1)-diag(ones(1,N-2),1));
K2=(1/(2*h))*(-diag(ones(1,N-2),1)+diag(ones(1,N-2),-1));
G=diag(1./rr);
K2=G*K2;
%Definicion de F y valores de contorno
F=zeros(1,N-1);
Ua=3;Ub=0;
F(1)=F(1)+Ua*((1/h^2)-(1/(2*h*rr(1))));
F(end)=F(end)+Ub*((1/h^2)+(1/(2*h*rr(end))));
%Discretizacion temporal
t0=0;tM=10;
k=h/4;
t=t0:k:tM;
M=length(t)-1;
%Definicion del valor inicial
U0=exp(-8.*((rr-2).^2))+3-rr;
Y=(K1+K2)\F';
Y=[Ua;Y;Ub];
sol(:,1)=U0;
%Solucion
for i=1:M
    sol(:,i+1)=(eye(size(K1))+k*(K1+K2))\(sol(:,i)+k*F');
end
%Grafico
UA=Ua*ones(1,length(t));
UB=Ub*ones(1,length(t));
sol=[UA;sol;UB];
Y1=sol(:,1);
Y2=sol(:,41);
Y3=sol(:,81);
Y4=sol(:,401);
figure(1)
hold on
plot(r,Y,'r');
plot(r,Y1,'g');
plot(r,Y2,'b');
plot(r,Y3,'m');
plot(r,Y4,'k');
legend('Estacionaria','Para t=0','Para t=1','Parat=2','Para t=10');
hold off
error1=abs(Y(2)-Y1(2));
error2=abs(Y(2)-Y2(2));
error3=abs(Y(2)-Y3(2));
error4=abs(Y(2)-Y4(2));
errores=[error1;error2;error3;error4];
}}

El gráfico resultante nos muestra que a medida que avanza el tiempo el error es cada vez mas pequeño 


[[Archivo:Comparacionerrores.jpg|marco|centro|Gráfico de errores]]


Tomando el punto central ρ=2, estos son los valores obtenidos de los errores
              
[[Archivo:Valorerrores.jpg|centro]]

==Planteamiento del problema con una pieza aislante en ρ=1==

En este caso el problema cambia ligeramente y obliga a que la condicion de contorno en el interior del disco sea:

<math> U_ρ(1,t)=0 </math>

Siguiendo el procedimiento y aproximando la derivada parcial:

<math> U_ρ = \frac{U_{n+1}(t)-U_{n-1}(t)}{2h} = 0 </math>

Donde para el radio interior, es decir n=0 la ecuación queda en:

<math> U_1(t)=0  </math>

Se ha utilizado el siguiente código para resolver este problema:

{{matlab|codigo= 
clear all
clc
%Discretización espacial
a=1;b=3;
h=0.1;
N=(b-a)/h;
r=a:h:b;
rr=r(2:N); %Elementos interiores del vector 
%Definicion de las variables K1 K2
K1=(1/(h^2))*(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1)-diag(ones(1,N-2),1));
K1(1,1)=-1;K1(2,1)=0;
K2=(1/(2*h))*(-diag(ones(1,N-2),1)+diag(ones(1,N-2),-1));
K2(2,1)=0;K2(1,1)=1;
G=diag(1./rr);
K2=G*K2;
%Definicion de F y valores de contorno
F=zeros(1,N-1);Ub=0;
%Discretizacion temporal
t0=0;tM=10;
k=h/4;
t=t0:k:tM;
M=length(t)-1;
%Definicion del valor inicial
U0=exp(-8.*((rr-2).^2))+3-rr;
sol(:,1)=U0;
%Solucion
for i=1:M
    sol(:,i+1)=(eye(size(K1))+k*(K1+K2))\(sol(:,i)+k*F');
end
A=zeros(1,length(t));
%Grafico
UB=Ub*ones(1,length(t));
sol=[A;sol;UB];
sol(1,:)=sol(2,:);
sol(2,:)=0;
[Mr,Mt]=meshgrid(r,t);
s=surf(Mr,Mt,sol');
set(s,'edgecolor','flat');
colormap summer
}}

En el grafico resultante se puede observar como la temperatura del radio interior sube a medida que avanza el tiempo debido al aislamiento térmico

[[Archivo:Graficaaislante.jpg|marco|centro|Gráfico para  ρ=1 aislado]]

En cuanto al valor estacionario se puede deducir analíticamente que es U=0.

==Disco de radio 3==

En este apartado consideramos un disco completo de radio 3 que cumple la misma ecuación que el disco anterior con la diferencia de que en este caso la condicion inicial es: <math> u(3,t)=0 </math>

===Problema de autovalores===

La condicion de que la funcion sea diferenciable en el punto ρ=0 nos lleva a la conclusion de que la derivada en ese mismo punto es igual a 0:

<math> u_ρ(0,t)=0 </math>

A partir de este resultado y con la ecuación inicial, tenemos el siguiente sistema:

<math> \left\{\begin{matrix} U_t - U_{ρρ} -\frac{1 }{ρ} U_ρ = 0  \\ U_ρ(0,t)=0 ;  U(3,t)=0 \\ t>0 \end{matrix}\right.</math>

Aplicando el método de separación de variables donde U(ρ,t)=ρT llegamos al desarrollo de dos problemas, uno de contorno en ρ y uno de valor inicial en T. El problema de contorno en ρ queda de la siguiente manera:

<math> \left\{\begin{matrix} ρ'' + \frac{1}{ρ}ρ'+λρ=0 \\ ρ'(0)=0;ρ(3)=0 \end{matrix}\right.</math>

===Bessel===

Tomando la ecuacion anterior y multiplicándola por <math>ρ^2</math> da como resultado una nueva ecuación igual que la que se presenta en el ejemplo del enunciado:


<math> \left\{\begin{matrix} {ρ^2}ρ'' + {ρ}ρ'+{ρ^2}λρ=0 \\ ρ'(0)=0;ρ(3)=0 \end{matrix}\right. </math>

Las autofunciones <math> φ_k = J_0 (ρ{λ}^{1/2}) </math> son soluciones de este problema y se puede demostrar derivando y sustituyendo en el problema que queremos resolver. Los ceros de la funcion de Bessel son los valores para los cuales la funcion es nula. Por lo tanto, en este caso tienen que ser los valores <math> ρ{λ}^{1/2} </math> para los cuales las autofunciones <math> φ_k = J_0 (ρ{λ}^{1/2}) </math> son cero. Tomamos <math> C_0  </math>  a los ceros de estas autofunciones:

<math> C_0 = ρ{λ}^{1/2} </math> 

Ahora, aplicamos la condicion de contorno <math> ρ(3)=0 </math> , de lo que resulta:
 
<math>φ_k = J_0 (3{λ}^{1/2}) = 0 </math> :
<math>C_0 = 3{λ}^{1/2}</math> 

Con el siguiente codigo se pueden obtener los autovalores a partir de los ceros indicados:

{{matlab|codigo=
C=[2.4048,5.5201,8.6537,11.7915,14.9309];
for i=1:length(C)
    L(i)=(C(i)/3)^2;
end
}}

Los resultados son:
<math>λ=[ 0.6426  ,  3.3857  ,  8.3207  , 15.4488  , 24.7702]</math>

===Autofunciones===

Para dibujar las autofunciones creamos el siguiente codigo:

{{matlab|codigo=
clear all
clc
%Cálculo de los autovalores
a=0;b=3;
C=[2.4048,5.5201,8.6537,11.7915,14.9309];
for i=1:length(C)
    L(i)=(C(i)/3)^2;
end
h=0.01;
%Radios
r=a:h:b;
N=(b-a)/h;
G=zeros(N,length(L));
f=zeros(N,length(L));
j0=besselj(0,0);
G(:,1)=j0;
f0=0;
f(:,1)=f0;
for k=1:length(L)
    for j=1:N
        f(j,k)=(r(j))'*sqrt(L(k));
        G(j,k)=besselj(0,f(j,k));
    end
end
%Representacion gráfica
plot(G)
}}

Se obtiene así el siguiente gráfico donde se observan las 5 autofunciones correspondientes a los 5 autovalores obtenidos:

[[Archivo:BESSEL1.jpg|marco|centro|Autofunciones]]

===Fourier===
Para desarrollar este método y sabiendo que la condicion inicial es <math> u(ρ,0)=-log \frac{|ρ+0.1|}{3.1} </math>, ensayaremos soluciones de la forma:

<math> u(ρ,t)= T_k(t)φ_k(ρ) </math>

Para esto se utiliza el código:

{{matlab|codigo=
clear all
clc
%Cálculo de los autovalores
a=0;b=3;
C=[2.4048,5.5201,8.6537,11.7915,14.9309];
for i=1:length(C)
    L(i)=(C(i)/3)^2;
end
h=0.1;
%Radios
r=a:h:b;
N=(b-a)/h;
G=zeros(N+1,length(L));
f=zeros(N+1,length(L));
g0=besselj(0,0);
G(:,1)=g0;
f0=0*r;
f(:,1)=f0;
W=0;
t0=0;tM=30;
m=0.1;
t=t0:m:tM;
[Mr,Mt]=meshgrid(r,t);
for k=1:length(L)
 
        f(:,k)=(r)'*sqrt(L(k));
        G(:,k)=besselj(0,f(:,k));
        h=-log((abs(r+0.1))/3.1);
        e=trapz(r,r'.*h'.*G(:,k))/trapz(r,r'.*G(:,k).*G(:,k));
        W=W+e*exp(-L(k)*t)'*G(:,k)';
 
end
%Representacion gráfica
s=surf(Mr,Mt,W);
set(s,'edgecolor','flat');
colormap summer
}}

El grafico resultante es:

[[Archivo:BESSEL2.jpg|marco|centro|Grafico método Fourier]]




[[Categoría:Ecuaciones Diferentiales]]
[[Categoría:ED14/15]]
[[Categoría:Trabajos 2014-15]]