Ecuación del calor RAJ

Trabajo realizado por estudiantes
Título	Ecuación del calor. Grupo LÁJ
Asignatura	EDP
Curso	2025-26
Autores	Luis García Suárez Álvaro Moreno Cisneros Juan Pérez Guerra
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

PDF del póster

Abajo se puede ver el código que se ha utilizado para conseguir las gráficas.

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from mpl_toolkits.mplot3d import Axes3D

#Visualización unidimensional

def nucleo_calor(x, t):
    return (1.0 / np.sqrt(4 * np.pi * t)) * np.exp(-(x**2) / (4 * t))

def condicion_inicial_escalon(x):
    return np.where((x >= 0) & (x <= 1), 1.0, 0.0)

def condicion_inicial_gaussiana(x):
    return np.exp(-x**2)

def resolver_calor_integral(x_eval, t, u0_func, y_min=-10, y_max=10, dy=0.002):
    y = np.arange(y_min, y_max, dy)
    u0_y = u0_func(y)

    u_vals = np.zeros_like(x_eval, dtype=float)

    for i, x in enumerate(x_eval):
        integrando = nucleo_calor(x - y, t) * u0_y
        u_vals[i] = np.trapezoid(integrando, y)

    return u_vals

x_dominio = np.linspace(-2, 3, 500)
tiempos = [0.001, 0.01, 0.1, 0.5, 1]


plt.figure(figsize=(10, 6))
plt.plot(x_dominio, condicion_inicial_escalon(x_dominio), 'k--', label='t = 0 (Dato Inicial)')

for t in tiempos:
    u_t = resolver_calor_integral(x_dominio, t, condicion_inicial_escalon)
    plt.plot(x_dominio, u_t, label=f't = {t}')

plt.title(r'Evolución de la Ecuación del Calor por Convolución en $\mathbb{R}$')
plt.xlabel('Posición x')
plt.ylabel('Temperatura u(x,t)')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
plt.close()

x_vals = np.linspace(-2, 3, 100)
t_vals = np.linspace(0.001, 2, 100)

X_mesh, T_mesh = np.meshgrid(x_vals, t_vals)

y_int = np.linspace(0, 1, 300)

diffs = X_mesh[:, :, np.newaxis] - y_int[np.newaxis, np.newaxis, :]
kernel_vals = nucleo_calor(diffs, T_mesh[:, :, np.newaxis])
u0_y = condicion_inicial_escalon(y_int)[np.newaxis, np.newaxis, :]


U_solution = np.trapezoid(kernel_vals * u0_y, x=y_int, axis=2)


fig = plt.figure(figsize=(9, 7))
ax = fig.add_subplot(111, projection='3d')

surf = ax.plot_surface(X_mesh, T_mesh, U_solution, cmap='viridis',
                       edgecolor='none', alpha=0.9, antialiased=True)


ax.set_title(r'Evolución Espacio-Temporal de la Ecuación del Calor por Convolución en $\mathbb{R}$')
ax.set_xlabel('Posición x')
ax.set_ylabel('Tiempo t')
ax.set_zlabel('Temperatura u(x,t)')


fig.colorbar(surf, ax=ax, shrink=0.5, aspect=10, label='Temperatura')


ax.view_init(elev=30, azim=-120)

plt.tight_layout()

plt.show()

#Velocidad Infinita de Propagación

x_lejos = np.linspace(1.5, 10, 300)

plt.figure(figsize=(10, 6))
for t in tiempos:
    u_t_lejos = resolver_calor_integral(x_lejos, t, condicion_inicial_escalon)
    plt.plot(x_lejos, u_t_lejos + 1e-300, label=f't = {t}')

plt.yscale('log')
plt.title('Velocidad Infinita de Propagación')
plt.xlabel('Posición x')
plt.ylabel('Temperatura u(x,t) - Escala Logarítmica')
plt.legend()
plt.grid(True, which="both", ls="--", alpha=0.5)
plt.tight_layout()
plt.show()
plt.close()

puntos_obs = [-1.0, -0.1, 0.5, 1.1, 2.0]
tabla_datos = {'Posición (x)': puntos_obs, 'u(x, t=0)': [0.0, 0.0, 1.0, 0.0, 0.0]}

for t in tiempos:
    valores = resolver_calor_integral(np.array(puntos_obs), t, condicion_inicial_escalon)
    tabla_datos[f't={t}'] = valores

df = pd.DataFrame(tabla_datos)
print("\n--- Tabla de Propagación ---")
print(df.to_string(index=False))
df.to_csv('tabla_propagacion.csv', index=False)

#Comparación del decaimiento

def max_R_t(t):
    y = np.linspace(-5, 6, 2000)
    integrando = nucleo_calor(0.5 - y, t) * condicion_inicial_escalon(y)
    return np.trapezoid(integrando, y)

def max_acotado_t(t, L=5.0, X_centro=0.5, terms=50):
    val = 0
    for n in range(1, terms + 1):
        bn = (2.0 / (n * np.pi)) * (1 - np.cos(n * np.pi * 1.0 / L))
        val += bn * np.exp(-((n * np.pi / L)**2) * t) * np.sin(n * np.pi * X_centro / L)
    return val

