Ecuación del calor (Grupo MAMBD))

Trabajo realizado por estudiantes
Título	Ecuación del calor. Grupo MAMBD
Asignatura	EDP
Curso	2024-25
Autores	Matilde Rubio Arranz, Antonio Lozano Fernández, Marcos Gil García, Bárbara Jiménez Pérez y Daniel Marcos Viña
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 Introducción

La ecuación del calor es un modelo matemático que describe la difusión térmica en un medio. En una dimensión, se expresa como:

[math]u_t-u_{xx}=0.[/math]

donde [math]u(x,t)[/math] representa la temperatura en función del tiempo y la posición.

Resolveremos un problema específico utilizando la solución fundamental:

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

permitiéndonos analizar la propagación del calor en una región no acotada e interpretar los resultados.

2 Mensaje secreto en la ecuación del calor

Imaginamos que escribimos un mensaje secreto en una barra de metal infinita usando calor (interpretando teóricamente dicha barra como una recta, ya que estamos en dimensión 1). Definimos la temperatura en el instante inicial

[math]u_0(x)=\left\{\begin{array}{ll} 5, & \text{$x$ está en una de las letras}\\ 0, & \text{en otro caso} \end{array}\right. ,[/math]

es decir, el mensaje inicial está formado por letras dibujadas en el eje [math]x[/math] con calor. Con el tiempo, el calor empieza a difundirse según la ecuación [math]u_t-u_{xx}=0[/math]. Su solución viene dada por la convolución

[math]u(x,t)=\int_{-\infty}^{\infty} G(x-y,t)u_0(y) dy=\frac{1}{4\sqrt{\pi t}}\int_{-\infty}^{\infty} u_0(y)e^{-\frac{(x-y)^2}{4t}}.[/math]

Dibujando la solución para diferentes tiempos, obtenemos lo siguiente:

Vídeo de la difusión del mensaje en la barra de metal con [math] t \in [0,1.60] [/math], dimensión 1 (1D).

Conforme avanza el tiempo, la difusión del calor hace que los bordes de las letras pierdan definición. Al inicio, las regiones calientes se distinguen claramente, pero poco a poco el calor se extiende a las áreas frías, suavizando los contornos. A medida que esta difusión continúa, las letras se vuelven borrosas hasta que el contraste desaparece por completo y el mensaje se vuelve ilegible, ya que la temperatura se ha distribuido de manera uniforme en la barra.

Difusión del mensaje en la barra de metal

import numpy as np
import matplotlib.pyplot as plt

# Parámetros
L = 10.0  # Longitud de la barra
T = 0.1   # Tiempo final
nx = 100  # Número de puntos en el espacio
nt = 1000 # Número de puntos en el tiempo
alpha = 1 # Coeficiente de difusión térmica

dx = L / (nx - 1)
dt = T / (nt - 1)
r = alpha * dt / dx**2

# Inicialización
x = np.linspace(0, L, nx)
u = np.zeros(nx)

# Condición inicial
def initial_condition(x):
    # Definir las regiones de las letras
    if 2 < x < 3 or 4 < x < 5:
        return 5
    else:
        return 0

u = np.array([initial_condition(xi) for xi in x])

# Evolución temporal
for t in range(nt):
    un = u.copy()
    for i in range(1, nx-1):
        u[i] = un[i] + r * (un[i+1] - 2*un[i] + un[i-1])

    # Visualización en tiempos específicos
    if t in [int(0.001/dt), int(0.01/dt), int(0.1/dt)]:
        plt.plot(x, u, label=f't={t*dt:.3f}')

plt.xlabel('Posición x')
plt.ylabel('Temperatura u(x,t)')
plt.title('Difusión térmica del mensaje')
plt.legend()
plt.show()

3 Regularidad y dimensiones

Vamos a explorar que ocurre al pasar a dos dimensiones (2D), analizando como la difusión del calor se dispersa más rápido en 2D debido a la propagación en múltiples direcciones.

En 1D, la ecuación del calor viene dada por [math]u_t-u_{xx}=0[/math]. Sabemos que dicha ecuación tiene un efecto regularizador, es decir, aunque la condición inicial tenga discontinuidades, la solución se vuelve suave para [math]t\gt0[/math]. Por otro lado, la difusión en 1D ocurre a lo largo de una sola dirección (el eje [math]x[/math]) y la solución se suaviza rápidamente con el tiempo.

Veamos que ocurre en 2D. Ahora, nuestra ecuación es [math]u_t-u_{xx}-u_{yy}=0[/math]. La temperatura inicial se define ahora como:

[math]u_0(x)=\left\{\begin{array}{ll} 5, & \text{si $(x,y)$ está dentro de una de las letras}\\ 0, & \text{en otro caso} \end{array}\right. ,[/math]

Continuemos con la generalización a [math]n[/math] dimensiones tal que la ecuación del calor vendrá dada por: [math]u_t - \Delta u = 0[/math]

Cuya solución fundamental (análoga a la hallada para una dimensión) se distribuye como una gaussiana en varias variables ([math]n[/math]): [math]u(x,t) = \frac{1}{(4\pi t)^{n/2}} e^{-\frac{|x|^2}{4 t}}[/math]

Aquí, el efecto regularizador del que hablábamos sigue presente. En cambio, en dimensiones superiores la difusión ocurre en múltiples direcciones, lo que supone que a mayor dimensión, más rápido se dispersa el calor, ya que hay mas direcciones hacia las que puede propagarse.

Vídeo de la difusión del mensaje en la barra de metal con [math] t \in [0,10] [/math], dimensión 2 (2D).

% Parámetros del espacio y tiempo
Lx = 10; Ly = 10; dx = 0.1; dy = 0.1;
x = -Lx:dx:Lx;
y = -Ly:dy:Ly;
[X, Y] = meshgrid(x, y);
t_values = logspace(-3, 1, 50); % Valores de tiempo en escala logarítmica

% Condición inicial: Definir "letras" con calor
U0 = zeros(size(X));
U0((X > -3) & (X < -2) & (Y > -1) & (Y < 1)) = 5; % Letra 1
U0((X > -1) & (X < 1) & (Y > -2) & (Y < 2)) = 5;  % Letra 2
U0((X > 2) & (X < 3) & (Y > -1) & (Y < 1)) = 5;  % Letra 3

% Configurar video
video = VideoWriter('difusion_termica_2D.avi');
video.FrameRate = 10;
open(video);

% Configurar la figura
figure;
cmap = hot;
clim = [0, 5]; % Fijar la escala de colores entre 0 y 5

for i = 1:length(t_values)
    t = t_values(i);
    sigma = sqrt(2 * t) / dx;
    U_t = gaussian_filter(U0, sigma);
    
    imagesc(x, y, U_t);
    axis([-Lx Lx -Ly Ly]);
    colorbar;
    colormap(cmap);
    caxis(clim); % Mantener la escala de colores fija
    title(sprintf('Difusión térmica en t = %.4f', t));
    xlabel('x');
    ylabel('y');
    drawnow;
    
    % Capturar frame para el video
    frame = getframe(gcf);
    writeVideo(video, frame);
end

close(video);

function G = gaussian_filter(U, sigma)
    tam = ceil(6 * sigma); % Tamaño del kernel basado en sigma
    [X, Y] = meshgrid(-tam:tam, -tam:tam);
    
    % Definir el filtro gaussiano
    G = exp(-(X.^2 + Y.^2) / (2 * sigma^2));
    G = G / sum(G(:)); % Normalizar
    
    % Aplicar la convolución con el filtro gaussiano
    G = conv2(U, G, 'same');
end