MateWiki - Contribuciones del usuario [es]

Ecuación del calor. (Grupo A8)

2015-05-17T15:52:45Z

Antonio Carrero Muñoz: /* Cambio de condiciones de contorno */

{{ TrabajoED |Ecuación del calor (Grupo A8)| [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED14/15|Curso 2014-15]] | Valentina Salazar; Antonio Carrero; José Francisco Aguilera }}
=Introducción y modelización del problema.=
En este artículo trataremos el comportamiento térmico de una varilla sometida a ciertas condiciones térmicas y físicas mediante diferentes métodos numéricos, daremos interpretación física a los resultados arrojados por estos métodos. Basaremos el estudio analítico y numérico necesario en la ecuación del calor propuesta por Jean-Baptiste Joseph Fourier en 1822.

Para estudiar el problema consideraremos una varilla delgada, de sección constante y de un material homogéneo, de longitud <math>L=4</math>. La situaremos en el intervalo <math>x\in{(0,4)}</math> de la recta real. Para llevar acabo esta modelización tendremos que asumir unas ciertas hipótesis y simplificaciones:

*Para el primer caso a plantear la varilla estará infinitamente aislada del entorno y por tanto no habrá flujo de calor sobre la superficie lateral de la misma.
*Al estar considerando el caso unidimensional de la ecuación del calor, asumiremos que al ser una varilla delgada la temperatura a lo largo de una sección ortogonal al eje <math>x</math> se mantiene constante en toda la sección.
*Consideraremos que el calor específico del material, <math>c</math>, es constante y no depende de la temperatura, y por tanto la difusividad térmica,<math>\alpha=\frac{k}{c\rho}</math> también lo será.

Si damos un valor a las constantes <math>c</math>, <math>k</math> y <math>\rho</math>, tal que <math>\alpha=2</math>; y no hay fuentes ni sumideros de calor dentro de la varilla, la ecuación en derivadas parciales que tendrá que cumplir la distribución de temperaturas a lo largo de la misma es::

<math>u_t(x,t)-2u_{xx}(x,t)=0</math>

Donde:
*<math>u(x,t)</math> nos define como evoluciona la temperatura a lo largo del eje <math>x</math>.
*<math>u_t(x,t)</math> es la derivada de la temperatura con respecto del tiempo.
*<math>u_{xx}(x,t)</math> es la segunda derivada de la temperatura con respecto de <math>x</math>.

Junto con está ecuación en derivadas parciales consideraremos una serie de condiciones: iniciales y de contorno. Como condición inicial para está modelización propondremos::

<math>u(x,0)=e^{-6(x-2)^2}+5-\frac{5}{4}x</math>

Donde <math>u(x,0)</math> no es otra cosa que la distribución de temperaturas inicial.

A lo largo del artículo variaremos las condiciones de contorno para dar diferentes ejemplos, pero en el bloque de casos prácticos trabajaremos con las siguientes::

<math>u(0,t)=5</math>:
<math>u(4,t)=0</math>

Éstas indican la temperatura (en este caso, otras pueden ser los flujos de calor o combinaciones de ambos) en cada uno de los extremos de la varilla.

=Casos prácticos=
En los apartados que siguen estudiaremos a través del cálculo numérico el caso anteriormente expuesto e interpretaremos los resultados obtenidos mediante dicho análisis. El lenguaje usado en el que se implementarán estos modelos será código Matlab u Octave.
==Método de diferencias finitas o método de lineas==
Trataremos de reducir el sistema de ecuaciones a una sola variable en una discretización de <math>n</math> valores (siendo <math>h</math> el tamaño de los subintervalos) haciendo la aproximación::

<math>u_{xx}(x,t)=\frac{-U^{N-1}(t)+2U^N(t)-U^{N+1}(t)}{h^2}</math>

Así podemos aplicar uno de los diferentes métodos iterativos en la variable <math>t</math> que se muestran en los subapartados que siguen.

===Método del trapecio===
Se trata de un método implícito, es decir que hay que despejar en la ecuación matricial que nos da la <math>U(t)</math>. A continuación se muestra el código que lo implementa en Matlab u Octave, dandonos una gráfica en 3D del comportamiento de la temperatura frente al tiempo en la varilla y otro de la evolución temporal de temperaturas del punto medio de la misma.

