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

Revisión del 16:49 7 mar 2024

1 Introducción

El objetivo de este trabajo es el estudio de la ecuación del calor, comenzando con un análisis en un dominio acotado. Se investigará el efecto de la variación de los distintos parámetros involucrados en esta ecuación. Se planteará un problema inicial que será resuelto, interpretado y modificado para continuar con el estudio.

Además, se abordará la solución de la ecuación del calor en diferentes dimensiones, condiciones frontera y condiciones iniciales, viendo de esta forma como varía en cada caso y que características comparte.

2 Ecuación del Calor en dimensión [math]n=1[/math]

2.1 Definición

Se considera una varilla metálica que ocupa el intervalo [math] I [/math] y que se encuentra aislada por su superficie lateral, de manera que la conducción de calor sólo se produce en la dirección longitudinal. Inicialmente se planteará el problema de manera que la temperatura inicial de la varilla es [math]u1(x)[/math]. Además, como condiciones frontera se considera en el extremo izquierdo una temperatura [math]u2(t)[/math] y en el extremo derecho [math]u3(t)[/math]. Dado esto, a partir del principio de conservación de la energía^[1], y la Ley de Fourier^[2], se plantea el problema de EDP que modeliza el comportamiento de la temperatura.

Donde la función [math]u(x,t)[/math] es la ecuación EDP que determina el calor dependiendo del espacio y del tiempo; [math]k[/math] es la conductividad térmica y [math]c[/math] el calor específico.

Este problema se ha planteado con condiciones frontera de Dirichlet^[3]. Sin embargo, en estos problemas se puede trabajar con condiciones frontera de Neumann^[4] y Robin^[5], tal y como se verá más adelante.

2.2 Ejemplo

A continuación, se estudiará el comportamiento de este problema a partir de un caso particular.

Se considerará el intervalo [math][0,1] [/math], temperatura inicial de la varilla de 0º, y en los extremos, izquierdo y derecho, 0º y 1º respectivamente.

Además, se considerará que tanto la conductividad térmica, [math]k[/math], como el calor específico [math]c[/math] son 1. Entonces se plantea el siguiente problema a modelizar.

Una vez definido el problema, vamos a resolverlo para así observar su comportamiento.

2.2.1 Homogeneización y solución estacionaria

Para resolver esta EDP, en primer lugar se deben tener las condiciones frontera homogeneizadas, de manera que, [math] u(0,t)= 0, u(1,t)= 0 [/math]. Por tanto, primero se debe estudiar la solución estacionaria del problema. Para ello, se supondá que [math]{t \to \infty} [/math], en dicho caso, [math] u_{t} \approx 0[/math] y [math] u(x,t) \approx v(x)[/math].

Reescribiendo el problema,

Se obtiene entonces una EDO homogénea de orden 2, cuya solución es la siguiente.

A continuación, se muestra el código utilizado para representar la gráfica de la solución estacionaria, y la respectiva gráfica.

Representación de la solución estacionaria

%Definimos la solución estacionaria.
v=@(x) x;
%definimos el intervalo donde vamos a representar la solución. 
x=linspace(0,1,1000);

%Comando para representar la solución en R^2.
hold on 
plot(x, v(x), 'b-')
xlabel('x');
ylabel('v(x)');
title('Solución extacionaria.');
grid on; 
hold off

2.2.2 Separación de Variables

Posterioremente se utilizará el método de separación de variables^[6], a partir del cual, se obtiene la solución de a siguiente EDP, donde se ha realizado el cambio [math] w(x,t)=v(x) - u(x,t)[/math] para la homogenización de las condiciones frontera.

Una vez aplicado el método, se han obtenido los coeficientes de Fourier correspondientes, tal y como se explicó en el anterior trabajo ^[7] utilizando la extensión de funciones impares, y finalmente la solución de [math] w(x,t)[/math] será la siguiente.

Un vez deshecho el cambio de variable para obtener la solución del problema original, donde [math] u(x,t)=v(x) - w(x,t)[/math], se llega a la expresión final del problema original.

2.2.3 Representación y análisis

Una vez obtenida la solución, se verá la propagación del calor en función del tiempo y el espacio. Para ello, se empleará el siguiente programa que representará la solución obtenida tomando los 10 primeros términos de la serie.

Representación de la solución u(x,t)

% Definimos los términos de la serie
n=10;

%inicializamos la función en 0. 
f=@(x,t) 0;

