Revisión actual del 12:57 28 abr 2024

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
- 4.1 Programa
5 Interpretación del flujo en los extremos
- 5.1 Programa
6 Solución al considerar otro coeficiente de conductividad
- 6.1 Programa
7 Solución al considerar otra condición inicial
8 Cambio de condiciones frontera
- 8.1 Principio del máximo
- 8.2 Programa
9 Soluciones de la ecuación del calor usando la solución fundamental
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

[math] q=-k\cdot T_{x} [/math]

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:

[math] e= c \cdot T [/math]

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:

[math]\left \{ \begin{array}{ll} \frac{\partial T}{\partial t}-\frac{k}{c} \frac{\partial ^2 T}{\partial x^2}=0 & \quad 0 \lt x \lt 1, 0 \lt t \\ T(x,0)=0, & \quad 0 \lt x \lt 1, \\ T(0,t)=0, & \quad 0 \lt t , \\ T(1,t)=1 & \quad 0 \lt t . \end{array} \right. [/math]

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:

[math] \begin{array}{ll} T_{est}(x)=x, & \quad x \in [0,1] \end{array}. [/math]

En la siguiente gráfica se representa esta solución:

Representación de la solución de la ecuación del estado estacionario, [math] T_{est}(x)=x [/math], para [math] x \in [0,1] [/math]

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

[math]\left \{ \begin{array}{ll} \frac{\partial T_{hom}}{\partial t}-\frac{\partial ^2 T_{hom}}{\partial x^2}=0 & \quad 0 \lt x \lt 1, 0 \lt t , \\ T_{hom}(0,t)=0, & \quad 0 \lt t , \\ T_{hom}(1,t)=0 & \quad 0 \lt t , \\ T_{hom}(x,0)=x, & \quad 0 \lt x \lt 1. \end{array} \right. [/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:

[math] T_{hom}(x,t)=u(x)v(t) [/math]

Haciendo las cuentas pertinentes se tiene que

[math] T_{hom}(x,t)=\sum_{k=1}^{\infty} \frac{2(-1)^{k+1}}{k \pi}\sin{(k \pi x)} e^{-k^2\pi^2t} [/math].

Finalmente, la solución del problema original es [math] T(x,t)=x-\sum_{k=1}^{\infty} \frac{2(-1)^{k+1}}{k \pi}\sin{(k \pi x)} e^{-k^2\pi^2t} [/math].

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:

Representación de la solución estacionaria [math] T_{est} [/math] y de [math] T(x,t) [/math] para [math] t \in [0,1] [/math], tomando [math] 1000 [/math] términos de la serie de Fourier. Se ha considerdo una malla con una división de [math] 10 ^{-2}[/math] para facilitar la visualización de la pendiente de la gráfica, ya que con una división más fina es menos visible en Matlab.

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:

Representación de la solución [math] T(x,t) [/math] sobre la varilla al variar el tiempo entre [math] 0[/math] y [math] 1[/math] , así como la solución estacionaria [math] T_{est}=x [/math] .

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.

Representación de la solución [math] T(x,t) [/math] para [math] t \in [0,1] [/math] sobre la varilla de metal, en función del 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:

[math] T_{x}(x,t)= 1- \sum_{k=1}^{\infty} 2(-1)^{k+1}\cos{(k \pi x)} e^{-k^2\pi^2t} [/math].

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:

[math] T_{x}(0,t)=1- \sum_{k=1}^{\infty} 2(-1)^{k+1} e^{-k^2\pi^2t} [/math]. [math] T_{x}(1,t)=1- \sum_{k=1}^{\infty} -2 e^{-k^2\pi^2t} [/math].

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.

[math] T_{x}(0,0)=1- \sum_{k=1}^{\infty} 2(-1)^{k+1} [/math].

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.

[math] T_{x}(1,0)=1- \sum_{k=1}^{\infty} -2 \rightarrow -\infty [/math].

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:

Evolución del flujo con respecto al tiempo en los extremos izquierdo ([math]x=0[/math]) y derecho ([math]x=1[/math]). Considerando [math]t \in [10^{-3}, 2][/math] y [math]1000[/math] 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:

[math] \left \{ \begin{array}{ll} \frac{\partial T}{\partial t}-\frac{1}{2} \frac{\partial ^2 T}{\partial x^2}=0 & \quad 0 \lt x \lt 1, 0 \lt t , \\ T(x,0)=0, & \quad 0 \lt x \lt 1, \\ T(0,t)=0, & \quad 0 \lt t , \\ T(1,t)=1 & \quad 0 \lt t . \end{array} \right. [/math]

Siguiendo los mismos pasos que el problema inicial, obtenemos que la solución es

[math] T(x,t)=\sum_{k=1}^{\infty} 2(-1)^{k+1}\sin{(k \pi x)} e^\frac{-k^2\pi^2t}{2} [/math].

En las siguientes gráficas, se muestran las soluciones de los problemas considerando [math]k=1[/math] y [math]k=1/2[/math]:

Representación de las soluciones [math] T(x,t) [/math] para [math] t \in [0,1] [/math], tomando [math] 1000 [/math] términos de la serie, para [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].

Representación de la diferencia entre la solución estacionaria y la solución en [math] x=\frac{1}{2}[/math] en escala logarítmica, para [math] k=1[/math] y [math] k=1/2[/math]. Hemos trabajado con [math] 1000[/math] términos de la serie de Fourier.

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:

Representación del flujo en [math] x=0[/math] y [math] x=1[/math], considerando la solución con [math] 1000[/math] términos de la serie de Fourier, tiempo [math]t \in [10^-3, 2] [/math], con [math] k=1[/math] y [math] k=\frac{1}{2}[/math].

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:

Comparación de las soluciones para [math]k=1[/math] y [math]k=\frac{1}{2}[/math], considerando [math]t \in [10^{-3}, 2][/math] y [math]1000[/math] términos de la serie.

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=1/2 %%% 

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:

[math]\left \{ \begin{array}{ll} \frac{\partial T}{\partial t}- \frac{\partial ^2 T}{\partial x^2}=0 & \quad 0 \lt x \lt 1, 0 \lt t , \\ T(x,0)=max\{0,1-4|x-\frac{1}{2}|\}, & \quad 0 \lt x \lt 1, \\ T(0,t)=T(1,t)=0 & \quad 0 \lt t \end{array} \right. [/math]

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:

[math] u(x,t)=\sum_{k=1}^{\infty} c_{k}\sin{(k \pi x)} e^{-k^2 \pi^2t} [/math].

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:

Representación de la solución en [math] t \in [0,1][/math], [math] x \in [0,1][/math] y considerando [math] 1000[/math] términos de la serie de Fourier.

A continuación, mostramos un vídeo en el que representamos la solución en función del tiempo:

Representación de la solución en [math] t \in [0,1][/math], [math] x \in [0,1][/math] y considerando [math] 1000[/math] términos de la serie de Fourier.

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:

[math] T_{x}(x,t)=\sum_{k=1}^{\infty} c_{k} \cos{(k \pi x)} e^{-k^2\pi^2t} [/math].

La gráfica que muestra la evolución del flujo en los extremos de la barra con respecto al tiempo es:

Representación del flujo en los extremos izquierdo ( [math] x=0 [/math]) y derecho ([math] x=1[/math]), en [math] t \in [0,1][/math] tomando [math] 1000 [/math] términos de la serie de Fourier en la solución.

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

[math]\left \{ \begin{array}{ll} \frac{\partial T}{\partial t}- \frac{\partial ^2 T}{\partial x^2}=0 & \quad 0 \lt x \lt 1, 0 \lt t , \\ T(0,t)=0, & \quad 0 \lt t , \\ T_{x}(1,t)=0 & \quad 0 \lt t , \\ T(x,0)=max \{0,1-4|x-1/2|\}, & \quad 0 \lt x \lt 1. \end{array} \right. [/math]

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{\left(\left(\frac{\pi}{2}+k\pi\right)x\right)} e^{-\left(\frac{\pi}{2}+k\pi\right)^2t}[/math] para [math]k=0, ~..., ~\infty[/math]. Con esta solución, no podemos emplear la base trigonométrica usual [math]\left\{ \frac{1}{2}, sen(k \pi x), cos(k \pi x) \right\}_{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]\left\{sen\left(\left(\frac{\pi}{2}+k\pi\right)x\right)\right\}_{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.

[math] \langle sen\left(\left(\frac{\pi}{2}+n\pi\right)x\right), sen\left(\left(\frac{\pi}{2}+m\pi\right)x\right)\rangle _{L^2(0,1)}= \int_{0}^{1}sen\left(\left(\frac{\pi}{2}+n\pi\right)x\right) \cdot sen\left(\left(\frac{\pi}{2}+m\pi\right)x\right) dx= \left[ \frac{sen(((n-m) \pi )x)}{\pi (2n-2m)}- \frac{sen((\pi(n+m+1)x)}{\pi(2n+2m+2)} \right]_{0}^{1}=0 [/math].

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:

[math] u(x,t)=\sum_{k=0}^{\infty} c_{k}\sin{\left(\left(\frac{\pi}{2}+k\pi\right)x\right)} e^{-\left(\frac{\pi}{2}+k\pi\right)^2t} [/math].

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:

Representación de la solución tomando [math] 1000 [/math] términos de la serie de Fourier en la solución, [math] t \in [0,1] [/math] y [math] x \in [0,1] [/math].

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:

Vídeo de la solución tomando [math] 1000 [/math] términos de la serie de Fourier en la solución, [math] t \in [0,1] [/math] y [math] x \in [0,1] [/math].

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:

[math] \Phi_D (x,t) = \frac{1}{(4 \pi Dt)^{n/2}}e^{-\frac{|x|^2}{4Dt}}[/math]

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

[math]\left \{ \begin{array}{ll} \frac{\partial \Phi_D}{\partial t}- D\Delta \Phi_D=0 & \quad x \in \mathbb{R}^n, 0 \lt t , \\ \Phi_D(x,t)=0, & \quad |x| \rightarrow \infty, 0 \lt t , \\ \Phi_D(x,0)=\delta(x), & \quad x \in \mathbb{R}^n, \\ \int_{\mathbb{R}^n}\Phi_D(x,t)dx = 1 & \quad t \gt0. \end{array} \right. [/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:

[math] \Phi_1 (x,t) = \frac{1}{(4 \pi t)^{1/2}}e^{-\frac{|x|^2}{4t}}[/math].

A continuación podemos ver la representación gráfica de dicha función:

Representación de la solución fundamental para [math]x \in [-1,1][/math], [math]t \in [10^{-2},1][/math] ,[math]D=1[/math] y [math]n=1[/math].

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:

[math]\left \{ \begin{array}{ll} \frac{\partial T}{\partial t}- \frac{\partial ^2 T}{\partial x^2}=0 & \quad 0 \lt x , 0 \lt t \\ T(0,t)=1 & \quad 0 \lt t \\ T(x,t)=0 & \quad 0 \lt t , x \rightarrow \infty \\ T(x,0)=0 & \quad x\gt0 \end{array} \right. [/math]

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

[math] 0=\frac{\partial T}{\partial t}- \frac{\partial ^2 T}{\partial x^2}=u_t\left(\frac{x}{\sqrt{t}}\right)-u_{xx}\left(\frac{x}{\sqrt{t}}\right)= -\frac{x}{2t\sqrt{t}}u'\left(\frac{x}{\sqrt{t}}\right) -\frac{1}{t}u''\left(\frac{x}{\sqrt{t}}\right) \Longleftrightarrow 0= \frac{1}{2}\frac{x}{\sqrt{t}}u'\left(\frac{x}{\sqrt{t}}\right) + u''\left(\frac{x}{\sqrt{t}}\right) [/math]

Llamamos [math] \zeta = \frac{x}{\sqrt{t}} [/math] y obtenemos que [math]u[/math] cumple la siguiente ecuación diferencial:

[math] 0=\frac{1}{2}\zeta u'(\zeta) +u''(\zeta) [/math]

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:

[math]\left \{ \begin{array}{ll} \frac{\partial T_1}{\partial t}- \frac{\partial ^2 T_1}{\partial x^2}=0 & \quad 0 \lt x , 0 \lt t \\ T(0,t)=0 & \quad 0 \lt t \\ T(x,t)= \sqrt{\pi} & \quad 0 \lt t , x \rightarrow \infty \\ T(x,t)=\sqrt{\pi} & \quad x\gt0 , x \rightarrow 0 \end{array} \right. [/math]

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:

[math] T(x,t) = 1 - \frac{2}{\sqrt{\pi}}\int^{\frac{x}{2\sqrt{t}}}_0e^{-z^2}dz = 1 - erf\left(\frac{x}{2\sqrt{t}}\right) [/math]

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.

[math] 10^{-\alpha} \gt T(x,t) =1 - erf\left(\frac{x}{2\sqrt{t}}\right)=1 - erf(z)= 1 - \frac{2}{\sqrt{\pi}}\int^{z}_0e^{-s^2}ds =\frac{2}{\sqrt{\pi}}\int^{\infty}_0e^{-s^2}ds - \frac{2}{\sqrt{\pi}}\int^{z}_0e^{-s^2}ds = \frac{2}{\sqrt{\pi}}\int^{\infty}_ze^{-s^2}ds ~~~~~~ \left(z=\frac{x}{2\sqrt{t}}\right) [/math]

siendo [math]10^{-\alpha}[/math] el valor de significación que consideremos.

[math] \left(\int^{\infty}_ze^{-s^2}ds \right)^2 = \int_z^{\infty}\int_z^{\infty}e^{-(s_1^2+s_2^2)}ds_1ds_2 \leq \int_0^{\frac{\pi}{2}}\int_z^{\infty}re^{-(r^2)}drd\theta [/math]

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.

[math] \int_0^{\frac{\pi}{2}}\int_z^{\infty}re^{-(r^2)}drd\theta = \frac{\pi}{2}\left[\frac{e^{-r^2}}{2}\right]^{\infty}_z=\frac{\pi}{4}e^{-z^2} [/math] [math] 10^{-2\alpha}\gt \frac{4}{\pi} \frac{\pi}{4}e^{-z^2} \Leftrightarrow -2\alpha \log (10) \gt -z^2 \Leftrightarrow z \gt \sqrt{ 2\alpha\log (10)} [/math]

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á:

[math] x \gt 2\sqrt{max(t)2\alpha\log (10)} [/math]

siendo [math]max(t)[/math] el máximo valor de tiempo que queramos considerar.

Obtenemos la siguiente la gráfica:

Representación de la solución de la ecuación del calor en el semiespacio [math]x \geq 0[/math] , [math]t \in [0,1][/math] ,[math]D=1[/math], dimensión 1 con una significación [math] 10^{-6} [/math] ([math] \alpha =6 [/math]) .

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:

Vídeo de la solución de la ecuación del calor en el semiespacio [math]x \geq 0[/math], con [math] x \in [0,5] [/math] , [math]t \in [0,500][/math] ,[math]D=1[/math], dimensión 1.

Observamos cómo, efectivamente, la temperatura parece tender a [math]1[/math] cuando el tiempo aumenta.

9.2.1 Programa

clc 
clear all 
close all 

%%% SOLUCIÓN EN EL SEMIESPACIO POSITIVO %%%
  
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(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 

%%% VIDEO SOLUCIÓN EN EL SEMIESPACIO POSITIVO %%%

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:

[math] U(x,t)= \Phi(x,t) \ast_{x} u_0(x) = \int_{-\infty}^{\infty} \Phi(x-y,t) \cdot u_0(y) ~ dy [/math]

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

[math] \int_{-\infty}^{\infty} \Phi(x-y,t) \cdot u_0(y) ~ dy = \int_{x-1}^{x+1} \Phi(z,t) ~ dz \leq 2 \Lambda [/math]

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

[math] 10^{-\alpha}\gt 2 \Phi(x,t) = 2 \frac{1}{(4 \pi t)^{1/2}}e^{-\frac{x^2}{4t}} ~ \Leftrightarrow ~ - \alpha \log(10)\gt \log(2) - \frac{1}{2} \log(4 \pi t) - \frac{x^2}{4t} ~ \Leftrightarrow ~ x \gt \sqrt{4 max(t) ( \log(2) - \frac{1}{2} \log(4 \pi max(t)) + \alpha \log(10))} [/math]

siendo [math] max(t) [/math] el mayor valor de [math] t[/math] en el intervalo considerado.

Representación de la solución fundamental con [math]x \in \mathbb{R}[/math] con un nivel de significación de [math]10^{-6}[/math], [math]D=1[/math], [math]n=1[/math] y considerando los tiempos [math] t= 0.001[/math], [math] t= 0.01[/math] y [math] t= 0.1[/math].

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:

Representación de la solución fundamental con [math]x \in \mathbb{R}[/math] con un nivel de significación de [math]10^{-6}[/math], [math]D=1[/math], [math]n=1[/math] y considerando [math] t\in [10^{-2}, 3][/math].

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.3.1 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.4 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:

[math] \Phi_1 (\vec{x},t) = \frac{1}{4 \pi t}e^{-\frac{|\vec{x}|^2}{4t}}[/math]

A continuación se representa para diferentes instantes de tiempo ([math] t=0.001,~t=0.01[/math] y [math] t=0.1[/math]):

Representación de la solución fundamental con [math](x_1,x_2)=[-1,1]^2[/math] [math]D=1[/math], [math]n=2[/math] y considerando [math] t=0.001[/math], [math] 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:

Representación de la solución fundamental con [math](x_1,x_2)=[-1,1]^2[/math] [math]D=1[/math], [math]n=2[/math] y considerando [math] t \in [0.02,0.3][/math].

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

%%%% VIDEO DE LA SOLUCIÓN FUNDAMENTAL DIM 2 %%%

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);

10 Referencias

¿Qué es la Ley de conducción térmica de Fourier? Definición. Nick Connor

↑ Delta de Dirac. Fernando Andino, Marlon Recarte, Michael Spilsbury [1]

[1]

Diferencia entre revisiones de «Ecuación del calor (GRwM)»