{{matlab|codigo=clear all
%datos del problema
L=4;
h=0.1;
T=8;
%discretizacion espacial
N=L/h;
%vector de puntos en el espacio
x=0:h:L;
xint=h:h:L-h;%nodos interiores
%aproximacion de -Uxx
K=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
K=(2/h^2)*K;
%calculamos F
F=zeros(N-1,1);
F(1)=F(1)+2*5/h^2;
%calculamos U0
U0=((exp((-6)*(xint-2).^2))+(5-5/4.*xint))';
%discretizacion temporal
dt=h; %paso en el espacio
t=0:dt:T;%vector en tiempos
%metodo del trapecio
I=eye(N-1);
M=I+(dt./2)*K;
sol(:,1)=U0;
for n=1:length(t)-1
sol(:,n+1)=M\(sol(:,n)+dt/2*(F)+dt/2*(-K*sol(:,n)+F));
end
UA=5*ones(1,length(t));
UB=0*ones(1,length(t));
sol=[UA;sol;UB];;
%dibujamos la solucion (apartado 3)

[xx,tt]=meshgrid(x,t);

figure(1)
surf(xx,tt,sol')
xlabel('espacio')
ylabel('tiempo')
zlabel('temperatura')
figure(2)
plot(t,sol(20,:))
xlabel('tiempo')
ylabel('temperatura')}}

Obteniendo los siguientes gráficos:

[[Archivo:Dftrap.png|marco|centro|Evolución de la temperatura en la varilla frente al tiempo y la posición de los puntos de la misma.]]

Como podemos ver la temperatura apenas tarda en estabilizarse en el tiempo a lo largo de los puntos de la varilla, esto nos sugiere que una buena aproximación será la solución estacionaria de la ecuación.

[[Archivo:Dftrappm.png|marco|centro|Sucesivas temperaturas del punto medio frente al tiempo.]]

La temperatura del punto medio descenderá hasta volverse casi asintótica hacia un valor aproximado de 2.625.

===Método de Euler explícito===
A diferencia del método del trapecio, como su nombre indica, este método es explicito, lo que es una ventaja a la hora de implementar ya que nos ahorramos el despejar; pero sin embargo se trata de un método bastante inestable, y para que arroje soluciones lógicas hemos de disminuir el tamaño del paso temporal frente al espacial del orden de <math>dt=\frac{h^2}{4}</math>, siendo <math>h</math>, <math>dt</math> dichos pasos espacial y temporal respectivamente. A continuación se muestra el código que lo implementa en Matlab u Octave, dandonos una gráfica en 3D del comportamiento de la temperatura frente al tiempo en la varilla y otro de la evolución temporal de temperaturas del punto medio de la misma.

{{matlab|codigo=clear all
%datos del problema
L=4;
h=0.1;
T=8;
%discretizacion espacial
N=L/h;
%vector de puntos en el espacio
x=0:h:L;
xint=h:h:L-h;%nodos interiores
%aproximacion de -Uxx
K=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
K=(2/h^2)*K;
%calculamos F
F=zeros(N-1,1);
F(1)=F(1)+2*5/h^2;
%calculamos U0
U0=((exp((-6)*(xint-2).^2))+(5-5/4.*xint))';
%discretizacion temporal
dt=h^2/4; %paso en el espacio
t=0:dt:T;%vector en tiempos
%metodo del trapecio
I=eye(N-1);
M=I-dt*K;
sol(:,1)=U0;
for n=1:length(t)-1
sol(:,n+1)=M*sol(:,n)+dt*F;
end
UA=5*ones(1,length(t));
UB=0*ones(1,length(t));
sol=[UA;sol;UB];
%dibujamos la solucion

[xx,tt]=meshgrid(x,t);

figure(1)
mesh(xx,tt,sol')
xlabel('espacio')
ylabel('tiempo')
zlabel('Temperatura')
figure(2)
plot(t,sol(20,:))
xlabel('tiempo')
ylabel('temperatura')}}

Obteniendo los gráficos:

[[Archivo:Dfeulerexp.png|marco|centro|Evolución de la temperatura en la varila frente al tiempo y la posición de los puntos de la misma.]]

Al igual que en el apartado anterior vemos que se estabiliza muy rápido. Debido al tener que disminuir bastante el tamaño del paso temporal hemos tenido que cambiar la función de Octave para representar este gráfico.

[[Archivo:Dfeulerexppm.png|marco|centro|Sucesivas temperaturas del punto medio frente al tiempo.]]

Vemos que en esta implementación los cambios de temperatura en el punto medio son más ''continuos'', pero también acaba siendo asintótica en un valor aproximado de 2.625.
===Método de Euler implícito===
Es más fácil hacerlo converger que su homónimo explícito, pero tiene la desventaja de tener que despejar en la ecuación para implementarlo. A continuación se muestra el código que lo implementa en Matlab u Octave, dándonos una gráfica en 3D del comportamiento de la temperatura frente al tiempo en la varilla y otro de la evolución temporal de temperaturas del punto medio de la misma.

{{matlab|codigo=clear all
%datos del problema
L=4;
h=0.1;
T=8;
%discretizacion espacial
N=L/h;
%vector de puntos en el espacio
x=0:h:L;
xint=h:h:L-h;%nodos interiores
%aproximacion de -Uxx
K=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
K=(2/h^2)*K;
%calculamos F
F=zeros(N-1,1);
F(1)=F(1)+2*5/h^2;
%calculamos U0
U0=((exp((-6)*(xint-2).^2))+(5-5/4.*xint))';
%discretizacion temporal
dt=h; %paso en el espacio
t=0:dt:T;%vector en tiempos
%metodo del trapecio
I=eye(N-1);
M=I+dt*K;
sol(:,1)=U0;
for n=1:length(t)-1
sol(:,n+1)=M\(sol(:,n)+dt*(F));
end
UA=5*ones(1,length(t));
UB=0*ones(1,length(t));
sol=[UA;sol;UB];;
%dibujamos la solucion

[xx,tt]=meshgrid(x,t);

figure(1)
surf(xx,tt,sol')
xlabel('espacio')
ylabel('tiempo')
zlabel('temperatura')
figure(2)
plot(t,sol(20,:))
xlabel('tiempo')
ylabel('temperatura')}}

Obteniendo los siguientes gráficos:

[[Archivo:Dfeulerimp.png|marco|centro|Evolución de la temperatura en la varilla frente al tiempo y la posición de los puntos de la misma.]]

Volvemos a ver que la estabilización de temperaturas en la varilla es bastante rápida.

[[Archivo:Dfeulerimppm.png|marco|centro|Sucesivas temperaturas del punto medio frente al tiempo.]]

Una vez más observamos un descenso asintótico hacia el valor de temperatura 2.625.
===Método de Runge-Kutta===
Se trata de un método explícito, por tanto necesitaremos aumentar el tamaño del paso una vez más hasta los mismos valores que para Euler explícito para asegurar su convergencia. A continuación se muestra el código que lo implementa en Matlab u Octave, dandonos una gráfica en 3D del comportamiento de la temperatura frente al tiempo en la varilla y otro de la evolución temporal de temperaturas del punto medio de la misma.

{{matlab|codigo=clear all
%datos del problema
L=4;
h=0.1;
T=8;
%discretizacion espacial
N=L/h;
%vector de puntos en el espacio
x=0:h:L;
xint=h:h:L-h;%nodos interiores
%aproximacion de -Uxx
K=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
K=(2/h^2)*K;
%calculamos F
F=zeros(N-1,1);
F(1)=F(1)+2*5/h^2;
%calculamos U0
U0=((exp((-6)*(xint-2).^2))+(5-5/4.*xint))';
%discretizacion temporal
dt=h^2/4; %paso en el espacio
t=0:dt:T;%vector en tiempos
%metodo del trapecio
I=eye(N-1);
M=I-dt.*K;
sol(:,1)=U0;
for n=1:length(t)-1
k1=-K*sol(:,n)+F;
k2=-K*(sol(:,n)+k1*dt/2)+F;
k3=-K*(sol(:,n)+k2*dt/2)+F;
k4=-K*(sol(:,n)+k3*dt)+F;
sol(:,n+1)=sol(:,n)+(dt/6)*(k1+2*k2+2*k3+k4);
end
UA=5*ones(1,length(t));
UB=0*ones(1,length(t));
sol=[UA;sol;UB];
%dibujamos la solucion

[xx,tt]=meshgrid(x,t);

figure(1)
mesh(xx,tt,sol')
xlabel('espacio')
ylabel('tiempo')
zlabel('temperatura')
figure(2)
plot(t,sol(20,:))
xlabel('tiempo')
ylabel('temperatura')}}

Obteniendo las siguientes gráficas:

[[Archivo:Dfrk.png|marco|centro|Evolución de la temperatura en la varilla frente al tiempo y la posición de los puntos de la misma.]]

Al igual que en todos los métodos anteriores se produce la ya esperada estabilización de temperaturas en la varilla de manera bastante acelerada.

[[Archivo:Dfrkpm.png|marco|centro|Sucesivas temperaturas del punto medio frente al tiempo.]]

De nuevo aparece un descenso asintótico hasta el valor 2.625.

==Método de Fourier==
Partiendo de una solución analítica previa obtenemos la solución en forma de sumatorio de infinitos términos, reduciendo este número de términos numéricamente calcularemos una solución aproximada. Los autovalores y autovectores obtenidos para este problema y con condiciones de contorno son::

<math>\mu_k=\frac{k^2\pi^2}{16}</math>:

<math>\varphi_k(x)=sen(\frac{k\pi x}{4})</math>:

<math>k=1,2,3,...</math>

Dados estos valores y la función de homogenización para las condiciones de contorno::

<math>v(x,t)=5-\frac{5}{4}x</math>

Mediante el siguiente código de Matlab implementamos el problema numéricamente, obteniendo así la gráfica en 3D de la solución para <math>k=20</math> (por ser la más representativa y exacta),y la comparación de las gráficas en el plano de <math>t=0,5</math> para <math>k=1,3,5,10 y 20</math>.

{{matlab|codigo=clear all
N=100;
a=0; b=4;
h=(b-a)/N; x=0:h:4; % space discretization
ht=h; t=0:ht:5; % time discretization
[xx,tt]=meshgrid(x,t); % time discretization
f=(exp((-6)*(x-2).^2)); % Initial data
aprox=0; % Compute the approximation
q=1;
for k=1:q
autovalor=((k*pi)^2)/16; % eigenvalue
autofuncion=sin(k*pi*x/4); % eigenfunction
c(k)=(trapz(x,f.*autofuncion))/trapz(x,autofuncion.^2); % Fourier coef.
aprox=aprox+c(k)*exp(-2*autovalor.*tt).*sin(k*pi*xx/4);
end
aprox2=aprox+5-(5/4)*xx;
figure(1)
plot(x,aprox2(0.5/ht-0.5,:)) % Draw the approximation
hold on
q=3;
for k=1:q
autovalor=((k*pi)^2)/16; % eigenvalue
autofuncion=sin(k*pi*x/4); % eigenfunction
c(k)=(trapz(x,f.*autofuncion))/trapz(x,autofuncion.^2); % Fourier coef.
aprox=aprox+c(k)*exp(-2*autovalor.*tt).*sin(k*pi*xx/4);
end
aprox2=aprox+5-(5/4)*xx;
figure(1)
plot(x,aprox2(0.5/ht-0.5,:),'r') % Draw the approximation
hold on
q=5;
for k=1:q
autovalor=((k*pi)^2)/16; % eigenvalue
autofuncion=sin(k*pi*x/4); % eigenfunction
c(k)=(trapz(x,f.*autofuncion))/trapz(x,autofuncion.^2); % Fourier coef.
aprox=aprox+c(k)*exp(-2*autovalor.*tt).*sin(k*pi*xx/4);
end
aprox2=aprox+5-(5/4)*xx;
figure(1)
plot(x,aprox2(0.5/ht-0.5,:),'m') % Draw the approximation
hold on
q=10;
for k=1:q
autovalor=((k*pi)^2)/16; % eigenvalue
autofuncion=sin(k*pi*x/4); % eigenfunction
c(k)=(trapz(x,f.*autofuncion))/trapz(x,autofuncion.^2); % Fourier coef.
aprox=aprox+c(k)*exp(-2*autovalor.*tt).*sin(k*pi*xx/4);
end
aprox2=aprox+5-(5/4)*xx;
figure(1)
plot(x,aprox2(0.5/ht-0.5,:),'g') % Draw the approximation
hold on
q=20;
for k=1:q
autovalor=((k*pi)^2)/16; % eigenvalue
autofuncion=sin(k*pi*x/4); % eigenfunction
c(k)=(trapz(x,f.*autofuncion))/trapz(x,autofuncion.^2); % Fourier coef.
aprox=aprox+c(k)*exp(-2*autovalor.*tt).*sin(k*pi*xx/4);
end
aprox2=aprox+5-(5/4)*xx;
figure(1)
plot(x,aprox2(0.5/ht-0.5,:),'y') % Draw the approximation
hold on
plot(x,5-5/4*x,'k')
figure(2)
mesh(xx,tt,aprox2)
hold on}}

Obteniendo así:

[[Archivo:MF.png|marco|centro|Gráfico que muestra de la temperatura de la varilla frente a al tiempo a lo largo de sus diferentes posiciones.]]

Se puede ver como a diferencia del método de lineas el error viene inducido por la aproximación mediante series de Fourier de la condición inicial. Por lo demás evoluciona de manera similar a los métodos aplicados anteriormente.

[[Archivo:MF20.png|marco|centro|Gráficas de las temperaturas a lo largo de la varilla con aproximaciones de diferentes índices del sumatorio en la solución. k=1(azul). k=3(rojo). k=5(magenta). k=10(verde). k=20(amarillo). solución estacionaria(negro)]]

En la gráfica anterior se puede apreciar el nivel de aproximación con relación al número de coeficientes que utilizamos de la serie. Vemos como cuantos menos coeficientes usamos más nos acercamos incluso en tiempo cortos a la solución estacionaria, lo que nos indica que desciende el nivel de precisión de estas aproximaciones.

=Cambio de condiciones de contorno=
Ahora manteniendo la ecuación que nos modelizaba el comportamiento térmico de la varilla, cambiamos las condiciones de contorno en los extremos. Las nuevas condiciones propuestas se pueden asimilar como el aislamiento térmico total del extremo izquierdo, lo que anula el flujo de calor en ese lado. Estas son:

<math>u_x(0,t)=0</math>:
<math>u(4,t)=0</math>

Para implementarlo en sus nuevas condiciones usaremos el método de lineas con el método del trapecio. El código que implementa esta modelización en lenguaje Matlab y nos proporciona un gráfico de la temperatura frente a las variables espacial y temporal es:

{{matlab|codigo=clear all
%datos del problema
L=4;
h=0.1;
T=8;
%discretizacion espacial
N=L/h;
%vector de puntos en el espacio
x=0:h:L;
xint=0:h:L-h;%nodos interiores
%aproximacion de -Uxx
K=2*diag(ones(1,N))-diag(ones(1,N-1),1)-diag(ones(1,N-1),-1);
K=(2/h^2)*K;
%calculamos F
F=zeros(N,1);
F(1)=F(1)+2*5/h^2;
%calculamos U0
U0=((exp((-6)*(xint-2).^2))+(5-5/4.*xint))';
%discretizacion temporal
dt=h; %paso en el espacio
t=0:dt:T;%vector en tiempos
%metodo del trapecio
I=eye(N);
M=I+(dt./2)*K;
sol(:,1)=U0;
for n=1:length(t)-1
sol(:,n+1)=M\(sol(:,n)+dt/2*(F)+dt/2*(-K*sol(:,n)+F));
end
UB=zeros(length(t),1);
sol=[sol;UB'];
%dibujamos la solucion (apartado 3)
[xx,tt]=meshgrid(x,t);
figure(1)
surf(xx,tt,sol')
xlabel('espacio')
ylabel('tiempo')
zlabel('Temperatura')}}

Del cual obtenemos:

[[Archivo:Dfneu.png|marco|centro|Progreso de la temperatura frente a tiempo y espacio con las nuevas condiciones.]]

Como se puede ver en la imagen anterior la temperatura de los puntos interiores de la varilla va descendiendo lenta y progresivamente, a diferencia de en el caso anterior donde esta tendía a homogeneizarse.

=Interpretaciones estacionarias=
Como hemos visto en los apartados anteriores de este artículo, las soluciones por norma general de esta ecuación tienden a volverse estacionarias en unos intervalos de tiempo bastante breves. Por ello en esta sección estudiaremos estas soluciones estacionarias y el error que se comete asumiéndolas como ciertas en tiempos relativamente cortos.

Para calcular analiticamente estas soluciones estacionarias basta con considerar nulo el término en <math>u_t(x,t)</math> de la ecuación en derivadas parciales, y resolver lo que ahora nos queda un problema de contorno con dichas condiciones y la ecuación diferencial ordinaria que nos queda de hacer la despreciación anterior.

Tratando con las condiciones de contorno dadas inicialmente::

<math>u(0,t)=5</math>:
<math>u(4,t)=0</math>

Analíticamente para estas condiciones la solución estacionaria hallada es::

<math>u_{est}(x)=5-\frac{5}{4}x</math>

Mediante cálculo numérico y el siguiente código de Octave, conseguimos una gráfica de los diferentes valores de la temperatura en ciertos tiempos que tendrán que ir acercándose a la solución estacionaria, y un gráfico del error cometido si aproximásemos en estos tiempos el valor numérico aproximado de la temperatura por el estacionario calculado analíticamente. Además, nos aporta una tabla con en la primera fila el tiempo para el que se calcula el error respecto de la estacionaria y en las sucesivas filar los errores correspondientes.

{{matlab|codigo=clear all
%datos del problema
L=4;
h=0.1;
T=10;
%discretizacion espacial
N=L/h;
%vector de puntos en el espacio
x=0:h:L;
xint=h:h:L-h;%nodos interiores
%aproximacion de -Uxx
K=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
K=(2/h^2)*K;
%calculamos F
F=zeros(N-1,1);
F(1)=F(1)+2*5/h^2;
%calculamos U0
U0=((exp((-6)*(xint-2).^2))+(5-5/4.*xint))';
%discretizacion temporal
dt=T/L*h; %paso en el espacio
t=0:dt:T;%vector en tiempos
%metodo del trapecio
I=eye(N-1);
M=I+(dt./2)*K;
sol(:,1)=U0;
for n=1:length(t)-1
sol(:,n+1)=M\(sol(:,n)+dt/2*(F)+dt/2*(-K*sol(:,n)+F));
end
UA=5*ones(1,length(t));
UB=0*ones(1,length(t));
sol=[UA;sol;UB];;
%dibujamos la solucion (apartado 3)

[xx,tt]=meshgrid(x,t);

%Posiciones
tiempos=[0,1,2,10];
posicion=tiempos/dt;
figure(2)
%dibujamos para ciertos tiempos
for m=1:length(posicion)
hold on
plot(x,sol(:,posicion(m)+1),'m');
end
xlabel('espacio')
ylabel('temperaturas en cada t')
%Estacionaria
hold on
est=5*(1-x/4);
plot(x,est);
figure(3)
%errores
for p=1:length(posicion)
error=sol(:,posicion(p)+1).-est';
Error(:,p)=error;
hold on
plot(x,error)
end
xlabel('espacio')
ylabel('error cometido')
[0,1,2,10;Error]}}

Obteniendo:

[[Archivo:Estacionaria1.png|marco|centro|Contraste de la gráfica de la estacionaria con las soluciones calculadas mediante numérico en diferentes tiempos.]]

Se puede observar como se reduce el error a medida que aumentamos el tiempo de modelización cor respecto de la estacionaria.

[[Archivo:Estacionaria2.png|marco|centro|Error cometido al aproximar en diferentes tiempos por la estacionaria.]]

En la imagen anterior se puede ver una comparativa del el error con respecto al eje x, donde se ve más claramente que arriba como se reduce.

[[Archivo:Tablaerror.png|marco|centro|Tabla que nos muestra el error cometido en los diferentes tiempos.]]

En los valores numéricos dispuestos en la enterior tabla ya se muestra con total seguridad como se reduce el error a medida que aumentamos el tiempo.

Si ahora por el contrario tratamos las condiciones propuestas en el apartado de cambio de condiciones de contorno, las cuales son:

<math>u_x(0,t)=0</math>:
<math>u(4,t)=0</math>

Analíticamente para estas condiciones la solución estacionaria hallada es::

<math>u_{est}(x)=0</math>

Lo que al principio nos puede chocar al ver la gráfica, ya que en las otras condiciones se llega al estado estacionario de manera mucho más rápida. Si consideramos como estado estacionario que la diferencia entre dos valores consecutivos sea menor de 0.05, mediante el siguiente programa do Octave podemos calcular este valor de tiempo:

{{matlab|codigo=clear all
%datos del problema
L=4;
h=0.1;
T=8;
%discretizacion espacial
N=L/h;
%vector de puntos en el espacio
x=0:h:L;
xint=0:h:L-h;%nodos interiores
%aproximacion de -Uxx
K=2*diag(ones(1,N))-diag(ones(1,N-1),1)-diag(ones(1,N-1),-1);
K=(2/h^2)*K;
%calculamos F
F=zeros(N,1);
F(1)=F(1)+2*5/h^2;
%calculamos U0
U0=((exp((-6)*(xint-2).^2))+(5-5/4.*xint))';
%discretizacion temporal
dt=h; %paso en el espacio
t=0:dt:T;%vector en tiempos
%metodo del trapecio
I=eye(N);
M=I+(dt./2)*K;
sol(:,1)=U0;
for n=1:length(t)-1
sol(:,n+1)=M\(sol(:,n)+dt/2*(F)+dt/2*(-K*sol(:,n)+F));
end
UB=zeros(length(t),1);
sol=[sol;UB'];for i=1:length(t)
er=abs(sol(:,i+1)-sol(:,i));
a=max(er);
if a<=0.05
cont=i;
tiempo=(cont-1)*dt;
break
end
end
tiempo}}

El programa nos dice que el tiempo en el que se cumple esta condición es <math>t=0.7</math>.

=Intercambio de calor con el entorno. Varilla no aislada longitudinalmente=
Si ahora panteamos el intercambio de calor con el exterior, suponiendo que este ambiente se encuentra a una temperatura de 16ºC, la ecuación diferencial se ve afectada de manera que nos queda::

<math>u_t(x,t)-2u_{xx}(x,t)+u(x,t)-16=0</math>

Implementando numéricamente esta ecuación con las condiciones iniciales propuestas en el primer apartado de este artículo, mediante el método de lineas y el método del trapecio, mediante el siguiente código, obteniendo así la gráfica de la temperatura en cada punto de la varilla con respecto al tiempo.

{{matlab|codigo=clear all
%datos del problema
L=4;
h=0.1;
T=10;
%discretizacion espacial
N=L/h;
%vector de puntos en el espacio
x=0:h:L;
xint=h:h:L-h;%nodos interiores
%aproximacion de -Uxx
K=(2+h^2)*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
K=(1/h^2)*K;
%calculamos F
F=zeros(N-1,1).+16;;
F(1)=F(1)+2*5/h^2;
%calculamos U0
U0=((exp((-6)*(xint-2).^2))+(5-5/4.*xint))';
%discretizacion temporal
dt=h; %paso en el espacio
t=0:dt:T;%vector en tiempos
%metodo del trapecio
I=eye(N-1);
M=I+(dt./2)*K;
sol(:,1)=U0;
for n=1:length(t)-1
sol(:,n+1)=M\(sol(:,n)+dt/2*(F)+dt/2*(-K*sol(:,n)+F));
end
UA=5*ones(1,length(t));
UB=0*ones(1,length(t));
sol=[UA;sol;UB];;
%dibujamos la solucion (apartado 3)

[xx,tt]=meshgrid(x,t);

figure(1)
surf(xx,tt,sol')
xlabel('espacio')
ylabel('tiempo')
zlabel('Temperatura')}}

Obteniendo de este código:

[[Archivo:Perdidasaire.png|marco|centro|Representación de la temperatura en cada punto de la varilla con respecto al tiempo]]

Donde se ve como el exterior aporta calor a la varilla elevando sus temperaturas frente a las iniciales.

Si considerásemos el caso estacionario, se puede calcular mediante la inclusión del siguiente fragmento de código en el anterior programa que modelizaba el fenómeno no estacionario. Además de darnos la situación de temperaturas estacionarias, nos da el tiempo en el que este se alcanza con un error relativo de <math>10^{-3}</math>.

{{matlab|codigo=%Estacionaria
figure(2)
xlabel('espacio')
ylabel('temperatura')
c1=det([-11,1;-16,exp(-4)])/det([1,1;exp(4),exp(-4)]);
c2=det([1,-11;exp(4),-16])/det([1,1;exp(4),exp(-4)]);
est=c1*exp(x)+c2*exp(-x)+16;
plot(x,est);
hold on
for i=1:length(t)-1
er=abs(sol(:,i+1)-sol(:,i));
a=max(er);
if a<=0.001
cont=i;
ti=(cont-1)*dt;
break
end
end
ti
cont}}

La gráfica conseguida del código anterior es:

[[Archivo:Perdidasaireest.png|marco|centro|Temperatura estacionaria de la varilla con respecto a la posición de cada punto de ella.]]

El valor de tiempo obtenido es <math>t=5.9</math>.

==Cambio de condiciones de contorno==
Si cambiamos las condiciones de contorno de la ecuación por estas:

<math>u_x(0,t)=2</math>:
<math>u(4,t)=-2</math>

Físicamente podrían interpretarse como la conexión de un deposito de fluido en el extremo derecho que obligara un flujo constante de calor en el mismo extremo, y que el extremo izquierdo estuviera posicionado sobre un objeto que le obligara a tener una temperatura constante.

Si implementamos numéricamente este nuevo modelo propuesto mediante el siguiente código Matlab, obtenemos la representación de las temperaturas de la varilla frente al tiempo.

{{matlab|codigo=clear all
%datos del problema
L=4;
h=0.1;
T=10;
%discretizacion espacial
N=L/h;
%vector de puntos en el espacio
x=0:h:L;
xint=h:h:L;%nodos interiores
%aproximacion de -Uxx
K=(2+h^2)*diag(ones(1,N))-diag(ones(1,N-1),1)-diag(ones(1,N-1),-1);
K(1,1)=K(1,1)-1;
K=(1/h^2)*K;
%calculamos F
F=zeros(N,1).+16;;
F(1)=F(1)+2/h;
F(N)=F(N)+2/h^2;
%calculamos U0
U0=((exp((-6)*(xint-2).^2))+(5-5/4.*xint))';
%discretizacion temporal
dt=h; %paso en el espacio
t=0:dt:T;%vector en tiempos
%metodo del trapecio
I=eye(N);
M=I+(dt./2)*K;
sol(:,1)=U0;
for n=1:length(t)-1
sol(:,n+1)=M\(sol(:,n)+dt/2*(F)+dt/2*(-K*sol(:,n)+F));
end
UB=2*ones(1,length(t));
sol=[sol;UB];;
%dibujamos la solucion (apartado 3)

[xx,tt]=meshgrid(x,t);

figure(1)
surf(xx,tt,sol')
xlabel('espacio')
ylabel('tiempo')
zlabel('Temperatura')

}}

Logrando así la gráfica:

[[Archivo:Perdidasneu111.png|marco|centro|Representación de las temperaturas de la varilla frente al tiempo]]

Se puede observar como la temperatura de los puntos de la varilla va aumentado debido al flujo entrante por el lado derecho de esta hasta que alcanza un estado estacionario.

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

Ecuación del calor. (Grupo A8)

2015-05-15T19:53:33Z

Antonio Carrero Muñoz: /* Interpretaciones estacionarias */

{{ TrabajoED |Ecuación del calor (Grupo A8)| [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED14/15|Curso 2014-15]] | Valentina Salazar; Antonio Carrero; José Francisco Aguilera }}
=Introducción y modelización del problema.=
En este artículo trataremos el comportamiento térmico de una varilla sometida a ciertas condiciones térmicas y físicas mediante diferentes métodos numéricos, daremos interpretación física a los resultados arrojados por estos métodos. Basaremos el estudio analítico y numérico necesario en la ecuación del calor propuesta por Jean-Baptiste Joseph Fourier en 1822.

Para estudiar el problema consideraremos una varilla delgada, de sección constante y de un material homogéneo, de longitud <math>L=4</math>. La situaremos en el intervalo <math>x\in{(0,4)}</math> de la recta real. Para llevar acabo esta modelización tendremos que asumir unas ciertas hipótesis y simplificaciones:

*Para el primer caso a plantear la varilla estará infinitamente aislada del entorno y por tanto no habrá flujo de calor sobre la superficie lateral de la misma.
*Al estar considerando el caso unidimensional de la ecuación del calor, asumiremos que al ser una varilla delgada la temperatura a lo largo de una sección ortogonal al eje <math>x</math> se mantiene constante en toda la sección.
*Consideraremos que el calor específico del material, <math>c</math>, es constante y no depende de la temperatura, y por tanto la difusividad térmica,<math>\alpha=\frac{k}{c\rho}</math> también lo será.

Si damos un valor a las constantes <math>c</math>, <math>k</math> y <math>\rho</math>, tal que <math>\alpha=2</math>; y no hay fuentes ni sumideros de calor dentro de la varilla, la ecuación en derivadas parciales que tendrá que cumplir la distribución de temperaturas a lo largo de la misma es::

<math>u_t(x,t)-2u_{xx}(x,t)=0</math>

Donde:
*<math>u(x,t)</math> nos define como evoluciona la temperatura a lo largo del eje <math>x</math>.
*<math>u_t(x,t)</math> es la derivada de la temperatura con respecto del tiempo.
*<math>u_{xx}(x,t)</math> es la segunda derivada de la temperatura con respecto de <math>x</math>.

Junto con está ecuación en derivadas parciales consideraremos una serie de condiciones: iniciales y de contorno. Como condición inicial para está modelización propondremos::

<math>u(x,0)=e^{-6(x-2)^2}+5-\frac{5}{4}x</math>

Donde <math>u(x,0)</math> no es otra cosa que la distribución de temperaturas inicial.

A lo largo del artículo variaremos las condiciones de contorno para dar diferentes ejemplos, pero en el bloque de casos prácticos trabajaremos con las siguientes::

<math>u(0,t)=5</math>:
<math>u(4,t)=0</math>

Éstas indican la temperatura (en este caso, otras pueden ser los flujos de calor o combinaciones de ambos) en cada uno de los extremos de la varilla.

=Casos prácticos=
En los apartados que siguen estudiaremos a través del cálculo numérico el caso anteriormente expuesto e interpretaremos los resultados obtenidos mediante dicho análisis. El lenguaje usado en el que se implementarán estos modelos será código Matlab u Octave.
==Método de diferencias finitas o método de lineas==
Trataremos de reducir el sistema de ecuaciones a una sola variable en una discretización de <math>n</math> valores (siendo <math>h</math> el tamaño de los subintervalos) haciendo la aproximación::

<math>u_{xx}(x,t)=\frac{-U^{N-1}(t)+2U^N(t)-U^{N+1}(t)}{h^2}</math>

Así podemos aplicar uno de los diferentes métodos iterativos en la variable <math>t</math> que se muestran en los subapartados que siguen.

===Método del trapecio===
Se trata de un método implícito, es decir que hay que despejar en la ecuación matricial que nos da la <math>U(t)</math>. A continuación se muestra el código que lo implementa en Matlab u Octave, dandonos una gráfica en 3D del comportamiento de la temperatura frente al tiempo en la varilla y otro de la evolución temporal de temperaturas del punto medio de la misma.

{{matlab|codigo=clear all
%datos del problema
L=4;
h=0.1;
T=8;
%discretizacion espacial
N=L/h;
%vector de puntos en el espacio
x=0:h:L;
xint=h:h:L-h;%nodos interiores
%aproximacion de -Uxx
K=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
K=(2/h^2)*K;
%calculamos F
F=zeros(N-1,1);
F(1)=F(1)+2*5/h^2;
%calculamos U0
U0=((exp((-6)*(xint-2).^2))+(5-5/4.*xint))';
%discretizacion temporal
dt=h; %paso en el espacio
t=0:dt:T;%vector en tiempos
%metodo del trapecio
I=eye(N-1);
M=I+(dt./2)*K;
sol(:,1)=U0;
for n=1:length(t)-1
sol(:,n+1)=M\(sol(:,n)+dt/2*(F)+dt/2*(-K*sol(:,n)+F));
end
UA=5*ones(1,length(t));
UB=0*ones(1,length(t));
sol=[UA;sol;UB];;
%dibujamos la solucion (apartado 3)

[xx,tt]=meshgrid(x,t);

figure(1)
surf(xx,tt,sol')
xlabel('espacio')
ylabel('tiempo')
zlabel('temperatura')
figure(2)
plot(t,sol(20,:))
xlabel('tiempo')
ylabel('temperatura')}}

Obteniendo los siguientes gráficos:

[[Archivo:Dftrap.png|marco|centro|Evolución de la temperatura en la varilla frente al tiempo y la posición de los puntos de la misma.]]

Como podemos ver la temperatura apenas tarda en estabilizarse en el tiempo a lo largo de los puntos de la varilla, esto nos sugiere que una buena aproximación será la solución estacionaria de la ecuación.

[[Archivo:Dftrappm.png|marco|centro|Sucesivas temperaturas del punto medio frente al tiempo.]]

La temperatura del punto medio descenderá hasta volverse casi asintótica hacia un valor aproximado de 2.625.

===Método de Euler explícito===
A diferencia del método del trapecio, como su nombre indica, este método es explicito, lo que es una ventaja a la hora de implementar ya que nos ahorramos el despejar; pero sin embargo se trata de un método bastante inestable, y para que arroje soluciones lógicas hemos de disminuir el tamaño del paso temporal frente al espacial del orden de <math>dt=\frac{h^2}{4}</math>, siendo <math>h</math>, <math>dt</math> dichos pasos espacial y temporal respectivamente. A continuación se muestra el código que lo implementa en Matlab u Octave, dandonos una gráfica en 3D del comportamiento de la temperatura frente al tiempo en la varilla y otro de la evolución temporal de temperaturas del punto medio de la misma.

{{matlab|codigo=clear all
%datos del problema
L=4;
h=0.1;
T=8;
%discretizacion espacial
N=L/h;
%vector de puntos en el espacio
x=0:h:L;
xint=h:h:L-h;%nodos interiores
%aproximacion de -Uxx
K=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
K=(2/h^2)*K;
%calculamos F
F=zeros(N-1,1);
F(1)=F(1)+2*5/h^2;
%calculamos U0
U0=((exp((-6)*(xint-2).^2))+(5-5/4.*xint))';
%discretizacion temporal
dt=h^2/4; %paso en el espacio
t=0:dt:T;%vector en tiempos
%metodo del trapecio
I=eye(N-1);
M=I-dt*K;
sol(:,1)=U0;
for n=1:length(t)-1
sol(:,n+1)=M*sol(:,n)+dt*F;
end
UA=5*ones(1,length(t));
UB=0*ones(1,length(t));
sol=[UA;sol;UB];
%dibujamos la solucion

[xx,tt]=meshgrid(x,t);

figure(1)
mesh(xx,tt,sol')
xlabel('espacio')
ylabel('tiempo')
zlabel('Temperatura')
figure(2)
plot(t,sol(20,:))
xlabel('tiempo')
ylabel('temperatura')}}

Obteniendo los gráficos:

[[Archivo:Dfeulerexp.png|marco|centro|Evolución de la temperatura en la varila frente al tiempo y la posición de los puntos de la misma.]]

Al igual que en el apartado anterior vemos que se estabiliza muy rápido. Debido al tener que disminuir bastante el tamaño del paso temporal hemos tenido que cambiar la función de Octave para representar este gráfico.

[[Archivo:Dfeulerexppm.png|marco|centro|Sucesivas temperaturas del punto medio frente al tiempo.]]

Vemos que en esta implementación los cambios de temperatura en el punto medio son más ''continuos'', pero también acaba siendo asintótica en un valor aproximado de 2.625.
===Método de Euler implícito===
Es más fácil hacerlo converger que su homónimo explícito, pero tiene la desventaja de tener que despejar en la ecuación para implementarlo. A continuación se muestra el código que lo implementa en Matlab u Octave, dándonos una gráfica en 3D del comportamiento de la temperatura frente al tiempo en la varilla y otro de la evolución temporal de temperaturas del punto medio de la misma.

{{matlab|codigo=clear all
%datos del problema
L=4;
h=0.1;
T=8;
%discretizacion espacial
N=L/h;
%vector de puntos en el espacio
x=0:h:L;
xint=h:h:L-h;%nodos interiores
%aproximacion de -Uxx
K=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
K=(2/h^2)*K;
%calculamos F
F=zeros(N-1,1);
F(1)=F(1)+2*5/h^2;
%calculamos U0
U0=((exp((-6)*(xint-2).^2))+(5-5/4.*xint))';
%discretizacion temporal
dt=h; %paso en el espacio
t=0:dt:T;%vector en tiempos
%metodo del trapecio
I=eye(N-1);
M=I+dt*K;
sol(:,1)=U0;
for n=1:length(t)-1
sol(:,n+1)=M\(sol(:,n)+dt*(F));
end
UA=5*ones(1,length(t));
UB=0*ones(1,length(t));
sol=[UA;sol;UB];;
%dibujamos la solucion

[xx,tt]=meshgrid(x,t);

figure(1)
surf(xx,tt,sol')
xlabel('espacio')
ylabel('tiempo')
zlabel('temperatura')
figure(2)
plot(t,sol(20,:))
xlabel('tiempo')
ylabel('temperatura')}}

Obteniendo los siguientes gráficos:

[[Archivo:Dfeulerimp.png|marco|centro|Evolución de la temperatura en la varilla frente al tiempo y la posición de los puntos de la misma.]]

Volvemos a ver que la estabilización de temperaturas en la varilla es bastante rápida.

[[Archivo:Dfeulerimppm.png|marco|centro|Sucesivas temperaturas del punto medio frente al tiempo.]]

Una vez más observamos un descenso asintótico hacia el valor de temperatura 2.625.
===Método de Runge-Kutta===
Se trata de un método explícito, por tanto necesitaremos aumentar el tamaño del paso una vez más hasta los mismos valores que para Euler explícito para asegurar su convergencia. A continuación se muestra el código que lo implementa en Matlab u Octave, dandonos una gráfica en 3D del comportamiento de la temperatura frente al tiempo en la varilla y otro de la evolución temporal de temperaturas del punto medio de la misma.

{{matlab|codigo=clear all
%datos del problema
L=4;
h=0.1;
T=8;
%discretizacion espacial
N=L/h;
%vector de puntos en el espacio
x=0:h:L;
xint=h:h:L-h;%nodos interiores
%aproximacion de -Uxx
K=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
K=(2/h^2)*K;
%calculamos F
F=zeros(N-1,1);
F(1)=F(1)+2*5/h^2;
%calculamos U0
U0=((exp((-6)*(xint-2).^2))+(5-5/4.*xint))';
%discretizacion temporal
dt=h^2/4; %paso en el espacio
t=0:dt:T;%vector en tiempos
%metodo del trapecio
I=eye(N-1);
M=I-dt.*K;
sol(:,1)=U0;
for n=1:length(t)-1
k1=-K*sol(:,n)+F;
k2=-K*(sol(:,n)+k1*dt/2)+F;
k3=-K*(sol(:,n)+k2*dt/2)+F;
k4=-K*(sol(:,n)+k3*dt)+F;
sol(:,n+1)=sol(:,n)+(dt/6)*(k1+2*k2+2*k3+k4);
end
UA=5*ones(1,length(t));
UB=0*ones(1,length(t));
sol=[UA;sol;UB];
%dibujamos la solucion

[xx,tt]=meshgrid(x,t);

figure(1)
mesh(xx,tt,sol')
xlabel('espacio')
ylabel('tiempo')
zlabel('temperatura')
figure(2)
plot(t,sol(20,:))
xlabel('tiempo')
ylabel('temperatura')}}

Obteniendo las siguientes gráficas:

[[Archivo:Dfrk.png|marco|centro|Evolución de la temperatura en la varilla frente al tiempo y la posición de los puntos de la misma.]]

Al igual que en todos los métodos anteriores se produce la ya esperada estabilización de temperaturas en la varilla de manera bastante acelerada.

[[Archivo:Dfrkpm.png|marco|centro|Sucesivas temperaturas del punto medio frente al tiempo.]]

De nuevo aparece un descenso asintótico hasta el valor 2.625.

==Método de Fourier==
Partiendo de una solución analítica previa obtenemos la solución en forma de sumatorio de infinitos términos, reduciendo este número de términos numéricamente calcularemos una solución aproximada. Los autovalores y autovectores obtenidos para este problema y con condiciones de contorno son::

<math>\mu_k=\frac{k^2\pi^2}{16}</math>:

<math>\varphi_k(x)=sen(\frac{k\pi x}{4})</math>:

<math>k=1,2,3,...</math>

Dados estos valores y la función de homogenización para las condiciones de contorno::

<math>v(x,t)=5-\frac{5}{4}x</math>

Mediante el siguiente código de Matlab implementamos el problema numéricamente, obteniendo así la gráfica en 3D de la solución para <math>k=20</math> (por ser la más representativa y exacta),y la comparación de las gráficas en el plano de <math>t=0,5</math> para <math>k=1,3,5,10 y 20</math>.

{{matlab|codigo=clear all
N=100;
a=0; b=4;
h=(b-a)/N; x=0:h:4; % space discretization
ht=h; t=0:ht:5; % time discretization
[xx,tt]=meshgrid(x,t); % time discretization
f=(exp((-6)*(x-2).^2)); % Initial data
aprox=0; % Compute the approximation
q=1;
for k=1:q
autovalor=((k*pi)^2)/16; % eigenvalue
autofuncion=sin(k*pi*x/4); % eigenfunction
c(k)=(trapz(x,f.*autofuncion))/trapz(x,autofuncion.^2); % Fourier coef.
aprox=aprox+c(k)*exp(-2*autovalor.*tt).*sin(k*pi*xx/4);
end
aprox2=aprox+5-(5/4)*xx;
figure(1)
plot(x,aprox2(0.5/ht-0.5,:)) % Draw the approximation
hold on
q=3;
for k=1:q
autovalor=((k*pi)^2)/16; % eigenvalue
autofuncion=sin(k*pi*x/4); % eigenfunction
c(k)=(trapz(x,f.*autofuncion))/trapz(x,autofuncion.^2); % Fourier coef.
aprox=aprox+c(k)*exp(-2*autovalor.*tt).*sin(k*pi*xx/4);
end
aprox2=aprox+5-(5/4)*xx;
figure(1)
plot(x,aprox2(0.5/ht-0.5,:),'r') % Draw the approximation
hold on
q=5;
for k=1:q
autovalor=((k*pi)^2)/16; % eigenvalue
autofuncion=sin(k*pi*x/4); % eigenfunction
c(k)=(trapz(x,f.*autofuncion))/trapz(x,autofuncion.^2); % Fourier coef.
aprox=aprox+c(k)*exp(-2*autovalor.*tt).*sin(k*pi*xx/4);
end
aprox2=aprox+5-(5/4)*xx;
figure(1)
plot(x,aprox2(0.5/ht-0.5,:),'m') % Draw the approximation
hold on
q=10;
for k=1:q
autovalor=((k*pi)^2)/16; % eigenvalue
autofuncion=sin(k*pi*x/4); % eigenfunction
c(k)=(trapz(x,f.*autofuncion))/trapz(x,autofuncion.^2); % Fourier coef.
aprox=aprox+c(k)*exp(-2*autovalor.*tt).*sin(k*pi*xx/4);
end
aprox2=aprox+5-(5/4)*xx;
figure(1)
plot(x,aprox2(0.5/ht-0.5,:),'g') % Draw the approximation
hold on
q=20;
for k=1:q
autovalor=((k*pi)^2)/16; % eigenvalue
autofuncion=sin(k*pi*x/4); % eigenfunction
c(k)=(trapz(x,f.*autofuncion))/trapz(x,autofuncion.^2); % Fourier coef.
aprox=aprox+c(k)*exp(-2*autovalor.*tt).*sin(k*pi*xx/4);
end
aprox2=aprox+5-(5/4)*xx;
figure(1)
plot(x,aprox2(0.5/ht-0.5,:),'y') % Draw the approximation
hold on
plot(x,5-5/4*x,'k')
figure(2)
mesh(xx,tt,aprox2)
hold on}}

Obteniendo así:

[[Archivo:MF.png|marco|centro|Gráfico que muestra de la temperatura de la varilla frente a al tiempo a lo largo de sus diferentes posiciones.]]

Se puede ver como a diferencia del método de lineas el error viene inducido por la aproximación mediante series de Fourier de la condición inicial. Por lo demás evoluciona de manera similar a los métodos aplicados anteriormente.

[[Archivo:MF20.png|marco|centro|Gráficas de las temperaturas a lo largo de la varilla con aproximaciones de diferentes índices del sumatorio en la solución. k=1(azul). k=3(rojo). k=5(magenta). k=10(verde). k=20(amarillo). solución estacionaria(negro)]]

En la gráfica anterior se puede apreciar el nivel de aproximación con relación al número de coeficientes que utilizamos de la serie. Vemos como cuantos menos coeficientes usamos más nos acercamos incluso en tiempo cortos a la solución estacionaria, lo que nos indica que desciende el nivel de precisión de estas aproximaciones.

=Cambio de condiciones de contorno=
Ahora manteniendo la ecuación que nos modelizaba el comportamiento térmico de la varilla, cambiamos las condiciones de contorno en los extremos. Las nuevas condiciones propuestas se pueden asimilar como el aislamiento térmico total del extremo izquierdo, lo que anula el flujo de calor en ese lado. Estas son:

<math>u_x(0,t)=0</math>:
<math>u(4,t)=0</math>

Para implementarlo en sus nuevas condiciones usaremos el método de lineas con el método del trapecio. El código que implementa esta modelización en lenguaje Matlab y nos proporciona un gráfico de la temperatura frente a las variables espacial y temporal es:

{{matlab|codigo=clear all
%datos del problema
L=4;
h=0.1;
T=8;
%discretizacion espacial
N=L/h;
%vector de puntos en el espacio
x=0:h:L;
xint=0:h:L-h;%nodos interiores
%aproximacion de -Uxx
K=2*diag(ones(1,N))-diag(ones(1,N-1),1)-diag(ones(1,N-1),-1);
K=(2/h^2)*K;
%calculamos F
F=zeros(N,1);
F(1)=F(1)+2*5/h^2;
%calculamos U0
U0=((exp((-6)*(xint-2).^2))+(5-5/4.*xint))';
%discretizacion temporal
dt=h; %paso en el espacio
t=0:dt:T;%vector en tiempos
%metodo del trapecio
I=eye(N);
M=I+(dt./2)*K;
sol(:,1)=U0;
for n=1:length(t)-1
sol(:,n+1)=M\(sol(:,n)+dt/2*(F)+dt/2*(-K*sol(:,n)+F));
end
UB=zeros(length(t),1);
sol=[sol;UB'];
%dibujamos la solucion (apartado 3)
[xx,tt]=meshgrid(x,t);
figure(1)
surf(xx,tt,sol')
xlabel('espacio')
ylabel('tiempo')
zlabel('Temperatura')}}

Del cual obtenemos:

[[Archivo:Dfneu.png|marco|centro|Progreso de la temperatura frente a tiempo y espacio con las nuevas condiciones.]]

Como se puede ver en la imagen anterior la temperatura de los puntos interiores de la varilla va descendiendo lenta y progresivamente, a diferencia de en el caso anterior donde esta tendía a homogeneizarse.

=Interpretaciones estacionarias=
Como hemos visto en los apartados anteriores de este artículo, las soluciones por norma general de esta ecuación tienden a volverse estacionarias en unos intervalos de tiempo bastante breves. Por ello en esta sección estudiaremos estas soluciones estacionarias y el error que se comete asumiéndolas como ciertas en tiempos relativamente cortos.

Para calcular analiticamente estas soluciones estacionarias basta con considerar nulo el término en <math>u_t(x,t)</math> de la ecuación en derivadas parciales, y resolver lo que ahora nos queda un problema de contorno con dichas condiciones y la ecuación diferencial ordinaria que nos queda de hacer la despreciación anterior.

Tratando con las condiciones de contorno dadas inicialmente::

<math>u(0,t)=5</math>:
<math>u(4,t)=0</math>

Analíticamente para estas condiciones la solución estacionaria hallada es::

<math>u_{est}(x)=5-\frac{5}{4}x</math>

Mediante cálculo numérico y el siguiente código de Octave, conseguimos una gráfica de los diferentes valores de la temperatura en ciertos tiempos que tendrán que ir acercándose a la solución estacionaria, y un gráfico del error cometido si aproximásemos en estos tiempos el valor numérico aproximado de la temperatura por el estacionario calculado analíticamente. Además, nos aporta una tabla con en la primera fila el tiempo para el que se calcula el error respecto de la estacionaria y en las sucesivas filar los errores correspondientes.

{{matlab|codigo=clear all
%datos del problema
L=4;
h=0.1;
T=10;
%discretizacion espacial
N=L/h;
%vector de puntos en el espacio
x=0:h:L;
xint=h:h:L-h;%nodos interiores
%aproximacion de -Uxx
K=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
K=(2/h^2)*K;
%calculamos F
F=zeros(N-1,1);
F(1)=F(1)+2*5/h^2;
%calculamos U0
U0=((exp((-6)*(xint-2).^2))+(5-5/4.*xint))';
%discretizacion temporal
dt=T/L*h; %paso en el espacio
t=0:dt:T;%vector en tiempos
%metodo del trapecio
I=eye(N-1);
M=I+(dt./2)*K;
sol(:,1)=U0;
for n=1:length(t)-1
sol(:,n+1)=M\(sol(:,n)+dt/2*(F)+dt/2*(-K*sol(:,n)+F));
end
UA=5*ones(1,length(t));
UB=0*ones(1,length(t));
sol=[UA;sol;UB];;
%dibujamos la solucion (apartado 3)

[xx,tt]=meshgrid(x,t);

%Posiciones
tiempos=[0,1,2,10];
posicion=tiempos/dt;
figure(2)
%dibujamos para ciertos tiempos
for m=1:length(posicion)
hold on
plot(x,sol(:,posicion(m)+1),'m');
end
xlabel('espacio')
ylabel('temperaturas en cada t')
%Estacionaria
hold on
est=5*(1-x/4);
plot(x,est);
figure(3)
%errores
for p=1:length(posicion)
error=sol(:,posicion(p)+1).-est';
Error(:,p)=error;
hold on
plot(x,error)
end
xlabel('espacio')
ylabel('error cometido')
[0,1,2,10;Error]}}

Obteniendo:

[[Archivo:Estacionaria1.png|marco|centro|Contraste de la gráfica de la estacionaria con las soluciones calculadas mediante numérico en diferentes tiempos.]]

Se puede observar como se reduce el error a medida que aumentamos el tiempo de modelización cor respecto de la estacionaria.

[[Archivo:Estacionaria2.png|marco|centro|Error cometido al aproximar en diferentes tiempos por la estacionaria.]]

En la imagen anterior se puede ver una comparativa del el error con respecto al eje x, donde se ve más claramente que arriba como se reduce.

[[Archivo:Tablaerror.png|marco|centro|Tabla que nos muestra el error cometido en los diferentes tiempos.]]

En los valores numéricos dispuestos en la enterior tabla ya se muestra con total seguridad como se reduce el error a medida que aumentamos el tiempo.

Si ahora por el contrario tratamos las condiciones propuestas en el apartado de cambio de condiciones de contorno, las cuales son:

<math>u_x(0,t)=0</math>:
<math>u(4,t)=0</math>

Analíticamente para estas condiciones la solución estacionaria hallada es::

<math>u_{est}(x)=0</math>

Lo que al principio nos puede chocar al ver la gráfica, ya que en las otras condiciones se llega al estado estacionario de manera mucho más rápida. Si consideramos como estado estacionario que la diferencia entre dos valores consecutivos sea menor de 0.05, mediante el siguiente programa do Octave podemos calcular este valor de tiempo:

{{matlab|codigo=clear all
%datos del problema
L=4;
h=0.1;
T=8;
%discretizacion espacial
N=L/h;
%vector de puntos en el espacio
x=0:h:L;
xint=0:h:L-h;%nodos interiores
%aproximacion de -Uxx
K=2*diag(ones(1,N))-diag(ones(1,N-1),1)-diag(ones(1,N-1),-1);
K=(2/h^2)*K;
%calculamos F
F=zeros(N,1);
F(1)=F(1)+2*5/h^2;
%calculamos U0
U0=((exp((-6)*(xint-2).^2))+(5-5/4.*xint))';
%discretizacion temporal
dt=h; %paso en el espacio
t=0:dt:T;%vector en tiempos
%metodo del trapecio
I=eye(N);
M=I+(dt./2)*K;
sol(:,1)=U0;
for n=1:length(t)-1
sol(:,n+1)=M\(sol(:,n)+dt/2*(F)+dt/2*(-K*sol(:,n)+F));
end
UB=zeros(length(t),1);
sol=[sol;UB'];for i=1:length(t)
er=abs(sol(:,i+1)-sol(:,i));
a=max(er);
if a<=0.05
cont=i;
tiempo=(cont-1)*dt;
break
end
end
tiempo}}

El programa nos dice que el tiempo en el que se cumple esta condición es <math>t=0.7</math>.

=Intercambio de calor con el entorno. Varilla no aislada longitudinalmente=
Si ahora panteamos el intercambio de calor con el exterior, suponiendo que este ambiente se encuentra a una temperatura de 16ºC, la ecuación diferencial se ve afectada de manera que nos queda::

<math>u_t(x,t)-2u_{xx}(x,t)+u(x,t)-16=0</math>

Implementando numéricamente esta ecuación con las condiciones iniciales propuestas en el primer apartado de este artículo, mediante el método de lineas y el método del trapecio, mediante el siguiente código, obteniendo así la gráfica de la temperatura en cada punto de la varilla con respecto al tiempo.

{{matlab|codigo=clear all
%datos del problema
L=4;
h=0.1;
T=10;
%discretizacion espacial
N=L/h;
%vector de puntos en el espacio
x=0:h:L;
xint=h:h:L-h;%nodos interiores
%aproximacion de -Uxx
K=(2+h^2)*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
K=(1/h^2)*K;
%calculamos F
F=zeros(N-1,1).+16;;
F(1)=F(1)+2*5/h^2;
%calculamos U0
U0=((exp((-6)*(xint-2).^2))+(5-5/4.*xint))';
%discretizacion temporal
dt=h; %paso en el espacio
t=0:dt:T;%vector en tiempos
%metodo del trapecio
I=eye(N-1);
M=I+(dt./2)*K;
sol(:,1)=U0;
for n=1:length(t)-1
sol(:,n+1)=M\(sol(:,n)+dt/2*(F)+dt/2*(-K*sol(:,n)+F));
end
UA=5*ones(1,length(t));
UB=0*ones(1,length(t));
sol=[UA;sol;UB];;
%dibujamos la solucion (apartado 3)

[xx,tt]=meshgrid(x,t);

figure(1)
surf(xx,tt,sol')
xlabel('espacio')
ylabel('tiempo')
zlabel('Temperatura')}}

Obteniendo de este código:

[[Archivo:Perdidasaire.png|marco|centro|Representación de la temperatura en cada punto de la varilla con respecto al tiempo]]

Donde se ve como el exterior aporta calor a la varilla elevando sus temperaturas frente a las iniciales.

Si considerásemos el caso estacionario, se puede calcular mediante la inclusión del siguiente fragmento de código en el anterior programa que modelizaba el fenómeno no estacionario. Además de darnos la situación de temperaturas estacionarias, nos da el tiempo en el que este se alcanza con un error relativo de <math>10^{-3}</math>.

{{matlab|codigo=%Estacionaria
figure(2)
xlabel('espacio')
ylabel('temperatura')
c1=det([-11,1;-16,exp(-4)])/det([1,1;exp(4),exp(-4)]);
c2=det([1,-11;exp(4),-16])/det([1,1;exp(4),exp(-4)]);
est=c1*exp(x)+c2*exp(-x)+16;
plot(x,est);
hold on
for i=1:length(t)-1
er=abs(sol(:,i+1)-sol(:,i));
a=max(er);
if a<=0.001
cont=i;
ti=(cont-1)*dt;
break
end
end
ti
cont}}

La gráfica conseguida del código anterior es:

[[Archivo:Perdidasaireest.png|marco|centro|Temperatura estacionaria de la varilla con respecto a la posición de cada punto de ella.]]

El valor de tiempo obtenido es <math>t=5.9</math>.

==Cambio de condiciones de contorno==
Si cambiamos las condiciones de contorno de la ecuación por estas:

<math>u_x(0,t)=2</math>:
<math>u(4,t)=-2</math>

Físicamente podrían interpretarse como la conexión de un deposito de fluido en el extremo derecho que obligara un flujo constante de calor en el mismo extremo, y que el extremo derecho estuviera posicionado sobre un objeto que le obligara a tener una temperatura constante.

Si implementamos numéricamente este nuevo modelo propuesto mediante el siguiente código Matlab, obtenemos la representación de las temperaturas de la varilla frente al tiempo.

{{matlab|codigo=clear all
%datos del problema
L=4;
h=0.1;
T=10;
%discretizacion espacial
N=L/h;
%vector de puntos en el espacio
x=0:h:L;
xint=h:h:L;%nodos interiores
%aproximacion de -Uxx
K=(2+h^2)*diag(ones(1,N))-diag(ones(1,N-1),1)-diag(ones(1,N-1),-1);
K(1,1)=K(1,1)-1;
K=(1/h^2)*K;
%calculamos F
F=zeros(N,1).+16;;
F(1)=F(1)+2/h;
F(N)=F(N)+2/h^2;
%calculamos U0
U0=((exp((-6)*(xint-2).^2))+(5-5/4.*xint))';
%discretizacion temporal
dt=h; %paso en el espacio
t=0:dt:T;%vector en tiempos
%metodo del trapecio
I=eye(N);
M=I+(dt./2)*K;
sol(:,1)=U0;
for n=1:length(t)-1
sol(:,n+1)=M\(sol(:,n)+dt/2*(F)+dt/2*(-K*sol(:,n)+F));
end
UB=2*ones(1,length(t));
sol=[sol;UB];;
%dibujamos la solucion (apartado 3)

[xx,tt]=meshgrid(x,t);

figure(1)
surf(xx,tt,sol')
xlabel('espacio')
ylabel('tiempo')
zlabel('Temperatura')

}}

Logrando así la gráfica:

[[Archivo:Perdidasneu111.png|marco|centro|Representación de las temperaturas de la varilla frente al tiempo]]

Se puede observar como la temperatura de los puntos de la varilla va aumentado debido al flujo entrante por el lado izquierdo de esta hasta que alcanza un estado estacionario.

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

Ecuación del calor. (Grupo A8)

2015-05-15T19:24:23Z

Antonio Carrero Muñoz: /* Método de Fourier */

{{ TrabajoED |Ecuación del calor (Grupo A8)| [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED14/15|Curso 2014-15]] | Valentina Salazar; Antonio Carrero; José Francisco Aguilera }}
=Introducción y modelización del problema.=
En este artículo trataremos el comportamiento térmico de una varilla sometida a ciertas condiciones térmicas y físicas mediante diferentes métodos numéricos, daremos interpretación física a los resultados arrojados por estos métodos. Basaremos el estudio analítico y numérico necesario en la ecuación del calor propuesta por Jean-Baptiste Joseph Fourier en 1822.

Para estudiar el problema consideraremos una varilla delgada, de sección constante y de un material homogéneo, de longitud <math>L=4</math>. La situaremos en el intervalo <math>x\in{(0,4)}</math> de la recta real. Para llevar acabo esta modelización tendremos que asumir unas ciertas hipótesis y simplificaciones:

*Para el primer caso a plantear la varilla estará infinitamente aislada del entorno y por tanto no habrá flujo de calor sobre la superficie lateral de la misma.
*Al estar considerando el caso unidimensional de la ecuación del calor, asumiremos que al ser una varilla delgada la temperatura a lo largo de una sección ortogonal al eje <math>x</math> se mantiene constante en toda la sección.
*Consideraremos que el calor específico del material, <math>c</math>, es constante y no depende de la temperatura, y por tanto la difusividad térmica,<math>\alpha=\frac{k}{c\rho}</math> también lo será.

Si damos un valor a las constantes <math>c</math>, <math>k</math> y <math>\rho</math>, tal que <math>\alpha=2</math>; y no hay fuentes ni sumideros de calor dentro de la varilla, la ecuación en derivadas parciales que tendrá que cumplir la distribución de temperaturas a lo largo de la misma es::

<math>u_t(x,t)-2u_{xx}(x,t)=0</math>

Donde:
*<math>u(x,t)</math> nos define como evoluciona la temperatura a lo largo del eje <math>x</math>.
*<math>u_t(x,t)</math> es la derivada de la temperatura con respecto del tiempo.
*<math>u_{xx}(x,t)</math> es la segunda derivada de la temperatura con respecto de <math>x</math>.

Junto con está ecuación en derivadas parciales consideraremos una serie de condiciones: iniciales y de contorno. Como condición inicial para está modelización propondremos::

<math>u(x,0)=e^{-6(x-2)^2}+5-\frac{5}{4}x</math>

Donde <math>u(x,0)</math> no es otra cosa que la distribución de temperaturas inicial.

A lo largo del artículo variaremos las condiciones de contorno para dar diferentes ejemplos, pero en el bloque de casos prácticos trabajaremos con las siguientes::

<math>u(0,t)=5</math>:
<math>u(4,t)=0</math>

Éstas indican la temperatura (en este caso, otras pueden ser los flujos de calor o combinaciones de ambos) en cada uno de los extremos de la varilla.

=Casos prácticos=
En los apartados que siguen estudiaremos a través del cálculo numérico el caso anteriormente expuesto e interpretaremos los resultados obtenidos mediante dicho análisis. El lenguaje usado en el que se implementarán estos modelos será código Matlab u Octave.
==Método de diferencias finitas o método de lineas==
Trataremos de reducir el sistema de ecuaciones a una sola variable en una discretización de <math>n</math> valores (siendo <math>h</math> el tamaño de los subintervalos) haciendo la aproximación::

<math>u_{xx}(x,t)=\frac{-U^{N-1}(t)+2U^N(t)-U^{N+1}(t)}{h^2}</math>

Así podemos aplicar uno de los diferentes métodos iterativos en la variable <math>t</math> que se muestran en los subapartados que siguen.

===Método del trapecio===
Se trata de un método implícito, es decir que hay que despejar en la ecuación matricial que nos da la <math>U(t)</math>. A continuación se muestra el código que lo implementa en Matlab u Octave, dandonos una gráfica en 3D del comportamiento de la temperatura frente al tiempo en la varilla y otro de la evolución temporal de temperaturas del punto medio de la misma.

{{matlab|codigo=clear all
%datos del problema
L=4;
h=0.1;
T=8;
%discretizacion espacial
N=L/h;
%vector de puntos en el espacio
x=0:h:L;
xint=h:h:L-h;%nodos interiores
%aproximacion de -Uxx
K=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
K=(2/h^2)*K;
%calculamos F
F=zeros(N-1,1);
F(1)=F(1)+2*5/h^2;
%calculamos U0
U0=((exp((-6)*(xint-2).^2))+(5-5/4.*xint))';
%discretizacion temporal
dt=h; %paso en el espacio
t=0:dt:T;%vector en tiempos
%metodo del trapecio
I=eye(N-1);
M=I+(dt./2)*K;
sol(:,1)=U0;
for n=1:length(t)-1
sol(:,n+1)=M\(sol(:,n)+dt/2*(F)+dt/2*(-K*sol(:,n)+F));
end
UA=5*ones(1,length(t));
UB=0*ones(1,length(t));
sol=[UA;sol;UB];;
%dibujamos la solucion (apartado 3)

[xx,tt]=meshgrid(x,t);

figure(1)
surf(xx,tt,sol')
xlabel('espacio')
ylabel('tiempo')
zlabel('temperatura')
figure(2)
plot(t,sol(20,:))
xlabel('tiempo')
ylabel('temperatura')}}

Obteniendo los siguientes gráficos:

[[Archivo:Dftrap.png|marco|centro|Evolución de la temperatura en la varilla frente al tiempo y la posición de los puntos de la misma.]]

Como podemos ver la temperatura apenas tarda en estabilizarse en el tiempo a lo largo de los puntos de la varilla, esto nos sugiere que una buena aproximación será la solución estacionaria de la ecuación.

[[Archivo:Dftrappm.png|marco|centro|Sucesivas temperaturas del punto medio frente al tiempo.]]

La temperatura del punto medio descenderá hasta volverse casi asintótica hacia un valor aproximado de 2.625.

===Método de Euler explícito===
A diferencia del método del trapecio, como su nombre indica, este método es explicito, lo que es una ventaja a la hora de implementar ya que nos ahorramos el despejar; pero sin embargo se trata de un método bastante inestable, y para que arroje soluciones lógicas hemos de disminuir el tamaño del paso temporal frente al espacial del orden de <math>dt=\frac{h^2}{4}</math>, siendo <math>h</math>, <math>dt</math> dichos pasos espacial y temporal respectivamente. A continuación se muestra el código que lo implementa en Matlab u Octave, dandonos una gráfica en 3D del comportamiento de la temperatura frente al tiempo en la varilla y otro de la evolución temporal de temperaturas del punto medio de la misma.

{{matlab|codigo=clear all
%datos del problema
L=4;
h=0.1;
T=8;
%discretizacion espacial
N=L/h;
%vector de puntos en el espacio
x=0:h:L;
xint=h:h:L-h;%nodos interiores
%aproximacion de -Uxx
K=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
K=(2/h^2)*K;
%calculamos F
F=zeros(N-1,1);
F(1)=F(1)+2*5/h^2;
%calculamos U0
U0=((exp((-6)*(xint-2).^2))+(5-5/4.*xint))';
%discretizacion temporal
dt=h^2/4; %paso en el espacio
t=0:dt:T;%vector en tiempos
%metodo del trapecio
I=eye(N-1);
M=I-dt*K;
sol(:,1)=U0;
for n=1:length(t)-1
sol(:,n+1)=M*sol(:,n)+dt*F;
end
UA=5*ones(1,length(t));
UB=0*ones(1,length(t));
sol=[UA;sol;UB];
%dibujamos la solucion

[xx,tt]=meshgrid(x,t);

figure(1)
mesh(xx,tt,sol')
xlabel('espacio')
ylabel('tiempo')
zlabel('Temperatura')
figure(2)
plot(t,sol(20,:))
xlabel('tiempo')
ylabel('temperatura')}}

Obteniendo los gráficos:

[[Archivo:Dfeulerexp.png|marco|centro|Evolución de la temperatura en la varila frente al tiempo y la posición de los puntos de la misma.]]

Al igual que en el apartado anterior vemos que se estabiliza muy rápido. Debido al tener que disminuir bastante el tamaño del paso temporal hemos tenido que cambiar la función de Octave para representar este gráfico.

[[Archivo:Dfeulerexppm.png|marco|centro|Sucesivas temperaturas del punto medio frente al tiempo.]]

Vemos que en esta implementación los cambios de temperatura en el punto medio son más ''continuos'', pero también acaba siendo asintótica en un valor aproximado de 2.625.
===Método de Euler implícito===
Es más fácil hacerlo converger que su homónimo explícito, pero tiene la desventaja de tener que despejar en la ecuación para implementarlo. A continuación se muestra el código que lo implementa en Matlab u Octave, dándonos una gráfica en 3D del comportamiento de la temperatura frente al tiempo en la varilla y otro de la evolución temporal de temperaturas del punto medio de la misma.

{{matlab|codigo=clear all
%datos del problema
L=4;
h=0.1;
T=8;
%discretizacion espacial
N=L/h;
%vector de puntos en el espacio
x=0:h:L;
xint=h:h:L-h;%nodos interiores
%aproximacion de -Uxx
K=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
K=(2/h^2)*K;
%calculamos F
F=zeros(N-1,1);
F(1)=F(1)+2*5/h^2;
%calculamos U0
U0=((exp((-6)*(xint-2).^2))+(5-5/4.*xint))';
%discretizacion temporal
dt=h; %paso en el espacio
t=0:dt:T;%vector en tiempos
%metodo del trapecio
I=eye(N-1);
M=I+dt*K;
sol(:,1)=U0;
for n=1:length(t)-1
sol(:,n+1)=M\(sol(:,n)+dt*(F));
end
UA=5*ones(1,length(t));
UB=0*ones(1,length(t));
sol=[UA;sol;UB];;
%dibujamos la solucion

[xx,tt]=meshgrid(x,t);

figure(1)
surf(xx,tt,sol')
xlabel('espacio')
ylabel('tiempo')
zlabel('temperatura')
figure(2)
plot(t,sol(20,:))
xlabel('tiempo')
ylabel('temperatura')}}

Obteniendo los siguientes gráficos:

[[Archivo:Dfeulerimp.png|marco|centro|Evolución de la temperatura en la varilla frente al tiempo y la posición de los puntos de la misma.]]

Volvemos a ver que la estabilización de temperaturas en la varilla es bastante rápida.

[[Archivo:Dfeulerimppm.png|marco|centro|Sucesivas temperaturas del punto medio frente al tiempo.]]

Una vez más observamos un descenso asintótico hacia el valor de temperatura 2.625.
===Método de Runge-Kutta===
Se trata de un método explícito, por tanto necesitaremos aumentar el tamaño del paso una vez más hasta los mismos valores que para Euler explícito para asegurar su convergencia. A continuación se muestra el código que lo implementa en Matlab u Octave, dandonos una gráfica en 3D del comportamiento de la temperatura frente al tiempo en la varilla y otro de la evolución temporal de temperaturas del punto medio de la misma.

{{matlab|codigo=clear all
%datos del problema
L=4;
h=0.1;
T=8;
%discretizacion espacial
N=L/h;
%vector de puntos en el espacio
x=0:h:L;
xint=h:h:L-h;%nodos interiores
%aproximacion de -Uxx
K=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
K=(2/h^2)*K;
%calculamos F
F=zeros(N-1,1);
F(1)=F(1)+2*5/h^2;
%calculamos U0
U0=((exp((-6)*(xint-2).^2))+(5-5/4.*xint))';
%discretizacion temporal
dt=h^2/4; %paso en el espacio
t=0:dt:T;%vector en tiempos
%metodo del trapecio
I=eye(N-1);
M=I-dt.*K;
sol(:,1)=U0;
for n=1:length(t)-1
k1=-K*sol(:,n)+F;
k2=-K*(sol(:,n)+k1*dt/2)+F;
k3=-K*(sol(:,n)+k2*dt/2)+F;
k4=-K*(sol(:,n)+k3*dt)+F;
sol(:,n+1)=sol(:,n)+(dt/6)*(k1+2*k2+2*k3+k4);
end
UA=5*ones(1,length(t));
UB=0*ones(1,length(t));
sol=[UA;sol;UB];
%dibujamos la solucion

[xx,tt]=meshgrid(x,t);

figure(1)
mesh(xx,tt,sol')
xlabel('espacio')
ylabel('tiempo')
zlabel('temperatura')
figure(2)
plot(t,sol(20,:))
xlabel('tiempo')
ylabel('temperatura')}}

Obteniendo las siguientes gráficas:

[[Archivo:Dfrk.png|marco|centro|Evolución de la temperatura en la varilla frente al tiempo y la posición de los puntos de la misma.]]

Al igual que en todos los métodos anteriores se produce la ya esperada estabilización de temperaturas en la varilla de manera bastante acelerada.

[[Archivo:Dfrkpm.png|marco|centro|Sucesivas temperaturas del punto medio frente al tiempo.]]

De nuevo aparece un descenso asintótico hasta el valor 2.625.

==Método de Fourier==
Partiendo de una solución analítica previa obtenemos la solución en forma de sumatorio de infinitos términos, reduciendo este número de términos numéricamente calcularemos una solución aproximada. Los autovalores y autovectores obtenidos para este problema y con condiciones de contorno son::

<math>\mu_k=\frac{k^2\pi^2}{16}</math>:

<math>\varphi_k(x)=sen(\frac{k\pi x}{4})</math>:

<math>k=1,2,3,...</math>

Dados estos valores y la función de homogenización para las condiciones de contorno::

<math>v(x,t)=5-\frac{5}{4}x</math>

Mediante el siguiente código de Matlab implementamos el problema numéricamente, obteniendo así la gráfica en 3D de la solución para <math>k=20</math> (por ser la más representativa y exacta),y la comparación de las gráficas en el plano de <math>t=0,5</math> para <math>k=1,3,5,10 y 20</math>.

{{matlab|codigo=clear all
N=100;
a=0; b=4;
h=(b-a)/N; x=0:h:4; % space discretization
ht=h; t=0:ht:5; % time discretization
[xx,tt]=meshgrid(x,t); % time discretization
f=(exp((-6)*(x-2).^2)); % Initial data
aprox=0; % Compute the approximation
q=1;
for k=1:q
autovalor=((k*pi)^2)/16; % eigenvalue
autofuncion=sin(k*pi*x/4); % eigenfunction
c(k)=(trapz(x,f.*autofuncion))/trapz(x,autofuncion.^2); % Fourier coef.
aprox=aprox+c(k)*exp(-2*autovalor.*tt).*sin(k*pi*xx/4);
end
aprox2=aprox+5-(5/4)*xx;
figure(1)
plot(x,aprox2(0.5/ht-0.5,:)) % Draw the approximation
hold on
q=3;
for k=1:q
autovalor=((k*pi)^2)/16; % eigenvalue
autofuncion=sin(k*pi*x/4); % eigenfunction
c(k)=(trapz(x,f.*autofuncion))/trapz(x,autofuncion.^2); % Fourier coef.
aprox=aprox+c(k)*exp(-2*autovalor.*tt).*sin(k*pi*xx/4);
end
aprox2=aprox+5-(5/4)*xx;
figure(1)
plot(x,aprox2(0.5/ht-0.5,:),'r') % Draw the approximation
hold on
q=5;
for k=1:q
autovalor=((k*pi)^2)/16; % eigenvalue
autofuncion=sin(k*pi*x/4); % eigenfunction
c(k)=(trapz(x,f.*autofuncion))/trapz(x,autofuncion.^2); % Fourier coef.
aprox=aprox+c(k)*exp(-2*autovalor.*tt).*sin(k*pi*xx/4);
end
aprox2=aprox+5-(5/4)*xx;
figure(1)
plot(x,aprox2(0.5/ht-0.5,:),'m') % Draw the approximation
hold on
q=10;
for k=1:q
autovalor=((k*pi)^2)/16; % eigenvalue
autofuncion=sin(k*pi*x/4); % eigenfunction
c(k)=(trapz(x,f.*autofuncion))/trapz(x,autofuncion.^2); % Fourier coef.
aprox=aprox+c(k)*exp(-2*autovalor.*tt).*sin(k*pi*xx/4);
end
aprox2=aprox+5-(5/4)*xx;
figure(1)
plot(x,aprox2(0.5/ht-0.5,:),'g') % Draw the approximation
hold on
q=20;
for k=1:q
autovalor=((k*pi)^2)/16; % eigenvalue
autofuncion=sin(k*pi*x/4); % eigenfunction
c(k)=(trapz(x,f.*autofuncion))/trapz(x,autofuncion.^2); % Fourier coef.
aprox=aprox+c(k)*exp(-2*autovalor.*tt).*sin(k*pi*xx/4);
end
aprox2=aprox+5-(5/4)*xx;
figure(1)
plot(x,aprox2(0.5/ht-0.5,:),'y') % Draw the approximation
hold on
plot(x,5-5/4*x,'k')
figure(2)
mesh(xx,tt,aprox2)
hold on}}

Obteniendo así:

[[Archivo:MF.png|marco|centro|Gráfico que muestra de la temperatura de la varilla frente a al tiempo a lo largo de sus diferentes posiciones.]]

Se puede ver como a diferencia del método de lineas el error viene inducido por la aproximación mediante series de Fourier de la condición inicial. Por lo demás evoluciona de manera similar a los métodos aplicados anteriormente.

[[Archivo:MF20.png|marco|centro|Gráficas de las temperaturas a lo largo de la varilla con aproximaciones de diferentes índices del sumatorio en la solución. k=1(azul). k=3(rojo). k=5(magenta). k=10(verde). k=20(amarillo). solución estacionaria(negro)]]

En la gráfica anterior se puede apreciar el nivel de aproximación con relación al número de coeficientes que utilizamos de la serie. Vemos como cuantos menos coeficientes usamos más nos acercamos incluso en tiempo cortos a la solución estacionaria, lo que nos indica que desciende el nivel de precisión de estas aproximaciones.

=Cambio de condiciones de contorno=
Ahora manteniendo la ecuación que nos modelizaba el comportamiento térmico de la varilla, cambiamos las condiciones de contorno en los extremos. Las nuevas condiciones propuestas se pueden asimilar como el aislamiento térmico total del extremo izquierdo, lo que anula el flujo de calor en ese lado. Estas son:

<math>u_x(0,t)=0</math>:
<math>u(4,t)=0</math>

Para implementarlo en sus nuevas condiciones usaremos el método de lineas con el método del trapecio. El código que implementa esta modelización en lenguaje Matlab y nos proporciona un gráfico de la temperatura frente a las variables espacial y temporal es:

{{matlab|codigo=clear all
%datos del problema
L=4;
h=0.1;
T=8;
%discretizacion espacial
N=L/h;
%vector de puntos en el espacio
x=0:h:L;
xint=0:h:L-h;%nodos interiores
%aproximacion de -Uxx
K=2*diag(ones(1,N))-diag(ones(1,N-1),1)-diag(ones(1,N-1),-1);
K=(2/h^2)*K;
%calculamos F
F=zeros(N,1);
F(1)=F(1)+2*5/h^2;
%calculamos U0
U0=((exp((-6)*(xint-2).^2))+(5-5/4.*xint))';
%discretizacion temporal
dt=h; %paso en el espacio
t=0:dt:T;%vector en tiempos
%metodo del trapecio
I=eye(N);
M=I+(dt./2)*K;
sol(:,1)=U0;
for n=1:length(t)-1
sol(:,n+1)=M\(sol(:,n)+dt/2*(F)+dt/2*(-K*sol(:,n)+F));
end
UB=zeros(length(t),1);
sol=[sol;UB'];
%dibujamos la solucion (apartado 3)
[xx,tt]=meshgrid(x,t);
figure(1)
surf(xx,tt,sol')
xlabel('espacio')
ylabel('tiempo')
zlabel('Temperatura')}}

Del cual obtenemos:

[[Archivo:Dfneu.png|marco|centro|Progreso de la temperatura frente a tiempo y espacio con las nuevas condiciones.]]

Como se puede ver en la imagen anterior la temperatura de los puntos interiores de la varilla va descendiendo lenta y progresivamente, a diferencia de en el caso anterior donde esta tendía a homogeneizarse.

=Interpretaciones estacionarias=
Como hemos visto en los apartados anteriores de este artículo, las soluciones por norma general de esta ecuación tienden a volverse estacionarias en unos intervalos de tiempo bastante breves. Por ello en esta sección estudiaremos estas soluciones estacionarias y el error que se comete asumiéndolas como ciertas en tiempos relativamente cortos.

Para calcular analiticamente estas soluciones estacionarias basta con considerar nulo el término en <math>u_t(x,t)</math> dela ecuación en derivadas parciales, y resolver lo que ahora nos queda un problema de contorno con dichas condiciones y la ecuación diferencial ordinaria que nos queda de hacer la despreciación anterior.

Tratando con las condiciones de contorno dadas inicialmente::

<math>u(0,t)=5</math>:
<math>u(4,t)=0</math>

Analíticamente para estas condiciones la solución estacionaria hallada es::

<math>u_{est}(x)=5-\frac{5}{4}x</math>

Mediante cálculo numérico y el siguiente código de Octave, conseguimos una gráfica de los diferentes valores de la temperatura en ciertos tiempos que tendrán que ir acercándose a la solución estacionaria, y un gráfico del error cometido si aproximásemos en estos tiempos el valor numérico aproximado de la temperatura por el estacionario calculado analíticamente. Además, nos aporta una tabla con en la primera fila el tiempo para el que se calcula el error respecto de la estacionaria y en las sucesivas filar los errores correspondientes.

{{matlab|codigo=clear all
%datos del problema
L=4;
h=0.1;
T=10;
%discretizacion espacial
N=L/h;
%vector de puntos en el espacio
x=0:h:L;
xint=h:h:L-h;%nodos interiores
%aproximacion de -Uxx
K=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
K=(2/h^2)*K;
%calculamos F
F=zeros(N-1,1);
F(1)=F(1)+2*5/h^2;
%calculamos U0
U0=((exp((-6)*(xint-2).^2))+(5-5/4.*xint))';
%discretizacion temporal
dt=T/L*h; %paso en el espacio
t=0:dt:T;%vector en tiempos
%metodo del trapecio
I=eye(N-1);
M=I+(dt./2)*K;
sol(:,1)=U0;
for n=1:length(t)-1
sol(:,n+1)=M\(sol(:,n)+dt/2*(F)+dt/2*(-K*sol(:,n)+F));
end
UA=5*ones(1,length(t));
UB=0*ones(1,length(t));
sol=[UA;sol;UB];;
%dibujamos la solucion (apartado 3)

[xx,tt]=meshgrid(x,t);

%Posiciones
tiempos=[0,1,2,10];
posicion=tiempos/dt;
figure(2)
%dibujamos para ciertos tiempos
for m=1:length(posicion)
hold on
plot(x,sol(:,posicion(m)+1),'m');
end
xlabel('espacio')
ylabel('temperaturas en cada t')
%Estacionaria
hold on
est=5*(1-x/4);
plot(x,est);
figure(3)
%errores
for p=1:length(posicion)
error=sol(:,posicion(p)+1).-est';
Error(:,p)=error;
hold on
plot(x,error)
end
xlabel('espacio')
ylabel('error cometido')
[0,1,2,10;Error]}}

Obteniendo:

[[Archivo:Estacionaria1.png|marco|centro|Contraste de la gráfica de la estacionaria con las soluciones calculadas mediante numérico en diferentes tiempos.]]

Se puede observar como se reduce el error a medida que aumentamos el tiempo de modelización cor respecto de la estacionaria.

[[Archivo:Estacionaria2.png|marco|centro|Error cometido al aproximar en diferentes tiempos por la estacionaria.]]

En la imagen anterior se puede ver una comparativa del el error con respecto al eje x, donde se ve más claramente que arriba como se reduce.

[[Archivo:Tablaerror.png|marco|centro|Tabla que nos muestra el error cometido en los diferentes tiempos.]]

En los valores numéricos dispuestos en la enterior tabla ya se muestra con total seguridad como se reduce el error a medida que aumentamos el tiempo.

Si ahora por el contrario tratamos las condiciones propuestas en el apartado de cambio de condiciones de contorno, las cuales son:

<math>u_x(0,t)=0</math>:
<math>u(4,t)=0</math>

Analíticamente para estas condiciones la solución estacionaria hallada es::

<math>u_{est}(x)=0</math>

Lo que al principio nos puede chocar al ver la gráfica, ya que en las otras condiciones se llega al estado estacionario de manera mucho más rápida. Si consideramos como estado estacionario que la diferencia entre dos valores consecutivos sea menor de 0.05, mediante el siguiente programa do Octave podemos calcular este valor de tiempo:

{{matlab|codigo=clear all
%datos del problema
L=4;
h=0.1;
T=8;
%discretizacion espacial
N=L/h;
%vector de puntos en el espacio
x=0:h:L;
xint=0:h:L-h;%nodos interiores
%aproximacion de -Uxx
K=2*diag(ones(1,N))-diag(ones(1,N-1),1)-diag(ones(1,N-1),-1);
K=(2/h^2)*K;
%calculamos F
F=zeros(N,1);
F(1)=F(1)+2*5/h^2;
%calculamos U0
U0=((exp((-6)*(xint-2).^2))+(5-5/4.*xint))';
%discretizacion temporal
dt=h; %paso en el espacio
t=0:dt:T;%vector en tiempos
%metodo del trapecio
I=eye(N);
M=I+(dt./2)*K;
sol(:,1)=U0;
for n=1:length(t)-1
sol(:,n+1)=M\(sol(:,n)+dt/2*(F)+dt/2*(-K*sol(:,n)+F));
end
UB=zeros(length(t),1);
sol=[sol;UB'];for i=1:length(t)
er=abs(sol(:,i+1)-sol(:,i));
a=max(er);
if a<=0.05
cont=i;
tiempo=(cont-1)*dt;
break
end
end
tiempo}}

El programa nos dice que el tiempo en el que se cumple esta condición es <math>t=0.7</math>.

=Intercambio de calor con el entorno. Varilla no aislada longitudinalmente=
Si ahora panteamos el intercambio de calor con el exterior, suponiendo que este ambiente se encuentra a una temperatura de 16ºC, la ecuación diferencial se ve afectada de manera que nos queda::

<math>u_t(x,t)-2u_{xx}(x,t)+u(x,t)-16=0</math>

Implementando numéricamente esta ecuación con las condiciones iniciales propuestas en el primer apartado de este artículo, mediante el método de lineas y el método del trapecio, mediante el siguiente código, obteniendo así la gráfica de la temperatura en cada punto de la varilla con respecto al tiempo.

{{matlab|codigo=clear all
%datos del problema
L=4;
h=0.1;
T=10;
%discretizacion espacial
N=L/h;
%vector de puntos en el espacio
x=0:h:L;
xint=h:h:L-h;%nodos interiores
%aproximacion de -Uxx
K=(2+h^2)*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
K=(1/h^2)*K;
%calculamos F
F=zeros(N-1,1).+16;;
F(1)=F(1)+2*5/h^2;
%calculamos U0
U0=((exp((-6)*(xint-2).^2))+(5-5/4.*xint))';
%discretizacion temporal
dt=h; %paso en el espacio
t=0:dt:T;%vector en tiempos
%metodo del trapecio
I=eye(N-1);
M=I+(dt./2)*K;
sol(:,1)=U0;
for n=1:length(t)-1
sol(:,n+1)=M\(sol(:,n)+dt/2*(F)+dt/2*(-K*sol(:,n)+F));
end
UA=5*ones(1,length(t));
UB=0*ones(1,length(t));
sol=[UA;sol;UB];;
%dibujamos la solucion (apartado 3)

[xx,tt]=meshgrid(x,t);

figure(1)
surf(xx,tt,sol')
xlabel('espacio')
ylabel('tiempo')
zlabel('Temperatura')}}

Obteniendo de este código:

[[Archivo:Perdidasaire.png|marco|centro|Representación de la temperatura en cada punto de la varilla con respecto al tiempo]]

Donde se ve como el exterior aporta calor a la varilla elevando sus temperaturas frente a las iniciales.

Si considerásemos el caso estacionario, se puede calcular mediante la inclusión del siguiente fragmento de código en el anterior programa que modelizaba el fenómeno no estacionario. Además de darnos la situación de temperaturas estacionarias, nos da el tiempo en el que este se alcanza con un error relativo de <math>10^{-3}</math>.

{{matlab|codigo=%Estacionaria
figure(2)
xlabel('espacio')
ylabel('temperatura')
c1=det([-11,1;-16,exp(-4)])/det([1,1;exp(4),exp(-4)]);
c2=det([1,-11;exp(4),-16])/det([1,1;exp(4),exp(-4)]);
est=c1*exp(x)+c2*exp(-x)+16;
plot(x,est);
hold on
for i=1:length(t)-1
er=abs(sol(:,i+1)-sol(:,i));
a=max(er);
if a<=0.001
cont=i;
ti=(cont-1)*dt;
break
end
end
ti
cont}}

La gráfica conseguida del código anterior es:

[[Archivo:Perdidasaireest.png|marco|centro|Temperatura estacionaria de la varilla con respecto a la posición de cada punto de ella.]]

El valor de tiempo obtenido es <math>t=5.9</math>.

==Cambio de condiciones de contorno==
Si cambiamos las condiciones de contorno de la ecuación por estas:

<math>u_x(0,t)=2</math>:
<math>u(4,t)=-2</math>

Físicamente podrían interpretarse como la conexión de un deposito de fluido en el extremo derecho que obligara un flujo constante de calor en el mismo extremo, y que el extremo derecho estuviera posicionado sobre un objeto que le obligara a tener una temperatura constante.

Si implementamos numéricamente este nuevo modelo propuesto mediante el siguiente código Matlab, obtenemos la representación de las temperaturas de la varilla frente al tiempo.

{{matlab|codigo=clear all
%datos del problema
L=4;
h=0.1;
T=10;
%discretizacion espacial
N=L/h;
%vector de puntos en el espacio
x=0:h:L;
xint=h:h:L;%nodos interiores
%aproximacion de -Uxx
K=(2+h^2)*diag(ones(1,N))-diag(ones(1,N-1),1)-diag(ones(1,N-1),-1);
K(1,1)=K(1,1)-1;
K=(1/h^2)*K;
%calculamos F
F=zeros(N,1).+16;;
F(1)=F(1)+2/h;
F(N)=F(N)+2/h^2;
%calculamos U0
U0=((exp((-6)*(xint-2).^2))+(5-5/4.*xint))';
%discretizacion temporal
dt=h; %paso en el espacio
t=0:dt:T;%vector en tiempos
%metodo del trapecio
I=eye(N);
M=I+(dt./2)*K;
sol(:,1)=U0;
for n=1:length(t)-1
sol(:,n+1)=M\(sol(:,n)+dt/2*(F)+dt/2*(-K*sol(:,n)+F));
end
UB=2*ones(1,length(t));
sol=[sol;UB];;
%dibujamos la solucion (apartado 3)

[xx,tt]=meshgrid(x,t);

figure(1)
surf(xx,tt,sol')
xlabel('espacio')
ylabel('tiempo')
zlabel('Temperatura')

}}

Logrando así la gráfica:

[[Archivo:Perdidasneu111.png|marco|centro|Representación de las temperaturas de la varilla frente al tiempo]]

Se puede observar como la temperatura de los puntos de la varilla va aumentado debido al flujo entrante por el lado izquierdo de esta hasta que alcanza un estado estacionario.

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

Ecuación del calor. (Grupo A8)

2015-05-15T19:23:55Z

Antonio Carrero Muñoz: /* Método de Fourier */

{{ TrabajoED |Ecuación del calor (Grupo A8)| [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED14/15|Curso 2014-15]] | Valentina Salazar; Antonio Carrero; José Francisco Aguilera }}
=Introducción y modelización del problema.=
En este artículo trataremos el comportamiento térmico de una varilla sometida a ciertas condiciones térmicas y físicas mediante diferentes métodos numéricos, daremos interpretación física a los resultados arrojados por estos métodos. Basaremos el estudio analítico y numérico necesario en la ecuación del calor propuesta por Jean-Baptiste Joseph Fourier en 1822.

Para estudiar el problema consideraremos una varilla delgada, de sección constante y de un material homogéneo, de longitud <math>L=4</math>. La situaremos en el intervalo <math>x\in{(0,4)}</math> de la recta real. Para llevar acabo esta modelización tendremos que asumir unas ciertas hipótesis y simplificaciones:

*Para el primer caso a plantear la varilla estará infinitamente aislada del entorno y por tanto no habrá flujo de calor sobre la superficie lateral de la misma.
*Al estar considerando el caso unidimensional de la ecuación del calor, asumiremos que al ser una varilla delgada la temperatura a lo largo de una sección ortogonal al eje <math>x</math> se mantiene constante en toda la sección.
*Consideraremos que el calor específico del material, <math>c</math>, es constante y no depende de la temperatura, y por tanto la difusividad térmica,<math>\alpha=\frac{k}{c\rho}</math> también lo será.

Si damos un valor a las constantes <math>c</math>, <math>k</math> y <math>\rho</math>, tal que <math>\alpha=2</math>; y no hay fuentes ni sumideros de calor dentro de la varilla, la ecuación en derivadas parciales que tendrá que cumplir la distribución de temperaturas a lo largo de la misma es::

<math>u_t(x,t)-2u_{xx}(x,t)=0</math>

Donde:
*<math>u(x,t)</math> nos define como evoluciona la temperatura a lo largo del eje <math>x</math>.
*<math>u_t(x,t)</math> es la derivada de la temperatura con respecto del tiempo.
*<math>u_{xx}(x,t)</math> es la segunda derivada de la temperatura con respecto de <math>x</math>.

Junto con está ecuación en derivadas parciales consideraremos una serie de condiciones: iniciales y de contorno. Como condición inicial para está modelización propondremos::

<math>u(x,0)=e^{-6(x-2)^2}+5-\frac{5}{4}x</math>

Donde <math>u(x,0)</math> no es otra cosa que la distribución de temperaturas inicial.

A lo largo del artículo variaremos las condiciones de contorno para dar diferentes ejemplos, pero en el bloque de casos prácticos trabajaremos con las siguientes::

<math>u(0,t)=5</math>:
<math>u(4,t)=0</math>

Éstas indican la temperatura (en este caso, otras pueden ser los flujos de calor o combinaciones de ambos) en cada uno de los extremos de la varilla.

=Casos prácticos=
En los apartados que siguen estudiaremos a través del cálculo numérico el caso anteriormente expuesto e interpretaremos los resultados obtenidos mediante dicho análisis. El lenguaje usado en el que se implementarán estos modelos será código Matlab u Octave.
==Método de diferencias finitas o método de lineas==
Trataremos de reducir el sistema de ecuaciones a una sola variable en una discretización de <math>n</math> valores (siendo <math>h</math> el tamaño de los subintervalos) haciendo la aproximación::

<math>u_{xx}(x,t)=\frac{-U^{N-1}(t)+2U^N(t)-U^{N+1}(t)}{h^2}</math>

Así podemos aplicar uno de los diferentes métodos iterativos en la variable <math>t</math> que se muestran en los subapartados que siguen.

===Método del trapecio===
Se trata de un método implícito, es decir que hay que despejar en la ecuación matricial que nos da la <math>U(t)</math>. A continuación se muestra el código que lo implementa en Matlab u Octave, dandonos una gráfica en 3D del comportamiento de la temperatura frente al tiempo en la varilla y otro de la evolución temporal de temperaturas del punto medio de la misma.

{{matlab|codigo=clear all
%datos del problema
L=4;
h=0.1;
T=8;
%discretizacion espacial
N=L/h;
%vector de puntos en el espacio
x=0:h:L;
xint=h:h:L-h;%nodos interiores
%aproximacion de -Uxx
K=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
K=(2/h^2)*K;
%calculamos F
F=zeros(N-1,1);
F(1)=F(1)+2*5/h^2;
%calculamos U0
U0=((exp((-6)*(xint-2).^2))+(5-5/4.*xint))';
%discretizacion temporal
dt=h; %paso en el espacio
t=0:dt:T;%vector en tiempos
%metodo del trapecio
I=eye(N-1);
M=I+(dt./2)*K;
sol(:,1)=U0;
for n=1:length(t)-1
sol(:,n+1)=M\(sol(:,n)+dt/2*(F)+dt/2*(-K*sol(:,n)+F));
end
UA=5*ones(1,length(t));
UB=0*ones(1,length(t));
sol=[UA;sol;UB];;
%dibujamos la solucion (apartado 3)

[xx,tt]=meshgrid(x,t);

figure(1)
surf(xx,tt,sol')
xlabel('espacio')
ylabel('tiempo')
zlabel('temperatura')
figure(2)
plot(t,sol(20,:))
xlabel('tiempo')
ylabel('temperatura')}}

Obteniendo los siguientes gráficos:

[[Archivo:Dftrap.png|marco|centro|Evolución de la temperatura en la varilla frente al tiempo y la posición de los puntos de la misma.]]

Como podemos ver la temperatura apenas tarda en estabilizarse en el tiempo a lo largo de los puntos de la varilla, esto nos sugiere que una buena aproximación será la solución estacionaria de la ecuación.

[[Archivo:Dftrappm.png|marco|centro|Sucesivas temperaturas del punto medio frente al tiempo.]]

La temperatura del punto medio descenderá hasta volverse casi asintótica hacia un valor aproximado de 2.625.

===Método de Euler explícito===
A diferencia del método del trapecio, como su nombre indica, este método es explicito, lo que es una ventaja a la hora de implementar ya que nos ahorramos el despejar; pero sin embargo se trata de un método bastante inestable, y para que arroje soluciones lógicas hemos de disminuir el tamaño del paso temporal frente al espacial del orden de <math>dt=\frac{h^2}{4}</math>, siendo <math>h</math>, <math>dt</math> dichos pasos espacial y temporal respectivamente. A continuación se muestra el código que lo implementa en Matlab u Octave, dandonos una gráfica en 3D del comportamiento de la temperatura frente al tiempo en la varilla y otro de la evolución temporal de temperaturas del punto medio de la misma.

{{matlab|codigo=clear all
%datos del problema
L=4;
h=0.1;
T=8;
%discretizacion espacial
N=L/h;
%vector de puntos en el espacio
x=0:h:L;
xint=h:h:L-h;%nodos interiores
%aproximacion de -Uxx
K=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
K=(2/h^2)*K;
%calculamos F
F=zeros(N-1,1);
F(1)=F(1)+2*5/h^2;
%calculamos U0
U0=((exp((-6)*(xint-2).^2))+(5-5/4.*xint))';
%discretizacion temporal
dt=h^2/4; %paso en el espacio
t=0:dt:T;%vector en tiempos
%metodo del trapecio
I=eye(N-1);
M=I-dt*K;
sol(:,1)=U0;
for n=1:length(t)-1
sol(:,n+1)=M*sol(:,n)+dt*F;
end
UA=5*ones(1,length(t));
UB=0*ones(1,length(t));
sol=[UA;sol;UB];
%dibujamos la solucion

[xx,tt]=meshgrid(x,t);

figure(1)
mesh(xx,tt,sol')
xlabel('espacio')
ylabel('tiempo')
zlabel('Temperatura')
figure(2)
plot(t,sol(20,:))
xlabel('tiempo')
ylabel('temperatura')}}

Obteniendo los gráficos:

[[Archivo:Dfeulerexp.png|marco|centro|Evolución de la temperatura en la varila frente al tiempo y la posición de los puntos de la misma.]]

Al igual que en el apartado anterior vemos que se estabiliza muy rápido. Debido al tener que disminuir bastante el tamaño del paso temporal hemos tenido que cambiar la función de Octave para representar este gráfico.

[[Archivo:Dfeulerexppm.png|marco|centro|Sucesivas temperaturas del punto medio frente al tiempo.]]

Vemos que en esta implementación los cambios de temperatura en el punto medio son más ''continuos'', pero también acaba siendo asintótica en un valor aproximado de 2.625.
===Método de Euler implícito===
Es más fácil hacerlo converger que su homónimo explícito, pero tiene la desventaja de tener que despejar en la ecuación para implementarlo. A continuación se muestra el código que lo implementa en Matlab u Octave, dándonos una gráfica en 3D del comportamiento de la temperatura frente al tiempo en la varilla y otro de la evolución temporal de temperaturas del punto medio de la misma.

{{matlab|codigo=clear all
%datos del problema
L=4;
h=0.1;
T=8;
%discretizacion espacial
N=L/h;
%vector de puntos en el espacio
x=0:h:L;
xint=h:h:L-h;%nodos interiores
%aproximacion de -Uxx
K=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
K=(2/h^2)*K;
%calculamos F
F=zeros(N-1,1);
F(1)=F(1)+2*5/h^2;
%calculamos U0
U0=((exp((-6)*(xint-2).^2))+(5-5/4.*xint))';
%discretizacion temporal
dt=h; %paso en el espacio
t=0:dt:T;%vector en tiempos
%metodo del trapecio
I=eye(N-1);
M=I+dt*K;
sol(:,1)=U0;
for n=1:length(t)-1
sol(:,n+1)=M\(sol(:,n)+dt*(F));
end
UA=5*ones(1,length(t));
UB=0*ones(1,length(t));
sol=[UA;sol;UB];;
%dibujamos la solucion

[xx,tt]=meshgrid(x,t);

figure(1)
surf(xx,tt,sol')
xlabel('espacio')
ylabel('tiempo')
zlabel('temperatura')
figure(2)
plot(t,sol(20,:))
xlabel('tiempo')
ylabel('temperatura')}}

Obteniendo los siguientes gráficos:

[[Archivo:Dfeulerimp.png|marco|centro|Evolución de la temperatura en la varilla frente al tiempo y la posición de los puntos de la misma.]]

Volvemos a ver que la estabilización de temperaturas en la varilla es bastante rápida.

[[Archivo:Dfeulerimppm.png|marco|centro|Sucesivas temperaturas del punto medio frente al tiempo.]]

Una vez más observamos un descenso asintótico hacia el valor de temperatura 2.625.
===Método de Runge-Kutta===
Se trata de un método explícito, por tanto necesitaremos aumentar el tamaño del paso una vez más hasta los mismos valores que para Euler explícito para asegurar su convergencia. A continuación se muestra el código que lo implementa en Matlab u Octave, dandonos una gráfica en 3D del comportamiento de la temperatura frente al tiempo en la varilla y otro de la evolución temporal de temperaturas del punto medio de la misma.

{{matlab|codigo=clear all
%datos del problema
L=4;
h=0.1;
T=8;
%discretizacion espacial
N=L/h;
%vector de puntos en el espacio
x=0:h:L;
xint=h:h:L-h;%nodos interiores
%aproximacion de -Uxx
K=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
K=(2/h^2)*K;
%calculamos F
F=zeros(N-1,1);
F(1)=F(1)+2*5/h^2;
%calculamos U0
U0=((exp((-6)*(xint-2).^2))+(5-5/4.*xint))';
%discretizacion temporal
dt=h^2/4; %paso en el espacio
t=0:dt:T;%vector en tiempos
%metodo del trapecio
I=eye(N-1);
M=I-dt.*K;
sol(:,1)=U0;
for n=1:length(t)-1
k1=-K*sol(:,n)+F;
k2=-K*(sol(:,n)+k1*dt/2)+F;
k3=-K*(sol(:,n)+k2*dt/2)+F;
k4=-K*(sol(:,n)+k3*dt)+F;
sol(:,n+1)=sol(:,n)+(dt/6)*(k1+2*k2+2*k3+k4);
end
UA=5*ones(1,length(t));
UB=0*ones(1,length(t));
sol=[UA;sol;UB];
%dibujamos la solucion

[xx,tt]=meshgrid(x,t);

figure(1)
mesh(xx,tt,sol')
xlabel('espacio')
ylabel('tiempo')
zlabel('temperatura')
figure(2)
plot(t,sol(20,:))
xlabel('tiempo')
ylabel('temperatura')}}

Obteniendo las siguientes gráficas:

[[Archivo:Dfrk.png|marco|centro|Evolución de la temperatura en la varilla frente al tiempo y la posición de los puntos de la misma.]]

Al igual que en todos los métodos anteriores se produce la ya esperada estabilización de temperaturas en la varilla de manera bastante acelerada.

[[Archivo:Dfrkpm.png|marco|centro|Sucesivas temperaturas del punto medio frente al tiempo.]]

De nuevo aparece un descenso asintótico hasta el valor 2.625.

==Método de Fourier==
Partiendo de una solución analítica previa obtenemos la solución en forma de sumatorio de infinitos términos, reduciendo este número de términos numéricamente calcularemos una solución aproximada. Los autovalores y autovectores obtenidos para este problema y con condiciones de contorno son::

<math>\mu_k=\frac{k^2\pi^2}{16}</math>:

<math>\varphi_k(x)=sen(\frac{k\pi x}{4})</math>:

<math>k=1,2,3,...</math>

Dados estos valores y la función de homogenización para las condiciones de contorno::

<math>v(x,t)=5-\frac{5}{4}x</math>

Mediante el siguiente código de Matlab implementamos el problema numéricamente, obteniendo así la gráfica en 3D de la solución para <math>k=20</math> (por ser la más representativa y exacta),y la comparaciónde las gráficas en el plano de <math>t=0,5</math> para <math>k=1,3,5,10 y 20</math>.

{{matlab|codigo=clear all
N=100;
a=0; b=4;
h=(b-a)/N; x=0:h:4; % space discretization
ht=h; t=0:ht:5; % time discretization
[xx,tt]=meshgrid(x,t); % time discretization
f=(exp((-6)*(x-2).^2)); % Initial data
aprox=0; % Compute the approximation
q=1;
for k=1:q
autovalor=((k*pi)^2)/16; % eigenvalue
autofuncion=sin(k*pi*x/4); % eigenfunction
c(k)=(trapz(x,f.*autofuncion))/trapz(x,autofuncion.^2); % Fourier coef.
aprox=aprox+c(k)*exp(-2*autovalor.*tt).*sin(k*pi*xx/4);
end
aprox2=aprox+5-(5/4)*xx;
figure(1)
plot(x,aprox2(0.5/ht-0.5,:)) % Draw the approximation
hold on
q=3;
for k=1:q
autovalor=((k*pi)^2)/16; % eigenvalue
autofuncion=sin(k*pi*x/4); % eigenfunction
c(k)=(trapz(x,f.*autofuncion))/trapz(x,autofuncion.^2); % Fourier coef.
aprox=aprox+c(k)*exp(-2*autovalor.*tt).*sin(k*pi*xx/4);
end
aprox2=aprox+5-(5/4)*xx;
figure(1)
plot(x,aprox2(0.5/ht-0.5,:),'r') % Draw the approximation
hold on
q=5;
for k=1:q
autovalor=((k*pi)^2)/16; % eigenvalue
autofuncion=sin(k*pi*x/4); % eigenfunction
c(k)=(trapz(x,f.*autofuncion))/trapz(x,autofuncion.^2); % Fourier coef.
aprox=aprox+c(k)*exp(-2*autovalor.*tt).*sin(k*pi*xx/4);
end
aprox2=aprox+5-(5/4)*xx;
figure(1)
plot(x,aprox2(0.5/ht-0.5,:),'m') % Draw the approximation
hold on
q=10;
for k=1:q
autovalor=((k*pi)^2)/16; % eigenvalue
autofuncion=sin(k*pi*x/4); % eigenfunction
c(k)=(trapz(x,f.*autofuncion))/trapz(x,autofuncion.^2); % Fourier coef.
aprox=aprox+c(k)*exp(-2*autovalor.*tt).*sin(k*pi*xx/4);
end
aprox2=aprox+5-(5/4)*xx;
figure(1)
plot(x,aprox2(0.5/ht-0.5,:),'g') % Draw the approximation
hold on
q=20;
for k=1:q
autovalor=((k*pi)^2)/16; % eigenvalue
autofuncion=sin(k*pi*x/4); % eigenfunction
c(k)=(trapz(x,f.*autofuncion))/trapz(x,autofuncion.^2); % Fourier coef.
aprox=aprox+c(k)*exp(-2*autovalor.*tt).*sin(k*pi*xx/4);
end
aprox2=aprox+5-(5/4)*xx;
figure(1)
plot(x,aprox2(0.5/ht-0.5,:),'y') % Draw the approximation
hold on
plot(x,5-5/4*x,'k')
figure(2)
mesh(xx,tt,aprox2)
hold on}}

Obteniendo así:

[[Archivo:MF.png|marco|centro|Gráfico que muestra de la temperatura de la varilla frente a al tiempo a lo largo de sus diferentes posiciones.]]

Se puede ver como a diferencia del método de lineas el error viene inducido por la aproximación mediante series de Fourier de la condición inicial. Por lo demás evoluciona de manera similar a los métodos aplicados anteriormente.

[[Archivo:MF20.png|marco|centro|Gráficas de las temperaturas a lo largo de la varilla con aproximaciones de diferentes índices del sumatorio en la solución. k=1(azul). k=3(rojo). k=5(magenta). k=10(verde). k=20(amarillo). solución estacionaria(negro)]]

En la gráfica anterior se puede apreciar el nivel de aproximación con relación al número de coeficientes que utilizamos de la serie. Vemos como cuantos menos coeficientes usamos más nos acercamos incluso en tiempo cortos a la solución estacionaria, lo que nos indica que desciende el nivel de precisión de estas aproximaciones.

=Cambio de condiciones de contorno=
Ahora manteniendo la ecuación que nos modelizaba el comportamiento térmico de la varilla, cambiamos las condiciones de contorno en los extremos. Las nuevas condiciones propuestas se pueden asimilar como el aislamiento térmico total del extremo izquierdo, lo que anula el flujo de calor en ese lado. Estas son:

<math>u_x(0,t)=0</math>:
<math>u(4,t)=0</math>

Para implementarlo en sus nuevas condiciones usaremos el método de lineas con el método del trapecio. El código que implementa esta modelización en lenguaje Matlab y nos proporciona un gráfico de la temperatura frente a las variables espacial y temporal es:

{{matlab|codigo=clear all
%datos del problema
L=4;
h=0.1;
T=8;
%discretizacion espacial
N=L/h;
%vector de puntos en el espacio
x=0:h:L;
xint=0:h:L-h;%nodos interiores
%aproximacion de -Uxx
K=2*diag(ones(1,N))-diag(ones(1,N-1),1)-diag(ones(1,N-1),-1);
K=(2/h^2)*K;
%calculamos F
F=zeros(N,1);
F(1)=F(1)+2*5/h^2;
%calculamos U0
U0=((exp((-6)*(xint-2).^2))+(5-5/4.*xint))';
%discretizacion temporal
dt=h; %paso en el espacio
t=0:dt:T;%vector en tiempos
%metodo del trapecio
I=eye(N);
M=I+(dt./2)*K;
sol(:,1)=U0;
for n=1:length(t)-1
sol(:,n+1)=M\(sol(:,n)+dt/2*(F)+dt/2*(-K*sol(:,n)+F));
end
UB=zeros(length(t),1);
sol=[sol;UB'];
%dibujamos la solucion (apartado 3)
[xx,tt]=meshgrid(x,t);
figure(1)
surf(xx,tt,sol')
xlabel('espacio')
ylabel('tiempo')
zlabel('Temperatura')}}

Del cual obtenemos:

[[Archivo:Dfneu.png|marco|centro|Progreso de la temperatura frente a tiempo y espacio con las nuevas condiciones.]]

Como se puede ver en la imagen anterior la temperatura de los puntos interiores de la varilla va descendiendo lenta y progresivamente, a diferencia de en el caso anterior donde esta tendía a homogeneizarse.

=Interpretaciones estacionarias=
Como hemos visto en los apartados anteriores de este artículo, las soluciones por norma general de esta ecuación tienden a volverse estacionarias en unos intervalos de tiempo bastante breves. Por ello en esta sección estudiaremos estas soluciones estacionarias y el error que se comete asumiéndolas como ciertas en tiempos relativamente cortos.

Para calcular analiticamente estas soluciones estacionarias basta con considerar nulo el término en <math>u_t(x,t)</math> dela ecuación en derivadas parciales, y resolver lo que ahora nos queda un problema de contorno con dichas condiciones y la ecuación diferencial ordinaria que nos queda de hacer la despreciación anterior.

Tratando con las condiciones de contorno dadas inicialmente::

<math>u(0,t)=5</math>:
<math>u(4,t)=0</math>

Analíticamente para estas condiciones la solución estacionaria hallada es::

<math>u_{est}(x)=5-\frac{5}{4}x</math>

Mediante cálculo numérico y el siguiente código de Octave, conseguimos una gráfica de los diferentes valores de la temperatura en ciertos tiempos que tendrán que ir acercándose a la solución estacionaria, y un gráfico del error cometido si aproximásemos en estos tiempos el valor numérico aproximado de la temperatura por el estacionario calculado analíticamente. Además, nos aporta una tabla con en la primera fila el tiempo para el que se calcula el error respecto de la estacionaria y en las sucesivas filar los errores correspondientes.

{{matlab|codigo=clear all
%datos del problema
L=4;
h=0.1;
T=10;
%discretizacion espacial
N=L/h;
%vector de puntos en el espacio
x=0:h:L;
xint=h:h:L-h;%nodos interiores
%aproximacion de -Uxx
K=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
K=(2/h^2)*K;
%calculamos F
F=zeros(N-1,1);
F(1)=F(1)+2*5/h^2;
%calculamos U0
U0=((exp((-6)*(xint-2).^2))+(5-5/4.*xint))';
%discretizacion temporal
dt=T/L*h; %paso en el espacio
t=0:dt:T;%vector en tiempos
%metodo del trapecio
I=eye(N-1);
M=I+(dt./2)*K;
sol(:,1)=U0;
for n=1:length(t)-1
sol(:,n+1)=M\(sol(:,n)+dt/2*(F)+dt/2*(-K*sol(:,n)+F));
end
UA=5*ones(1,length(t));
UB=0*ones(1,length(t));
sol=[UA;sol;UB];;
%dibujamos la solucion (apartado 3)

[xx,tt]=meshgrid(x,t);

%Posiciones
tiempos=[0,1,2,10];
posicion=tiempos/dt;
figure(2)
%dibujamos para ciertos tiempos
for m=1:length(posicion)
hold on
plot(x,sol(:,posicion(m)+1),'m');
end
xlabel('espacio')
ylabel('temperaturas en cada t')
%Estacionaria
hold on
est=5*(1-x/4);
plot(x,est);
figure(3)
%errores
for p=1:length(posicion)
error=sol(:,posicion(p)+1).-est';
Error(:,p)=error;
hold on
plot(x,error)
end
xlabel('espacio')
ylabel('error cometido')
[0,1,2,10;Error]}}

Obteniendo:

[[Archivo:Estacionaria1.png|marco|centro|Contraste de la gráfica de la estacionaria con las soluciones calculadas mediante numérico en diferentes tiempos.]]

Se puede observar como se reduce el error a medida que aumentamos el tiempo de modelización cor respecto de la estacionaria.

[[Archivo:Estacionaria2.png|marco|centro|Error cometido al aproximar en diferentes tiempos por la estacionaria.]]

En la imagen anterior se puede ver una comparativa del el error con respecto al eje x, donde se ve más claramente que arriba como se reduce.

[[Archivo:Tablaerror.png|marco|centro|Tabla que nos muestra el error cometido en los diferentes tiempos.]]

En los valores numéricos dispuestos en la enterior tabla ya se muestra con total seguridad como se reduce el error a medida que aumentamos el tiempo.

Si ahora por el contrario tratamos las condiciones propuestas en el apartado de cambio de condiciones de contorno, las cuales son:

<math>u_x(0,t)=0</math>:
<math>u(4,t)=0</math>

Analíticamente para estas condiciones la solución estacionaria hallada es::

<math>u_{est}(x)=0</math>

Lo que al principio nos puede chocar al ver la gráfica, ya que en las otras condiciones se llega al estado estacionario de manera mucho más rápida. Si consideramos como estado estacionario que la diferencia entre dos valores consecutivos sea menor de 0.05, mediante el siguiente programa do Octave podemos calcular este valor de tiempo:

{{matlab|codigo=clear all
%datos del problema
L=4;
h=0.1;
T=8;
%discretizacion espacial
N=L/h;
%vector de puntos en el espacio
x=0:h:L;
xint=0:h:L-h;%nodos interiores
%aproximacion de -Uxx
K=2*diag(ones(1,N))-diag(ones(1,N-1),1)-diag(ones(1,N-1),-1);
K=(2/h^2)*K;
%calculamos F
F=zeros(N,1);
F(1)=F(1)+2*5/h^2;
%calculamos U0
U0=((exp((-6)*(xint-2).^2))+(5-5/4.*xint))';
%discretizacion temporal
dt=h; %paso en el espacio
t=0:dt:T;%vector en tiempos
%metodo del trapecio
I=eye(N);
M=I+(dt./2)*K;
sol(:,1)=U0;
for n=1:length(t)-1
sol(:,n+1)=M\(sol(:,n)+dt/2*(F)+dt/2*(-K*sol(:,n)+F));
end
UB=zeros(length(t),1);
sol=[sol;UB'];for i=1:length(t)
er=abs(sol(:,i+1)-sol(:,i));
a=max(er);
if a<=0.05
cont=i;
tiempo=(cont-1)*dt;
break
end
end
tiempo}}

El programa nos dice que el tiempo en el que se cumple esta condición es <math>t=0.7</math>.

=Intercambio de calor con el entorno. Varilla no aislada longitudinalmente=
Si ahora panteamos el intercambio de calor con el exterior, suponiendo que este ambiente se encuentra a una temperatura de 16ºC, la ecuación diferencial se ve afectada de manera que nos queda::

<math>u_t(x,t)-2u_{xx}(x,t)+u(x,t)-16=0</math>

Implementando numéricamente esta ecuación con las condiciones iniciales propuestas en el primer apartado de este artículo, mediante el método de lineas y el método del trapecio, mediante el siguiente código, obteniendo así la gráfica de la temperatura en cada punto de la varilla con respecto al tiempo.

{{matlab|codigo=clear all
%datos del problema
L=4;
h=0.1;
T=10;
%discretizacion espacial
N=L/h;
%vector de puntos en el espacio
x=0:h:L;
xint=h:h:L-h;%nodos interiores
%aproximacion de -Uxx
K=(2+h^2)*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
K=(1/h^2)*K;
%calculamos F
F=zeros(N-1,1).+16;;
F(1)=F(1)+2*5/h^2;
%calculamos U0
U0=((exp((-6)*(xint-2).^2))+(5-5/4.*xint))';
%discretizacion temporal
dt=h; %paso en el espacio
t=0:dt:T;%vector en tiempos
%metodo del trapecio
I=eye(N-1);
M=I+(dt./2)*K;
sol(:,1)=U0;
for n=1:length(t)-1
sol(:,n+1)=M\(sol(:,n)+dt/2*(F)+dt/2*(-K*sol(:,n)+F));
end
UA=5*ones(1,length(t));
UB=0*ones(1,length(t));
sol=[UA;sol;UB];;
%dibujamos la solucion (apartado 3)

[xx,tt]=meshgrid(x,t);

figure(1)
surf(xx,tt,sol')
xlabel('espacio')
ylabel('tiempo')
zlabel('Temperatura')}}

Obteniendo de este código:

[[Archivo:Perdidasaire.png|marco|centro|Representación de la temperatura en cada punto de la varilla con respecto al tiempo]]

Donde se ve como el exterior aporta calor a la varilla elevando sus temperaturas frente a las iniciales.

Si considerásemos el caso estacionario, se puede calcular mediante la inclusión del siguiente fragmento de código en el anterior programa que modelizaba el fenómeno no estacionario. Además de darnos la situación de temperaturas estacionarias, nos da el tiempo en el que este se alcanza con un error relativo de <math>10^{-3}</math>.

{{matlab|codigo=%Estacionaria
figure(2)
xlabel('espacio')
ylabel('temperatura')
c1=det([-11,1;-16,exp(-4)])/det([1,1;exp(4),exp(-4)]);
c2=det([1,-11;exp(4),-16])/det([1,1;exp(4),exp(-4)]);
est=c1*exp(x)+c2*exp(-x)+16;
plot(x,est);
hold on
for i=1:length(t)-1
er=abs(sol(:,i+1)-sol(:,i));
a=max(er);
if a<=0.001
cont=i;
ti=(cont-1)*dt;
break
end
end
ti
cont}}

La gráfica conseguida del código anterior es:

[[Archivo:Perdidasaireest.png|marco|centro|Temperatura estacionaria de la varilla con respecto a la posición de cada punto de ella.]]

El valor de tiempo obtenido es <math>t=5.9</math>.

==Cambio de condiciones de contorno==
Si cambiamos las condiciones de contorno de la ecuación por estas:

<math>u_x(0,t)=2</math>:
<math>u(4,t)=-2</math>

Físicamente podrían interpretarse como la conexión de un deposito de fluido en el extremo derecho que obligara un flujo constante de calor en el mismo extremo, y que el extremo derecho estuviera posicionado sobre un objeto que le obligara a tener una temperatura constante.

Si implementamos numéricamente este nuevo modelo propuesto mediante el siguiente código Matlab, obtenemos la representación de las temperaturas de la varilla frente al tiempo.

{{matlab|codigo=clear all
%datos del problema
L=4;
h=0.1;
T=10;
%discretizacion espacial
N=L/h;
%vector de puntos en el espacio
x=0:h:L;
xint=h:h:L;%nodos interiores
%aproximacion de -Uxx
K=(2+h^2)*diag(ones(1,N))-diag(ones(1,N-1),1)-diag(ones(1,N-1),-1);
K(1,1)=K(1,1)-1;
K=(1/h^2)*K;
%calculamos F
F=zeros(N,1).+16;;
F(1)=F(1)+2/h;
F(N)=F(N)+2/h^2;
%calculamos U0
U0=((exp((-6)*(xint-2).^2))+(5-5/4.*xint))';
%discretizacion temporal
dt=h; %paso en el espacio
t=0:dt:T;%vector en tiempos
%metodo del trapecio
I=eye(N);
M=I+(dt./2)*K;
sol(:,1)=U0;
for n=1:length(t)-1
sol(:,n+1)=M\(sol(:,n)+dt/2*(F)+dt/2*(-K*sol(:,n)+F));
end
UB=2*ones(1,length(t));
sol=[sol;UB];;
%dibujamos la solucion (apartado 3)

[xx,tt]=meshgrid(x,t);

figure(1)
surf(xx,tt,sol')
xlabel('espacio')
ylabel('tiempo')
zlabel('Temperatura')

}}

Logrando así la gráfica:

[[Archivo:Perdidasneu111.png|marco|centro|Representación de las temperaturas de la varilla frente al tiempo]]

Se puede observar como la temperatura de los puntos de la varilla va aumentado debido al flujo entrante por el lado izquierdo de esta hasta que alcanza un estado estacionario.

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

Ecuación del calor. (Grupo A8)

2015-05-15T19:21:55Z

Antonio Carrero Muñoz: /* Introducción y modelización del problema. */

{{ TrabajoED |Ecuación del calor (Grupo A8)| [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED14/15|Curso 2014-15]] | Valentina Salazar; Antonio Carrero; José Francisco Aguilera }}
=Introducción y modelización del problema.=
En este artículo trataremos el comportamiento térmico de una varilla sometida a ciertas condiciones térmicas y físicas mediante diferentes métodos numéricos, daremos interpretación física a los resultados arrojados por estos métodos. Basaremos el estudio analítico y numérico necesario en la ecuación del calor propuesta por Jean-Baptiste Joseph Fourier en 1822.

Para estudiar el problema consideraremos una varilla delgada, de sección constante y de un material homogéneo, de longitud <math>L=4</math>. La situaremos en el intervalo <math>x\in{(0,4)}</math> de la recta real. Para llevar acabo esta modelización tendremos que asumir unas ciertas hipótesis y simplificaciones:

*Para el primer caso a plantear la varilla estará infinitamente aislada del entorno y por tanto no habrá flujo de calor sobre la superficie lateral de la misma.
*Al estar considerando el caso unidimensional de la ecuación del calor, asumiremos que al ser una varilla delgada la temperatura a lo largo de una sección ortogonal al eje <math>x</math> se mantiene constante en toda la sección.
*Consideraremos que el calor específico del material, <math>c</math>, es constante y no depende de la temperatura, y por tanto la difusividad térmica,<math>\alpha=\frac{k}{c\rho}</math> también lo será.

Si damos un valor a las constantes <math>c</math>, <math>k</math> y <math>\rho</math>, tal que <math>\alpha=2</math>; y no hay fuentes ni sumideros de calor dentro de la varilla, la ecuación en derivadas parciales que tendrá que cumplir la distribución de temperaturas a lo largo de la misma es::

<math>u_t(x,t)-2u_{xx}(x,t)=0</math>

Donde:
*<math>u(x,t)</math> nos define como evoluciona la temperatura a lo largo del eje <math>x</math>.
*<math>u_t(x,t)</math> es la derivada de la temperatura con respecto del tiempo.
*<math>u_{xx}(x,t)</math> es la segunda derivada de la temperatura con respecto de <math>x</math>.

Junto con está ecuación en derivadas parciales consideraremos una serie de condiciones: iniciales y de contorno. Como condición inicial para está modelización propondremos::

<math>u(x,0)=e^{-6(x-2)^2}+5-\frac{5}{4}x</math>

Donde <math>u(x,0)</math> no es otra cosa que la distribución de temperaturas inicial.

A lo largo del artículo variaremos las condiciones de contorno para dar diferentes ejemplos, pero en el bloque de casos prácticos trabajaremos con las siguientes::

<math>u(0,t)=5</math>:
<math>u(4,t)=0</math>

Éstas indican la temperatura (en este caso, otras pueden ser los flujos de calor o combinaciones de ambos) en cada uno de los extremos de la varilla.

=Casos prácticos=
En los apartados que siguen estudiaremos a través del cálculo numérico el caso anteriormente expuesto e interpretaremos los resultados obtenidos mediante dicho análisis. El lenguaje usado en el que se implementarán estos modelos será código Matlab u Octave.
==Método de diferencias finitas o método de lineas==
Trataremos de reducir el sistema de ecuaciones a una sola variable en una discretización de <math>n</math> valores (siendo <math>h</math> el tamaño de los subintervalos) haciendo la aproximación::

<math>u_{xx}(x,t)=\frac{-U^{N-1}(t)+2U^N(t)-U^{N+1}(t)}{h^2}</math>

Así podemos aplicar uno de los diferentes métodos iterativos en la variable <math>t</math> que se muestran en los subapartados que siguen.

===Método del trapecio===
Se trata de un método implícito, es decir que hay que despejar en la ecuación matricial que nos da la <math>U(t)</math>. A continuación se muestra el código que lo implementa en Matlab u Octave, dandonos una gráfica en 3D del comportamiento de la temperatura frente al tiempo en la varilla y otro de la evolución temporal de temperaturas del punto medio de la misma.

{{matlab|codigo=clear all
%datos del problema
L=4;
h=0.1;
T=8;
%discretizacion espacial
N=L/h;
%vector de puntos en el espacio
x=0:h:L;
xint=h:h:L-h;%nodos interiores
%aproximacion de -Uxx
K=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
K=(2/h^2)*K;
%calculamos F
F=zeros(N-1,1);
F(1)=F(1)+2*5/h^2;
%calculamos U0
U0=((exp((-6)*(xint-2).^2))+(5-5/4.*xint))';
%discretizacion temporal
dt=h; %paso en el espacio
t=0:dt:T;%vector en tiempos
%metodo del trapecio
I=eye(N-1);
M=I+(dt./2)*K;
sol(:,1)=U0;
for n=1:length(t)-1
sol(:,n+1)=M\(sol(:,n)+dt/2*(F)+dt/2*(-K*sol(:,n)+F));
end
UA=5*ones(1,length(t));
UB=0*ones(1,length(t));
sol=[UA;sol;UB];;
%dibujamos la solucion (apartado 3)

[xx,tt]=meshgrid(x,t);

figure(1)
surf(xx,tt,sol')
xlabel('espacio')
ylabel('tiempo')
zlabel('temperatura')
figure(2)
plot(t,sol(20,:))
xlabel('tiempo')
ylabel('temperatura')}}

Obteniendo los siguientes gráficos:

[[Archivo:Dftrap.png|marco|centro|Evolución de la temperatura en la varilla frente al tiempo y la posición de los puntos de la misma.]]

Como podemos ver la temperatura apenas tarda en estabilizarse en el tiempo a lo largo de los puntos de la varilla, esto nos sugiere que una buena aproximación será la solución estacionaria de la ecuación.

[[Archivo:Dftrappm.png|marco|centro|Sucesivas temperaturas del punto medio frente al tiempo.]]

La temperatura del punto medio descenderá hasta volverse casi asintótica hacia un valor aproximado de 2.625.

===Método de Euler explícito===
A diferencia del método del trapecio, como su nombre indica, este método es explicito, lo que es una ventaja a la hora de implementar ya que nos ahorramos el despejar; pero sin embargo se trata de un método bastante inestable, y para que arroje soluciones lógicas hemos de disminuir el tamaño del paso temporal frente al espacial del orden de <math>dt=\frac{h^2}{4}</math>, siendo <math>h</math>, <math>dt</math> dichos pasos espacial y temporal respectivamente. A continuación se muestra el código que lo implementa en Matlab u Octave, dandonos una gráfica en 3D del comportamiento de la temperatura frente al tiempo en la varilla y otro de la evolución temporal de temperaturas del punto medio de la misma.

{{matlab|codigo=clear all
%datos del problema
L=4;
h=0.1;
T=8;
%discretizacion espacial
N=L/h;
%vector de puntos en el espacio
x=0:h:L;
xint=h:h:L-h;%nodos interiores
%aproximacion de -Uxx
K=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
K=(2/h^2)*K;
%calculamos F
F=zeros(N-1,1);
F(1)=F(1)+2*5/h^2;
%calculamos U0
U0=((exp((-6)*(xint-2).^2))+(5-5/4.*xint))';
%discretizacion temporal
dt=h^2/4; %paso en el espacio
t=0:dt:T;%vector en tiempos
%metodo del trapecio
I=eye(N-1);
M=I-dt*K;
sol(:,1)=U0;
for n=1:length(t)-1
sol(:,n+1)=M*sol(:,n)+dt*F;
end
UA=5*ones(1,length(t));
UB=0*ones(1,length(t));
sol=[UA;sol;UB];
%dibujamos la solucion

[xx,tt]=meshgrid(x,t);

figure(1)
mesh(xx,tt,sol')
xlabel('espacio')
ylabel('tiempo')
zlabel('Temperatura')
figure(2)
plot(t,sol(20,:))
xlabel('tiempo')
ylabel('temperatura')}}

Obteniendo los gráficos:

[[Archivo:Dfeulerexp.png|marco|centro|Evolución de la temperatura en la varila frente al tiempo y la posición de los puntos de la misma.]]

Al igual que en el apartado anterior vemos que se estabiliza muy rápido. Debido al tener que disminuir bastante el tamaño del paso temporal hemos tenido que cambiar la función de Octave para representar este gráfico.

[[Archivo:Dfeulerexppm.png|marco|centro|Sucesivas temperaturas del punto medio frente al tiempo.]]

Vemos que en esta implementación los cambios de temperatura en el punto medio son más ''continuos'', pero también acaba siendo asintótica en un valor aproximado de 2.625.
===Método de Euler implícito===
Es más fácil hacerlo converger que su homónimo explícito, pero tiene la desventaja de tener que despejar en la ecuación para implementarlo. A continuación se muestra el código que lo implementa en Matlab u Octave, dándonos una gráfica en 3D del comportamiento de la temperatura frente al tiempo en la varilla y otro de la evolución temporal de temperaturas del punto medio de la misma.

{{matlab|codigo=clear all
%datos del problema
L=4;
h=0.1;
T=8;
%discretizacion espacial
N=L/h;
%vector de puntos en el espacio
x=0:h:L;
xint=h:h:L-h;%nodos interiores
%aproximacion de -Uxx
K=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
K=(2/h^2)*K;
%calculamos F
F=zeros(N-1,1);
F(1)=F(1)+2*5/h^2;
%calculamos U0
U0=((exp((-6)*(xint-2).^2))+(5-5/4.*xint))';
%discretizacion temporal
dt=h; %paso en el espacio
t=0:dt:T;%vector en tiempos
%metodo del trapecio
I=eye(N-1);
M=I+dt*K;
sol(:,1)=U0;
for n=1:length(t)-1
sol(:,n+1)=M\(sol(:,n)+dt*(F));
end
UA=5*ones(1,length(t));
UB=0*ones(1,length(t));
sol=[UA;sol;UB];;
%dibujamos la solucion

[xx,tt]=meshgrid(x,t);

figure(1)
surf(xx,tt,sol')
xlabel('espacio')
ylabel('tiempo')
zlabel('temperatura')
figure(2)
plot(t,sol(20,:))
xlabel('tiempo')
ylabel('temperatura')}}

Obteniendo los siguientes gráficos:

[[Archivo:Dfeulerimp.png|marco|centro|Evolución de la temperatura en la varilla frente al tiempo y la posición de los puntos de la misma.]]

Volvemos a ver que la estabilización de temperaturas en la varilla es bastante rápida.

[[Archivo:Dfeulerimppm.png|marco|centro|Sucesivas temperaturas del punto medio frente al tiempo.]]

Una vez más observamos un descenso asintótico hacia el valor de temperatura 2.625.
===Método de Runge-Kutta===
Se trata de un método explícito, por tanto necesitaremos aumentar el tamaño del paso una vez más hasta los mismos valores que para Euler explícito para asegurar su convergencia. A continuación se muestra el código que lo implementa en Matlab u Octave, dandonos una gráfica en 3D del comportamiento de la temperatura frente al tiempo en la varilla y otro de la evolución temporal de temperaturas del punto medio de la misma.

{{matlab|codigo=clear all
%datos del problema
L=4;
h=0.1;
T=8;
%discretizacion espacial
N=L/h;
%vector de puntos en el espacio
x=0:h:L;
xint=h:h:L-h;%nodos interiores
%aproximacion de -Uxx
K=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
K=(2/h^2)*K;
%calculamos F
F=zeros(N-1,1);
F(1)=F(1)+2*5/h^2;
%calculamos U0
U0=((exp((-6)*(xint-2).^2))+(5-5/4.*xint))';
%discretizacion temporal
dt=h^2/4; %paso en el espacio
t=0:dt:T;%vector en tiempos
%metodo del trapecio
I=eye(N-1);
M=I-dt.*K;
sol(:,1)=U0;
for n=1:length(t)-1
k1=-K*sol(:,n)+F;
k2=-K*(sol(:,n)+k1*dt/2)+F;
k3=-K*(sol(:,n)+k2*dt/2)+F;
k4=-K*(sol(:,n)+k3*dt)+F;
sol(:,n+1)=sol(:,n)+(dt/6)*(k1+2*k2+2*k3+k4);
end
UA=5*ones(1,length(t));
UB=0*ones(1,length(t));
sol=[UA;sol;UB];
%dibujamos la solucion

[xx,tt]=meshgrid(x,t);

figure(1)
mesh(xx,tt,sol')
xlabel('espacio')
ylabel('tiempo')
zlabel('temperatura')
figure(2)
plot(t,sol(20,:))
xlabel('tiempo')
ylabel('temperatura')}}

Obteniendo las siguientes gráficas:

[[Archivo:Dfrk.png|marco|centro|Evolución de la temperatura en la varilla frente al tiempo y la posición de los puntos de la misma.]]

Al igual que en todos los métodos anteriores se produce la ya esperada estabilización de temperaturas en la varilla de manera bastante acelerada.

[[Archivo:Dfrkpm.png|marco|centro|Sucesivas temperaturas del punto medio frente al tiempo.]]

De nuevo aparece un descenso asintótico hasta el valor 2.625.

==Método de Fourier==
Partiendo de una solución analítica previa obtenemos la solución en forma de sumatorio de infinitos términos, reduciendo este número de términos numéricamente calcularemos una solución aproximada. Los autovalores y autovectores obtenidos para este problema y con condiciones de contorno son::

<math>\mu_k=\frac{k^2\pi^2}{16}</math>:

<math>\varphi_k(x)=sen(\frac{k\pi x}{4})</math>:

<math>k=1,2,3,...</math>

Dados estos valores y la función de homogenización para las condiciones de contorno::

<math>v(x,t)=5-\frac{5}{4}x</math>

Mediante el siguiente código de Matlab implementamos el problema numéricamente, obteniendo así la gráfica en 3D de la solución para <math>k=20</math> (por ser la más representativa y exacta),y la comparaciónde las gráficas en el plano de <math>t=0,5</math> para <math>k=1,3,5,10 y 20</math>.

{{matlab|codigo=clear all
N=100;
a=0; b=4;
h=(b-a)/N; x=0:h:4; % space discretization
ht=h; t=0:ht:5; % time discretization
[xx,tt]=meshgrid(x,t); % time discretization
f=(exp((-6)*(x-2).^2)); % Initial data
aprox=0; % Compute the approximation
q=1;
for k=1:q
autovalor=((k*pi)^2)/16; % eigenvalue
autofuncion=sin(k*pi*x/4); % eigenfunction
c(k)=(trapz(x,f.*autofuncion))/trapz(x,autofuncion.^2); % Fourier coef.
aprox=aprox+c(k)*exp(-2*autovalor.*tt).*sin(k*pi*xx/4);
end
aprox2=aprox+5-(5/4)*xx;
figure(1)
plot(x,aprox2(0.5/ht-0.5,:)) % Draw the approximation
hold on
q=3;
for k=1:q
autovalor=((k*pi)^2)/16; % eigenvalue
autofuncion=sin(k*pi*x/4); % eigenfunction
c(k)=(trapz(x,f.*autofuncion))/trapz(x,autofuncion.^2); % Fourier coef.
aprox=aprox+c(k)*exp(-2*autovalor.*tt).*sin(k*pi*xx/4);
end
aprox2=aprox+5-(5/4)*xx;
figure(1)
plot(x,aprox2(0.5/ht-0.5,:),'r') % Draw the approximation
hold on
q=5;
for k=1:q
autovalor=((k*pi)^2)/16; % eigenvalue
autofuncion=sin(k*pi*x/4); % eigenfunction
c(k)=(trapz(x,f.*autofuncion))/trapz(x,autofuncion.^2); % Fourier coef.
aprox=aprox+c(k)*exp(-2*autovalor.*tt).*sin(k*pi*xx/4);
end
aprox2=aprox+5-(5/4)*xx;
figure(1)
plot(x,aprox2(0.5/ht-0.5,:),'m') % Draw the approximation
hold on
q=10;
for k=1:q
autovalor=((k*pi)^2)/16; % eigenvalue
autofuncion=sin(k*pi*x/4); % eigenfunction
c(k)=(trapz(x,f.*autofuncion))/trapz(x,autofuncion.^2); % Fourier coef.
aprox=aprox+c(k)*exp(-2*autovalor.*tt).*sin(k*pi*xx/4);
end
aprox2=aprox+5-(5/4)*xx;
figure(1)
plot(x,aprox2(0.5/ht-0.5,:),'g') % Draw the approximation
hold on
q=20;
for k=1:q
autovalor=((k*pi)^2)/16; % eigenvalue
autofuncion=sin(k*pi*x/4); % eigenfunction
c(k)=(trapz(x,f.*autofuncion))/trapz(x,autofuncion.^2); % Fourier coef.
aprox=aprox+c(k)*exp(-2*autovalor.*tt).*sin(k*pi*xx/4);
end
aprox2=aprox+5-(5/4)*xx;
figure(1)
plot(x,aprox2(0.5/ht-0.5,:),'y') % Draw the approximation
hold on
plot(x,5-5/4*x,'k')
figure(2)
mesh(xx,tt,aprox2)
hold on}}

Obteniendo así:

[[Archivo:MF.png|marco|centro|Gráfico que muestra de la temperatura de la varilla frente a al tiempo a lo largo de sus diferentes posiciones.]]

Se puede ver como a diferencia del método de lineas el error viene inducido por la aproximación mediante series de Fourier de la condición inicial. Por lo demás evoluciona de manera similar a los métodos aplicados anteriormente.

[[Archivo:MF20.png|marco|centro|Gráficas de las temperaturas a lo largo de la varilla con aproximaciones de diferentes índices del sumatorio en la solución. k=1(azul). k=3(rojo). k=5(magenta). k=10(verde). k=20(amarillo). solución estacionaria(negro)]]

En la gráfica anterior se puede apreciar el nivel de aproximación con relación al número de coeficientes que utilizamos de la serie. Vemos como cuantos menos coeficientes usamos más nos acercamos incluso en tiempo cortos a la solución estacionaria, lo que nos indica que desciende el nivel de precisión de estas aproximaciones.

=Cambio de condiciones de contorno=
Ahora manteniendo la ecuación que nos modelizaba el comportamiento térmico de la varilla, cambiamos las condiciones de contorno en los extremos. Las nuevas condiciones propuestas se pueden asimilar como el aislamiento térmico total del extremo izquierdo, lo que anula el flujo de calor en ese lado. Estas son:

<math>u_x(0,t)=0</math>:
<math>u(4,t)=0</math>

Para implementarlo en sus nuevas condiciones usaremos el método de lineas con el método del trapecio. El código que implementa esta modelización en lenguaje Matlab y nos proporciona un gráfico de la temperatura frente a las variables espacial y temporal es:

{{matlab|codigo=clear all
%datos del problema
L=4;
h=0.1;
T=8;
%discretizacion espacial
N=L/h;
%vector de puntos en el espacio
x=0:h:L;
xint=0:h:L-h;%nodos interiores
%aproximacion de -Uxx
K=2*diag(ones(1,N))-diag(ones(1,N-1),1)-diag(ones(1,N-1),-1);
K=(2/h^2)*K;
%calculamos F
F=zeros(N,1);
F(1)=F(1)+2*5/h^2;
%calculamos U0
U0=((exp((-6)*(xint-2).^2))+(5-5/4.*xint))';
%discretizacion temporal
dt=h; %paso en el espacio
t=0:dt:T;%vector en tiempos
%metodo del trapecio
I=eye(N);
M=I+(dt./2)*K;
sol(:,1)=U0;
for n=1:length(t)-1
sol(:,n+1)=M\(sol(:,n)+dt/2*(F)+dt/2*(-K*sol(:,n)+F));
end
UB=zeros(length(t),1);
sol=[sol;UB'];
%dibujamos la solucion (apartado 3)
[xx,tt]=meshgrid(x,t);
figure(1)
surf(xx,tt,sol')
xlabel('espacio')
ylabel('tiempo')
zlabel('Temperatura')}}

Del cual obtenemos:

[[Archivo:Dfneu.png|marco|centro|Progreso de la temperatura frente a tiempo y espacio con las nuevas condiciones.]]

Como se puede ver en la imagen anterior la temperatura de los puntos interiores de la varilla va descendiendo lenta y progresivamente, a diferencia de en el caso anterior donde esta tendía a homogeneizarse.

=Interpretaciones estacionarias=
Como hemos visto en los apartados anteriores de este artículo, las soluciones por norma general de esta ecuación tienden a volverse estacionarias en unos intervalos de tiempo bastante breves. Por ello en esta sección estudiaremos estas soluciones estacionarias y el error que se comete asumiéndolas como ciertas en tiempos relativamente cortos.

Para calcular analiticamente estas soluciones estacionarias basta con considerar nulo el término en <math>u_t(x,t)</math> dela ecuación en derivadas parciales, y resolver lo que ahora nos queda un problema de contorno con dichas condiciones y la ecuación diferencial ordinaria que nos queda de hacer la despreciación anterior.

Tratando con las condiciones de contorno dadas inicialmente::

<math>u(0,t)=5</math>:
<math>u(4,t)=0</math>

Analíticamente para estas condiciones la solución estacionaria hallada es::

<math>u_{est}(x)=5-\frac{5}{4}x</math>

Mediante cálculo numérico y el siguiente código de Octave, conseguimos una gráfica de los diferentes valores de la temperatura en ciertos tiempos que tendrán que ir acercándose a la solución estacionaria, y un gráfico del error cometido si aproximásemos en estos tiempos el valor numérico aproximado de la temperatura por el estacionario calculado analíticamente. Además, nos aporta una tabla con en la primera fila el tiempo para el que se calcula el error respecto de la estacionaria y en las sucesivas filar los errores correspondientes.

{{matlab|codigo=clear all
%datos del problema
L=4;
h=0.1;
T=10;
%discretizacion espacial
N=L/h;
%vector de puntos en el espacio
x=0:h:L;
xint=h:h:L-h;%nodos interiores
%aproximacion de -Uxx
K=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
K=(2/h^2)*K;
%calculamos F
F=zeros(N-1,1);
F(1)=F(1)+2*5/h^2;
%calculamos U0
U0=((exp((-6)*(xint-2).^2))+(5-5/4.*xint))';
%discretizacion temporal
dt=T/L*h; %paso en el espacio
t=0:dt:T;%vector en tiempos
%metodo del trapecio
I=eye(N-1);
M=I+(dt./2)*K;
sol(:,1)=U0;
for n=1:length(t)-1
sol(:,n+1)=M\(sol(:,n)+dt/2*(F)+dt/2*(-K*sol(:,n)+F));
end
UA=5*ones(1,length(t));
UB=0*ones(1,length(t));
sol=[UA;sol;UB];;
%dibujamos la solucion (apartado 3)

[xx,tt]=meshgrid(x,t);

%Posiciones
tiempos=[0,1,2,10];
posicion=tiempos/dt;
figure(2)
%dibujamos para ciertos tiempos
for m=1:length(posicion)
hold on
plot(x,sol(:,posicion(m)+1),'m');
end
xlabel('espacio')
ylabel('temperaturas en cada t')
%Estacionaria
hold on
est=5*(1-x/4);
plot(x,est);
figure(3)
%errores
for p=1:length(posicion)
error=sol(:,posicion(p)+1).-est';
Error(:,p)=error;
hold on
plot(x,error)
end
xlabel('espacio')
ylabel('error cometido')
[0,1,2,10;Error]}}

Obteniendo:

[[Archivo:Estacionaria1.png|marco|centro|Contraste de la gráfica de la estacionaria con las soluciones calculadas mediante numérico en diferentes tiempos.]]

Se puede observar como se reduce el error a medida que aumentamos el tiempo de modelización cor respecto de la estacionaria.

[[Archivo:Estacionaria2.png|marco|centro|Error cometido al aproximar en diferentes tiempos por la estacionaria.]]

En la imagen anterior se puede ver una comparativa del el error con respecto al eje x, donde se ve más claramente que arriba como se reduce.

[[Archivo:Tablaerror.png|marco|centro|Tabla que nos muestra el error cometido en los diferentes tiempos.]]

En los valores numéricos dispuestos en la enterior tabla ya se muestra con total seguridad como se reduce el error a medida que aumentamos el tiempo.

Si ahora por el contrario tratamos las condiciones propuestas en el apartado de cambio de condiciones de contorno, las cuales son:

<math>u_x(0,t)=0</math>:
<math>u(4,t)=0</math>

Analíticamente para estas condiciones la solución estacionaria hallada es::

<math>u_{est}(x)=0</math>

Lo que al principio nos puede chocar al ver la gráfica, ya que en las otras condiciones se llega al estado estacionario de manera mucho más rápida. Si consideramos como estado estacionario que la diferencia entre dos valores consecutivos sea menor de 0.05, mediante el siguiente programa do Octave podemos calcular este valor de tiempo:

{{matlab|codigo=clear all
%datos del problema
L=4;
h=0.1;
T=8;
%discretizacion espacial
N=L/h;
%vector de puntos en el espacio
x=0:h:L;
xint=0:h:L-h;%nodos interiores
%aproximacion de -Uxx
K=2*diag(ones(1,N))-diag(ones(1,N-1),1)-diag(ones(1,N-1),-1);
K=(2/h^2)*K;
%calculamos F
F=zeros(N,1);
F(1)=F(1)+2*5/h^2;
%calculamos U0
U0=((exp((-6)*(xint-2).^2))+(5-5/4.*xint))';
%discretizacion temporal
dt=h; %paso en el espacio
t=0:dt:T;%vector en tiempos
%metodo del trapecio
I=eye(N);
M=I+(dt./2)*K;
sol(:,1)=U0;
for n=1:length(t)-1
sol(:,n+1)=M\(sol(:,n)+dt/2*(F)+dt/2*(-K*sol(:,n)+F));
end
UB=zeros(length(t),1);
sol=[sol;UB'];for i=1:length(t)
er=abs(sol(:,i+1)-sol(:,i));
a=max(er);
if a<=0.05
cont=i;
tiempo=(cont-1)*dt;
break
end
end
tiempo}}

El programa nos dice que el tiempo en el que se cumple esta condición es <math>t=0.7</math>.

=Intercambio de calor con el entorno. Varilla no aislada longitudinalmente=
Si ahora panteamos el intercambio de calor con el exterior, suponiendo que este ambiente se encuentra a una temperatura de 16ºC, la ecuación diferencial se ve afectada de manera que nos queda::

<math>u_t(x,t)-2u_{xx}(x,t)+u(x,t)-16=0</math>

Implementando numéricamente esta ecuación con las condiciones iniciales propuestas en el primer apartado de este artículo, mediante el método de lineas y el método del trapecio, mediante el siguiente código, obteniendo así la gráfica de la temperatura en cada punto de la varilla con respecto al tiempo.

{{matlab|codigo=clear all
%datos del problema
L=4;
h=0.1;
T=10;
%discretizacion espacial
N=L/h;
%vector de puntos en el espacio
x=0:h:L;
xint=h:h:L-h;%nodos interiores
%aproximacion de -Uxx
K=(2+h^2)*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
K=(1/h^2)*K;
%calculamos F
F=zeros(N-1,1).+16;;
F(1)=F(1)+2*5/h^2;
%calculamos U0
U0=((exp((-6)*(xint-2).^2))+(5-5/4.*xint))';
%discretizacion temporal
dt=h; %paso en el espacio
t=0:dt:T;%vector en tiempos
%metodo del trapecio
I=eye(N-1);
M=I+(dt./2)*K;
sol(:,1)=U0;
for n=1:length(t)-1
sol(:,n+1)=M\(sol(:,n)+dt/2*(F)+dt/2*(-K*sol(:,n)+F));
end
UA=5*ones(1,length(t));
UB=0*ones(1,length(t));
sol=[UA;sol;UB];;
%dibujamos la solucion (apartado 3)

[xx,tt]=meshgrid(x,t);

figure(1)
surf(xx,tt,sol')
xlabel('espacio')
ylabel('tiempo')
zlabel('Temperatura')}}

Obteniendo de este código:

[[Archivo:Perdidasaire.png|marco|centro|Representación de la temperatura en cada punto de la varilla con respecto al tiempo]]

Donde se ve como el exterior aporta calor a la varilla elevando sus temperaturas frente a las iniciales.

Si considerásemos el caso estacionario, se puede calcular mediante la inclusión del siguiente fragmento de código en el anterior programa que modelizaba el fenómeno no estacionario. Además de darnos la situación de temperaturas estacionarias, nos da el tiempo en el que este se alcanza con un error relativo de <math>10^{-3}</math>.

{{matlab|codigo=%Estacionaria
figure(2)
xlabel('espacio')
ylabel('temperatura')
c1=det([-11,1;-16,exp(-4)])/det([1,1;exp(4),exp(-4)]);
c2=det([1,-11;exp(4),-16])/det([1,1;exp(4),exp(-4)]);
est=c1*exp(x)+c2*exp(-x)+16;
plot(x,est);
hold on
for i=1:length(t)-1
er=abs(sol(:,i+1)-sol(:,i));
a=max(er);
if a<=0.001
cont=i;
ti=(cont-1)*dt;
break
end
end
ti
cont}}

La gráfica conseguida del código anterior es:

[[Archivo:Perdidasaireest.png|marco|centro|Temperatura estacionaria de la varilla con respecto a la posición de cada punto de ella.]]

El valor de tiempo obtenido es <math>t=5.9</math>.

==Cambio de condiciones de contorno==
Si cambiamos las condiciones de contorno de la ecuación por estas:

<math>u_x(0,t)=2</math>:
<math>u(4,t)=-2</math>

Físicamente podrían interpretarse como la conexión de un deposito de fluido en el extremo derecho que obligara un flujo constante de calor en el mismo extremo, y que el extremo derecho estuviera posicionado sobre un objeto que le obligara a tener una temperatura constante.

Si implementamos numéricamente este nuevo modelo propuesto mediante el siguiente código Matlab, obtenemos la representación de las temperaturas de la varilla frente al tiempo.

{{matlab|codigo=clear all
%datos del problema
L=4;
h=0.1;
T=10;
%discretizacion espacial
N=L/h;
%vector de puntos en el espacio
x=0:h:L;
xint=h:h:L;%nodos interiores
%aproximacion de -Uxx
K=(2+h^2)*diag(ones(1,N))-diag(ones(1,N-1),1)-diag(ones(1,N-1),-1);
K(1,1)=K(1,1)-1;
K=(1/h^2)*K;
%calculamos F
F=zeros(N,1).+16;;
F(1)=F(1)+2/h;
F(N)=F(N)+2/h^2;
%calculamos U0
U0=((exp((-6)*(xint-2).^2))+(5-5/4.*xint))';
%discretizacion temporal
dt=h; %paso en el espacio
t=0:dt:T;%vector en tiempos
%metodo del trapecio
I=eye(N);
M=I+(dt./2)*K;
sol(:,1)=U0;
for n=1:length(t)-1
sol(:,n+1)=M\(sol(:,n)+dt/2*(F)+dt/2*(-K*sol(:,n)+F));
end
UB=2*ones(1,length(t));
sol=[sol;UB];;
%dibujamos la solucion (apartado 3)

[xx,tt]=meshgrid(x,t);

figure(1)
surf(xx,tt,sol')
xlabel('espacio')
ylabel('tiempo')
zlabel('Temperatura')

}}

Logrando así la gráfica:

[[Archivo:Perdidasneu111.png|marco|centro|Representación de las temperaturas de la varilla frente al tiempo]]

Se puede observar como la temperatura de los puntos de la varilla va aumentado debido al flujo entrante por el lado izquierdo de esta hasta que alcanza un estado estacionario.

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