%Tomamos un bucle para poder obtener la solución w con el nº de términos
%seleccionado.
for k=1:n
    f = @(x,t) exp(-k^2*pi^2*t).*sin(pi * k * x).*(2/(k*pi)).*((-1)^(k+1)) + f(x,t);
end

% Queremos representar u entonces definimos la solución deshaciendo el cambio
% w(x,t)=v(x)-u(x,t),donde x(x)=x.
g = @(x,t) x-f(x,t)

% Dibujar la solución en el intervalode [0,1].
[x,t] = meshgrid(0:.001:1);

%utilizamos el comando mesh para representar la solución en R^3.
mesh(x,t,g(x,t)) 
xlabel('x');
ylabel('t');
zlabel('u(x,t)');
title('Solución para k = 10');

Se puede observar como la representación de la solución final es muy similar a la que sería la representación en [math]R^{3} [/math] de la solución estacionaria, en el intervalo[math][0.4,1] [/math]. Fuera de este, la solución no se asemeja a la solución estacionaria. Debido a esto, debemos señalar que la solución obtenida del problema modificado, la [math]w(x,t) [/math], será similar al plano generado en [math]x=0[/math], fuera del entorno de [math]t=0 [/math], ya que como se ha visto, en el intervalo [math][0,0.4] [/math] la gráfica de [math]u(x,t) [/math] no se asemeja a la solución estacionaria.

Representación de la comparativa entre u(x,t) y v(x)

En este gráfica se puede observar lo anteriormente mencionado, donde la solución estacionaria y la original en gran parte del intervalo son prácticamente idénticas, formando prácticamente el mismo plano, tal y como se observa en la foto de la derecha.

2.2.4 Flujo

El flujo de calor en los extremos de una ecuación de calor describe cómo la energía térmica se transfiere entre dos puntos.

Para estudiar el flujo en este caso, hay que calcular [math] u_{x}(x,t)[/math] y estudiar su comportamiento en los extremos, es decir en [math] u_{x}(0,t)[/math] y [math] u_{x}(1,t)[/math]. A partir de este estudio, se puede determinar si el flujo es entrante o saliente en cada uno de los extremos. Esto se debe a que el flujo es [math]q(x_{i},t)= -ku_{x}(x_{i},t)[/math], donde [math]x_{i} \in [0,1][/math], así que el flujo, de izquierda a derecha, tiene signo contrario a la derivada.

Primero, definiremos el flujo de calor en los extremos.

Para la interpretación, se ha implementado un código que nos representa dichas funciones en función del tiempo.

Representación del flujo en función del tiempo.

% Definimos los términos de la serie
n=10;

%Inicializamos la función en 0. 
dif_f=@(x,t) 0;

%Tomamos un bucle para poder obtener w_{x} con el nº de términos
%seleccionado.
for k=1:n
    dif_f = @(x,t) 2*exp(-k^2*pi^2*t).*cos(pi * k * x).*((-1)^(k+1)) + dif_f(x,t);
end

%Expresión de la derivada de u(x,t).
g=@(x,t)(1-dif_f(x,t));

%Evaluamos la derivada en x=0 y x=1 
l1=@(t) g(0,t);
h1=@(t) g(1,t);

%calculamos el flujo en los dos puntos tomando k=1.
l=@(t) -l1(t);
h=@(t) -h1(t);

%definimos el intervalo donde vamos a representar la solución. 
t=linspace(0,1,1000);

%Comando para representar la solución en R^2.
hold on 
plot(t, h(t), 'r--')
plot(t, l(t), 'b--')
line([0 1], [0.1 -0.1], 'Color', 'black') %Delimitamos en el eje x positivo. 
xlabel('t');
ylabel('Flujo');
legend('Flujo en x=0', 'Flujo en x=1')
title('Representación del flujo en los 2 extremos');
grid on; 
hold off

En ambas gráficas se observa un flujo negativo constante, indicando que el calor se desplaza del extremo derecho al izquierdo de la barra.

En el extremo izquierdo, el flujo inicial es nulo debido a la constancia de la temperatura en su entorno, pero a medida que avanza el tiempo, el flujo aumenta hasta estabilizarse en [math] q(0,t) = -1 [/math] . Este valor es igual al flujo en el otro extremo de la barra, sin embargo, su comportamiento en función del tiempo es contrario, aumenta hasta llegar a [math] q(1,t) = -1 [/math] .

