Ecuación del calor (GRwM)
| Trabajo realizado por estudiantes | |
|---|---|
| Título | Ecuación del calor. Grupo GRwM |
| Asignatura | EDP |
| Curso | 2023-24 |
| Autores | Guillermo Gómez Tejedor, Marina Jiménez Barrantes y Rocío Tajuelo Díaz |
| Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura | |
Contenido
- 1 Introducción
- 2 Conceptos previos
- 3 Planteamiento del problema
- 4 Resolución del sistema EDP
- 5 Interpretación del flujo en los extremos
- 6 Solución al considerar otro coeficiente de conductividad
- 7 Solución al considerar otra condición inicial
- 8 Cambio de condiciones frontera
- 9 Soluciones de la ecuación del calor usando la solución fundamental
- 9.1 Solución fundamental de la ecuación del calor en dimensión [math]1[/math]
- 9.2 Ecuación del calor en una dimensión en el semiespacio [math] x\gt0[/math].
- 9.3 Solución del calor con dato inicial [math] u_{0}(x) = 1_{[-1,1]} [/math]
- 9.4 Programa
- 9.5 Solución fundamental de la ecuación del calor en dimensión 2
- 9.6 Programa
- 10 Referencias
1 Introducción
En el trabajo que se presenta a continuación, vamos a estudiar la ecuación del calor en una dimensión. Para ello, vamos a considerar una varilla metálica que se encuentra aislada por su superficie lateral, de modo que la conducción de calor se produzca en la dirección longitudinal.
Partiendo de esto y añadiendo ciertas condiciones de frontera que iremos modificando a lo largo del trabajo, vamos a calcular la solución de la ecuación del calor y la vamos a representar en los distintos casos.
Posteriormente y empleando la solución fundamental de la ecuación del calor, estudiaremos la solución para dimensión 2.
Cabe destacar que los cálculos y visualizaciones han sido programados con MatLab.
2 Conceptos previos
En esta sección, vamos a presentar algunos conceptos esenciales para comprender la obtención del sistema EDP que permite modelizar el problema.
Ley de Fourier: Establece que el flujo de transferencia de calor por conducción es proporcional y de sentido contrario al gradiente de temperatura ([math] T [/math]) en esa dirección. En nuestro caso, al trabajar en una dimensión, [math] \nabla T(\textbf{x},t)= T_{x} [/math] y la ley de Fourier es
siendo [math] q [/math] el flujo de transferencia de calor y [math] k [/math] el coeficiente de conductividad térmica, que es un valor que indica la capacidad de un material para transferir calor a través de él cuando existe una diferencia de temperatura.
Energía térmica: La energía térmica [math] e [/math] es una forma de energía que se produce como resultado del movimiento de átomos en un objeto y es una medida de la cantidad total de energía cinética de estos. La expresión que relaciona la temperatura del objeto con la energía calorífica es:
siendo [math] T [/math] la temperatura y [math] c [/math] el calor específico, que se define como la cantidad de calor que hay que suministrar a la unidad de masa de una sustancia para elevar su temperatura en una unidad.
Principio de conservación de la energía térmica: Establece que la cantidad total de energía térmica para cierto volumen de control permanece constante cuando se tiene en cuenta el flujo de calor entrante y saliente con respecto al tiempo, siempre y cuando no haya ninguna fuente externa que intervenga.
Además, en este trabajo vamos a emplear series de Fourier.
Por otro lado, también vamos a comprobar que el principio del máximo se verifica en este problema.
Principio del máximo: sea [math] u \in C^{2,1} (Q_T) \cap C(\overline{Q_T})[/math] tal que [math]u_t - \Delta u \leq 0 [/math] en [math] Q_T = \Omega \times (0,T)[/math]. Entonces [math]u[/math] alcanza un máximo en la frontera parabólica, es decir, [math] \max \limits_{(x,t) \in \overline{Q_T}} u = \max \limits_{(x,t) \in \partial _P Q_T} u [/math].
Además, en el último ejercicio vamos a emplear la Delta de Dirac. [1]
3 Planteamiento del problema
Para comenzar con el estudio de la ecuación del calor, primero debemos plantear el problema a resolver, que involucra esta ecuación junto a ciertas condiciones de frontera y condición inicial.
Como ya se ha mencionado, en nuestro estudio vamos a considerar una varilla metálica de un metro de largo que se encuentra aislada por su superficie lateral. De esta manera, la conducción de calor se produce únicamente en la dirección longitudinal.
Además, vamos a considerar que la temperatura inicial de la varilla es 0 ºC. También vamos a fijar la temperatura en el extremo izquierdo en 0ºC y en el derecho a 1ºC.
Teniendo en cuenta el principio de conservación la energía y la definición de energía térmica en función de la temperatura, así como las condición inicial y de frontera, se obtiene el siguiente sistema de EDP:
En este primer ejemplo suponemos que tanto la conductividad térmica [math]k[/math] como el calor específico [math]c[/math] toman el valor constante [math]1[/math].
4 Resolución del sistema EDP
Una vez planteado el sistema, procedemos a resolverlo. Para ello, con el objetivo de homogeneizar las condiciones de frontera, vamos a comenzar obteniendo la solución estacionaria.
Para calcular la ecuación del estado estacionario, vamos a tomar el tiempo [math] t \rightarrow \infty [/math]. Esto hace que la variación de la temperatura con respecto al tiempo desaparezca, de modo que la ecuación del sistema es ahora [math] T_{xx}=\frac{\partial ^2 T}{\partial x^2}=0[/math].
Considerando además las condiciones [math] T(0)=0[/math] y [math] T(1)=0[/math], provenientes de las condiciones frontera del problema original, se tiene que la solución de la ecuación del estado estacionario es:
En la siguiente gráfica se representa esta solución:
Consideramos ahora el problema equivalente con condiciones de frontera homogéneas, donde se define [math] T_{hom}(x,t):= T_{est}(x)-T(x,t)[/math]:
Obtenemos la solución de este problema aplicando separación de variables. Para ello, vamos a suponer que la solución es de la forma:
Haciendo las cuentas pertinentes se tiene que
A continuación, se muestra una gráfica con su representación en el intervalo de tiempo [math] t \in [0,1] [/math], tomando los [math] 1000[/math] primeros términos de la serie:
Como podemos observar, al aumentar el tiempo, la solución [math] T(x,t) [/math] tiende a la estacionaria [math] T_{est} =x[/math].
Para ver esto de forma más clara, presentamos a continuación un vídeo que ilustra la evolución de la temperatura sobre la varilla a lo largo del tiempo:
Como se puede apreciar, la variación de la solución es muy notoria al aumentar inicialmente el tiempo, aspecto que se observaba también en la gráfica tridimensional anterior. Sin embargo, para valores de tiempo superiores a [math] 0.5 [/math] esta variación a penas se puede observar, siendo la solución como ya hemos mencionado, muy próxima a la estacionaria.
Además, cabe destacar que, como la solución involucra series de Fourier, el número de términos que consideramos de estas tendrán consecuencias sobre el resultado, concretamente sobre la condición inicial.
Todas las sumas parciales o términos de la serie cumplen la ecuación del calor junto a las condiciones frontera. Sin embargo, solo la serie cumple la condición incial, pero la cumple en el sentido de [math] L^2 [/math], siguiendo la teoría de series de Fourier.
Para poner en manifiesto esto, se presenta a continuación un vídeo en el que se ilustra la solución de la ecuación al variar el número de términos de la serie.
Como podemos observar, al considerar pocos términos de la serie de Fourier es más notorio que no se cumple la condición inicial. Al aumentar el número de términos,
se observa cómo en el instante inicial la solución se aproxima mejor a la condición inicial requerida.
Por otro lado, podemos apreciar que se producen grandes oscilaciones en [math] t=0 [/math], conocidas como fenómeno de Gibbs. Esto se debe a que, la condición inicial que se aproxima mediante las series de Fourier para la solución homogénea, no es continua entendiéndola como una función periódica.
4.1 Programa
A continuación se dejan los diferentes programas de los que obtenemos las gráficas.
1- La solución estacionaria
clear all
close all
%%% SOLUCIÓN ESTACIONARIA
x=linspace(0,1,1000);
plot(x,x, 'm', LineWidth=2)
xlabel('longitud (m)')
ylabel('temperatura (ºC)')
title('Solución de la ecuación del estado estacionario')
2- La solución al problema planteado.
En este programa y en los siguientes la solución se obtiene de forma numérica. Discretizaremos el espacio mediante una malla en la cual obtenemos un valor aproximado de la solución. En los siguientes programas los comentarios indicarán los valores que se pueden modificar y qué cambian.
clc
clear all
close all
%%% REPRESENTACIÓN DE LA ECUACIÓN DEL CALOR (4) %%%
a=0; b=1; % Intervalo de definición de la x ([a,b])
t0=0; tf=1; % Intervalo de definición de la t ([t0,tf])
N=1000; % Número de términos de la serie que aproximan la solución
divx=10^-2; % División del intervalo en x (a más pequeño el valor la malla será más fina en x)
X=a:divx:b; % Vector de x
divt=10^-2; % División del intervalo en t (a más pequeño el valor la malla será más fina en t)
tt=t0:divt:tf; % Vector de t
[Xx,Tt]=meshgrid(X,tt); % Define la malla
T=zeros(length(tt),length(X)); % Inicializamos los valores de la ecuación del calor
for k=1:N % Primero se obtiene la solución homogénea
c=2*(-1)^(k+1)/(k*pi); % Coeficiente de Fourier de la solución
sen = sin(k*pi*X); % Evaluamos sin(kpix) en el vector x
expo = exp(-k^2*pi^2*tt); % Evaluamos exp(-k^2*pi^2*t) en el vector tiempo
T = T + c*expo'*sen; % Añadimos el término k de la serie a la solución
end
Tfinal = ones(length(tt),1)*X - T; % Deshomogeneizamos para obtener la solución requerida
surf(Xx',Tt',Tfinal'); %Visualización de la gráfica
%shading interp (en caso de poner una division más fina se debe de utilizar este comando para que se pueda ver bien)
colormap jet % Hacemos mapa de colores
hold on
plot3(X,tf*ones(length(X)),X,"m","LineWidth",2) % Representamos solución estacionaria
%otros detalles
xlabel('longitud (m)')
ylabel('tiempo (s)')
ylim([t0,tf])
zlabel('T(x,t) (ºC)')
legend('T(x,t)','T_{est}')
title('Solución de la ecuación del calor')
Código de los vídeos:
clc
clear all
close all
%%% VÍDEO DE LA ECUACIÓN DEL CALOR EN FUNCIÓN DEL TIEMPO (Ej 4) %%%
a=0; b=1; % Intervalo de definición de la x ([a,b])
t0=0; tf=1; % Intervalo de definición de la t ([t0,tf])
N=1000; % Número de términos de la serie que aproximan la solución
divx=10^-2; % División del intervalo en x (a más pequeño el valor la malla será más fina en x)
X=a:divx:b; % Vector de x
divt=10^-2; % División del intervalo en t (a más pequeño el valor la malla será más fina en t)
tt=t0:divt:tf; % Vector de t
[Xx,Tt]=meshgrid(X,tt); % Define la malla
T=zeros(length(tt),length(X)); % Inicializamos los valores de la ecuación del calor
for k=1:N
c=2*(-1)^(k+1)/(k*pi); % Coeficiente de Fourier
sen = sin(k*pi*X); % Evaluamos sin(kpix) en el vector x
expo = exp(-k^2*pi^2*tt); % Evaluamos exp(-k^2*pi^2*t) en el vector t
T = T + c*expo'*sen; % Añadimos el término k de la serie
end
Tfinal = ones(length(tt),1)*X - T; % Deshomogeneizamos
pelicula=VideoWriter('Video_Ej_1.4'); % Creamos el video
pelicula.FrameRate=10; % Numero de frames por segundo
open(pelicula);
for j=1:length(tt)
figura = figure(1);
plot(X,Tfinal(j,:),"r","LineWidth",1.5)
hold on
plot(X',X',"k","LineWidth",1)
ylim([0,1])
hold off
xlabel("longitud (m)")
ylabel("temperatura (ºC)")
legend("Solución original","Solución estacionaria")
title("Solución de la Ecuación del Calor")
subtitle("t= "+num2str(tt(j)+" s"))
%legend("solución de la ecuación del calor","solución estacionaria")
imagen=getframe(figura);
writeVideo(pelicula,imagen); % Insertamos la imagen
end
close(pelicula);
clc
clear all
close all
%%% VÍDEO DE LA ECUACIÓN DEL CALOR EN FUNCIÓN DEL NÚMERO DE TÉRMINOS (Ej 4) %%%
a=0; b=1; % Intervalo de definición de la x ([a,b])
t0=0; tf=1; % Intervalo de definición de la t ([t0,tf])
N=100; % Número de términos de la serie que aproximan la solución
divx=10^-2; % División del intervalo en x (a más pequeño el valor la malla será más fina en x)
X=a:divx:b; % Vector de x
divt=10^-2; % División del intervalo en t (a más pequeño el valor la malla será más fina en t)
tt=t0:divt:tf; % Vector de t
[Xx,Tt]=meshgrid(X,tt); % Define la malla
T=zeros(length(tt),length(X)); % Inicializamos los valores de la ecuación del calor
pelicula=VideoWriter('Video 3 Ej 1.4.2'); % Creamos el video
pelicula.FrameRate=10;
open(pelicula);
for k=1:N
c=2*(-1)^(k+1)/(k*pi); % Coeficiente de Fourier
sen = sin(k*pi*X); % Evaluamos sin(kpix) en el vector x
expo = exp(-k^2*pi^2*tt); % Evaluamos exp(-k^2*pi^2*t) en el vector t
T = T + c*expo'*sen; % Añadimos el término k de la serie
Taux=ones(length(tt),1)*X - T; % Dehomogeneizamos
figura = figure(1);
surf(Xx',Tt',Taux');
colormap jet
%shading interp
zlim([-0.2,1])
ylim([t0,tf])
xlabel("longitud (m)")
ylabel("tiempo (s)")
zlabel("temperatura (ºC)")
title("Solución de la Ec del Calor")
subtitle("Aproximación hasta término "+ num2str(k))
imagen=getframe(figura);
writeVideo(pelicula,imagen); % Insertamos la imagen
end
close(pelicula);
5 Interpretación del flujo en los extremos
Para finalizar el estudio de la solución del sistema EDP, vamos obtener e interpretar el flujo de calor saliente y entrante en ambos extremos a lo largo del tiempo.
Para ello, por la propia definición de flujo, debemos obtener la derivada de la solución con respecto a la variable x, que es:
Además, para su interpretación, debemos tener en cuenta que el flujo es positivo cuando su sentido es el del eje longitudinal, es decir, de izquierda a derecha. De modo que será entrante o saliente dependiendo de su signo en los extremos. Calculamos así la derivada en los extremos:
Es importante recalcar que en el instante [math] t=0 [/math] la derivada no está definida, ya que, como hemos visto, para este valor del tiempo hay una discontinuidad. También se puede observar analíticamente pues la serie no converge.
En este caso, oscila entre [math] -1 [/math] y [math] 1 [/math] dependiendo de si consideramos un número par o impar de términos, respectivamente.
Y en este caso diverge. Por esta razón, para el estudio del flujo, consideraremos tiempos estrictamente positivos.
En la siguiente gráfica se ilustra el flujo para 1000 términos de la serie:
Como podemos observar, el flujo en ambos extremos converge con el tiempo al valor [math]-1 ~ \frac{Cal}{m^2s}[/math]. Este resultado no es arbitario, es el flujo de la solución estacionaria.
Además, que sea un valor negativo implica que, en el extremo izquierdo el flujo sea saliente, y en el derecho sea entrante. Este resultado es coherente con el problema de estudio pues, por definición, el flujo va del punto con mayor temperatura al de menor. Por tanto, como uno de los extremos, el derecho, no está homogeneizado y su temperatura es superior, el flujo irá de ese punto al extremo izquierdo.
5.1 Programa
clc
clear all
close all
%%% REPRESENTACIÓN DEL FLUJO EN LA FRONTERA (5)%%%
a=0; b=1; % Intervalo de definición de la x ([a,b])
t0=10^-3; tf=2; % Intervalo de definición de la t ([t0,tf])
N=1000; % Número de términos de la serie que aproximan la solución
X=[a b]; % Vector de extremos del intervalo en x
div=10^-3; % Número de divisiones del intervalo t
tt=t0:div:tf; % Vector de tiempo
Tx=zeros(length(tt),length(X)); %Inicializamos los valores de la ecuación del calor
for k=1:N % Calculamos la solución homogénea
c=2*(-1)^(k+1)/(k*pi); % Coeficiente de Fourier
kpcose = (k*pi )*cos(k*pi*X); % Derivada con respecto a x de sin(kpix) evaluado
expo = exp(-k^2*pi^2*tt); % Calculamos exp(-k^2*pi^2*t) evaluada en el vector de t
Tx = Tx + c*expo'*cose ; % Añadimos el término k de la serie
end
Flujo = Tx - ones(length(tt),2); % Deshomogeneizamos
% Dibujamos
figure(1)
plot(tt,Flujo(:,1),"b",LineWidth=1.5)
hold on
plot(tt,Flujo(:,2),"r",LineWidth=1.5)
hold off
legend("x=0","x=1")
xlabel("tiempo (s)")
ylabel("Flujo (Cal/(m^2s))")
ylim([-6 1])
title("Flujo en los extremos de la barra")
6 Solución al considerar otro coeficiente de conductividad
Si consideramos ahora el coeficiente de conductividad térmica igual a [math]\frac{1}{2}[/math], se tiene el siguiente problema:
Siguiendo los mismos pasos que el problema inicial, obtenemos que la solución es
En las siguientes gráficas, se muestran las soluciones de los problemas considerando [math]k=1[/math] y [math]k=1/2[/math]:
A partir de estas gráficas, el único aspecto destacable es la diferencia en la pendiente de las soluciones para los instantes iniciales. Además, para el tiempo considerado, visualmente ambas soluciones convergen a la estacionaria. Sin embargo, no podemos distinguir la velocidad con la cual lo hacen. Para poder observar esto, vamos a dibujar la diferencia entre la solución y el estado estacionario en [math]x= \frac{1}{2} [/math] para [math] t \in [0, 1][/math], es decir, [math] u(1/2, t) - u_{est}(1/2)[/math].
En esta gráficas se observa que para [math] k=1[/math] la diferencia entre la solución estacionaria y la original en [math] x=\frac{1}{2}[/math] es menor que al considerar [math] k=\frac{1}{2}[/math]. Por lo que podemos intuir que la solución con [math] k=1[/math] converge más rápido a la estacionaria.
Para terminar de comparar ambos problemas y confirmar esta idea, vamos a estudiar el flujo. Lo representamos en la siguiente gráfica:
Vemos que, al considerar un coeficiente de conductividad menor, en este caso [math]k=\frac{1}{2}[/math], aumenta el tiempo necesario para que se alcance el valor [math] -1[/math], que es el flujo al considerar la solución estacionaria. Esto indica que efectivamente, la solución para [math] k=1[/math] converge más rápido a la estacionaria.
Este resultado se debe a que, como se mencionaba en la sección de conceptos previos, la constante [math] k[/math] indica la capacidad del material para transferir calor a través de él cuando existe una diferencia de temperatura. Por tanto, si este valor es menor, el tiempo en alcanzar el valor del flujo de la solución estacionaria es mayor.
Para visualizar la convergencia, se presenta a continuación el siguiente vídeo:
6.1 Programa
Gráfica de solución y diferencia:
clc
clear all
close all
%%% REPRESENTACIÓN DE LA ECUACIÓN DEL CALOR con k=1/2 y k=1 y diferencias(6) %%%
a=0; b=1; % Intervalo de definición de la x ([a,b])
t0=0; tf=1; % Intervalo de definición de la t ([t0,tf])
N=1000; % Número de términos de la serie que aproximan la solución
divx=10^-2; % División del intervalo en x (a más pequeño el valor la malla será más fina en x)
X=a:divx:b; % Vector de x
divt=10^-2; % División del intervalo en t (a más pequeño el valor la malla será más fina en t)
tt=t0:divt:tf; % Vector de t
[Xx,Tt]=meshgrid(X,tt); % Define la malla
T=zeros(length(tt),length(X)); % Inicializamos los valores de la ecuación del calor con k=1
T2=zeros(length(tt),length(X)); % Inicializamos los valores de la ecuación del calor con k=1/2
for k=1:N
c=2*(-1)^(k+1)/(k*pi); % Calculamos los coeficientes de Fourier
sen = sin(k*pi*X); % Evaluamos sin(kpix) en el vector de X
expo2 = exp(-k^2*pi^2*tt/2); % Evaluamos en exp(-k^2*pi^2*t/2) (para el coef k=1/2)
expo = exp(-k^2*pi^2*tt); % Evaluamos en exp(-k^2*pi^2*t) (para el coef k=1)
% Añadimos el término k de la serie a ambas soluciones
T2 = T2 + c*expo2'*sen;
T = T + c*expo'*sen;
end
% Deshomogeneizamos
T = ones(length(tt),1)*X - T;
T2 = ones(length(tt),1)*X - T2;
N=(b-a)/divx/2 +1; %Valor en 0.5 (si se pone N impar dará fallos)
medio=1/2*ones(length(tt),1); % Vector de 1/2
error=abs(T(:,N)- medio); % Error respecto a la solución estacionaria en x=1/2 con k=1
error2=abs(T2(:,N)- medio); % Error respecto a la solución estacionaria en x=1/2 con k=1/2
% Dibujamos la solución
figure(1)
subplot(1,2,1)
surf(Xx',Tt',T');
%shading interp % Poner si se disminuye div
colormap jet
hold on
plot3(X,tf*ones(length(X)),X,"m","LineWidth",2)
zlim([0 1])
hold off
xlabel('longitud (m)')
ylabel('tiempo (s)')
zlabel('T(x,t) (ºC)')
title('Solución para k=1 en t \in [0,1]')
legend('T(x,t)','T_{est}')
subplot(1,2,2)
surf(Xx',Tt',T2');
colormap jet
hold on
plot3(X,tf*ones(length(X)),X,"m","LineWidth",2)
zlim([0 1])
hold off
xlabel('longitud (m)')
ylabel('tiempo (s)')
zlabel('T(x,t) (ºC)')
title('Solución para k=1/2 en t \in [0,1]')
legend('T(x,t)','T_{est}')
%Dibujamos el error en x=1/2 en escala logarítmica
figure(2)
semilogy(tt,error,"r",LineWidth=2)
hold on
semilogy(tt,error2,"b",LineWidth=2)
ylabel("error (ºC)")
xlabel("tiempo (s)")
title('
Gráfica de comparación de flujos
clc
clear all
close all
%%% COMPARACIÓN DEL FLUJO EN LA FRONTERA (6) %%%
a=0; b=1; % Intervalo de definición en x
t0=10^-3; tf=2; % Intervalo de definición en t
N=1000; % Término hasta el que se aproxima la serie
X=[a b]; % Vector de extremos
div=10^-3; % División del intervalo en t
tt=t0:div:tf; % Vector de t
Tx1=zeros(length(tt),length(X)); % Inicializamos los valores de la ecuación del calor para k=1
Tx2=zeros(length(tt),length(X)); % Inicializamos los valores de la ecuación del calor para k=1/2
for k=1:N
c=2*(-1)^(k+1)/(k*pi); % Coeficiente de Fourier
cose = k*pi*cos(k*pi*X); % Evaluamos (d/dx)sin(kpix) en X
expo = exp(-k^2*pi^2*tt); % Evaluamos exp(-k^2*pi^2*t) en t (k=1)
expo2 = exp(-k^2*pi^2*tt/2); % Evaluamos exp(-k^2*pi^2*t) en t (k=1/2)
% Añadimos el término k
Tx1 = Tx1 + c*expo'*cose;
Tx2 = Tx2 + c*expo2'*cose;
end
% Deshomogeneizamos
Flujo = Tx1 - ones(length(tt),2);
Flujo2 = Tx2 - ones(length(tt),2);
% Dibujamos
figure(1)
plot(tt,Flujo(:,1),"b","LineWidth",1.5)
hold on
plot(tt,Flujo(:,2),"r","LineWidth",1.5)
hold on
plot(tt,Flujo2(:,1),"g","LineWidth",1.5)
hold on
plot(tt,Flujo2(:,2),"m","LineWidth",1.5)
hold off
legend("x=0,k=1","x=1,k=1","x=0,k=1/2","x=1,k=1/2")
xlabel("tiempo (s)")
ylabel("flujo (Cal/(m^2s))")
ylim([-3 1])
title("Flujo en los extremos")
Vídeo de comparación de soluciones para [math]k=1[/math] y [math]k=\frac{1}{2}[/math].
clc
clear all
close all
%%% VÍDEO DE COMPARACIÓN DE SOLUCIONES PARA K=1 Y K=\FRAC %%%
a=0; b=1; % Intervalo de definición en x
t0=0; tf=1; % Intervalo de definición en t
N=1000; % Número de términos de la serie
divx=10^-2; % División del intervalo en x
X=a:divx:b; % Vector de x
div=10^-2; % División del intervalo en t
tt=t0:div:tf; % Vector de t
[Xx,Tt]=meshgrid(X,tt);
T1=zeros(length(tt),length(X));
T2=zeros(length(tt),length(X)); % Valores de la ecuación del calor
for k=1:N
c=2*(-1)^(k+1)/(k*pi); % Coeficiente de Fourier
sen = sin(k*pi*X); % Evaluamos el vector de x en sin(kpix)
expo1 = exp(-k^2*pi^2*tt); % Evaluamos el vector de t en exp(-k^2*pi^2*t)
expo2 = exp(-k^2*pi^2*tt/2); % Evaluamos el vector de t en exp(-k^2*pi^2*t/2)
T1 = T1 + c*expo1'*sen; % Añadimos el término k
T2 = T2 + c*expo2'*sen;
end
% Dehomogeneizamos
Tfinal1 = ones(length(tt),1)*X - T1;
Tfinal2 = ones(length(tt),1)*X - T2;
pelicula=VideoWriter('Video_Ej_1.6'); % Creamos el video
pelicula.FrameRate=10;
open(pelicula);
for j=1:length(tt)
figura = figure(1);
plot(X,Tfinal1(j,:),"r","LineWidth",1.5)
hold on
plot(X,Tfinal2(j,:),"b","LineWidth",1.5)
plot(X',X',"k","LineWidth",1)
ylim([0,1])
hold off
xlabel("longitud (m)")
ylabel("temperatura (ºC)")
legend("Solución con k=1", "Solución con k=1/2","Solución estacionaria")
title("Solución de la Ecuación del Calor")
subtitle("t= "+num2str(tt(j)+" s"))
%legend("solución de la ecuación del calor","solución estacionaria")
imagen=getframe(figura);
writeVideo(pelicula,imagen); %inserto la imagen
end
close(pelicula);
7 Solución al considerar otra condición inicial
En esta sección, vamos a resolver el problema de la ecuación del calor considerando otra condición inicial.
7.1 Replanteamos el problema
Ahora vamos a suponer que la temperatura en el extremo derecho es también [math]0[/math]ºC, es decir, tenemos condiciones de frontera homogéneas. Además, vamos a cambiar la condición inicial por una nueva función. Mantendremos tanto la conductividad térmica [math]k[/math] como el calor específico [math]c[/math] con el valor constante [math]1[/math], siendo así el sistema EDP:
7.2 Resolución del sistema EDP
Al haber considerado condiciones de frontera nulas se tiene que la solución estacionaria del problema es [math]0[/math]. Por tanto, no es necesario el proceso de homogeneización.
Aplicando el método de separación de variables obtenemos que la solución de la ecuación del calor es de la forma:
En este caso, en lugar de obtener los coeficientes de Fourier de forma analítica, vamos a aproximar las integrales por la fórmula del trapecio.
Teniendo esto en cuenta, la gráfica de la solución es:
A continuación, mostramos un vídeo en el que representamos la solución en función del tiempo:
Como podemos observar, la solución a medida que avanza el tiempo se aproxima al estado estacionario, que en este caso es la función nula.
Vemos que hay tramos en el instante [math] t=0[/math] donde la temperatura es nula, concretamente en [math] \left[0,\frac{1}{4} \right][/math] y [math] \left[ \frac{3}{4},1 \right][/math]. Además, se tiene que, en cuanto [math] t[/math] aumenta, la temperatura a lo largo de toda la varilla pasa a ser estrictamente positiva. Es decir, el calor que se concentraba en el intervalo [math] \left[ \frac{1}{4},\frac{3}{4} \right][/math] en el instante inicial, se traslada al resto de la varilla de forma instantánea. Esto significa que el calor ha recorrido dos intervalos de longitud [math] 0.25 ~ m[/math] instantaneamente. Por tanto, podemos interpretar que la difusión ’transporta’ el calor con velocidad infinita.
7.2.1 Programa
Hay que destacar que en este caso no hemos dado valores concretos en función de [math] k[/math] de los coeficientes de Fourier, esto se debe a que el valor es más rápido y cómodo de obtener numéricamente.
clc
clear all
close all
%%% REPRESENTACIÓN DE LA ECUACIÓN DEL CALOR (7) %%%
a=0; b=1; % Intervalo de definición de x ([a,b])
t0=0; tf=0.5; % Intervalo de definición de t ([t0,tf])
N=1000; % Número de términos considerados en la serie
divx=10^-3; % División del intervalo en x
X=a:divx:b; % Vector de x
divt=10^-3; % División del intervalo en t
tt=t0:divt:tf; % Vector de t
[Xx,Tt]=meshgrid(X,tt); % Creamos la malla
f=@(x)(1-4*abs(x-1/2)).*((1-4*abs(x-1/2))>0); % Condición inicial
Y=f(X); % Evaluamos la condición inicial en x
T=zeros(length(tt),length(X)); % Inicializamos los valores de la ecuación del calor
for k=1:N
sen = sin(k*pi*X); % Evaluamos sin(kpix) en el vector de x
c=2*trapz(X,Y.*sen); % Coeficiente de Fourier aproximados por la regla del trapecio divididos por la norma del seno que es 1/2.
expo = exp(-k^2*pi^2*tt); % Evaluamos exp(-k^2*pi^2*t) en el vector de t
T = T + c*expo'*sen; % Añadimos el término k de la serie
end
% Dibujar las gráficas
surf(Xx',Tt',T');
colormap jet
shading interp
zlim([0,1])
ylim([0,tf])
xlabel('longitud (m)')
ylabel('tiempo (s)')
zlabel('T(x,t) (ºC)')
title('Solución en t \in [0,1]')
legend('T(x,t)')
Código del vídeo:
clc
clear all
close all
%%% VÍDEO DE LA ECUACIÓN DEL CALOR (7) %%%
a=0; b=1; % Intervalo de definición de x ([a,b])
t0=0; tf=1; % Intervalo de definición de t ([t0,tf])
N=1000; % Número de términos de la serie
divx=10^-3; % División del intervalo en x y t
X=a:divx:b; % Vector de x
div=10^-2; % División del intervalo en x y t
tt=t0:div:tf; % Vector de t
[Xx,Tt]=meshgrid(X,tt); % Hacemos la malla
f=@(x)(1-4*abs(x-1/2)).*((1-4*abs(x-1/2))>0); % Condición inicial
Y=f(X); % Evaluamos x en la condición inicial
T=zeros(length(tt),length(X)); % Valores de la ecuación del calor
for k=1:N
sen = sin(k*pi*X); % Evaluamos el vector x en sin(kpix)
c=2*trapz(X,Y.*sen); % Coeficiente de Fourier
expo = exp(-k^2*pi^2*tt); % Evaluamos t en exp(-k^2*pi^2*t)
T = T + c*expo'*sen; % Añadimos el término k
end
pelicula=VideoWriter('Video Ej 1.7'); % Creamos el video
pelicula.FrameRate=10;
open(pelicula);
for j=1:length(tt)
figura = figure(1);
plot(X,T(j,:),"r","LineWidth",1.5)
ylim([0,1])
hold off
xlabel("longitud (m)")
ylabel("temperatura (ºC)")
title("Solución de la Ecuación del Calor")
subtitle("t="+num2str(tt(j))+" s")
%legend("solución de la ecuación del calor","solución estacionaria")
imagen=getframe(figura);
writeVideo(pelicula,imagen); % Insertamos la imagen
end
close(pelicula);
7.3 Interpretación del flujo en los extremos
Para finalizar el estudio de la solución del sistema, vamos obtener e interpretar el flujo de calor saliente y entrante en ambos extremos a lo largo del tiempo.
Para ello obtenemos la derivada de la solución con respecto a la variable x, que es:
La gráfica que muestra la evolución del flujo en los extremos de la barra con respecto al tiempo es:
Como podemos observar, el valor del flujo en ambos extremos tiende a cero, que es el valor del flujo para la solución estacionaria y visualmente parece que converge a este valor en torno a [math] 0.7 ~ s [/math]. También hay que destacar que en el extremo derecho el flujo es positivo mientras que en el izquierdo es negativo. Como ya hemos explicado, esto se interpreta como que en ambos extremos es saliente, aunque tiende a ser nulo.
7.3.1 Programa
Código de representación del flujo en la frontera:
clc
clear all
close all
%%% REPRESENTACIÓN DEL FLUJO EN LA FRONTERA (7)%%%
a=0; b=1; % Intervalo de definición de x ([a,b])
t0=0; tf=1; % Intervalo de definición de t ([t0,tf])
N=1000; % Número de términos de la serie
X=[a b]; % Vector de extremos del intervalo de definición de x
div=10^-3; % División del intervalo en t
tt=t0:div:tf; % Vector de t
f=@(x)(1-4*abs(x-1/2)).*((1-4*abs(x-1/2))>0); % Condición inicial
xx=a:div:b; % Vector en X para los coeficientes de Fourier
Y=f(xx); % Evaluamos el vector x en la condición incial
Tx=zeros(length(tt),length(X)); % Valores de la ecuación del calor
for k=1:N
sen = sin(k*pi*xx); % Evaluamos sin(kpix) en a y b
c=2*trapz(xx,Y.*sen); % Calculamos los coeficientes de Fourier
cose = k*pi*cos(k*pi*X); % Evaluamos (d/dx)sin(kpix) en a y b
expo = exp(-k^2*pi^2*tt); % Evaluamos exp(-k^2*pi^2*t) en t
Tx = Tx + c*expo'*cose; % Añadimos el término k a la solución
end
% Dibujamos
figure(1)
plot(tt,-Tx(:,1),"b",LineWidth=1.5)
hold on
plot(tt,-Tx(:,2),"r",LineWidth=1.5)
hold off
title("extremo derecho x=1")
legend("x=0","x=1")
xlabel("tiempo (s)")
ylabel("flujo (Cal/(m^2s))")
title("Flujo en los extremos")
8 Cambio de condiciones frontera
Ahora suponemos que el extremo derecho de la barra está aislada térmicamente, en lugar de que la temperatura sea de [math]0 ~ ºC[/math]. Entonces, tenemos el siguiente problema:
que tiene condiciones de frontera nulas y por tanto, no es necesario el proceso de homogeneización.
A continuación, obtenemos la solución de este problema aplicando separación de variables. En este caso, obtenemos que la solución es una combinación lineal de funciones de la forma [math]\sin{((\frac{\pi}{2}+k\pi)x)} e^{-(\frac{\pi}{2}+k\pi)^2t}[/math] para [math]k=0, ~..., ~\infty[/math]. Con esta solución, no podemos emplear la base trigonométrica usual [math]\{ \frac{1}{2}, sen(k \pi x), cos(k \pi x) \}_{k=1}^{\infty} [/math]. Para probar que existe una solución de la ecuación del calor con la condición inicial pedida, debemos probar que [math]\{sen((\frac{\pi}{2}+k\pi)x)\}_{k=1}^{\infty} [/math] es una base en [math]L^2([0,1])[/math].
Supondremos que se trata de una familia completa , faltaría comprobar que son ortogonales.
Teniendo esto en cuenta y siendo [math]c_k[/math] los coeficientes de Fourier en dicha base de la condición inicial, la solución al nuevo problema es:
Al igual que en el caso anterior obtendremos los coeficientes de forma numérica mediante el método del trapecio.
Podemos visualizar la gráfica:
Lo más destacable de esta solución respecto a la anterior es cómo varía la temperatura en el extremo derecho, pudiéndose apreciar que incide perpendicularmente al plano vertical en la gráfica. Esto es el significado visual de tener la derivada en x igual a 0. Para verlo mejor observemos el siguiente video:
Y en efecto, podemos observar como la gráfica incide perpendicularmente en el tramo vertical de [math]x=1[/math].
8.1 Principio del máximo
También es importante observar cómo, en efecto, la solución de la ecuación del calor verifica el principio del máximo, pues sabemos que la solución [math] T \in C^{2,1} (Q_T) \cap C(\overline{Q_T})[/math].
Tras observar las gráficas en los dos casos anteriores que el máximo se alcanza en la frontera parabólica. En estos dos últimos problemas, se alcanza concretamente en [math](x,t)=\left( \frac{1}{2},0 \right)[/math].
8.2 Programa
Código de la representación de la solución de la ecuación del calor
clc
clear all
close all
%%% REPRESENTACIÓN DE LA ECUACIÓN DEL CALOR (8) %%%
a=0; b=1; % Intervalo de definición de x ([a,b])
t0=0; tf=0.5; % Intervalo de definición de t ([t0,tf])
N=1000; % Término hasta el que se aproxima la serie
divx=10^-3; % División del intervalo en x
X=a:divx:b; % Vector de x
divt=10^-3; % División del intervalo en t
tt=t0:divt:tf; % Vector de t
[Xx,Tt]=meshgrid(X,tt); % Creamos la malla de puntos
f=@(x)(1-4*abs(x-1/2)).*((1-4*abs(x-1/2))>0); % Condición incial
Y=f(X); % Evaluamos la condición inicial en X
T=zeros(length(tt),length(X)); % Inicializamos los valores de la ecuación del calor
for k=0:N
sen = sin((pi/2+k*pi)*X); % Evaluamos el seno en x
c=trapz(X,Y.*sen)*2; % Calculamos los coeficientes de Fourier
expo = exp(-(pi/2+k*pi)^2*tt); % Evaluamos la exponencial en t
T = T + c*expo'*sen; % Añadimos el término k de la serie
end
% Dibujamos
surf(Xx',Tt',T');
colormap jet
shading interp
zlim([0,1])
ylim([0,tf])
xlabel('longitud (m)')
ylabel('tiempo (s)')
zlabel('T(x,t) (ºC)')
title('Solución de la ecuación del calor')
legend('T(x,t)')
Código del vídeo:
clc
clear all
close all
%%% VIDEO DE LA ECUACIÓN DEL CALOR (8) %%%
a=0; b=1; % Intervalo de definición de x ([a,b])
t0=0; tf=1; % Intervalo de definición de t ([t0,tf])
N=1000; % Término hasta el que se aproxima la serie
divx=10^-3; % División del intervalo en x
X=a:divx:b; % Vector de x
divt=10^-2; % División del intervalo en t
tt=t0:divt:tf; % Vector de t
[Xx,Tt]=meshgrid(X,tt); % Creamos la malla de puntos
f=@(x)(1-4*abs(x-1/2)).*((1-4*abs(x-1/2))>0); % Condición incial
Y=f(X); % Evaluamos la condición inicial en X
T=zeros(length(tt),length(X)); % Inicializamos los valores de la ecuación del calor
for k=0:N
sen = sin((pi/2+k*pi)*X); % Evaluamos el seno en x
c=trapz(X,Y.*sen)*2; % Calculamos los coeficientes de Fourier
expo = exp(-(pi/2+k*pi)^2*tt); % Evaluamos la exponencial en t
T = T + c*expo'*sen; % Añadimos el término k de la serie
end
pelicula=VideoWriter('Video Ej 1.8'); % Creamos el video
pelicula.FrameRate=10;
open(pelicula);
for j=1:length(tt)
figura = figure(1);
plot(X,T(j,:),"r","LineWidth",1.5)
ylim([0,1])
hold off
xlabel("longitud (m)")
ylabel("temperatura (ºC)")
title("Solución de la Ecuación del Calor")
subtitle("t="+num2str(tt(j))+" s")
%legend("solución de la ecuación del calor","solución estacionaria")
imagen=getframe(figura);
writeVideo(pelicula,imagen); %Insertamos la imagen
end
close(pelicula);
9 Soluciones de la ecuación del calor usando la solución fundamental
En esta sección vamos a estudiar la solución de la ecuación del calor en dimensión [math] n[/math] empleando la solución fundamental, que viene dada por la siguiente expresión:
donde [math]D=\frac{k}{c}[/math].
Esta se trata de la solución de la ecuación de difusión en todo [math]\mathbb{R}^n[/math]:
Siendo [math] \delta(x) [/math] la 'Delta de Dirac' en [math]\mathbb{R}^n[/math].
La importancia de esta solución viene a partir de la siguiente propiedad de la convolución junto a la delta de Dirac: [math] f \ast \delta = f [/math], para una [math] f [/math] suficientemente regular. Pues se verifica que [math] \Phi_D \ast_x f [/math] cumple la ecuación anterior junto a que la condición inicial sea: [math] \Phi_D(x,0) \ast_x f = \delta \ast f = f[/math]. Es decir, una solución de la ecuación de la difusión en [math]\mathbb{R}^n[/math] con una condición inicial cualquiera [math]f[/math] se obtiene haciendo la convolución con la solución fundamental.
9.1 Solución fundamental de la ecuación del calor en dimensión [math]1[/math]
Como ejemplo en dimensión [math]n=1[/math], vamos a considerar [math]x \in [-1,1][/math], [math]t \in [10^{-2},1][/math] y [math]k=c=1[/math], es decir, [math]D=1[/math]. Por tanto, la expresión de la solución fundamental queda de la siguiente manera:
A continuación podemos ver la representación gráfica de dicha función:
Nota: La singularidad de la delta de Dirac no nos permite representar la función en el instante [math]t=0[/math].
9.1.1 Programa
clc
clear all
close all
%%%% SOLUCIÓN FUNDAMENTAL DIM 1
a=-1; b=1; % Intervalo de definición de x ([a,b])
t0=10^-2; tf=1; % Intervalo de definición de t ([t0,tf])
divx=10^-3; divt=10^-3; % División del intervalo en x y t
X=a:divx:b; T=t0:divt:tf; % Vector de x y de t
Phi = @(x,t)(1./((4*pi*t).^(1/2)).*exp(-abs(x).^2./(4*t))); % Solución fundamental para dimensión 1
[XX,TT]=meshgrid(X,T); % Hacemos la malla
PHI = Phi(XX,TT); % Evaluamos la solución fundamental en la malla
% Dibujamos
surf(XX,TT,PHI);
shading interp
colormap jet
xlabel('x')
ylabel('tiempo t')
zlabel('\Phi_1(x,t)')
title('Solución fundamental de la ecuación del calor en dim=1')
legend('\Phi_1(x,t)')
9.2 Ecuación del calor en una dimensión en el semiespacio [math] x\gt0[/math].
En esta sección, vamos a considerar el problema:
La forma de resolver esta ecuación diferencial es muy diferente, y consiste en suponer que [math] T(x,t)= u\left(\frac{x}{\sqrt{t}}\right) [/math].
Recalculamos la ecuación diferencial para [math]u[/math]:
Llamamos [math] \zeta = \frac{x}{\sqrt{t}} [/math] y obtenemos que [math]u[/math] cumple la siguiente ecuación diferencial:
Resolviendo la ecuación diferencial, obtenemos que [math] u(\zeta)=\int^{\zeta}_0e^{-\frac{s^2}{4}}ds = \sqrt{\pi}\left(\frac{2}{\sqrt{\pi}}\int^{\frac{\zeta}{2}}_0e^{-z^2}dz \right) = \sqrt{\pi}~erf \left(\frac{\zeta}{2} \right)[/math] luego [math] T_1(x,t)=u\left(\frac{x}{\sqrt{t}}\right)= \sqrt{\pi} ~ erf\left(\frac{x}{2\sqrt{t}}\right) [/math] (siendo [math]erf(x)=\frac{2}{\sqrt{\pi}}\int^x_0e^{-z^2}dz[/math]). Esta solución cumple:
Por tanto, para obtener la solución que cumpla las condiciones propuestas en el problema original, debemos hacer el siguiente cambio de variable: [math]T(x,t) = 1-\frac{1}{\sqrt{\pi}}T_1(x,t)[/math] . Luego la solución que buscamos es:
Para representar la función, dado que la [math]x[/math] varía en todo [math]\mathbb{R}^+[/math] debemos elegir un subintervalo concreto. Para ello observaremos, a partir de qué valores en la [math]x[/math] fijado el intervalo de tiempo , la función toma valores despreciables.
Siendo [math]10^{-\alpha}[/math] el valor de significación que consideremos.
Nótese que se ha hecho un cambio de variable a polares. La última desigualdad se debe a que la región de integración de la última integral es mayor que la de la anterior.
Por tanto con [math]z[/math] mayor a ese valor, la solución ha de ser menor que [math]10^{-\alpha}[/math]. Como [math]z=\frac{x}{2\sqrt{t}}[/math] la cota correcta para [math]x[/math] será:
Siendo [math]max(t)[/math] el máximo valor de tiempo que queramos considerar.
Obtenemos la siguiente la gráfica:
Se puede observar cómo la función cumple las condiciones del problema, y en efecto, en el extremo derecho la solución parece tener un valor muy próximo a 0.
Es también importante observar qué ocurre con la solución cuando [math]t\rightarrow \infty[/math]. Analíticamente, observamos que [math]\lim\limits_{t\rightarrow \infty} T(x,t) = 1 – erf( 0 ) = 1 [/math]. Es decir, parece que en todos los puntos del semieje positivo la temperatura acabará teniendo el valor [math]1[/math].
Veamos si esto se cumple visualmente con el siguiente vídeo:
Observamos como, efectivamente, la temperatura parece tender a [math]1[/math] cuando el tiempo aumenta.
9.2.1 Programa
clc
clear all
close all
divt=10^-3; % División del intervalo en t
t0=10^-3; tf=1; % Intervalo de definición en t ([t0,tf])
T=t0:divt:tf; % Vector de t
imp=6;
divx=10^-2; % División del intervalo en x
a=0; b=2*sqrt(tf)*sqrt(log(sqrt(pi)/2) + imp*log(10)); % Intervalo de definición en x ([a,b])
X=a:divx:b; % Vector de x
[XX,TT]=meshgrid(X,T); % Hacemos la malla
w=@(x,t)(1-erf(x./(sqrt(t)))); % Definimos la solución
W=w(XX,TT); % Evaluamos la solución en la malla
% Dibujamos
surf(XX,TT,W)
shading interp
colormap jet
xlim([0,b])
ylim([t0,tf])
xlabel("longitud (m)")
ylabel("tiempo (s)")
zlabel("temperatura (ªC)")
title("Solución de la ecuación de calor en el semiespacio x>=0")
clc
clear all
close all
divt=1; % División del intervalo de tiempo
t0=0; tf=500; % Intervalo de definición de t ([t0,tf])
T=t0:divt:tf; % Vector de t
divx=10^-2; % División del intervalo de x
a=0; b=5; % Intervalo de definición de x ([a,b])
X=a:divx:b; % Vector de x
[XX,TT]=meshgrid(X,T); % Definimos la malla
w=@(x,t)(1-erf(x./(sqrt(t)))); % Definimos la solución
W=w(XX,TT); % Evaluamos la solución en la malla
pelicula=VideoWriter('Video 3.2'); % Creo el video
pelicula.FrameRate=50; %Velocidad del vídeo
open(pelicula);
for j=1:length(T)
figura = figure(1);
plot(X,W(j,:),"r","LineWidth",1.5)
hold off
ylim([0,1.1])
xlabel("longitud (m)")
ylabel("temperatura (ºC)")
title("Solución de la Ecuación del Calor")
subtitle("t= "+num2str(T(j)+" s"))
imagen=getframe(figura);
writeVideo(pelicula,imagen); %inserto la imagen
end
close(pelicula);
9.3 Solución del calor con dato inicial [math] u_{0}(x) = 1_{[-1,1]} [/math]
En este caso, buscamos la solución en dimensión [math] 1[/math] con condición inicial [math] u_{0}(x) = 1_{[-1,1]}[/math]. Para ello, como hemos explicado anteriormente, podemos emplear la solución fundamental. De modo que para este problema, la solución será de la forma:
Para poder interpretarla, vamos a representarla en diferentes tiempos: [math] t= 0.1[/math], [math] t= 0.01[/math] y [math] t= 0.001[/math]. Pero antes de poder hacerlo, debemos considerar un intervalo adecuado para [math] x[/math], ya que no podemos representar todo [math] \mathbb{R}[/math].
Para calcular este intervalo, vamos a buscar el valor de [math] x [/math] tal que, fijado un intervalo de tiempo, la solución evaluada en dicho valor tome valores menores que cierto nivel de significación ([math] 10^{-\alpha}[/math]).
Primero observamos que
siendo [math]\Lambda [/math] el máximo valor que alcanza [math]\Phi(x,t) [/math] a partir de la cota que buscamos de [math]x[/math].
siendo [math] max(t) [/math] el mayor valor de [math] t[/math] en el intervalo considerado.
En un inicio, por la condición inicial del problema, tenemos una función meseta y podemos observar cómo, a media que el tiempo aumenta, el calor se va esparciendo por el intervalo de [math] x[/math] considerdo.
Para poder apreciar mejor este suceso, se presenta el siguiente vídeo:
Se puede apreciar cómo la solución se va aplanando a medida que aumenta el tiempo, parece tender a la solución estacionaria nula (como teóricamente debe ocurrir).
Nota: En el vídeo el intervalo de [math]x[/math] considerado es muy superior pues el tiempo máximo es [math]t=3[/math] en lugar de [math]t=0.1[/math].
9.4 Programa
clc
clear all
close all
%%%% SOLUCIÓN DE LA ECUACIÓN DEL CALOR MEDIANTE LA SOLUCIÓN FUNDAMENTAL CON CONDICIÓN UNICIAL U0 %%%
imp = 6; % Exponente del nivel de significación
T=[0.001,0.01,0.1]; % Tiempos concretos de visualización
a=-1; b=1; tf=max(T); % Intervalo donde u0 vale 1 y tf valor máximo del tiempo
div=10^-3; % Division de los intervalos
Y=a:div:b; % Vector de x donde u0 vale 1
c=-sqrt((log(2)-1/2*log(4*pi*tf)+imp*log(10))*4*tf); % Extremo izquierdo de la x para que en los extremos para todos los tiempos valga menos de 10^-imp
d=-c; % Extremo derecho
X=c:div:d; % Vector de x
Phi = @(x,t)(1./((4*pi*t).^(1/2)).*exp(-abs(x).^2./(4*t))); % Solución fundamental en 1 dimensión
U=zeros(length(X),length(T)); % Inicializamos los valores de u en X para cada tiempo
for j=1:length(X)
for i=1:length(T)
U(j,i)=trapz(Y,Phi(X(j)*ones(1,length(Y))-Y,T(i)*ones(1,length(Y)))); % Para cada tiempo evaluamos en t y X
end
end
% Dibujamos
for i=1:length(T)
subplot(1,3,i)
plot(X,U(:,i))
title("Solución de la ecuación del calor para t="+num2str(T(i))+' s')
ylim([0 1.2])
xlabel('longitud (m)')
ylabel('U(x,t) (ºC)')
%title('Solución de la ecuación del calor para t='+num2str(T(i)))
legend('U(x,t)')
end
clc
clear all
close all
%%%% VÍDEO DE LA SOLUCIÓN DE LA ECUACIÓN DEL CALOR MEDIANTE LA SOLUCIÓN FUNDAMENTAL CON CONDICIÓN UNICIAL U0 %%%
imp = 6; % Exponente del nivel de significación
divt=10^-2; % Division del intervalo de tiempo
t0=10^-2; tf=3; % Intervalo de definición del tiempo
T=t0:divt:tf; % Vector de tiempo
a=-1; b=1; % Intervalo donde u0 vale 1
div=10^-2; % Division del intervalo de x
Y=a:div:b; % Vector de x donde u0 vale 1
c=-sqrt(4*tf*(log(2)+ imp*log(10)-log(4*tf*pi)/2)); % Extremo izquierdo de la x para que en los extremos para todos los tiempos valga menos de 10^-imp
d=-c; % Extremo derecho
X=c:div:d; % Vector de x
Phi = @(x,t)(1./((4*pi*t).^(1/2)).*exp(-abs(x).^2./(4*t))); % Solución fundamental en 1 dimensión
U=zeros(length(X),length(T)); % Inicializamos los valores de u en X para cada tiempo
for j=1:length(X)
for i=1:length(T)
U(j,i)=trapz(Y,Phi(X(j)*ones(1,length(Y))-Y,T(i)*ones(1,length(Y)))); % Para cada tiempo evaluamos en t y X
end
end
pelicula=VideoWriter('Video 3.3'); % Creamos el video
pelicula.FrameRate=10;
open(pelicula);
for j=1:length(T)
figura = figure(1);
plot(X,U(:,j),"r","LineWidth",1.5)
ylim([0,1])
hold off
xlabel("longitud (m)")
xlim([c, d])
ylabel("temperatura (ºC)")
title("Solución de la Ecuación del Calor con condición inicial u_0")
subtitle("t = "+num2str(T(j))+" s")
%legend("solución de la ecuación del calor","solución estacionaria")
imagen=getframe(figura);
writeVideo(pelicula,imagen);
end
close(pelicula);
9.5 Solución fundamental de la ecuación del calor en dimensión 2
En último lugar, vamos a estudiar la solución fundamental de la ecuación del calor en dimensión [math]2[/math]. En este caso, tomamos el intervalo espacial [math](x_1,x_2)=[-1,1]^2[/math] y [math]D=1[/math]. La expresión de la solución fundamental es la siguiente:
A continuación se representa para diferentes instantes de tiempo ([math] t=0.001,~t=0.01[/math] y [math] t=0.1[/math]):
Nótese que los valores en el eje [math]z[/math] son distintos en cada una de las gráficas. Lo más destacable de esta representación es que, cuanto menor sea el valor del tiempo considerado, mayor es el valor máximo que alcanza la solución fundamental y el calor se va concentrando en un entorno del [math] (0,0)[/math]. Esto se debe a que, en el instante [math] t=0[/math], se tiene la Delta de Dirac.
Por último, presentamos el siguiente vídeo para observar cómo varía la solución fundamental respecto del tiempo para más valores:
Como podemos observar, a medida que aumenta el tiempo, el calor se va dispersando por un área mayor y el máximo valor de la solución fundamental disminuye.
9.6 Programa
clc
clear all
close all
%%%% SOLUCIÓN FUNDAMENTAL DIM 2 %%%
T=[0.001,0.01,0.1]; % Evaluaremos en tiempos concretos
a=-1; b=1; % Intervalo de definición de la x ([a,b]^2)
div=10^-3; % División del vector en x
X=a:div:b; % Vector de X
[X1,X2]=meshgrid(X,X); % Creación de la malla en x
norma2=@(x1,x2)sqrt(x1.^2+x2.^2); %Función norma
Phi2 = @(x1,x2,t)(1./((4*pi*t)).*exp(-norma2(x1,x2).^2./(4*t))); % Solución fundamental en 2 dim
PHI2=zeros(length(X),length(X),length(T)); % Inicializamos la solución en la malla
for i=1:length(T)
PHI2(:,:,i) = Phi2(X1,X2,T(i)); % Evaluamos la solución en la malla para cada tiempo
end
% Dibujamos
for i=1:length(T)
subplot(1,3,i)
surf(X1,X2,PHI2(:,:,i))
shading interp
colormap jet
xlabel('longitud (m)')
ylabel('tiempo (s)')
zlabel('T(x,t) (ºC)')
legend('\Phi_1(x,t)')
title("Solución fundamental en dim=2 para t="+num2str(T(i))+' s')
end
clc
clear all
close all
divt=2*10^-3; % División del intervalo de tiempo
t0=0.02; tf=0.3; % Intervalo de definición del tiempo
T=t0:divt:tf; % Vector de tiempo
a=-1; b=1; % Intervalo de definición de x1 y x2
div=10^-2; % División del intervalo de x
X=a:div:b; % Vector de X
[X1,X2]=meshgrid(X,X); % Creación de la malla
norma2=@(x1,x2)sqrt(x1.^2+x2.^2); % Definición de la norma
Phi2 = @(x1,x2,t)(1./((4*pi*t)).*exp(-norma2(x1,x2).^2./(4*t))); % Solución fundamental de dimensión 2
PHI2=zeros(length(X),length(X),length(T)); % Inicializamos la solución fundamental los valores de la malla
for i=1:length(T)
PHI2(:,:,i) = Phi2(X1,X2,T(i)); % Evaluamos la solución en la malla para cada tiempo
end
pelicula=VideoWriter('Video 3.4'); % Creamos el vídeo
pelicula.FrameRate=10;
open(pelicula);
for j=1:length(T)
figura = figure(1);
surf(X1,X2,PHI2(:,:,j))
shading interp
colormap jet
zlim([0,4])
hold off
zlabel("\Phi(x,t)")
title("Solución Fundamental en dimensión 2")
subtitle("t = "+num2str(T(j))+ " s")
imagen=getframe(figura);
writeVideo(pelicula,imagen);
end
close(pelicula);
