Revisión del 00:58 17 mar 2025

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:

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:

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

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

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:

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 pues el término [math](4\pi t)^{n/2}[/math], que muestra como en [math]u(x ,t)[/math] con [math]x \in \mathb{R}^n[/math] este efecto se vuelve más pronunciado en dimensiones superiores, pues el denominador crece más rápido. Con esto nos surge la siguiente pregunta, cómo podemos comparar el tiempo de dispersión en 1D, 2D y \(n\)D. Sabemos que la varianza de la solución fundamental en una dimensión es \(\sigma^2 = 2t\) con lo que para el análisis cualitativo de la dispersión se tiene que \(\sigma^2 \sim t\).

En 2D, la dispersión es más rápida porque la normalización de la solución implica un decaimiento más rápido \( u \sim t^{-1}\) frente a \(u \sim t^{-1/2}\). En \( n \) dimensiones, el calor decae como \(u(x,t) \sim t^{-n/2}\) lo que reafirma que el calor se disipa más rápido a medida que aumenta la dimensión del espacio.

Si queremos comparar el tiempo característico de disipación \( t_d \) en distintas dimensiones, fácilmente llegamos a la relación \( t_d(n) \sim t_d(1) n/2 \) que nos dice que el tiempo de disipación en dimensión \( n \) es inversamente proporcional a \( n/2 \).

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