Esta igualdad de flujos en ambos extremos implica que no hay una diferencia neta de temperatura entre ellos, evitando así la presencia de fuentes o sumideros de calor en la barra. Además, la presencia de un flujo no nulo en los extremos confirma que la solución estacionaria implica una variación de temperatura a lo largo de la barra.

Por tanto, en el extremo izquierdo entra calor, sin embargo, en el extremo derecho, el calor saldrá.

2.2.5 Coeficiente de conductividad térmica

Una vez completado el ejemplo, se procede a realizar la primera modificación, en este caso será en el coeficiente de conductividad térmica para ver como afecta este a la solución del problema.

En este caso se establece [math] k = \frac{1}{2}[/math]; por tanto el nuevo problema a resolver será el siguiente.

Tras un proceso de adimensionalización, repitiendo el proceso explicado anteriormente se llega a que las solución estacionaria de este "nuevo" problema es igual que en el anterior caso, es decir,

Por otro lado, la solución final será de la siguiente forma.

Nuevamente se procede a la representación de las soluciones para así facilitar su estudio. Para ello se ha implementado el siguiente código.

Representación del flujo en función del tiempo.

% Definimos los términos de la serie
n=10;

%inicializamos la función en 0. 
f=@(x,t) 0;

%Tomamos un bucle para poder obetener la solución w con el nº de términos
%seleccionado.
for k=1:n
    f = @(x,t) exp(-k^2*pi^2*(t/2)).*sin(pi * k * x).*(2/(k*pi)).*((-1)^(k+1)) + f(x,t);
end
% Queremos representar u entonces definimos la solución deshaciendo el cambio
% w(x,t)=v(x)-u(x,t), donde v(x)=x.
g = @(x,t) x-f(x,t)

% Dibujar la solución en el intervalode tiempo t∈[0,1].
[x,t] = meshgrid(0:.001:1);

%utilizamos el comando mesh para representar la solución en R^3.
mesh(x,t,g(x,t)) 
xlabel('x');
ylabel('t');
zlabel('u(x,t)');

Aparentemente esta solución es muy similar a la obtenida cuando [math] k = 1[/math], sin embargo, en este caso podemos observar una mayor curvatura en el mallado generado por la solución de lo que se observaba en el anterior caso.

Para ver esto visualmente se va a representar nuevamente la solución estacionaria con la solución del problema.

Representación del flujo en función del tiempo.

Es muy obvio que en este caso la solución de la EDP no se aproxima de la misma manera que antes ya tiene forma aplanada únicamente en el intervalo [math] [0.75,1][/math] aproximadamente.

Por lo tanto, se observa un cambio respecto a la solución con [math] k = 1[/math].

A continuación, se va a representar la diferencia entre la solución y el estado estacionario en el punto medio de nuestro intervalo, [math] x = \frac{1}{2}[/math], para ver como se disipan en función del tiempo. Para esta representación se toma el intervalo [math] t \in [0,1] [/math] para ambos coeficientes de conductividad.

Representación de [math] u(1/2,t)- u_{est}(1/2)[/math] para [math] k = 1[/math] y [math] k = \frac{1}{2}[/math].

% Definimos los términos de la serie
n=10;

%Inicializamos ambas funciones en 0. 
f1=@(x,t) 0;
g1=@(x,t) 0;

%Solución de la EDP homogeneizada con k=1/2
for k=1:n
    f1 = @(x,t) exp(-k^2*pi^2*(t/2)).*sin(pi * k * x).*(2/(k*pi)).*((-1)^(k+1)) + f1(x,t);
end

%Solución de la EDP homogeneizada con k=1
for k=1:n
    g1 = @(x,t) exp(-k^2*pi^2*t).*sin(pi * k * x).*(2/(k*pi)).*((-1)^(k+1)) + g1(x,t);
end

%Solución estacionaria.
v=@(x) x;

%Solución de la EDP original con k=1/2 y k=1 resp.
f = @(x,t) v(x)- f1(x,t);
g = @(x,t) v(x)- g1(x,t);

%definimos el intervalo donde vamos a representar la solución.
t=linspace(0,1,1000);

%Comando para representar la solución en R^2.
hold on 
plot(t, g(1/2,t)-v(1/2), 'b-')
plot(t, f(1/2,t)-v(1/2), 'r-')
xlabel('x');
ylabel('y');
legend('k=1', 'k=1/2')
grid on; 
hold off