def norm_L2_R_t(t):
    x_vals = np.linspace(-15, 15, 1000)
    y_vals = np.linspace(-1, 2, 400)

    X, Y = np.meshgrid(x_vals, y_vals)
    K = nucleo_calor(X - Y, t)
    u0 = condicion_inicial_escalon(Y)

    u_vals = np.trapezoid(K * u0, x=y_vals, axis=0)

    integral_u2 = np.trapezoid(u_vals**2, x=x_vals)
    return np.sqrt(integral_u2)

def norm_L2_acotado_t(t, L=5.0, terms=50):
    norm_cuadrado = 0
    for n in range(1, terms + 1):
        bn = (2.0 / (n * np.pi)) * (1 - np.cos(n * np.pi * 1.0 / L))
        lambda_n = (n * np.pi / L)**2
        norm_cuadrado += (bn**2) * np.exp(-2 * lambda_n * t)
    return np.sqrt((L / 2.0) * norm_cuadrado)

tiempos = np.linspace(0.1, 50, 100)
max_R = [max_R_t(t) for t in tiempos]
max_acotado = [max_acotado_t(t) for t in tiempos]


plt.figure(figsize=(9, 6))


plt.plot(tiempos, max_R, label=r'Dominio $\mathbb{R}$ (Decaimiento $t^{-1/2}$)', color='blue', lw=2)
plt.plot(tiempos, max_acotado, label=r'Dominio Acotado (Decaimiento $e^{-\lambda_1 t}$)', color='red', lw=2)


plt.yscale('log')


plt.title('Comparativa de Decaimiento de la Temperatura Máxima')
plt.xlabel('Tiempo t')
plt.ylabel('Temperatura Máxima - Escala Logarítmica')
plt.grid(True, which="both", ls="--", alpha=0.5)
plt.legend()
plt.tight_layout()

plt.show()

tiempos = np.linspace(0.1, 50, 100)
norma_R = [norm_L2_R_t(t) for t in tiempos]
norma_acotado = [norm_L2_acotado_t(t) for t in tiempos]


plt.figure(figsize=(9, 6))

plt.plot(tiempos, norma_R, label=r'Dominio $\mathbb{R}$ (Decaimiento $t^{-1/4}$)', color='blue', lw=2)
plt.plot(tiempos, norma_acotado, label=r'Dominio Acotado (Decaimiento $e^{-\lambda_1 t}$)', color='red', lw=2)

plt.yscale('log')
plt.title('Comparativa de Decaimiento en Norma $L^2$')
plt.xlabel('Tiempo t')
plt.ylabel('Norma $L^2$ - Escala Logarítmica')
plt.grid(True, which="both", ls="--", alpha=0.5)
plt.legend()
plt.tight_layout()

plt.show()

#Visualización del Caso bidimensional

def nucleo_calor_2d(x, y, t):
    return (1.0 / (4 * np.pi * t)) * np.exp(-(x**2 + y**2) / (4 * t))


def u0_gaussiana_2d(x, y):
    return np.exp(-(x**2 + y**2))

def u0_escalon_2d(x, y):
    return np.where((np.abs(x) <= 0.5) & (np.abs(y) <= 0.5), 1.0, 0.0)


def resolver_calor_convolucion_2d_trapecio(X, Y, t, u0_func,
                                           xi_min=-2, xi_max=2,
                                           eta_min=-2, eta_max=2,
                                           dxi=0.05, deta=0.05):

    xi = np.arange(xi_min, xi_max, dxi)
    eta = np.arange(eta_min, eta_max, deta)
    XI, ETA = np.meshgrid(xi, eta)

    u0_vals = u0_func(XI, ETA)

    Z = np.zeros_like(X)

    for i in range(X.shape[0]):
        for j in range(X.shape[1]):

            kernel = nucleo_calor_2d(X[i, j] - XI, Y[i, j] - ETA, t)
            integrando = kernel * u0_vals

            integral_eta = np.trapezoid(integrando, eta, axis=0)
            integral_total = np.trapezoid(integral_eta, xi)

            Z[i, j] = integral_total

    return Z

x_vals = np.linspace(-1, 1, 100)
y_vals = np.linspace(-1, 1, 100)
X, Y = np.meshgrid(x_vals, y_vals)

tiempos = [0.01, 0.1, 1]


fig = plt.figure(figsize=(18, 8))

for i, t in enumerate(tiempos):

    L = 5 * np.sqrt(t)

    Z = resolver_calor_convolucion_2d_trapecio(
        X, Y, t,
        u0_gaussiana_2d,
        xi_min=-L, xi_max=L,
        eta_min=-L, eta_max=L,
        dxi=0.05, deta=0.05
    )

    ax = fig.add_subplot(1, 3, i+1, projection='3d')

    surf = ax.plot_surface(X, Y, Z, cmap='inferno', edgecolor='none', alpha=0.9)

    ax.set_title(f'Tiempo t = {t}', fontsize=12)
    ax.set_xlabel('$x_1$', fontsize=14)
    ax.set_ylabel('$x_2$', fontsize=14)
    ax.set_zlabel('$u(x,t)$', fontsize=14)

    ax.tick_params(axis='x', labelsize=12)
    ax.tick_params(axis='y', labelsize=12)
    ax.tick_params(axis='z', labelsize=12)

    ax.view_init(elev=25, azim=-45)

plt.suptitle(r'Solución por Convolución de la Ecuación del Calor en $\mathbb{R}^2$', fontsize=20, y=0.98)
plt.tight_layout(rect=[0, 0.03, 1, 0.95], w_pad=4.0)

plt.show()