Diferencia entre revisiones de «Ecuación del Calor LÁJ»

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import trapezoid
from matplotlib import cm

# Ecuación del calor acotada-----------------------------------------------------------------------
def Eca_Acotada(L,alpha,n_terms,f):

    """
    Aproxima la solución de la ecuación del calor u_t -alpha * u_xx = 0, con 
    condición inicial f(x) y condiciones frontera dirichlet homogéneas

    Inputs: 
        L: longitud de la barra
        alpha: constante
        n_terms: número de términos para la serie de Fourier
        f: condición inicial
    """

    # 1. Preparar la malla con T en el eje "X" y X en el eje "Y"
    t_plot = np.linspace(0, 0.1, 100)
    x_plot = np.linspace(0, L, 100)
    T, X = np.meshgrid(t_plot, x_plot)

    # 2. Cálculo de coeficientes a_n con la fórmula del trapecio
    x_fine = np.linspace(0, L, 1000)
    y_init = f(x_fine)
    a = np.zeros(n_terms + 1)
    for n in range(1, n_terms + 1):
        integrando = y_init * np.sin(n * np.pi * x_fine / L)
        a[n] = (2/L) * trapezoid(integrando, x=x_fine)

    # 3. Calcular la solución U(T, X)
    U = np.zeros_like(T)
    for n in range(1, n_terms + 1):
        term_exp = np.exp(-alpha * (n * np.pi / L)**2 * T)
        term_sin = np.sin(n * np.pi * X / L)
        U += a[n] * term_exp * term_sin

    # 4. Graficar
    fig = plt.figure(figsize=(12, 8))
    ax = fig.add_subplot(111, projection='3d')

    # Ahora graficamos T en el primer eje y X en el segundo
    surf = ax.plot_surface(T, X, U, cmap=cm.inferno, linewidth=0, antialiased=True)

    ax.set_xlabel('Tiempo (t)')
    ax.set_ylabel('Posición (x)')
    ax.set_zlabel('Temperatura (u)')
    ax.set_title('Visualización: Tiempo vs Posición')
    fig.colorbar(surf, shrink=0.5, aspect=5)

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

    plt.show()


# Parámetros 
L = 1.0        
alpha = 1.0    
n_terms = 50   

# Dato inicial
def f(x):
    return np.where((x >= 0.3) & (x <= 0.7), 10.0, 0.0)
# def f(x):
#     return 1.0 - 2.0 * np.abs(0.5 - x)


# LLamamos a la función
Eca_Acotada(L,alpha,n_terms,f)


# Ecuación del calor no acotada-----------------------------------------------------------------------
# Visualizar como la solución fundamental tiende a la delta de Dirac------------------------
# Parámetros
alpha = 1.0  # Difusividad térmica

# Definición de la solución fundamental
def solucion_fundamental(x, t, alpha):
    """
    Calcula el valor de la solución fundamental en (x, t).
    Representa la distribución de calor partiendo de una masa puntual en x=0.
    """
    # Evitamos división por cero si t=0
    if t <= 0:
        return np.zeros_like(x)
    
    coeficiente = 1.0 / np.sqrt(4 * np.pi * alpha * t)
    exponente = - (x**2) / (4 * alpha * t)
    return coeficiente * np.exp(exponente)

# Configuración de la visualización
x = np.linspace(-2, 2, 1000)
# Tiempos que se acercan a cero para ver la Delta de Dirac
tiempos = [0.1, 0.05, 0.01, 0.005, 0.001] 

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

for t in tiempos:
    u = solucion_fundamental(x, t, alpha)
    plt.plot(x, u, label=f't = {t}')

plt.title('Convergencia de la Solución Fundamental a la Delta de Dirac')
plt.xlabel('Posición (x)')
plt.ylabel('Temperatura (u)')
plt.legend()
plt.grid(True, linestyle=':')
plt.show()


#Visualizar la solución fundamental en 3D----------------
# Parámetros
alpha = 1.0

# 1. Crear malla de puntos (T en el primer eje, X en el segundo)
t_plot = np.linspace(0.005, 0.1, 100)
x_plot = np.linspace(-2, 2, 100)
T, X = np.meshgrid(t_plot, x_plot) # Invertimos el orden para el intercambio

# 2. Calcular la solución fundamental
U = (1.0 / np.sqrt(4 * np.pi * alpha * T)) * np.exp(-(X**2) / (4 * alpha * T))

# 3. Graficar
fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(111, projection='3d')

# Graficamos T en el eje X del gráfico, X en el eje Y
surf = ax.plot_surface(T, X, U, cmap=cm.viridis, linewidth=0, antialiased=True)

# Etiquetas corregidas
ax.set_xlabel('Tiempo (t)')
ax.set_ylabel('Posición (x)')
ax.set_zlabel('Temperatura (u)')
ax.set_title('Núcleo de Calor: Evolución del Pico (t -> 0)')

fig.colorbar(surf, shrink=0.5, aspect=5)

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

plt.show()


#Visualizar solución de la ecuación del calor no acotada----------------

# 1. Definición del núcleo y dato inicial
def kernel(x, t, alpha=1.0):
    return (1.0 / np.sqrt(4 * np.pi * alpha * t)) * np.exp(-(x**2) / (4 * alpha * t))

def u0(y):
    # Dato inicial
    return np.where((y >= 0) & (y <= 1), 1.0, 0.0)

# 2. Configuración de la malla
t_range = np.linspace(0.001, 1, 50)  # Eje T
x_range = np.linspace(-2, 3, 60)      # Eje X
T, X = np.meshgrid(t_range, x_range)

# Espacio de integración para la convolución (eje y)
y_int = np.linspace(-4, 5, 500)
u_initial_values = u0(y_int)

# 3. Cálculo de la superficie U(T, X)
U = np.zeros_like(X)

for i in range(len(x_range)):
    for j in range(len(t_range)):
        curr_x = X[i, j]
        curr_t = T[i, j]
        # Realizamos la convolución numérica para cada punto (x, t)
        integrando = kernel(curr_x - y_int, curr_t) * u_initial_values
        U[i, j] = trapezoid(integrando, x=y_int)

# 4. Graficar en 3D
fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(111, projection='3d')

# Graficamos con Tiempo en el primer eje y Posición en el segundo
surf = ax.plot_surface(T, X, U, cmap=cm.plasma, linewidth=0, antialiased=True)

ax.set_xlabel('Tiempo (t)')
ax.set_ylabel('Posición (x)')
ax.set_zlabel('Temperatura (u)')
ax.set_title('Solución por Convolución: Disipación de un bloque de calor')

fig.colorbar(surf, shrink=0.5, aspect=5)
ax.view_init(elev=30, azim=-130)

plt.show()

#Ver como el máximo va disminuyendo a lo largo del tiempo
maximos = [np.max(U[:, j]) for j in range(len(t_range))]

plt.figure(figsize=(8, 4))
plt.plot(t_range, maximos, color='red', linewidth=2)
plt.axhline(y=1.0, color='black', linestyle='--', label='Máximo Inicial (t=0)')
plt.title('Validación del Principio del Máximo')
plt.xlabel('Tiempo (t)')
plt.ylabel('Temperatura Máxima alcanzada')
plt.legend()
plt.grid(True)
plt.show()