En esta última gráfica es fácil de observar lo mencionado anteriormente donde se razonaba la principal diferencia entre ambas soluciones en función de la [math] k [/math]. Si [math] k = \frac{1}{2}[/math] aparentemente la diferencia entre la solución y la solución estacionaria no se aproxima a 0 hasta muy próximo al tiempo final. Sin embargo, en [math] k = 1 [/math] se comprueba que la diferencia es 0, es decir, la solución estacionaria y la del problema original son iguales.

Concluimos que el coeficiente de conductividad térmica afecta a la solución final, cuanto más pequeño sea este, menor será la aproximación de la solución de la EDP a la del problema estacionario.

2.2.6 Condición inicial

Ahora se repetirá el ejemplo que se acaba de ver pero cambiando ciertas condiciones, en primer lugar se cambiará la condición frontera del lado derecho de tal forma que su valor será ahora 0, y en segundo lugar se cambiará también la condición frontera, de tal forma que ahora será la función:

De tal forma que se obtiene una EDP de la forma:

Lo primero que se observa, es como en este caso las condiciones frontera son 0, por lo que no hay que homogeneizar el problema, de tal forma que se pasa directamente al método de separación de variables. Resolviéndolo, la solución que se obtiene es:

Donde [math]C_{k}[/math] son los coeficientes de la serie de Fourier. En este caso, en la condición inicial se tiene una función cuyos coeficientes son difíciles de obtener, de tal forma que se obtendrán de forma numérica con la fórmula del trapecio. Hay que decir que en la base trigonométrica del problema se ha obtenido un seno, es por eso que para obtener los coeficientes se realizará una extensión impar.

Vamos a ver ahora cómo quedaría la solución representada, y cómo es el flujo en los extremos de la barra metálica:

Representación de la solución u(x,t)

Representación de la solución u(x,t)

clc
clear all
format long
close all

% Definir la función a representar tomando los 100 primeros términos de la
% serie
% Definir la función a integrar
f = @(x, n) sin(pi * n * x) .* (1-4*abs(x-1/2));
% Definir los límites de integración
x1=linspace(1/4, 3/4, 100);

%numero de terminos de la serie de Fourier
i=100;

% Inicializar un vector para almacenar los resultados
resultados_trapecio_adaptativoa = zeros(i, 1);

% Calcular las integrales para n desde 1 hasta 100 usando el método del trapecio
for n = 1:i
    resultadoa = 2*trapz(x1,f(x1,n));
    % Almacenar el resultado
    resultados_trapecio_adaptativoa(n) = resultadoa;    
end
%Serie de fourier para n=100, definimos las funciones del calor, y su
%diferencial, y las metemos en un bucle para obtenerlas
ecalor = @(x,t) 0;
diffcalor=@(x,t) 0;
for n=1:i
    ecalor = @(x,t) resultados_trapecio_adaptativoa(n).* sin(n*pi*x).*exp(-n^2*pi^2*t) + ecalor(x,t);
    diffcalor=@(x,t) resultados_trapecio_adaptativoa(n).* cos(n*pi*x).*exp(-n^2*pi^2*t)*n*pi + diffcalor(x,t);

end

% Dibujar la solución en el intervalode tiempo t∈[0,1].
[x,t] = meshgrid(0:.001:1);
%Dibujamos la ecuación del calor en el tiempo y en x, obteniendo una
%gráfica de dimensión 3

figure 
mesh(x,t,ecalor(x,t)) 
xlabel('x');
ylabel('t');
zlabel('u(x,t)');
title('Solución para k = 10');
disp(resultados_trapecio_adaptativoa)
hold on
%Dibujamos ahora la gráfica del flujo de entrada y salida del calor, para
%ello definimos las funciones h y l, las cualesindican el flujo del calor
%en ambos extremos
figure
h=@(t) diffcalor(1,t);
l=@(t) diffcalor(0,t);
t=linspace(0,4,1000);
hold on 
plot(t, -h(t), 'r--')
plot(t, -l(t), 'b--')
xlabel('t');
ylabel('Flujo');
legend('Flujo en x=0', 'Flujo en x=1')
title('Representación del flujo en los 2 extremos');
grid on; % Activar la cuadrícula
hold off

El Método del Trapecio^[8] es una función de Matlab que te aproxima integrales de manera numérica.

En la gráfica se puede apreciar cómo en el instante [math] t=0[/math], entre [math] x=\left[\frac{1}{4},\frac{3}{4}\right][/math] la temperatura sube hasta llegar a un pico y luego decrece, tal y como marca la condición frontera, pero ha medida que avanza el tiempo, la temperatura se queda estable en 0, alcanzando un estado estacionario.

En cuanto al flujo, se aprecia como se va escapando por ambos extremos. También se ve cómo a medida que pasa el tiempo el flujo llega hasta cero en ambos casos, por lo que no hay transferencia de calor en ningún punto.

2.2.7 Condiciones Frontera

Ahora se cambiará la condición frontera del lado derecho, de tal forma que se considerará que está aislado térmicamente. Para representar esta condición se hará considerando:

De tal forma que se obtiene una EDP de la forma:

En este caso tendremos una condición frontera de Dirichlet, sin embargo, la otra condición frontera será de Neumann (la del extremos derecho de la varilla).

Nuevamente las condiciones frontera son 0, por lo que no es necesario homogeneizarlo. Se pasa directamente a la separación de variables, pero en este caso no se obtiene la base trigonométrica usual, en su lugar se obtiene:

Finalmente la solución de este problema será:

En este caso [math]C_{k}[/math] se calcula de la misma forma que en el apartado anterior, pero habrá que tener en cuenta la base trigonométrica en la que estamos. Ahora se mostrarán los gráficos de la solución obtenida, y el flujo en ambos extremos de la vara metálica.

Representación de la solución u(x,t)

Representación de la solución u(x,t)

clc
clear all
format long
close all

% Definir la función a representar tomando los 100 primeros términos de la
% serie
% Definir la función a integrar
f = @(x, n) sin((pi * n + pi/2 ) * x) .*max(0,(1-4*abs(x-1/2))) ;
% Definir los límites de integración
x1=linspace(0, 1, 100);
%numero de terminos de la serie de fourier
i=100;
% Inicializar un vector para almacenar los resultados
Coefs_Fourier = zeros(i, 1);
% Calcular los coeficientes de Fourier para n desde 1 hasta 100 usando el método del trapecio adaptativo
for n = 0:i
    resultadoa = 2*trapz(x1,f(x1,n));
    % Almacenar el resultado
    Coefs_Fourier(n+1) = resultadoa;    
end
%Serie de fourier para n=100, definimos las funciones del calor, y su
%diferencial, y las metemos en un bucle para obtenerlas
ecalor = @(x,t) 0;
diffcalor=@(x,t) 0;
for n=0:i
    ecalor = @(x,t) Coefs_Fourier(n+1).* sin((n*pi+pi/2)*x).*exp(-(n*pi+pi/2)^2*t) + ecalor(x,t);
    diffcalor=@(x,t) Coefs_Fourier(n+1).* cos((n*pi+pi/2)*x).*exp(-(n*pi+pi/2)*t)*(n*pi+pi/2) + diffcalor(x,t);

end

% Dibujar la solución en el intervalode tiempo t∈[0,1].
[x,t] = meshgrid(0:.001:1);

%Dibujamos la ecuación del calor en el tiempo y en x, obteniendo una
%gráfica de dimensión 3
figure 
mesh(x,t,ecalor(x,t)) 
xlabel('x');
ylabel('t');
zlabel('u(x,t)');
title('Solución para k = 10');
disp(Coefs_Fourier)
hold on
%Dibujamos ahora la gráfica del flujo de entrada y salida del calor, para
%ello definimos las funciones h y l, las cualesindican el flujo del calor
%en ambos extremos
figure
h=@(t) diffcalor(1,t);
l=@(t) diffcalor(0,t);
t=linspace(0,4,1000);
hold on 
plot(t, -h(t), 'r--')
plot(t, -l(t), 'b--')
xlabel('t');
ylabel('Flujo');
legend('Flujo en x=0', 'Flujo en x=1')
title('Representación del flujo en los 2 extremos');
grid on; % Activar la cuadrícula
hold off

Este caso es parecido al anterior, en t=0 alcanza un pico, pero posteriormente se estabiliza. Se ve como en en el extremo derecho de la barra la temperatura baja más drásticamente. Como en el caso anterior, el flujo entra por el extremo izquierdo y se acaba yendo a 0, y por el extremo derecho se mantiene en 0, tal como se ha establecido en la condición.

2.2.8 Principio del Máximo

El principio del máximo es un teorema que establece que, sea [math]U\in C^{2,1}(Q_{T}) \cap C(\overline{Q_t})[/math] que cumple

3 Ecuación del Calor en diferentes dimensiones

En este apartado, se verá que aspecto toma la función solución del calor en diferentes dimensiones y diferentes condiciones.

3.1 Solución fundamental de la ecuación del calor en dimensión 1

Considerando el problema siguiente,

Se obtiene la solución fundamental del calor en una dimensión:

Se representará para [math] x \in [-1,1], t\in [10^{-2},1] [/math]

% Definimos el recinto en el que vamos a representar la función
[X, T] = meshgrid(-1:0.05:1, 10^(-2):0.05:1);
%Definimos la solución fundamental del calor en una dimensión
u= 1./sqrt(4*pi*T).*exp(-X.^(2)./(4*T));
% Graficamos la función
figure;
surf(X, T, u);
title('Gráfico de la función');
xlabel('X');
ylabel('T');
zlabel('u(x,t)');

Como se puede observar se trata se una función normalizada, tal que a medida que se acerca al punto [math] (x,t)=(0,0) [/math], tiende a infinito.

Considerando las siguientes condiciones iniciales y frontera:

ecuaci´on del calor... en una dimensión en el semiespacio x > 0 con la condición inicial u(x, 0) = 0 y la condici´on frontera u(0, t) = 1. ...

3.2 Solución de la ecuación del calor en diferentes situaciones

Sea el siguiente problema:

Su solución es,

SOLUCION

En el siguiente caso, se considera el siguiente problema:

La solución, viene dada por la convolución respecto la solución fundamental del calor y la condición inicial, es decir,

Se debe notar que la función anterior aplicada a la condición inicial, en este caso resulta

Para obtener una mayor visión de dicha solución, se representará en diferentes instantes de tiempo, [math] t \in (0.001, 0.01, 0.1) [/math] y en el espacio [math] x \in [-1,1] [/math]

Solución en diferentes instantes de tiempo

clc
clear 
format long
close all

t_vector = [0.001, 0.01, 0.1]; % Parámetros de t
X= -1:10^(-3):1; %Intervalo para representar
Y= -1:10^(-3):1; %Intervalo para integración
[fila,columna]= size(X);

%Hacemos el bucle para los 3 valores de t
for j=1:3
    t= t_vector(j);
    F = ones(size(X));
    for i=1:columna
        %Definimos cada elemento dado por la solución
        F(i) = F(i)* 1/sqrt(4 * pi * t) * trapz(Y,exp(-(X(i)-Y).^2 ./ (4*t)));
    end
    %Graficamos la función
    subplot(1,3,j);
    plot(X,F);
    xlabel('x');
    ylabel('u(x,t)');
    title(['t = ',num2str(t)]);
end

Se puede ver en la gráfica que a menor [math]t [/math], la función se asemeja a la función [math] u(x)_{0} = 1_{[-1,1]} [/math]. Además a medida que [math]t [/math] crece, la función va tomando forma de parábola y decrece.

Por otro lado, la solución del calor en dimensión [math] 2[/math] es la siguiente:

De igual forma, se representará la solución para instantes de tiempo [math] t \in (0.001, 0.01, 0.1) [/math], y para [math] (x,y) \in [-1,1]\times[-1,1] [/math]

Como se aprecia en las gráficas, para valores de [math] t[/math] cercanos a cero, la función toma valores altos en el punto [math] (x,y)=(0,0)[/math]. A medida que [math] t[/math] va aumentando, la función va teniendo un aspecto más en forma de paraboloide. Esto se debe a que para [math] t=0[/math] la función no está definida.

clc
clear 
format long
close all

% Definimos el recinto en el que vamos a representar la función
[X1, X2] = meshgrid(-1:0.05:1, -1:0.05:1);
% Definimos el vector de tiempos
t_vector=[0.001,0.01,0.1];
%Preparamos la representación
tiledlayout(1,3);
%Hacemos un bucle para los tiempos
for i=1:3
t= t_vector(i);
%Definimos la solución fundamental del calor en dos dimensiones
u= 1./sqrt(4*pi*t).*exp((-X1.^(2)-X2.^(2))./(4*t));
%Graficamos la solución
nexttile
surf(X1, X2, u);
title(['t = ',num2str(t)]);
xlabel('X');
ylabel('X');
zlabel('u(x,t)');
end

4 Referencias

1. ↑ Conservación de la energía
2. ↑ Ley de Fourier
3. ↑ Condición de frontera de Dirichlet
4. ↑ Condición de frontera de Neumann
5. ↑ Condición de frontera de Robin
6. ↑ Método de separación de variables
7. ↑ Series de Fourier
8. ↑ Integración numérica trapezoidal